diff --git a/.gitignore b/.gitignore index 22c568409..2a118ef29 100644 --- a/.gitignore +++ b/.gitignore @@ -1,6 +1,7 @@ *.pyc *.mod *.o +*.hdf5 *.exe *.bak *~ diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index 2464ae7d7..f1af6259f 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -343,6 +343,15 @@ Spectral_MPI: - master - release +SpectralAll_restartMPI: + stage: spectral + script: + - module load $IntelCompiler $MPICH_Intel $PETSc_MPICH_Intel + - SpectralAll_restartMPI/test.py + except: + - master + - release + Plasticity_DetectChanges: stage: spectral script: Plasticity_DetectChanges/test.py diff --git a/installation/mods_Abaqus/abaqus_v6.env b/installation/mods_Abaqus/abaqus_v6.env index d09257a9d..0b4a7fd43 100644 --- a/installation/mods_Abaqus/abaqus_v6.env +++ b/installation/mods_Abaqus/abaqus_v6.env @@ -12,9 +12,17 @@ # import os, re, glob, driverUtils from damask import version as DAMASKVERSION +from damask import Environment +myEnv = Environment() -# Use the version in $PATH -fortCmd = "ifort" +if myEnv.options['DAMASK_HDF5'] == 'ON': + # 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" +else: + # Use the version in $PATH + fortCmd = "ifort" # -free to use free-format FORTRAN 90 syntax # -O <0-3> optimization level @@ -50,4 +58,6 @@ ask_delete=OFF # Remove the temporary names from the namespace del fortCmd +del Environment +del myEnv del DAMASKVERSION diff --git a/installation/mods_Abaqus/abaqus_v6_debug.env b/installation/mods_Abaqus/abaqus_v6_debug.env index de5189a52..c967c1e65 100644 --- a/installation/mods_Abaqus/abaqus_v6_debug.env +++ b/installation/mods_Abaqus/abaqus_v6_debug.env @@ -12,9 +12,17 @@ # import os, re, glob, driverUtils from damask import version as DAMASKVERSION +from damask import Environment +myEnv = Environment() -# Use the version in $PATH -fortCmd = "ifort" +if myEnv.options['DAMASK_HDF5'] == 'ON': + # 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" +else: + # Use the version in $PATH + fortCmd = "ifort" # -free to use free-format FORTRAN 90 syntax # -O <0-3> optimization level @@ -55,4 +63,6 @@ ask_delete=OFF # Remove the temporary names from the namespace del fortCmd +del Environment +del myEnv del DAMASKVERSION diff --git a/installation/mods_MarcMentat/2018.1/Marc_tools/include_linux64 b/installation/mods_MarcMentat/2018.1/Marc_tools/include_linux64 index 22b306cb1..53eef9d83 100644 --- a/installation/mods_MarcMentat/2018.1/Marc_tools/include_linux64 +++ b/installation/mods_MarcMentat/2018.1/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 -DDAMASKHDF5" + echo $FCOMP +else + FCOMP=ifort +fi + # AEM if test "$MARCDLLOUTDIR" = ""; then DLLOUTDIR="$MARC_LIB" @@ -535,23 +544,17 @@ else DAMASKVERSION="'N/A'" fi -if test "$DAMASK_HDF5" = "ON";then - DFCOMP="$(h5fc -show) -DDAMASKHDF5" -else - DFCOMP=$FCOMP -fi -# # DAMASK compiler calls: additional flags are in line 2 OpenMP flags in line 3 -DFORTLOWMP="$DFCOMP -c -implicitnone -stand f08 -standard-semantics -assume nostd_mod_proc_name -safe_cray_ptr $PROFILE -zero -mp1 -WB -O0 $I8FFLAGS -I$MARC_SOURCE/common \ +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.1 -DDAMASKVERSION=$DAMASKVERSION \ -qopenmp -qopenmp-threadprivate=compat\ $MUMPS_INCLUDE $I8DEFINES -DLinux -DLINUX -DLinux_intel $FDEFINES $DDM $SOLVERFLAGS -I$KDTREE2_MOD" -DFORTRANMP="$DFCOMP -c -implicitnone -stand f08 -standard-semantics -assume nostd_mod_proc_name -safe_cray_ptr $PROFILE -zero -mp1 -WB -O1 $I8FFLAGS -I$MARC_SOURCE/common \ +DFORTRANMP="$FCOMP -c -implicitnone -stand f08 -standard-semantics -assume nostd_mod_proc_name -safe_cray_ptr $PROFILE -zero -mp1 -WB -O1 $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.1 -DDAMASKVERSION=$DAMASKVERSION \ -qopenmp -qopenmp-threadprivate=compat\ $MUMPS_INCLUDE $I8DEFINES -DLinux -DLINUX -DLinux_intel $FDEFINES $DDM $SOLVERFLAGS -I$KDTREE2_MOD" -DFORTHIGHMP="$DFCOMP -c -implicitnone -stand f08 -standard-semantics -assume nostd_mod_proc_name -safe_cray_ptr $PROFILE -zero -mp1 -WB -fno-alias -O2 $I8FFLAGS -I$MARC_SOURCE/common \ +DFORTHIGHMP="$FCOMP -c -implicitnone -stand f08 -standard-semantics -assume nostd_mod_proc_name -safe_cray_ptr $PROFILE -zero -mp1 -WB -fno-alias -O2 $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.1 -DDAMASKVERSION=$DAMASKVERSION \ -qopenmp -qopenmp-threadprivate=compat\ $MUMPS_INCLUDE $I8DEFINES -DLinux -DLINUX -DLinux_intel $FDEFINES $DDM $SOLVERFLAGS -I$KDTREE2_MOD" @@ -570,15 +573,15 @@ then fi # DAMASK compiler calls: additional flags are in line 2 OpenMP flags in line 3 - DFORTLOWMP="$DFCOMP -c -implicitnone -stand f08 -standard-semantics -assume nostd_mod_proc_name -safe_cray_ptr $PROFILE -zero -mp1 -WB $I8FFLAGS -I$MARC_SOURCE/common \ + DFORTLOWMP="$FCOMP -c -implicitnone -stand f08 -standard-semantics -assume nostd_mod_proc_name -safe_cray_ptr $PROFILE -zero -mp1 -WB $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.1 -DDAMASKVERSION=$DAMASKVERSION \ -qopenmp -qopenmp-threadprivate=compat\ $MUMPS_INCLUDE $I8DEFINES -DLinux -DLINUX -DLinux_intel $FDEFINES $DDM $SOLVERFLAGS -I$KDTREE2_MOD" - DFORTRANMP="$DFCOMP -c -implicitnone -stand f08 -standard-semantics -assume nostd_mod_proc_name -safe_cray_ptr $PROFILE -zero -mp1 -WB $I8FFLAGS -I$MARC_SOURCE/common \ + DFORTRANMP="$FCOMP -c -implicitnone -stand f08 -standard-semantics -assume nostd_mod_proc_name -safe_cray_ptr $PROFILE -zero -mp1 -WB $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.1 -DDAMASKVERSION=$DAMASKVERSION \ -qopenmp -qopenmp-threadprivate=compat\ $MUMPS_INCLUDE $I8DEFINES -DLinux -DLINUX -DLinux_intel $FDEFINES $DDM $SOLVERFLAGS -I$KDTREE2_MOD" - DFORTHIGHMP="$DFCOMP -c -implicitnone -stand f08 -standard-semantics -assume nostd_mod_proc_name -safe_cray_ptr $PROFILE -zero -mp1 -WB -fno-alias $I8FFLAGS -I$MARC_SOURCE/common \ + DFORTHIGHMP="$FCOMP -c -implicitnone -stand f08 -standard-semantics -assume nostd_mod_proc_name -safe_cray_ptr $PROFILE -zero -mp1 -WB -fno-alias $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.1 -DDAMASKVERSION=$DAMASKVERSION \ -qopenmp -qopenmp-threadprivate=compat\ $MUMPS_INCLUDE $I8DEFINES -DLinux -DLINUX -DLinux_intel $FDEFINES $DDM $SOLVERFLAGS -I$KDTREE2_MOD" @@ -744,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/processing/post/ascii2hdf5.py b/processing/post/ascii2hdf5.py deleted file mode 100755 index effac981b..000000000 --- a/processing/post/ascii2hdf5.py +++ /dev/null @@ -1,142 +0,0 @@ -#!/usr/bin/env python2.7 -# -*- coding: UTF-8 no BOM -*- - -# ------------------------------------------------------------------- # -# NOTE: # -# 1. Not all output is defined in the DS_HDF5.xml, please add new # -# new one to the system wide definition file # -# /lib/damask/DS_HDF5.xml # -# or specify your own when initializing HDF5 class # -# 2. Somehow the point cloud structure cannot be properly handled # -# by Xdmf, which is a descriptive wrapper for visualizing HDF5 # -# using Paraview. The current solution is using cell structured # -# HDF5 so that Xdmf can describe the data shape as a rectangular # -# mesh rather than polyvertex. # -# TODO: # -# 1. remove the ._tmp file, basically need a way to # -# just load data from ASCII table. # -# 2. a progress monitor when transferring data from ASCII table # -# to HDF5. # -# 3. a more flexible way handle the data structure rather than a # -# xml file. # -# ------------------------------------------------------------------- # - -import os -import damask -import numpy as np -from optparse import OptionParser - - -scriptName = os.path.splitext(os.path.basename(__file__))[0] -scriptID = ' '.join([scriptName, damask.version]) - - -# ----- helper function ----- # -def get_rectMshVectors(xyz_array, posNum): - """Get Vx, Vy, Vz for rectLinear grid""" - # need some improvement, and only works for rectangular grid - v = sorted(list(set(xyz_array[:, posNum]))) - v_interval = (v[2]+v[1])/2.0 - (v[1]+v[0])/2.0 - v_start = (v[1]+v[0])/2.0 - v_interval - v_end = (v[-1]+v[-2])/2.0 + v_interval - V = np.linspace(v_start, v_end, len(v)+1) - return V - - -# ----- MAIN ---- # -desp_msg = "Convert DAMASK ascii table to HDF5 file" -parser = OptionParser(option_class=damask.extendableOption, - usage='%prog options [file[s]]', - description=desp_msg, - version=scriptID) -parser.add_option('-D', '--DefinitionFile', - dest='storage definition file', - type='string', - metavar='string', - help='definition file for H5 data storage') -parser.add_option('-p', '--pos', '--position', - dest='pos', - type='string', metavar='string', - help='label of coordinates [%default]') - -parser.set_defaults(DefinitionFile='default', - pos='pos') - -(options, filenames) = parser.parse_args() - -filename = filenames[0] - -if options.DefinitionFile == 'default': - defFile = None -else: - defFile = options.DefinitionFile - -# ----- read in data using DAMASK ASCII table class ----- # -asciiTable = damask.ASCIItable(name=filename, buffered=False) -asciiTable.head_read() -asciiTable.data_readArray() -incNum = int(asciiTable.data[asciiTable.label_index('inc'), 0]) -fullTable = np.copy(asciiTable.data) # deep copy all data, just to be safe -labels = asciiTable.labels() -labels_idx = [asciiTable.label_index(label) for label in labels] -featuresDim = [labels_idx[i+1] - labels_idx[i] for i in range(len(labels)-1)] -featuresDim.append(fullTable.shape[1] - labels_idx[-1]) - -# ----- figure out size and grid ----- # -pos_idx = asciiTable.label_index('pos') -xyz_array = asciiTable.data[:, pos_idx:pos_idx+3] -Vx = get_rectMshVectors(xyz_array, 0) -Vy = get_rectMshVectors(xyz_array, 1) -Vz = get_rectMshVectors(xyz_array, 2) -# use the dimension of the rectangular grid to reshape all other data -mshGridDim = [len(Vx)-1, len(Vy)-1, len(Vz)-1] - -# ----- compose cmd log ----- # -cmd_log = " ".join([scriptID, filename]) - -# ----- create a new HDF5 file and save the data -----# -# force remove existing HDF5 file -h5fName = filename.replace(".txt", ".h5") -try: - os.remove(h5fName) -except OSError: - pass -h5f = damask.H5Table(h5fName, - new_file=True, - dsXMLFile=defFile) -# adding increment number as root level attributes -h5f.add_attr('inc', incNum) -# add the mesh grid data now -h5f.add_data("Vx", Vx, cmd_log=cmd_log) -h5f.add_data("Vy", Vy, cmd_log=cmd_log) -h5f.add_data("Vz", Vz, cmd_log=cmd_log) - -# add the rest of data from table -labelsProcessed = ['inc'] -for fi in range(len(labels)): - featureName = labels[fi] - # remove trouble maker "("" and ")" from label/feature name - if "(" in featureName: - featureName = featureName.replace("(", "") - if ")" in featureName: - featureName = featureName.replace(")", "") - # skip increment and duplicated columns in the ASCII table - if featureName in labelsProcessed: - continue - - featureIdx = labels_idx[fi] - featureDim = featuresDim[fi] - # grab the data hook - dataset = fullTable[:, featureIdx:featureIdx+featureDim] - # mapping 2D data onto a 3D rectangular mesh to get 4D data - # WARNING: In paraview, the data for a recmesh is mapped as: - # --> len(z), len(y), len(x), size(data) - # dataset = dataset.reshape((mshGridDim[0], - # mshGridDim[1], - # mshGridDim[2], - # dataset.shape[1])) - # write out data - print("adding {}...".format(featureName)) - h5f.add_data(featureName, dataset, cmd_log=cmd_log) - # write down the processed label - labelsProcessed.append(featureName) diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index a84609087..23c7a5643 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -29,24 +29,28 @@ add_library(IO OBJECT "IO.f90") add_dependencies(IO DAMASK_INTERFACE) list(APPEND OBJECTFILES $) -add_library(HDF5_UTILITIES OBJECT "HDF5_utilities.f90") -add_dependencies(HDF5_UTILITIES IO) -list(APPEND OBJECTFILES $) - add_library(NUMERICS OBJECT "numerics.f90") -add_dependencies(NUMERICS HDF5_UTILITIES) +add_dependencies(NUMERICS IO) list(APPEND OBJECTFILES $) add_library(DEBUG OBJECT "debug.f90") add_dependencies(DEBUG NUMERICS) list(APPEND OBJECTFILES $) -add_library(CONFIG OBJECT "config.f90") -add_dependencies(CONFIG DEBUG) -list(APPEND OBJECTFILES $) +add_library(DAMASK_CONFIG OBJECT "config.f90") +add_dependencies(DAMASK_CONFIG DEBUG) +list(APPEND OBJECTFILES $) + +add_library(HDF5_UTILITIES OBJECT "HDF5_utilities.f90") +add_dependencies(HDF5_UTILITIES DAMASK_CONFIG) +list(APPEND OBJECTFILES $) + +add_library(RESULTS OBJECT "results.f90") +add_dependencies(RESULTS HDF5_UTILITIES) +list(APPEND OBJECTFILES $) add_library(FEsolving OBJECT "FEsolving.f90") -add_dependencies(FEsolving DEBUG) +add_dependencies(FEsolving RESULTS) list(APPEND OBJECTFILES $) add_library(DAMASK_MATH OBJECT "math.f90") @@ -68,7 +72,7 @@ elseif (PROJECT_NAME STREQUAL "DAMASK_FEM") endif() add_library(MATERIAL OBJECT "material.f90") -add_dependencies(MATERIAL MESH CONFIG) +add_dependencies(MATERIAL MESH DAMASK_CONFIG) list(APPEND OBJECTFILES $) add_library(DAMASK_HELPERS OBJECT "lattice.f90") diff --git a/src/CPFEM2.f90 b/src/CPFEM2.f90 index 29e1ac744..75f57f4c2 100644 --- a/src/CPFEM2.f90 +++ b/src/CPFEM2.f90 @@ -10,8 +10,8 @@ module CPFEM2 public :: & CPFEM_age, & - CPFEM_initAll - + CPFEM_initAll, & + CPFEM_results contains @@ -20,8 +20,7 @@ contains !-------------------------------------------------------------------------------------------------- subroutine CPFEM_initAll() use prec, only: & - pInt - use prec, only: & + pInt, & prec_init use numerics, only: & numerics_init @@ -39,6 +38,8 @@ subroutine CPFEM_initAll() material_init use HDF5_utilities, only: & HDF5_utilities_init + use results, only: & + results_init use lattice, only: & lattice_init use constitutive, only: & @@ -73,6 +74,7 @@ subroutine CPFEM_initAll() call lattice_init call material_init call HDF5_utilities_init + call results_init call constitutive_init call crystallite_init call homogenization_init @@ -105,8 +107,7 @@ subroutine CPFEM_init debug_levelBasic, & debug_levelExtensive use FEsolving, only: & - restartRead, & - modelName + restartRead use material, only: & material_phase, & homogState, & @@ -125,23 +126,22 @@ subroutine CPFEM_init use hdf5 use HDF5_utilities, only: & HDF5_openFile, & - HDF5_openGroup2, & + HDF5_closeFile, & + HDF5_openGroup, & + HDF5_closeGroup, & HDF5_read use DAMASK_interface, only: & getSolverJobName implicit none - integer(pInt) :: k,l,m,ph,homog + integer(pInt) :: ph,homog character(len=1024) :: rankStr, PlasticItem, HomogItem - integer(HID_T) :: fileReadID, groupPlasticID, groupHomogID - integer :: hdferr + integer(HID_T) :: fileHandle, groupPlasticID, groupHomogID - mainProcess: if (worldrank == 0) then - write(6,'(/,a)') ' <<<+- CPFEM init -+>>>' - write(6,'(a15,a)') ' Current time: ',IO_timeStamp() + write(6,'(/,a)') ' <<<+- CPFEM init -+>>>' + write(6,'(a15,a)') ' Current time: ',IO_timeStamp() #include "compilation_info.f90" - flush(6) - endif mainProcess + flush(6) ! *** restore the last converged values of each essential variable from the binary file if (restartRead) then @@ -152,34 +152,39 @@ subroutine CPFEM_init write(rankStr,'(a1,i0)')'_',worldrank - fileReadID = HDF5_openFile(trim(getSolverJobName())//trim(rankStr)//'.hdf5') + fileHandle = HDF5_openFile(trim(getSolverJobName())//trim(rankStr)//'.hdf5') - call HDF5_read(material_phase, fileReadID,'recordedPhase') - call HDF5_read(crystallite_F0, fileReadID,'convergedF') - call HDF5_read(crystallite_Fp0, fileReadID,'convergedFp') - call HDF5_read(crystallite_Fi0, fileReadID,'convergedFi') - call HDF5_read(crystallite_Lp0, fileReadID,'convergedLp') - call HDF5_read(crystallite_Li0, fileReadID,'convergedLi') - call HDF5_read(crystallite_dPdF0, fileReadID,'convergeddPdF') - call HDF5_read(crystallite_Tstar0_v,fileReadID,'convergedTstar') + call HDF5_read(fileHandle,material_phase, 'recordedPhase') + call HDF5_read(fileHandle,crystallite_F0, 'convergedF') + call HDF5_read(fileHandle,crystallite_Fp0, 'convergedFp') + call HDF5_read(fileHandle,crystallite_Fi0, 'convergedFi') + call HDF5_read(fileHandle,crystallite_Lp0, 'convergedLp') + call HDF5_read(fileHandle,crystallite_Li0, 'convergedLi') + call HDF5_read(fileHandle,crystallite_dPdF0, 'convergeddPdF') + call HDF5_read(fileHandle,crystallite_Tstar0_v,'convergedTstar') - groupPlasticID = HDF5_openGroup2(fileReadID,'PlasticPhases') + groupPlasticID = HDF5_openGroup(fileHandle,'PlasticPhases') do ph = 1_pInt,size(phase_plasticity) write(PlasticItem,*) ph,'_' - call HDF5_read(plasticState(ph)%state0,groupPlasticID,trim(PlasticItem)//'convergedStateConst') + call HDF5_read(groupPlasticID,plasticState(ph)%state0,trim(PlasticItem)//'convergedStateConst') enddo + call HDF5_closeGroup(groupPlasticID) - groupHomogID = HDF5_openGroup2(fileReadID,'HomogStates') + groupHomogID = HDF5_openGroup(fileHandle,'HomogStates') do homog = 1_pInt, material_Nhomogenization write(HomogItem,*) homog,'_' - call HDF5_read(homogState(homog)%state0, groupHomogID,trim(HomogItem)//'convergedStateHomog') + call HDF5_read(groupHomogID,homogState(homog)%state0, trim(HomogItem)//'convergedStateHomog') enddo + call HDF5_closeGroup(groupHomogID) + + call HDF5_closeFile(fileHandle) restartRead = .false. endif end subroutine CPFEM_init + !-------------------------------------------------------------------------------------------------- !> @brief forwards data after successful increment !-------------------------------------------------------------------------------------------------- @@ -229,82 +234,102 @@ subroutine CPFEM_age() use HDF5_utilities, only: & HDF5_openFile, & HDF5_closeFile, & + HDF5_addGroup, & HDF5_closeGroup, & - HDF5_addGroup2, & HDF5_write use hdf5 use DAMASK_interface, only: & getSolverJobName implicit none - - integer(pInt) :: i, k, l, m, ph, homog, mySource + integer(pInt) :: i, ph, homog, mySource character(len=32) :: rankStr, PlasticItem, HomogItem integer(HID_T) :: fileHandle, groupPlastic, groupHomog - integer :: hdferr - integer(HSIZE_T) :: hdfsize -if (iand(debug_level(debug_CPFEM), debug_levelBasic) /= 0_pInt) & - write(6,'(a)') '<< CPFEM >> aging states' + if (iand(debug_level(debug_CPFEM), debug_levelBasic) /= 0_pInt) & + write(6,'(a)') '<< CPFEM >> aging states' - crystallite_F0 = crystallite_partionedF ! crystallite deformation (_subF is perturbed...) - crystallite_Fp0 = crystallite_Fp ! crystallite plastic deformation - crystallite_Lp0 = crystallite_Lp ! crystallite plastic velocity - crystallite_Fi0 = crystallite_Fi ! crystallite intermediate deformation - crystallite_Li0 = crystallite_Li ! crystallite intermediate velocity - crystallite_dPdF0 = crystallite_dPdF ! crystallite stiffness - crystallite_Tstar0_v = crystallite_Tstar_v ! crystallite 2nd Piola Kirchhoff stress - - forall (i = 1:size(plasticState)) plasticState(i)%state0 = plasticState(i)%state ! copy state in this lengthy way because: A component cannot be an array if the encompassing structure is an array - - do i = 1, size(sourceState) - do mySource = 1,phase_Nsources(i) - sourceState(i)%p(mySource)%state0 = sourceState(i)%p(mySource)%state ! copy state in this lengthy way because: A component cannot be an array if the encompassing structure is an array - enddo; enddo - - do homog = 1_pInt, material_Nhomogenization - homogState (homog)%state0 = homogState (homog)%state - thermalState (homog)%state0 = thermalState (homog)%state - damageState (homog)%state0 = damageState (homog)%state - enddo + crystallite_F0 = crystallite_partionedF + crystallite_Fp0 = crystallite_Fp + crystallite_Lp0 = crystallite_Lp + crystallite_Fi0 = crystallite_Fi + crystallite_Li0 = crystallite_Li + crystallite_dPdF0 = crystallite_dPdF + crystallite_Tstar0_v = crystallite_Tstar_v + + forall (i = 1:size(plasticState)) plasticState(i)%state0 = plasticState(i)%state ! copy state in this lengthy way because: A component cannot be an array if the encompassing structure is an array + + do i = 1, size(sourceState) + do mySource = 1,phase_Nsources(i) + sourceState(i)%p(mySource)%state0 = sourceState(i)%p(mySource)%state ! copy state in this lengthy way because: A component cannot be an array if the encompassing structure is an array + enddo; enddo + + do homog = 1_pInt, material_Nhomogenization + homogState (homog)%state0 = homogState (homog)%state + thermalState (homog)%state0 = thermalState (homog)%state + damageState (homog)%state0 = damageState (homog)%state + enddo -if (restartWrite) then - if (iand(debug_level(debug_CPFEM), debug_levelBasic) /= 0_pInt) & - write(6,'(a)') '<< CPFEM >> writing restart variables of last converged step to hdf5 file' - write(rankStr,'(a1,i0)')'_',worldrank + if (restartWrite) then + if (iand(debug_level(debug_CPFEM), debug_levelBasic) /= 0_pInt) & + write(6,'(a)') '<< CPFEM >> writing restart variables of last converged step to hdf5 file' + + write(rankStr,'(a1,i0)')'_',worldrank + fileHandle = HDF5_openFile(trim(getSolverJobName())//trim(rankStr)//'.hdf5','w') + + call HDF5_write(fileHandle,material_phase, 'recordedPhase') + call HDF5_write(fileHandle,crystallite_F0, 'convergedF') + call HDF5_write(fileHandle,crystallite_Fp0, 'convergedFp') + call HDF5_write(fileHandle,crystallite_Fi0, 'convergedFi') + call HDF5_write(fileHandle,crystallite_Lp0, 'convergedLp') + call HDF5_write(fileHandle,crystallite_Li0, 'convergedLi') + call HDF5_write(fileHandle,crystallite_dPdF0, 'convergeddPdF') + call HDF5_write(fileHandle,crystallite_Tstar0_v,'convergedTstar') + + groupPlastic = HDF5_addGroup(fileHandle,'PlasticPhases') + do ph = 1_pInt,size(phase_plasticity) + write(PlasticItem,*) ph,'_' + call HDF5_write(groupPlastic,plasticState(ph)%state0,trim(PlasticItem)//'convergedStateConst') + enddo + call HDF5_closeGroup(groupPlastic) - fileHandle = HDF5_openFile(trim(getSolverJobName())//trim(rankStr)//'.hdf5','w') - - call HDF5_write(material_phase, fileHandle,'recordedPhase') - call HDF5_write(crystallite_F0, fileHandle,'convergedF') - call HDF5_write(crystallite_Fp0, fileHandle,'convergedFp') - call HDF5_write(crystallite_Fi0, fileHandle,'convergedFi') - call HDF5_write(crystallite_Lp0, fileHandle,'convergedLp') - call HDF5_write(crystallite_Li0, fileHandle,'convergedLi') - call HDF5_write(crystallite_dPdF0, fileHandle,'convergeddPdF') - call HDF5_write(crystallite_Tstar0_v,fileHandle,'convergedTstar') - - groupPlastic = HDF5_addGroup2(fileHandle,'PlasticPhases') - do ph = 1_pInt,size(phase_plasticity) - write(PlasticItem,*) ph,'_' - call HDF5_write(plasticState(ph)%state0,groupPlastic,trim(PlasticItem)//'convergedStateConst') - enddo - call HDF5_closeGroup(groupPlastic) + groupHomog = HDF5_addGroup(fileHandle,'HomogStates') + do homog = 1_pInt, material_Nhomogenization + write(HomogItem,*) homog,'_' + call HDF5_write(groupHomog,homogState(homog)%state0,trim(HomogItem)//'convergedStateHomog') + enddo + call HDF5_closeGroup(groupHomog) + + call HDF5_closeFile(fileHandle) + restartWrite = .false. + endif - groupHomog = HDF5_addGroup2(fileHandle,'HomogStates') - do homog = 1_pInt, material_Nhomogenization - write(HomogItem,*) homog,'_' - call HDF5_write(homogState(homog)%state0,groupHomog,trim(HomogItem)//'convergedStateHomog') - enddo - call HDF5_closeGroup(groupHomog) - - call HDF5_closeFile(fileHandle) - restartWrite = .false. -endif - -if (iand(debug_level(debug_CPFEM), debug_levelBasic) /= 0_pInt) & - write(6,'(a)') '<< CPFEM >> done aging states' + if (iand(debug_level(debug_CPFEM), debug_levelBasic) /= 0_pInt) & + write(6,'(a)') '<< CPFEM >> done aging states' end subroutine CPFEM_age +!-------------------------------------------------------------------------------------------------- +!> @brief triggers writing of the results +!-------------------------------------------------------------------------------------------------- +subroutine CPFEM_results(inc,time) + use prec, only: & + pInt + use results + use HDF5_utilities + use constitutive, only: & + constitutive_results + + implicit none + integer(pInt), intent(in) :: inc + real(pReal), intent(in) :: time + + call results_openJobFile + call results_addIncrement(inc,time) + call constitutive_results() + call results_removeLink('current') ! ToDo: put this into closeJobFile + call results_closeJobFile + +end subroutine CPFEM_results + end module CPFEM2 diff --git a/src/DAMASK_spectral.f90 b/src/DAMASK_spectral.f90 index d6827543a..fca67c97d 100644 --- a/src/DAMASK_spectral.f90 +++ b/src/DAMASK_spectral.f90 @@ -46,7 +46,8 @@ program DAMASK_spectral grid, & geomSize use CPFEM2, only: & - CPFEM_initAll + CPFEM_initAll, & + CPFEM_results use FEsolving, only: & restartWrite, & restartInc @@ -80,6 +81,7 @@ program DAMASK_spectral use spectral_mech_Polarisation use spectral_damage use spectral_thermal + use results implicit none @@ -157,6 +159,9 @@ program DAMASK_spectral write(6,'(a15,a)') ' Current time: ',IO_timeStamp() #include "compilation_info.f90" + + call results_openJobFile() + call results_closeJobFile() !-------------------------------------------------------------------------------------------------- ! initialize field solver information nActiveFields = 1 @@ -420,6 +425,7 @@ program DAMASK_spectral writeUndeformed: if (interface_restartInc < 1_pInt) then write(6,'(1/,a)') ' ... writing initial configuration to file ........................' + call CPFEM_results(0_pInt,0.0_pReal) do i = 1, size(materialpoint_results,3)/(maxByteOut/(materialpoint_sizeResults*pReal))+1 ! slice the output of my process in chunks not exceeding the limit for one output outputIndex = int([(i-1_pInt)*((maxRealOut)/materialpoint_sizeResults)+1_pInt, & ! QUESTION: why not starting i at 0 instead of murky 1? min(i*((maxRealOut)/materialpoint_sizeResults),size(materialpoint_results,3))],pLongInt) @@ -596,6 +602,7 @@ program DAMASK_spectral if(ierr /=0_pInt) call IO_error(894_pInt, ext_msg='MPI_file_write') enddo fileOffset = fileOffset + sum(outputSize) ! forward to current file position + call CPFEM_results(totalIncsCounter,time) endif if ( loadCases(currentLoadCase)%restartFrequency > 0_pInt & ! writing of restart info requested ... .and. mod(inc,loadCases(currentLoadCase)%restartFrequency) == 0_pInt) then ! ... and at frequency of writing restart information diff --git a/src/HDF5_utilities.f90 b/src/HDF5_utilities.f90 index 43a7a26e8..2a05f101c 100644 --- a/src/HDF5_utilities.f90 +++ b/src/HDF5_utilities.f90 @@ -4,39 +4,38 @@ !> @author Jennifer Nastola, Max-Planck-Institut für Eisenforschung GmbH !> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH !-------------------------------------------------------------------------------------------------- -module HDF5_Utilities - use prec - use IO - use HDF5 +module HDF5_utilities + use prec + use IO + use HDF5 #ifdef PETSc - use PETSC + use PETSC #endif implicit none - private - integer(HID_T), public, protected :: tempCoordinates, tempResults - integer(HID_T), private :: resultsFile, currentIncID, plist_id - integer(pInt), private :: currentInc + public + integer(pInt), parameter, private :: & + HDF5_ERR_TYPE = 4_pInt !< kind of the integer return in the HDF5 library !-------------------------------------------------------------------------------------------------- !> @brief reads pInt or pReal data of defined shape from file !-------------------------------------------------------------------------------------------------- interface HDF5_read - module procedure HDF5_read_pReal_1 - module procedure HDF5_read_pReal_2 - module procedure HDF5_read_pReal_3 - module procedure HDF5_read_pReal_4 - module procedure HDF5_read_pReal_5 - module procedure HDF5_read_pReal_6 - module procedure HDF5_read_pReal_7 + 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_pInt_1 - module procedure HDF5_read_pInt_2 - module procedure HDF5_read_pInt_3 - module procedure HDF5_read_pInt_4 - module procedure HDF5_read_pInt_5 - module procedure HDF5_read_pInt_6 - module procedure HDF5_read_pInt_7 !ABOVE 8 DIMENSIONS IT GIVES ERROR: THE CALL TO H5DREAD_F DOESNT WORK + 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 @@ -58,51 +57,46 @@ module HDF5_Utilities module procedure HDF5_write_pInt4 module procedure HDF5_write_pInt5 module procedure HDF5_write_pInt6 - module procedure HDF5_write_pInt7 !ABOVE 8 DIMENSIONS IT GIVES ERROR: THE CALL TO H5DREAD_F DOESNT WORK + module procedure HDF5_write_pInt7 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 + + +!-------------------------------------------------------------------------------------------------- public :: & - HDF5_Utilities_init, & - HDF5_mappingPhase, & - HDF5_mappingHomog, & - HDF5_mappingCrystallite, & - HDF5_backwardMappingPhase, & - HDF5_backwardMappingHomog, & - HDF5_backwardMappingCrystallite, & - HDF5_mappingCells, & - HDF5_addGroup ,& + HDF5_utilities_init, & + HDF5_openFile, & + HDF5_closeFile, & + HDF5_addAttribute, & HDF5_closeGroup ,& HDF5_openGroup, & - HDF5_openGroup2, & - HDF5_forwardResults, & - HDF5_writeVectorDataset, & - HDF5_writeScalarDataset, & - HDF5_writeTensorDataset, & - HDF5_closeJobFile, & - HDF5_removeLink, & - HDF5_createFile, & - HDF5_closeFile, & - HDF5_addGroup2, & - HDF5_openFile, & + HDF5_addGroup, & HDF5_read, & - HDF5_write + HDF5_write, & + HDF5_setLink, & + HDF5_objectExists contains -subroutine HDF5_Utilities_init +subroutine HDF5_utilities_init use, intrinsic :: & iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment) implicit none - integer :: hdferr - integer(SIZE_T) :: typeSize + integer(HDF5_ERR_TYPE) :: hdferr + integer(SIZE_T) :: typeSize write(6,'(/,a)') ' <<<+- HDF5_Utilities init -+>>>' #include "compilation_info.f90" - !currentInc = -1_pInt ToDo - !call HDF5_createJobFile ToDo - !-------------------------------------------------------------------------------------------------- !initialize HDF5 library and check if integer and float type size match call h5open_f(hdferr) @@ -114,97 +108,22 @@ subroutine HDF5_Utilities_init if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_Utilities_init: h5tget_size_f (double)') if (int(pReal,SIZE_T)/=typeSize) call IO_error(0_pInt,ext_msg='pReal does not match H5T_NATIVE_DOUBLE') -end subroutine HDF5_Utilities_init +end subroutine HDF5_utilities_init -!-------------------------------------------------------------------------------------------------- -!> @brief creates and initializes HDF5 output files -!-------------------------------------------------------------------------------------------------- -subroutine HDF5_createJobFile - use hdf5 - use DAMASK_interface, only: & - getSolverJobName - - implicit none - integer :: hdferr - character(len=1024) :: path -#ifdef PETSc -#include - - call h5pcreate_f(H5P_FILE_ACCESS_F, plist_id, hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_Utilities_init: h5pcreate_f') - 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_Utilities_init: h5pset_fapl_mpio_f') -#endif - -!-------------------------------------------------------------------------------------------------- -! open file - path = trim(getSolverJobName())//'.'//'hdf5' - !call h5fcreate_f(path,H5F_ACC_TRUNC_F,resultsFile,hdferr) - call h5fcreate_f(path,H5F_ACC_TRUNC_F,resultsFile,hdferr,access_prp = plist_id) - if (hdferr < 0) call IO_error(100_pInt,ext_msg=path) - call HDF5_addStringAttribute(resultsFile,'createdBy',DAMASKVERSION) - call h5pclose_f(plist_id, hdferr) !neu - -end subroutine HDF5_createJobFile - - -!-------------------------------------------------------------------------------------------------- -!> @brief creates and initializes HDF5 output files -!-------------------------------------------------------------------------------------------------- - integer(HID_T) function HDF5_createFile(path) - use hdf5 - use DAMASK_interface, only: & - getSolverJobName - - implicit none - integer :: hdferr - integer(SIZE_T) :: typeSize - character(len=*), intent(in) :: path -#ifdef PETSc -#include -#endif - call h5open_f(hdferr) !############################################################ DANGEROUS -#ifdef PETSc - call h5pcreate_f(H5P_FILE_ACCESS_F, plist_id, hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_Utilities_init: h5pcreate_f') - 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_Utilities_init: h5pset_fapl_mpio_f') -#endif -!-------------------------------------------------------------------------------------------------- -! create a file - !call h5fcreate_f(path,H5F_ACC_TRUNC_F,resultsFile,hdferr) - call h5fcreate_f(path,H5F_ACC_TRUNC_F,HDF5_createFile,hdferr,access_prp = plist_id) - if (hdferr < 0) call IO_error(100_pInt,ext_msg=path) - !call HDF5_addStringAttribute(HDF5_createFile,'createdBy',DAMASKVERSION) - call h5pclose_f(plist_id, hdferr) !neu - -end function HDF5_createFile - -!-------------------------------------------------------------------------------------------------- -!> @brief close the opened HDF5 output file -!-------------------------------------------------------------------------------------------------- -subroutine HDF5_closeJobFile() - use hdf5 - - implicit none - integer :: hdferr - call HDF5_removeLink('current') - call h5fclose_f(resultsFile,hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_closeJobFile: h5fclose_f',el=hdferr) -! call h5close_f(hdferr) - -end subroutine HDF5_closeJobFile !-------------------------------------------------------------------------------------------------- !> @brief open and initializes HDF5 output file !-------------------------------------------------------------------------------------------------- -integer(HID_T) function HDF5_openFile(fileName,mode) +integer(HID_T) function HDF5_openFile(fileName,mode,parallel) implicit none character(len=*), intent(in) :: fileName character, intent(in), optional :: mode + logical, intent(in), optional :: parallel + character :: m - integer :: hdferr + integer(HID_T) :: plist_id + integer(HDF5_ERR_TYPE) :: hdferr if (present(mode)) then m = mode @@ -212,2130 +131,2938 @@ integer(HID_T) function HDF5_openFile(fileName,mode) m = 'r' endif - if (m == 'w') then - call h5fcreate_f(fileName,H5F_ACC_TRUNC_F,HDF5_openFile,hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_openFile: h5fcreate_f',el=hdferr) + 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 +#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') elseif(m == 'a') then - call h5fopen_f(fileName,H5F_ACC_RDWR_F,HDF5_openFile,hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_openFile: h5fopen_f (a)',el=hdferr) + 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) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_openFile: h5fopen_f (r)',el=hdferr) + 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',el=hdferr) + call IO_error(1_pInt,ext_msg='HDF5_openFile: h5fopen_f unknown access mode') 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 + !-------------------------------------------------------------------------------------------------- !> @brief close the opened HDF5 output file !-------------------------------------------------------------------------------------------------- subroutine HDF5_closeFile(fileHandle) implicit none - integer :: hdferr integer(HID_T), intent(in) :: fileHandle + + integer(HDF5_ERR_TYPE) :: hdferr + call h5fclose_f(fileHandle,hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_closeFile: h5fclose_f',el=hdferr) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_closeFile: h5fclose_f') end subroutine HDF5_closeFile + !-------------------------------------------------------------------------------------------------- -!> @brief adds a new group to the results file +!> @brief adds a new group to the fileHandle !-------------------------------------------------------------------------------------------------- -integer(HID_T) function HDF5_addGroup(groupName) - use hdf5 +integer(HID_T) function HDF5_addGroup(fileHandle,groupName) implicit none + integer(HID_T), intent(in) :: fileHandle character(len=*), intent(in) :: groupName - integer :: hdferr - call h5gcreate_f(resultsFile, trim(groupName), HDF5_addGroup, hdferr) + integer(HDF5_ERR_TYPE) :: 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)//')') + + !------------------------------------------------------------------------------------------------- + ! setting I/O mode to collective + 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)//')') + + !------------------------------------------------------------------------------------------------- + ! 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)//')') end function HDF5_addGroup !-------------------------------------------------------------------------------------------------- -!> @brief adds a new group to the fileHandle (additional to addGroup2) +!> @brief open an existing group of a file !-------------------------------------------------------------------------------------------------- -integer(HID_T) function HDF5_addGroup2(fileHandle,groupName) - use hdf5 +integer(HID_T) function HDF5_openGroup(fileHandle,groupName) implicit none + integer(HID_T), intent(in) :: fileHandle character(len=*), intent(in) :: groupName - integer(HID_T), intent(in) :: fileHandle - integer :: hdferr - - call h5gcreate_f(fileHandle, trim(groupName), HDF5_addGroup2, hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg = 'HDF5_addGroup2: h5gcreate_f ('//trim(groupName)//')') - -end function HDF5_addGroup2 -!-------------------------------------------------------------------------------------------------- -!> @brief open a group from the results file -!-------------------------------------------------------------------------------------------------- -integer(HID_T) function HDF5_openGroup(groupName) - use hdf5 + integer(HDF5_ERR_TYPE) :: 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)//')') - implicit none - character(len=*), intent(in) :: groupName - integer :: hdferr - - call h5gopen_f(resultsFile, trim(groupName), HDF5_openGroup, hdferr) + !------------------------------------------------------------------------------------------------- + ! setting I/O mode to collective + 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)//')') + + !------------------------------------------------------------------------------------------------- + ! 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)//')') end function HDF5_openGroup -!-------------------------------------------------------------------------------------------------- -!> @brief open an existing group of a file -!-------------------------------------------------------------------------------------------------- -integer(HID_T) function HDF5_openGroup2(FileReadID,groupName) - use hdf5 - - implicit none - character(len=*), intent(in) :: groupName - integer :: hdferr - integer(HID_T), intent(in) :: FileReadID - - call h5gopen_f(FileReadID, trim(groupName), HDF5_openGroup2, hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg = 'HDF5_openGroup2: h5gopen_f ('//trim(groupName)//')') - -end function HDF5_openGroup2 - -!-------------------------------------------------------------------------------------------------- -!> @brief set link to object in results file -!-------------------------------------------------------------------------------------------------- -subroutine HDF5_setLink(path,link) - use hdf5 - - implicit none - character(len=*), intent(in) :: path, link - integer :: hdferr - logical :: linkExists - - call h5lexists_f(resultsFile, link,linkExists, hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg = 'HDF5_setLink: h5lexists_soft_f ('//trim(link)//')') - if (linkExists) then - call h5ldelete_f(resultsFile,link, hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg = 'HDF5_setLink: h5ldelete_soft_f ('//trim(link)//')') - endif - call h5lcreate_soft_f(path, resultsFile, link, hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg = 'HDF5_setLink: h5lcreate_soft_f ('//trim(path)//' '//trim(link)//')') - -end subroutine HDF5_setLink - -!-------------------------------------------------------------------------------------------------- -!> @brief remove link to an object -!-------------------------------------------------------------------------------------------------- -subroutine HDF5_removeLink(link) - use hdf5 - - implicit none - character(len=*), intent(in) :: link - integer :: hdferr - - call h5ldelete_f(resultsFile,link, hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg = 'HDF5_removeLink: h5ldelete_soft_f ('//trim(link)//')') - -end subroutine HDF5_removeLink !-------------------------------------------------------------------------------------------------- !> @brief close a group !-------------------------------------------------------------------------------------------------- -subroutine HDF5_closeGroup(ID) - use hdf5 +subroutine HDF5_closeGroup(group_id) implicit none - integer(HID_T), intent(in) :: ID - integer :: hdferr + integer(HID_T), intent(in) :: group_id + integer(HDF5_ERR_TYPE) :: hdferr - call h5gclose_f(ID, hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg = 'HDF5_closeGroup: h5gclose_f (el is ID)', el = int(ID,pInt)) + 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 + !-------------------------------------------------------------------------------------------------- -!> @brief adds a StringAttribute to the results file +!> @brief check whether a group or a dataset exists !-------------------------------------------------------------------------------------------------- -subroutine HDF5_addStringAttribute(entity,attrLabel,attrValue) - use hdf5 +logical function HDF5_objectExists(loc_id,path) implicit none - integer(HID_T), intent(in) :: entity + integer(HID_T), intent(in) :: loc_id + character(len=*), intent(in), optional :: path + integer(HDF5_ERR_TYPE) :: hdferr + character(len=256) :: p + + 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 + +end function HDF5_objectExists + + +!-------------------------------------------------------------------------------------------------- +!> @brief adds a string attribute to the path given relative to the location +!-------------------------------------------------------------------------------------------------- +subroutine HDF5_addAttribute_str(loc_id,attrLabel,attrValue,path) + + implicit none + integer(HID_T), intent(in) :: loc_id character(len=*), intent(in) :: attrLabel, attrValue - integer :: hdferr + character(len=*), intent(in), optional :: path + integer(HDF5_ERR_TYPE) :: 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 call h5screate_f(H5S_SCALAR_F,space_id,hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_addStringAttribute: h5screate_f') + 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='IO_addStringAttribute: h5tcopy_f') + 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='IO_addStringAttribute: h5tset_size_f') - call h5acreate_f(entity, trim(attrLabel),type_id,space_id,attr_id,hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_addStringAttribute: h5acreate_f') + 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='IO_addStringAttribute: h5awrite_f') + 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='IO_addStringAttribute: h5aclose_f') + 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='IO_addStringAttribute: h5tclose_f') + 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='IO_addStringAttribute: h5sclose_f') + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_addAttribute_str: h5sclose_f') + +end subroutine HDF5_addAttribute_str -end subroutine HDF5_addStringAttribute !-------------------------------------------------------------------------------------------------- -!> @brief adds the unique mapping from spatial position and constituent ID to results +!> @brief adds a integer attribute to the path given relative to the location !-------------------------------------------------------------------------------------------------- -subroutine HDF5_mappingPhase(mapping,mapping2,Nconstituents,material_phase,phase_name,dataspace_size,mpiOffset,mpiOffset_phase) +subroutine HDF5_addAttribute_pInt(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(HDF5_ERR_TYPE) :: 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 + + 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 + + +!-------------------------------------------------------------------------------------------------- +!> @brief adds a integer attribute to the path given relative to the location +!-------------------------------------------------------------------------------------------------- +subroutine HDF5_addAttribute_pReal(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(HDF5_ERR_TYPE) :: 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 + + 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 + + +!-------------------------------------------------------------------------------------------------- +!> @brief set link to object in results file +!-------------------------------------------------------------------------------------------------- +subroutine HDF5_setLink(loc_id,target_name,link_name) use hdf5 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=*), intent(in) :: target_name, link_name + integer(HID_T), intent(in) :: loc_id + integer(HDF5_ERR_TYPE) :: hdferr + logical :: linkExists - character(len=len(phase_name(1))), dimension(:), allocatable :: namesNA - character(len=len(phase_name(1))) :: a - character(len=*), parameter :: n = "NULL" + 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(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 = HDF5_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 +end subroutine HDF5_setLink !-------------------------------------------------------------------------------------------------- -! 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') - +!> @brief subroutine for reading dataset of type pReal with 1 dimensions !-------------------------------------------------------------------------------------------------- -! 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') +subroutine HDF5_read_pReal1(loc_id,dataset,datasetName,parallel) + use numerics, only: & + worldrank, & + worldsize - 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') + 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 -!-------------------------------------------------------------------------------------------------- -! 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') + integer(pInt), dimension(:), allocatable :: & + globalShape, & !< shape of the dataset (all processes) + localShape, & !< shape of the dataset (this process) + readSize !< contribution of all processes + integer :: ierr + integer(HDF5_ERR_TYPE) :: hdferr + integer(HID_T) :: dset_id, filespace_id, memspace_id, plist_id, aplist_id + integer(HSIZE_T), dimension(1) :: myStart - 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') +!------------------------------------------------------------------------------------------------- +! determine shape of dataset + localShape = shape(dataset) + if (any(localShape(1:size(localShape)) == 0)) return -!-------------------------------------------------------------------------------------------------- -! 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 -#ifdef PETSc +!------------------------------------------------------------------------------------------------- +! creating a property list for transfer properties 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') -#endif +!------------------------------------------------------------------------------------------------- +! creating a property list for data access properties + call h5pcreate_f(H5P_DATASET_ACCESS_F, aplist_id, hdferr) !-------------------------------------------------------------------------------------------------- -! 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') - - 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') - -!-------------------------------------------------------------------------------------------------- -! 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 - -!-------------------------------------------------------------------------------------------------- -!> @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 - - 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 - - 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 = HDF5_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 + allocate(readSize(worldsize), source = 0_pInt) + readSize(worldrank+1) = localShape(1) #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') + if (present(parallel)) then; 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='IO_backwardMappingPhase: h5pset_dxpl_mpio_f') + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_read_pReal1: h5pset_dxpl_mpio_f') + call MPI_allreduce(MPI_IN_PLACE,readSize,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='HDF5_read_pReal1: MPI_allreduce') + endif; endif #endif + myStart = int([sum(readSize(1:worldrank))],HSIZE_T) + globalShape = [localShape(1:0),sum(readSize)] + !-------------------------------------------------------------------------------------------------- - ! 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') +! create dataspace in memory (local shape) + call h5screate_simple_f(size(localShape), int(localShape,HSIZE_T), memspace_id, hdferr, & + int(localShape,HSIZE_T)) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_read_pReal{}: h5screate_simple_f/memspace_id') +!-------------------------------------------------------------------------------------------------- +! set I/O mode for read operations to collective + call h5pset_all_coll_metadata_ops_f(aplist_id, .true., hdferr) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_read_pReal1: h5pset_all_coll_metadata_ops_f') +!-------------------------------------------------------------------------------------------------- +! open the dataset in the file + call h5dopen_f(loc_id,datasetName,dset_id,hdferr, dapl_id = aplist_id) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_read_pReal1: h5dopen_f') !-------------------------------------------------------------------------------------------------- - !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 +! get the space_id of dataset in the file + call h5dget_space_f(dset_id, filespace_id, hdferr) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_read_pReal1: h5dget_space_f') !-------------------------------------------------------------------------------------------------- -!> @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 - - 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 = HDF5_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 +! select a hyperslab (the portion of the current process) in the file + call h5sselect_hyperslab_f(filespace_id, H5S_SELECT_SET_F, myStart, int(localShape,HSIZE_T), hdferr) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_read_pReal1: h5sselect_hyperslab_f') !-------------------------------------------------------------------------------------------------- -! 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(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 -#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') -#endif - -!-------------------------------------------------------------------------------------------------- -! 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') - - 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') +! read + call h5dread_f(dset_id, H5T_NATIVE_DOUBLE,dataset,int(globalShape,HSIZE_T), 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') !-------------------------------------------------------------------------------------------------- !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) + call h5pclose_f(plist_id, hdferr) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_read_pReal1: plist_id') + call h5dclose_f(dset_id, hdferr) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_read_pReal1: h5dclose_f') + call h5sclose_f(filespace_id, hdferr) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_write_pReal1: h5sclose_f/filespace_id') + call h5sclose_f(memspace_id, hdferr) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_write_pReal1: h5sclose_f/memspace_id') + +end subroutine HDF5_read_pReal1 -end subroutine HDF5_mappingHomog !-------------------------------------------------------------------------------------------------- -!> @brief adds the backward mapping from spatial position and constituent ID to results +!> @brief subroutine for reading dataset of type pReal with 2 dimensions !-------------------------------------------------------------------------------------------------- -subroutine HDF5_backwardMappingHomog(material_homog,homogmemberat,homogenization_name,dataspace_size,mpiOffset,mpiOffset_homog) - use hdf5 +subroutine HDF5_read_pReal2(loc_id,dataset,datasetName,parallel) + use numerics, only: & + worldrank, & + worldsize 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 + 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(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(pInt), dimension(:), allocatable :: & + globalShape, & !< shape of the dataset (all processes) + localShape, & !< shape of the dataset (this process) + readSize !< contribution of all processes + integer :: ierr + integer(HDF5_ERR_TYPE) :: hdferr + integer(HID_T) :: dset_id, filespace_id, memspace_id, plist_id, aplist_id + integer(HSIZE_T), dimension(2) :: myStart - integer(HSIZE_T), dimension(1) :: counter - integer(HSSIZE_T), dimension(1) :: fileOffset +!------------------------------------------------------------------------------------------------- +! determine shape of dataset + localShape = shape(dataset) + if (any(localShape(1:size(localShape)) == 0)) return - 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 = HDF5_openGroup('/current/homogenization/'//trim(homogID)//'_'//homogenization_name(i)) +!------------------------------------------------------------------------------------------------- +! creating a property list for transfer properties + call h5pcreate_f(H5P_DATASET_XFER_F, plist_id, hdferr) +!------------------------------------------------------------------------------------------------- +! creating a property list for data access properties + call h5pcreate_f(H5P_DATASET_ACCESS_F, aplist_id, hdferr) !-------------------------------------------------------------------------------------------------- - ! 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 + allocate(readSize(worldsize), source = 0_pInt) + readSize(worldrank+1) = localShape(2) #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') + if (present(parallel)) then; 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='IO_backwardMappingHomog: h5pset_dxpl_mpio_f') + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_read_pReal2: h5pset_dxpl_mpio_f') + call MPI_allreduce(MPI_IN_PLACE,readSize,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='HDF5_read_pReal2: MPI_allreduce') + endif; endif #endif + myStart = int([0,sum(readSize(1:worldrank))],HSIZE_T) + globalShape = [localShape(1:1),sum(readSize)] + !-------------------------------------------------------------------------------------------------- - ! 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') +! create dataspace in memory (local shape) + call h5screate_simple_f(size(localShape), int(localShape,HSIZE_T), memspace_id, hdferr, & + int(localShape,HSIZE_T)) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_read_pReal{}: h5screate_simple_f/memspace_id') +!-------------------------------------------------------------------------------------------------- +! set I/O mode for read operations to collective + call h5pset_all_coll_metadata_ops_f(aplist_id, .true., hdferr) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_read_pReal2: h5pset_all_coll_metadata_ops_f') +!-------------------------------------------------------------------------------------------------- +! open the dataset in the file + call h5dopen_f(loc_id,datasetName,dset_id,hdferr, dapl_id = aplist_id) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_read_pReal2: h5dopen_f') !-------------------------------------------------------------------------------------------------- - !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 +! get the space_id of dataset in the file + call h5dget_space_f(dset_id, filespace_id, hdferr) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_read_pReal2: h5dget_space_f') !-------------------------------------------------------------------------------------------------- -!> @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 = HDF5_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 +! select a hyperslab (the portion of the current process) in the file + call h5sselect_hyperslab_f(filespace_id, H5S_SELECT_SET_F, myStart, int(localShape,HSIZE_T), hdferr) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_read_pReal2: h5sselect_hyperslab_f') !-------------------------------------------------------------------------------------------------- -! 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 +! read + call h5dread_f(dset_id, H5T_NATIVE_DOUBLE,dataset,int(globalShape,HSIZE_T), 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') !-------------------------------------------------------------------------------------------------- !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) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_read_pReal2: plist_id') + call h5dclose_f(dset_id, hdferr) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_read_pReal2: h5dclose_f') + call h5sclose_f(filespace_id, hdferr) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_write_pReal2: h5sclose_f/filespace_id') + call h5sclose_f(memspace_id, hdferr) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_write_pReal2: h5sclose_f/memspace_id') -end subroutine HDF5_mappingCrystallite +end subroutine HDF5_read_pReal2 !-------------------------------------------------------------------------------------------------- -!> @brief adds the backward mapping from spatial position and constituent ID to results +!> @brief subroutine for reading dataset of type pReal with 3 dimensions !-------------------------------------------------------------------------------------------------- -subroutine HDF5_backwardMappingCrystallite(crystalliteAt,crystmemberAt,crystallite_name,dataspace_size,mpiOffset,mpiOffset_cryst) - use hdf5 +subroutine HDF5_read_pReal3(loc_id,dataset,datasetName,parallel) + use numerics, only: & + worldrank, & + worldsize 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 + 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(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(pInt), dimension(:), allocatable :: & + globalShape, & !< shape of the dataset (all processes) + localShape, & !< shape of the dataset (this process) + readSize !< contribution of all processes + integer :: ierr + integer(HDF5_ERR_TYPE) :: hdferr + integer(HID_T) :: dset_id, filespace_id, memspace_id, plist_id, aplist_id + integer(HSIZE_T), dimension(3) :: myStart - integer(HSIZE_T), dimension(1) :: counter - integer(HSSIZE_T), dimension(1) :: fileOffset +!------------------------------------------------------------------------------------------------- +! determine shape of dataset + localShape = shape(dataset) + if (any(localShape(1:size(localShape)) == 0)) return - 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 = HDF5_openGroup('/current/crystallite/'//trim(crystallID)//'_'//crystallite_name(i)) - NmatPoints = count(crystalliteAt == i) +!------------------------------------------------------------------------------------------------- +! creating a property list for transfer properties + call h5pcreate_f(H5P_DATASET_XFER_F, plist_id, hdferr) +!------------------------------------------------------------------------------------------------- +! creating a property list for data access properties + call h5pcreate_f(H5P_DATASET_ACCESS_F, aplist_id, hdferr) !-------------------------------------------------------------------------------------------------- - ! 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 + allocate(readSize(worldsize), source = 0_pInt) + readSize(worldrank+1) = localShape(3) #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') + if (present(parallel)) then; 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='IO_backwardMappingCrystallite: h5pset_dxpl_mpio_f') + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_read_pReal3: h5pset_dxpl_mpio_f') + call MPI_allreduce(MPI_IN_PLACE,readSize,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='HDF5_read_pReal3: MPI_allreduce') + endif; endif #endif + myStart = int([0,0,sum(readSize(1:worldrank))],HSIZE_T) + globalShape = [localShape(1:2),sum(readSize)] + !-------------------------------------------------------------------------------------------------- - ! 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') +! create dataspace in memory (local shape) + call h5screate_simple_f(size(localShape), int(localShape,HSIZE_T), memspace_id, hdferr, & + int(localShape,HSIZE_T)) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_read_pReal{}: h5screate_simple_f/memspace_id') +!-------------------------------------------------------------------------------------------------- +! set I/O mode for read operations to collective + call h5pset_all_coll_metadata_ops_f(aplist_id, .true., hdferr) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_read_pReal3: h5pset_all_coll_metadata_ops_f') +!-------------------------------------------------------------------------------------------------- +! open the dataset in the file + call h5dopen_f(loc_id,datasetName,dset_id,hdferr, dapl_id = aplist_id) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_read_pReal3: h5dopen_f') !-------------------------------------------------------------------------------------------------- - !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 +! get the space_id of dataset in the file + call h5dget_space_f(dset_id, filespace_id, hdferr) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_read_pReal3: h5dget_space_f') !-------------------------------------------------------------------------------------------------- -!> @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 = HDF5_openGroup("mapping") +! select a hyperslab (the portion of the current process) in the file + call h5sselect_hyperslab_f(filespace_id, H5S_SELECT_SET_F, myStart, int(localShape,HSIZE_T), hdferr) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_read_pReal3: h5sselect_hyperslab_f') !-------------------------------------------------------------------------------------------------- -! 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') +! read + call h5dread_f(dset_id, H5T_NATIVE_DOUBLE,dataset,int(globalShape,HSIZE_T), 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') !-------------------------------------------------------------------------------------------------- !close types, dataspaces + call h5pclose_f(plist_id, hdferr) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_read_pReal3: plist_id') 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) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_read_pReal3: h5dclose_f') + call h5sclose_f(filespace_id, hdferr) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_write_pReal3: h5sclose_f/filespace_id') + call h5sclose_f(memspace_id, hdferr) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_write_pReal3: h5sclose_f/memspace_id') + +end subroutine HDF5_read_pReal3 -end subroutine HDF5_mappingCells !-------------------------------------------------------------------------------------------------- -!> @brief creates a new 3D Tensor dataset in the given group location !!!TODO: really necessary? +!> @brief subroutine for reading dataset of type pReal with 4 dimensions !-------------------------------------------------------------------------------------------------- -subroutine HDF5_addTensor3DDataset(group,Nnodes,tensorSize,label,SIunit) - use hdf5 +subroutine HDF5_read_pReal4(loc_id,dataset,datasetName,parallel) + use numerics, only: & + worldrank, & + worldsize implicit none - integer(HID_T), intent(in) :: group - integer(pInt), intent(in) :: Nnodes, tensorSize - character(len=*), intent(in) :: SIunit, label + 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) :: space_id, dset_id - integer(HSIZE_T), dimension(3) :: dataShape - dataShape = int([tensorSize,tensorSize,Nnodes], HSIZE_T) + integer(pInt), dimension(:), allocatable :: & + globalShape, & !< shape of the dataset (all processes) + localShape, & !< shape of the dataset (this process) + readSize !< contribution of all processes + integer :: ierr + integer(HDF5_ERR_TYPE) :: hdferr + integer(HID_T) :: dset_id, filespace_id, memspace_id, plist_id, aplist_id + integer(HSIZE_T), dimension(4) :: myStart + +!------------------------------------------------------------------------------------------------- +! determine shape of dataset + localShape = shape(dataset) + if (any(localShape(1:size(localShape)) == 0)) return + +!------------------------------------------------------------------------------------------------- +! creating a property list for transfer properties + call h5pcreate_f(H5P_DATASET_XFER_F, plist_id, hdferr) + +!------------------------------------------------------------------------------------------------- +! creating a property list for data access properties + call h5pcreate_f(H5P_DATASET_ACCESS_F, aplist_id, hdferr) +!-------------------------------------------------------------------------------------------------- + allocate(readSize(worldsize), source = 0_pInt) + readSize(worldrank+1) = localShape(4) +#ifdef PETSc + if (present(parallel)) then; 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='HDF5_read_pReal4: h5pset_dxpl_mpio_f') + call MPI_allreduce(MPI_IN_PLACE,readSize,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='HDF5_read_pReal4: MPI_allreduce') + endif; endif +#endif + myStart = int([0,0,0,sum(readSize(1:worldrank))],HSIZE_T) + globalShape = [localShape(1:3),sum(readSize)] + !-------------------------------------------------------------------------------------------------- -! create dataspace - call h5screate_simple_f(3, dataShape, space_id, hdferr, dataShape) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_addTensor3DDataset: h5screate_simple_f') +! create dataspace in memory (local shape) + call h5screate_simple_f(size(localShape), int(localShape,HSIZE_T), memspace_id, hdferr, & + int(localShape,HSIZE_T)) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_read_pReal{}: h5screate_simple_f/memspace_id') +!-------------------------------------------------------------------------------------------------- +! set I/O mode for read operations to collective + call h5pset_all_coll_metadata_ops_f(aplist_id, .true., hdferr) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_read_pReal4: h5pset_all_coll_metadata_ops_f') +!-------------------------------------------------------------------------------------------------- +! open the dataset in the file + call h5dopen_f(loc_id,datasetName,dset_id,hdferr, dapl_id = aplist_id) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_read_pReal4: h5dopen_f') !-------------------------------------------------------------------------------------------------- -! create Dataset - call h5dcreate_f(group, trim(label),H5T_NATIVE_DOUBLE, space_id, dset_id, hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_addTensor3DDataset: h5dcreate_f') - call HDF5_addStringAttribute(dset_id,'unit',trim(SIunit)) +! get the space_id of dataset in the file + call h5dget_space_f(dset_id, filespace_id, hdferr) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_read_pReal4: h5dget_space_f') + +!-------------------------------------------------------------------------------------------------- +! select a hyperslab (the portion of the current process) in the file + call h5sselect_hyperslab_f(filespace_id, H5S_SELECT_SET_F, myStart, int(localShape,HSIZE_T), hdferr) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_read_pReal4: h5sselect_hyperslab_f') + +!-------------------------------------------------------------------------------------------------- +! read + call h5dread_f(dset_id, H5T_NATIVE_DOUBLE,dataset,int(globalShape,HSIZE_T), 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') !-------------------------------------------------------------------------------------------------- !close types, dataspaces - call h5dclose_f(dset_id, hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_addTensor3DDataset: h5dclose_f') - call h5sclose_f(space_id, hdferr) - -end subroutine HDF5_addTensor3DDataset - -!-------------------------------------------------------------------------------------------------- -!> @brief creates a new vector dataset in the given group location !!!TODO: really necessary? -!-------------------------------------------------------------------------------------------------- -subroutine HDF5_writeVectorDataset(group,dataset,label,SIunit,dataspace_size,mpiOffset) - use hdf5 - - implicit none - integer(HID_T), intent(in) :: group - character(len=*), intent(in) :: SIunit,label - integer(pInt), intent(in) :: dataspace_size, mpiOffset - real(pReal), intent(in), dimension(:,:) :: dataset - - integer :: hdferr, vectorSize - integer(HID_T) :: dset_id, space_id, memspace, plist_id - - integer(HSIZE_T), dimension(2) :: counter - integer(HSSIZE_T), dimension(2) :: fileOffset - - if(any(shape(dataset) == 0)) return - - vectorSize = size(dataset,1) - - call HDF5_addVectorDataset(group,dataspace_size,vectorSize,label,SIunit) ! here nNodes need to be global - call h5dopen_f(group, label, dset_id, hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_writeVectorDataset: h5dopen_f') - - ! Define and select hyperslabs - counter(1) = vectorSize ! how big i am - counter(2) = size(dataset,2) - 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='HDF5_writeVectorDataset: h5screate_simple_f') - call h5dget_space_f(dset_id, space_id, hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_writeVectorDataset: 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='HDF5_writeVectorDataset: 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='HDF5_writeVectorDataset: h5pcreate_f') - call h5pset_dxpl_mpio_f(plist_id, H5FD_MPIO_COLLECTIVE_F, hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_writeVectorDataset: h5pset_dxpl_mpio_f') -#endif - - ! Write the dataset collectively - call h5dwrite_f(dset_id, H5T_NATIVE_DOUBLE, dataset, int([vectorSize, 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='HDF5_writeVectorDataset: h5dwrite_f') - - call h5sclose_f(space_id, hdferr) - call h5sclose_f(memspace, hdferr) - call h5dclose_f(dset_id, hdferr) call h5pclose_f(plist_id, hdferr) - -end subroutine HDF5_writeVectorDataset - -!-------------------------------------------------------------------------------------------------- -!> @brief creates a new tensor dataset in the given group location -! by default, a 3x3 tensor is assumed !!!TODO: really necessary? -!-------------------------------------------------------------------------------------------------- -subroutine HDF5_writeTensorDataset(group,dataset,label,SIunit,dataspace_size,mpiOffset) - use hdf5 - - implicit none - integer(HID_T), intent(in) :: group - character(len=*), intent(in) :: SIunit,label - integer(pInt), intent(in) :: dataspace_size, mpiOffset - real(pReal), intent(in), dimension(:,:,:) :: dataset - - integer :: hdferr, tensorSize - integer(HID_T) :: dset_id, space_id, memspace, plist_id - - integer(HSIZE_T), dimension(3) :: counter - integer(HSSIZE_T), dimension(3) :: fileOffset - - if(any(shape(dataset) == 0)) return - - tensorSize = size(dataset,1) - - call HDF5_addTensor3DDataset(group,dataspace_size,tensorSize,label,SIunit) - call h5dopen_f(group, label, dset_id, hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_writeTensorDataset: h5dopen_f') - - ! Define and select hyperslabs - counter(1) = tensorSize ! how big i am - counter(2) = tensorSize - counter(3) = size(dataset,3) - fileOffset(1) = 0 ! where i start to write my data - fileOffset(2) = 0 - fileOffset(3) = mpiOffset - - call h5screate_simple_f(3, counter, memspace, hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_writeTensorDataset: h5screate_simple_f') - call h5dget_space_f(dset_id, space_id, hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_writeTensorDataset: 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='HDF5_writeTensorDataset: 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='HDF5_writeTensorDataset: h5pcreate_f') - call h5pset_dxpl_mpio_f(plist_id, H5FD_MPIO_COLLECTIVE_F, hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_writeTensorDataset: h5pset_dxpl_mpio_f') -#endif - - ! Write the dataset collectively - call h5dwrite_f(dset_id, H5T_NATIVE_DOUBLE, dataset, int([tensorSize, 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='HDF5_writeTensorDataset: h5dwrite_f') - - call h5sclose_f(space_id, hdferr) - call h5sclose_f(memspace, hdferr) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_read_pReal4: plist_id') call h5dclose_f(dset_id, hdferr) - call h5pclose_f(plist_id, hdferr) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_read_pReal4: h5dclose_f') + call h5sclose_f(filespace_id, hdferr) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_write_pReal4: h5sclose_f/filespace_id') + call h5sclose_f(memspace_id, hdferr) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_write_pReal4: h5sclose_f/memspace_id') + +end subroutine HDF5_read_pReal4 - end subroutine HDF5_writeTensorDataset !-------------------------------------------------------------------------------------------------- -!> @brief adds a new vector dataset to the given group location +!> @brief subroutine for reading dataset of type pReal with 5 dimensions !-------------------------------------------------------------------------------------------------- -subroutine HDF5_addVectorDataset(group,nnodes,vectorSize,label,SIunit) - use hdf5 +subroutine HDF5_read_pReal5(loc_id,dataset,datasetName,parallel) + use numerics, only: & + worldrank, & + worldsize implicit none - integer(HID_T), intent(in) :: group - integer(pInt), intent(in) :: nnodes,vectorSize - character(len=*), intent(in) :: SIunit,label + 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(pInt), dimension(:), allocatable :: & + globalShape, & !< shape of the dataset (all processes) + localShape, & !< shape of the dataset (this process) + readSize !< contribution of all processes + integer :: ierr + integer(HDF5_ERR_TYPE) :: hdferr + integer(HID_T) :: dset_id, filespace_id, memspace_id, plist_id, aplist_id + integer(HSIZE_T), dimension(5) :: myStart + +!------------------------------------------------------------------------------------------------- +! determine shape of dataset + localShape = shape(dataset) + if (any(localShape(1:size(localShape)) == 0)) return + +!------------------------------------------------------------------------------------------------- +! creating a property list for transfer properties + call h5pcreate_f(H5P_DATASET_XFER_F, plist_id, hdferr) + +!------------------------------------------------------------------------------------------------- +! creating a property list for data access properties + call h5pcreate_f(H5P_DATASET_ACCESS_F, aplist_id, hdferr) +!-------------------------------------------------------------------------------------------------- + allocate(readSize(worldsize), source = 0_pInt) + readSize(worldrank+1) = localShape(5) +#ifdef PETSc + if (present(parallel)) then; 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='HDF5_read_pReal5: h5pset_dxpl_mpio_f') + call MPI_allreduce(MPI_IN_PLACE,readSize,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='HDF5_read_pReal5: MPI_allreduce') + endif; endif +#endif + myStart = int([0,0,0,0,sum(readSize(1:worldrank))],HSIZE_T) + globalShape = [localShape(1:4),sum(readSize)] - integer :: hdferr - integer(HID_T) :: space_id, dset_id !-------------------------------------------------------------------------------------------------- -! create dataspace - call h5screate_simple_f(2, int([vectorSize,Nnodes],HSIZE_T), space_id, hdferr, & - int([vectorSize,Nnodes],HSIZE_T)) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_addVectorDataset: h5screate_simple_f') +! create dataspace in memory (local shape) + call h5screate_simple_f(size(localShape), int(localShape,HSIZE_T), memspace_id, hdferr, & + int(localShape,HSIZE_T)) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_read_pReal{}: h5screate_simple_f/memspace_id') +!-------------------------------------------------------------------------------------------------- +! set I/O mode for read operations to collective + call h5pset_all_coll_metadata_ops_f(aplist_id, .true., hdferr) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_read_pReal5: h5pset_all_coll_metadata_ops_f') +!-------------------------------------------------------------------------------------------------- +! open the dataset in the file + call h5dopen_f(loc_id,datasetName,dset_id,hdferr, dapl_id = aplist_id) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_read_pReal5: h5dopen_f') !-------------------------------------------------------------------------------------------------- -! create Dataset - call h5dcreate_f(group, trim(label), H5T_NATIVE_DOUBLE, space_id, dset_id, hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_addVectorDataset: h5dcreate_f') - call HDF5_addStringAttribute(dset_id,'unit',trim(SIunit)) +! get the space_id of dataset in the file + call h5dget_space_f(dset_id, filespace_id, hdferr) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_read_pReal5: h5dget_space_f') + +!-------------------------------------------------------------------------------------------------- +! select a hyperslab (the portion of the current process) in the file + call h5sselect_hyperslab_f(filespace_id, H5S_SELECT_SET_F, myStart, int(localShape,HSIZE_T), hdferr) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_read_pReal5: h5sselect_hyperslab_f') + +!-------------------------------------------------------------------------------------------------- +! read + call h5dread_f(dset_id, H5T_NATIVE_DOUBLE,dataset,int(globalShape,HSIZE_T), 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') !-------------------------------------------------------------------------------------------------- !close types, dataspaces + call h5pclose_f(plist_id, hdferr) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_read_pReal5: plist_id') call h5dclose_f(dset_id, hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_addVectorDataset: h5dclose_f') - call h5sclose_f(space_id, hdferr) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_read_pReal5: h5dclose_f') + call h5sclose_f(filespace_id, hdferr) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_write_pReal5: h5sclose_f/filespace_id') + call h5sclose_f(memspace_id, hdferr) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_write_pReal5: h5sclose_f/memspace_id') + +end subroutine HDF5_read_pReal5 -end subroutine HDF5_addVectorDataset !-------------------------------------------------------------------------------------------------- -!> @brief writes to a new scalar dataset in the given group location +!> @brief subroutine for reading dataset of type pReal with 6 dimensions !-------------------------------------------------------------------------------------------------- -subroutine HDF5_writeScalarDataset(group,dataset,label,SIunit,dataspace_size,mpiOffset) - use hdf5 +subroutine HDF5_read_pReal6(loc_id,dataset,datasetName,parallel) + use numerics, only: & + worldrank, & + worldsize implicit none - integer(HID_T), intent(in) :: group - character(len=*), intent(in) :: SIunit,label - integer(pInt), intent(in) :: dataspace_size, mpiOffset - real(pReal), intent(in), dimension(:) :: dataset + 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, nNodes - integer(HID_T) :: dset_id, space_id, memspace, plist_id - integer(HSIZE_T), dimension(1) :: counter - integer(HSIZE_T), dimension(1) :: fileOffset + integer(pInt), dimension(:), allocatable :: & + globalShape, & !< shape of the dataset (all processes) + localShape, & !< shape of the dataset (this process) + readSize !< contribution of all processes + integer :: ierr + integer(HDF5_ERR_TYPE) :: hdferr + integer(HID_T) :: dset_id, filespace_id, memspace_id, plist_id, aplist_id + integer(HSIZE_T), dimension(6) :: myStart - nNodes = size(dataset) - if (nNodes < 1) return +!------------------------------------------------------------------------------------------------- +! determine shape of dataset + localShape = shape(dataset) + if (any(localShape(1:size(localShape)) == 0)) return - call HDF5_addScalarDataset(group,dataspace_size,label,SIunit) - call h5dopen_f(group, label, dset_id, hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_writeScalarDataset: h5dopen_f') - - ! Define and select hyperslabs - counter = size(dataset) ! 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='HDF5_writeScalarDataset: h5screate_simple_f') - call h5dget_space_f(dset_id, space_id, hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_writeScalarDataset: 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='HDF5_writeScalarDataset: h5sselect_hyperslab_f') - - ! Create property list for collective dataset write -#ifdef PETSc +!------------------------------------------------------------------------------------------------- +! creating a property list for transfer properties call h5pcreate_f(H5P_DATASET_XFER_F, plist_id, hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_writeScalarDataset: h5pcreate_f') - call h5pset_dxpl_mpio_f(plist_id, H5FD_MPIO_COLLECTIVE_F, hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_writeScalarDataset: h5pset_dxpl_mpio_f') + +!------------------------------------------------------------------------------------------------- +! creating a property list for data access properties + call h5pcreate_f(H5P_DATASET_ACCESS_F, aplist_id, hdferr) +!-------------------------------------------------------------------------------------------------- + allocate(readSize(worldsize), source = 0_pInt) + readSize(worldrank+1) = localShape(6) +#ifdef PETSc + if (present(parallel)) then; 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='HDF5_read_pReal6: h5pset_dxpl_mpio_f') + call MPI_allreduce(MPI_IN_PLACE,readSize,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='HDF5_read_pReal6: MPI_allreduce') + endif; endif +#endif + myStart = int([0,0,0,0,0,sum(readSize(1:worldrank))],HSIZE_T) + globalShape = [localShape(1:5),sum(readSize)] + + +!-------------------------------------------------------------------------------------------------- +! create dataspace in memory (local shape) + call h5screate_simple_f(size(localShape), int(localShape,HSIZE_T), memspace_id, hdferr, & + int(localShape,HSIZE_T)) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_read_pReal{}: h5screate_simple_f/memspace_id') +!-------------------------------------------------------------------------------------------------- +! set I/O mode for read operations to collective + call h5pset_all_coll_metadata_ops_f(aplist_id, .true., hdferr) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_read_pReal6: h5pset_all_coll_metadata_ops_f') +!-------------------------------------------------------------------------------------------------- +! open the dataset in the file + call h5dopen_f(loc_id,datasetName,dset_id,hdferr, dapl_id = aplist_id) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_read_pReal6: h5dopen_f') + +!-------------------------------------------------------------------------------------------------- +! get the space_id of dataset in the file + call h5dget_space_f(dset_id, filespace_id, hdferr) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_read_pReal6: h5dget_space_f') + +!-------------------------------------------------------------------------------------------------- +! select a hyperslab (the portion of the current process) in the file + call h5sselect_hyperslab_f(filespace_id, H5S_SELECT_SET_F, myStart, int(localShape,HSIZE_T), hdferr) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_read_pReal6: h5sselect_hyperslab_f') + +!-------------------------------------------------------------------------------------------------- +! read + call h5dread_f(dset_id, H5T_NATIVE_DOUBLE,dataset,int(globalShape,HSIZE_T), 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') + +!-------------------------------------------------------------------------------------------------- +!close types, dataspaces + call h5pclose_f(plist_id, hdferr) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_read_pReal6: plist_id') + call h5dclose_f(dset_id, hdferr) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_read_pReal6: h5dclose_f') + call h5sclose_f(filespace_id, hdferr) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_write_pReal6: h5sclose_f/filespace_id') + call h5sclose_f(memspace_id, hdferr) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_write_pReal6: h5sclose_f/memspace_id') + +end subroutine HDF5_read_pReal6 + + +!-------------------------------------------------------------------------------------------------- +!> @brief subroutine for reading dataset of type pReal with 7 dimensions +!-------------------------------------------------------------------------------------------------- +subroutine HDF5_read_pReal7(loc_id,dataset,datasetName,parallel) + use numerics, only: & + worldrank, & + worldsize + + 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(pInt), dimension(:), allocatable :: & + globalShape, & !< shape of the dataset (all processes) + localShape, & !< shape of the dataset (this process) + readSize !< contribution of all processes + integer :: ierr + integer(HDF5_ERR_TYPE) :: hdferr + integer(HID_T) :: dset_id, filespace_id, memspace_id, plist_id, aplist_id + integer(HSIZE_T), dimension(7) :: myStart + +!------------------------------------------------------------------------------------------------- +! determine shape of dataset + localShape = shape(dataset) + if (any(localShape(1:size(localShape)) == 0)) return + +!------------------------------------------------------------------------------------------------- +! creating a property list for transfer properties + call h5pcreate_f(H5P_DATASET_XFER_F, plist_id, hdferr) + +!------------------------------------------------------------------------------------------------- +! creating a property list for data access properties + call h5pcreate_f(H5P_DATASET_ACCESS_F, aplist_id, hdferr) +!-------------------------------------------------------------------------------------------------- + allocate(readSize(worldsize), source = 0_pInt) + readSize(worldrank+1) = localShape(7) +#ifdef PETSc + if (present(parallel)) then; 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='HDF5_read_pReal7: h5pset_dxpl_mpio_f') + call MPI_allreduce(MPI_IN_PLACE,readSize,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='HDF5_read_pReal7: MPI_allreduce') + endif; endif +#endif + myStart = int([0,0,0,0,0,0,sum(readSize(1:worldrank))],HSIZE_T) + globalShape = [localShape(1:6),sum(readSize)] + + +!-------------------------------------------------------------------------------------------------- +! create dataspace in memory (local shape) + call h5screate_simple_f(size(localShape), int(localShape,HSIZE_T), memspace_id, hdferr, & + int(localShape,HSIZE_T)) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_read_pReal{}: h5screate_simple_f/memspace_id') +!-------------------------------------------------------------------------------------------------- +! set I/O mode for read operations to collective + call h5pset_all_coll_metadata_ops_f(aplist_id, .true., hdferr) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_read_pReal7: h5pset_all_coll_metadata_ops_f') +!-------------------------------------------------------------------------------------------------- +! open the dataset in the file + call h5dopen_f(loc_id,datasetName,dset_id,hdferr, dapl_id = aplist_id) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_read_pReal7: h5dopen_f') + +!-------------------------------------------------------------------------------------------------- +! get the space_id of dataset in the file + call h5dget_space_f(dset_id, filespace_id, hdferr) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_read_pReal7: h5dget_space_f') + +!-------------------------------------------------------------------------------------------------- +! select a hyperslab (the portion of the current process) in the file + call h5sselect_hyperslab_f(filespace_id, H5S_SELECT_SET_F, myStart, int(localShape,HSIZE_T), hdferr) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_read_pReal7: h5sselect_hyperslab_f') + +!-------------------------------------------------------------------------------------------------- +! read + call h5dread_f(dset_id, H5T_NATIVE_DOUBLE,dataset,int(globalShape,HSIZE_T), 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') + +!-------------------------------------------------------------------------------------------------- +!close types, dataspaces + call h5pclose_f(plist_id, hdferr) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_read_pReal7: plist_id') + call h5dclose_f(dset_id, hdferr) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_read_pReal7: h5dclose_f') + call h5sclose_f(filespace_id, hdferr) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_write_pReal7: h5sclose_f/filespace_id') + call h5sclose_f(memspace_id, hdferr) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_write_pReal7: h5sclose_f/memspace_id') + +end subroutine HDF5_read_pReal7 + + +!-------------------------------------------------------------------------------------------------- +!> @brief subroutine for reading dataset of type pInt with 1 dimensions +!-------------------------------------------------------------------------------------------------- +subroutine HDF5_read_pInt1(loc_id,dataset,datasetName,parallel) + use numerics, only: & + worldrank, & + worldsize + + 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(pInt), dimension(:), allocatable :: & + globalShape, & !< shape of the dataset (all processes) + localShape, & !< shape of the dataset (this process) + readSize !< contribution of all processes + integer :: ierr + integer(HDF5_ERR_TYPE) :: hdferr + integer(HID_T) :: dset_id, filespace_id, memspace_id, plist_id, aplist_id + integer(HSIZE_T), dimension(1) :: myStart + +!------------------------------------------------------------------------------------------------- +! determine shape of dataset + localShape = shape(dataset) + if (any(localShape(1:size(localShape)) == 0)) return + +!------------------------------------------------------------------------------------------------- +! creating a property list for transfer properties + call h5pcreate_f(H5P_DATASET_XFER_F, plist_id, hdferr) + +!-------------------------------------------------------------------------------------------------- +! creating a property list for data access properties + call h5pcreate_f(H5P_DATASET_ACCESS_F, aplist_id, hdferr) +!-------------------------------------------------------------------------------------------------- + allocate(readSize(worldsize), source = 0_pInt) + readSize(worldrank+1) = localShape(1) + +#ifdef PETSc + if (present(parallel)) then; 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='HDF5_read_pInt1: h5pset_dxpl_mpio_f') + call MPI_allreduce(MPI_IN_PLACE,readSize,worldsize,MPI_LONG,MPI_SUM,PETSC_COMM_WORLD,ierr) ! get total output size over each process + if (ierr /= 0) call IO_error(894_pInt,ext_msg='HDF5_read_pInt1: MPI_allreduce') + endif; endif #endif - ! Write the dataset collectively - call h5dwrite_f(dset_id, H5T_NATIVE_DOUBLE, dataset, 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='HDF5_writeScalarDataset: h5dwrite_f') + myStart = int([sum(readSize(1:worldrank))],HSIZE_T) + globalShape = [localShape(1:0),sum(readSize)] - call h5sclose_f(space_id, hdferr) - call h5sclose_f(memspace, hdferr) - call h5dclose_f(dset_id, hdferr) +!-------------------------------------------------------------------------------------------------- +! create dataspace in memory (local shape) + call h5screate_simple_f(size(localShape), int(localShape,HSIZE_T), memspace_id, hdferr, & + int(localShape,HSIZE_T)) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_read_pInt{}: h5screate_simple_f/memspace_id') +!-------------------------------------------------------------------------------------------------- +! set I/O mode for read operations to collective + call h5pset_all_coll_metadata_ops_f(aplist_id, .true., hdferr) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_read_pInt1: h5pset_all_coll_metadata_ops_f') + +!-------------------------------------------------------------------------------------------------- +! open the dataset in the file + call h5dopen_f(loc_id,datasetName,dset_id,hdferr,dapl_id = aplist_id) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_read_pInt1: h5dopen_f') + +!-------------------------------------------------------------------------------------------------- +! get the space_id of dataset in the file + call h5dget_space_f(dset_id, filespace_id, hdferr) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_read_pInt1: h5dget_space_f') + +!-------------------------------------------------------------------------------------------------- +! select a hyperslab (the portion of the current process) in the file + call h5sselect_hyperslab_f(filespace_id, H5S_SELECT_SET_F, myStart, int(localShape,HSIZE_T), hdferr) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_read_pInt1: h5sselect_hyperslab_f') + +!-------------------------------------------------------------------------------------------------- +! read + call h5dread_f(dset_id, H5T_NATIVE_INTEGER,dataset,int(globalShape,HSIZE_T), 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') + +!-------------------------------------------------------------------------------------------------- +!close types, dataspaces call h5pclose_f(plist_id, hdferr) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_read_pInt1: plist_id') + call h5dclose_f(dset_id, hdferr) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_read_pInt1: h5dclose_f') + call h5sclose_f(filespace_id, hdferr) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_read_pInt1: h5sclose_f/filespace_id') + call h5sclose_f(memspace_id, hdferr) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_read_pInt1: h5sclose_f/memspace_id') + +end subroutine HDF5_read_pInt1 -end subroutine HDF5_writeScalarDataset !-------------------------------------------------------------------------------------------------- -!> @brief subroutine for reading dataset of the type pReal with 1 dimension +!> @brief subroutine for reading dataset of type pInt with 2 dimensions !-------------------------------------------------------------------------------------------------- -subroutine HDF5_read_pReal_1(dataset,loc_id,datasetName) +subroutine HDF5_read_pInt2(loc_id,dataset,datasetName,parallel) + use numerics, only: & + worldrank, & + worldsize implicit none - real(pReal), intent(out), dimension(:) :: dataset + 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 - integer(pInt),dimension(:), allocatable :: myShape - integer :: hdferr - integer(HID_T) :: dset_id - myShape = shape(dataset) + logical, intent(in), optional :: parallel - call h5dopen_f(loc_id,datasetName,dset_id,hdferr) - if (hdferr /= 0) call IO_error(0_pInt,ext_msg='HDF5_read_pReal_shape1: h5dopen_f') - call h5dread_f(dset_id,H5T_NATIVE_DOUBLE,dataset,int(myShape,HSIZE_T),hdferr) - if (hdferr /= 0) call IO_error(0_pInt,ext_msg='HDF5_read_pReal_shape1: h5dread_f') - call h5dclose_f(dset_id,hdferr) - if (hdferr /= 0) call IO_error(0_pInt,ext_msg='HDF5_read_pReal_shape1: h5dclose_f') -end subroutine HDF5_read_pReal_1 + integer(pInt), dimension(:), allocatable :: & + globalShape, & !< shape of the dataset (all processes) + localShape, & !< shape of the dataset (this process) + readSize !< contribution of all processes + integer :: ierr + integer(HDF5_ERR_TYPE) :: hdferr + integer(HID_T) :: dset_id, filespace_id, memspace_id, plist_id, aplist_id + integer(HSIZE_T), dimension(2) :: myStart + +!------------------------------------------------------------------------------------------------- +! determine shape of dataset + localShape = shape(dataset) + if (any(localShape(1:size(localShape)) == 0)) return + +!------------------------------------------------------------------------------------------------- +! creating a property list for transfer properties + call h5pcreate_f(H5P_DATASET_XFER_F, plist_id, hdferr) !-------------------------------------------------------------------------------------------------- -!> @brief subroutine for reading dataset of the type pReal with 2 dimensions +! creating a property list for data access properties + call h5pcreate_f(H5P_DATASET_ACCESS_F, aplist_id, hdferr) !-------------------------------------------------------------------------------------------------- -subroutine HDF5_read_pReal_2(dataset,loc_id,datasetName) + allocate(readSize(worldsize), source = 0_pInt) + readSize(worldrank+1) = localShape(2) + +#ifdef PETSc + if (present(parallel)) then; 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='HDF5_read_pInt2: h5pset_dxpl_mpio_f') + call MPI_allreduce(MPI_IN_PLACE,readSize,worldsize,MPI_LONG,MPI_SUM,PETSC_COMM_WORLD,ierr) ! get total output size over each process + if (ierr /= 0) call IO_error(894_pInt,ext_msg='HDF5_read_pInt2: MPI_allreduce') + endif; endif +#endif + + myStart = int([0,sum(readSize(1:worldrank))],HSIZE_T) + globalShape = [localShape(1:1),sum(readSize)] + +!-------------------------------------------------------------------------------------------------- +! create dataspace in memory (local shape) + call h5screate_simple_f(size(localShape), int(localShape,HSIZE_T), memspace_id, hdferr, & + int(localShape,HSIZE_T)) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_read_pInt{}: h5screate_simple_f/memspace_id') +!-------------------------------------------------------------------------------------------------- +! set I/O mode for read operations to collective + call h5pset_all_coll_metadata_ops_f(aplist_id, .true., hdferr) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_read_pInt2: h5pset_all_coll_metadata_ops_f') + +!-------------------------------------------------------------------------------------------------- +! open the dataset in the file + call h5dopen_f(loc_id,datasetName,dset_id,hdferr,dapl_id = aplist_id) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_read_pInt2: h5dopen_f') + +!-------------------------------------------------------------------------------------------------- +! get the space_id of dataset in the file + call h5dget_space_f(dset_id, filespace_id, hdferr) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_read_pInt2: h5dget_space_f') + +!-------------------------------------------------------------------------------------------------- +! select a hyperslab (the portion of the current process) in the file + call h5sselect_hyperslab_f(filespace_id, H5S_SELECT_SET_F, myStart, int(localShape,HSIZE_T), hdferr) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_read_pInt2: h5sselect_hyperslab_f') + +!-------------------------------------------------------------------------------------------------- +! read + call h5dread_f(dset_id, H5T_NATIVE_INTEGER,dataset,int(globalShape,HSIZE_T), 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') + +!-------------------------------------------------------------------------------------------------- +!close types, dataspaces + call h5pclose_f(plist_id, hdferr) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_read_pInt2: plist_id') + call h5dclose_f(dset_id, hdferr) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_read_pInt2: h5dclose_f') + call h5sclose_f(filespace_id, hdferr) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_read_pInt2: h5sclose_f/filespace_id') + call h5sclose_f(memspace_id, hdferr) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_read_pInt2: h5sclose_f/memspace_id') + +end subroutine HDF5_read_pInt2 + + +!-------------------------------------------------------------------------------------------------- +!> @brief subroutine for reading dataset of type pInt with 3 dimensions +!-------------------------------------------------------------------------------------------------- +subroutine HDF5_read_pInt3(loc_id,dataset,datasetName,parallel) + use numerics, only: & + worldrank, & + worldsize implicit none - real(pReal), intent(out), dimension(:,:) :: dataset + 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 - integer(pInt),dimension(:), allocatable :: myShape + logical, intent(in), optional :: parallel - integer :: hdferr - integer(HID_T) :: dset_id - myShape = shape(dataset) - call h5dopen_f(loc_id,datasetName,dset_id,hdferr) - if (hdferr /= 0) call IO_error(0_pInt,ext_msg='HDF5_read_pReal_shape2: h5dopen_f') - call h5dread_f(dset_id,H5T_NATIVE_DOUBLE,dataset,int(myShape,HSIZE_T),hdferr) - if (hdferr /= 0) call IO_error(0_pInt,ext_msg='HDF5_read_pReal_shape2: h5dread_f') - call h5dclose_f(dset_id,hdferr) - if (hdferr /= 0) call IO_error(0_pInt,ext_msg='HDF5_read_pReal_shape2: h5dclose_f') + integer(pInt), dimension(:), allocatable :: & + globalShape, & !< shape of the dataset (all processes) + localShape, & !< shape of the dataset (this process) + readSize !< contribution of all processes + integer :: ierr + integer(HDF5_ERR_TYPE) :: hdferr + integer(HID_T) :: dset_id, filespace_id, memspace_id, plist_id, aplist_id + integer(HSIZE_T), dimension(3) :: myStart -end subroutine HDF5_read_pReal_2 +!------------------------------------------------------------------------------------------------- +! determine shape of dataset + localShape = shape(dataset) + if (any(localShape(1:size(localShape)) == 0)) return + +!------------------------------------------------------------------------------------------------- +! creating a property list for transfer properties + call h5pcreate_f(H5P_DATASET_XFER_F, plist_id, hdferr) !-------------------------------------------------------------------------------------------------- -!> @brief subroutine for reading dataset of the type pReal with 3 dimensions +! creating a property list for data access properties + call h5pcreate_f(H5P_DATASET_ACCESS_F, aplist_id, hdferr) !-------------------------------------------------------------------------------------------------- -subroutine HDF5_read_pReal_3(dataset,loc_id,datasetName) + allocate(readSize(worldsize), source = 0_pInt) + readSize(worldrank+1) = localShape(3) + +#ifdef PETSc + if (present(parallel)) then; 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='HDF5_read_pInt3: h5pset_dxpl_mpio_f') + call MPI_allreduce(MPI_IN_PLACE,readSize,worldsize,MPI_LONG,MPI_SUM,PETSC_COMM_WORLD,ierr) ! get total output size over each process + if (ierr /= 0) call IO_error(894_pInt,ext_msg='HDF5_read_pInt3: MPI_allreduce') + endif; endif +#endif + + myStart = int([0,0,sum(readSize(1:worldrank))],HSIZE_T) + globalShape = [localShape(1:2),sum(readSize)] + +!-------------------------------------------------------------------------------------------------- +! create dataspace in memory (local shape) + call h5screate_simple_f(size(localShape), int(localShape,HSIZE_T), memspace_id, hdferr, & + int(localShape,HSIZE_T)) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_read_pInt{}: h5screate_simple_f/memspace_id') +!-------------------------------------------------------------------------------------------------- +! set I/O mode for read operations to collective + call h5pset_all_coll_metadata_ops_f(aplist_id, .true., hdferr) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_read_pInt3: h5pset_all_coll_metadata_ops_f') + +!-------------------------------------------------------------------------------------------------- +! open the dataset in the file + call h5dopen_f(loc_id,datasetName,dset_id,hdferr,dapl_id = aplist_id) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_read_pInt3: h5dopen_f') + +!-------------------------------------------------------------------------------------------------- +! get the space_id of dataset in the file + call h5dget_space_f(dset_id, filespace_id, hdferr) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_read_pInt3: h5dget_space_f') + +!-------------------------------------------------------------------------------------------------- +! select a hyperslab (the portion of the current process) in the file + call h5sselect_hyperslab_f(filespace_id, H5S_SELECT_SET_F, myStart, int(localShape,HSIZE_T), hdferr) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_read_pInt3: h5sselect_hyperslab_f') + +!-------------------------------------------------------------------------------------------------- +! read + call h5dread_f(dset_id, H5T_NATIVE_INTEGER,dataset,int(globalShape,HSIZE_T), 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') + +!-------------------------------------------------------------------------------------------------- +!close types, dataspaces + call h5pclose_f(plist_id, hdferr) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_read_pInt3: plist_id') + call h5dclose_f(dset_id, hdferr) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_read_pInt3: h5dclose_f') + call h5sclose_f(filespace_id, hdferr) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_read_pInt3: h5sclose_f/filespace_id') + call h5sclose_f(memspace_id, hdferr) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_read_pInt3: h5sclose_f/memspace_id') + +end subroutine HDF5_read_pInt3 + + +!-------------------------------------------------------------------------------------------------- +!> @brief subroutine for reading dataset of type pInt with 4 dimensions +!-------------------------------------------------------------------------------------------------- +subroutine HDF5_read_pInt4(loc_id,dataset,datasetName,parallel) + use numerics, only: & + worldrank, & + worldsize implicit none - real(pReal), intent(out), dimension(:,:,:) :: dataset + 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 - integer(pInt),dimension(:), allocatable :: myShape - integer :: hdferr - integer(HID_T) :: dset_id - myShape = shape(dataset) + logical, intent(in), optional :: parallel - call h5dopen_f(loc_id,datasetName,dset_id,hdferr) - if (hdferr /= 0) call IO_error(0_pInt,ext_msg='HDF5_read_pReal_shape3: h5dopen_f') - call h5dread_f(dset_id,H5T_NATIVE_DOUBLE,dataset,int(myShape,HSIZE_T),hdferr) - if (hdferr /= 0) call IO_error(0_pInt,ext_msg='HDF5_read_pReal_shape3: h5dread_f') - call h5dclose_f(dset_id,hdferr) - if (hdferr /= 0) call IO_error(0_pInt,ext_msg='HDF5_read_pReal_shape3: h5dclose_f') -end subroutine HDF5_read_pReal_3 + integer(pInt), dimension(:), allocatable :: & + globalShape, & !< shape of the dataset (all processes) + localShape, & !< shape of the dataset (this process) + readSize !< contribution of all processes + integer :: ierr + integer(HDF5_ERR_TYPE) :: hdferr + integer(HID_T) :: dset_id, filespace_id, memspace_id, plist_id, aplist_id + integer(HSIZE_T), dimension(4) :: myStart + +!------------------------------------------------------------------------------------------------- +! determine shape of dataset + localShape = shape(dataset) + if (any(localShape(1:size(localShape)) == 0)) return + +!------------------------------------------------------------------------------------------------- +! creating a property list for transfer properties + call h5pcreate_f(H5P_DATASET_XFER_F, plist_id, hdferr) !-------------------------------------------------------------------------------------------------- -!> @brief subroutine for reading dataset of the type pReal with 4 dimensions +! creating a property list for data access properties + call h5pcreate_f(H5P_DATASET_ACCESS_F, aplist_id, hdferr) !-------------------------------------------------------------------------------------------------- -subroutine HDF5_read_pReal_4(dataset,loc_id,datasetName) + allocate(readSize(worldsize), source = 0_pInt) + readSize(worldrank+1) = localShape(4) + +#ifdef PETSc + if (present(parallel)) then; 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='HDF5_read_pInt4: h5pset_dxpl_mpio_f') + call MPI_allreduce(MPI_IN_PLACE,readSize,worldsize,MPI_LONG,MPI_SUM,PETSC_COMM_WORLD,ierr) ! get total output size over each process + if (ierr /= 0) call IO_error(894_pInt,ext_msg='HDF5_read_pInt4: MPI_allreduce') + endif; endif +#endif + + myStart = int([0,0,0,sum(readSize(1:worldrank))],HSIZE_T) + globalShape = [localShape(1:3),sum(readSize)] + +!-------------------------------------------------------------------------------------------------- +! create dataspace in memory (local shape) + call h5screate_simple_f(size(localShape), int(localShape,HSIZE_T), memspace_id, hdferr, & + int(localShape,HSIZE_T)) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_read_pInt{}: h5screate_simple_f/memspace_id') +!-------------------------------------------------------------------------------------------------- +! set I/O mode for read operations to collective + call h5pset_all_coll_metadata_ops_f(aplist_id, .true., hdferr) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_read_pInt4: h5pset_all_coll_metadata_ops_f') + +!-------------------------------------------------------------------------------------------------- +! open the dataset in the file + call h5dopen_f(loc_id,datasetName,dset_id,hdferr,dapl_id = aplist_id) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_read_pInt4: h5dopen_f') + +!-------------------------------------------------------------------------------------------------- +! get the space_id of dataset in the file + call h5dget_space_f(dset_id, filespace_id, hdferr) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_read_pInt4: h5dget_space_f') + +!-------------------------------------------------------------------------------------------------- +! select a hyperslab (the portion of the current process) in the file + call h5sselect_hyperslab_f(filespace_id, H5S_SELECT_SET_F, myStart, int(localShape,HSIZE_T), hdferr) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_read_pInt4: h5sselect_hyperslab_f') + +!-------------------------------------------------------------------------------------------------- +! read + call h5dread_f(dset_id, H5T_NATIVE_INTEGER,dataset,int(globalShape,HSIZE_T), 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') + +!-------------------------------------------------------------------------------------------------- +!close types, dataspaces + call h5pclose_f(plist_id, hdferr) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_read_pInt4: plist_id') + call h5dclose_f(dset_id, hdferr) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_read_pInt4: h5dclose_f') + call h5sclose_f(filespace_id, hdferr) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_read_pInt4: h5sclose_f/filespace_id') + call h5sclose_f(memspace_id, hdferr) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_read_pInt4: h5sclose_f/memspace_id') + +end subroutine HDF5_read_pInt4 + + +!-------------------------------------------------------------------------------------------------- +!> @brief subroutine for reading dataset of type pInt with 5 dimensions +!-------------------------------------------------------------------------------------------------- +subroutine HDF5_read_pInt5(loc_id,dataset,datasetName,parallel) + use numerics, only: & + worldrank, & + worldsize implicit none - real(pReal), intent(out), dimension(:,:,:,:) :: dataset + 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 - integer(pInt),dimension(:), allocatable :: myShape + logical, intent(in), optional :: parallel - integer :: hdferr - integer(HID_T) :: dset_id - myShape = shape(dataset) - call h5dopen_f(loc_id,datasetName,dset_id,hdferr) - if (hdferr /= 0) call IO_error(0_pInt,ext_msg='HDF5_read_pReal_shape4: h5dopen_f') - call h5dread_f(dset_id,H5T_NATIVE_DOUBLE,dataset,int(myShape,HSIZE_T),hdferr) - if (hdferr /= 0) call IO_error(0_pInt,ext_msg='HDF5_read_pReal_shape4: h5dread_f') - call h5dclose_f(dset_id,hdferr) - if (hdferr /= 0) call IO_error(0_pInt,ext_msg='HDF5_read_pReal_shape4: h5dclose_f') + integer(pInt), dimension(:), allocatable :: & + globalShape, & !< shape of the dataset (all processes) + localShape, & !< shape of the dataset (this process) + readSize !< contribution of all processes + integer :: ierr + integer(HDF5_ERR_TYPE) :: hdferr + integer(HID_T) :: dset_id, filespace_id, memspace_id, plist_id, aplist_id + integer(HSIZE_T), dimension(5) :: myStart -end subroutine HDF5_read_pReal_4 +!------------------------------------------------------------------------------------------------- +! determine shape of dataset + localShape = shape(dataset) + if (any(localShape(1:size(localShape)) == 0)) return + +!------------------------------------------------------------------------------------------------- +! creating a property list for transfer properties + call h5pcreate_f(H5P_DATASET_XFER_F, plist_id, hdferr) !-------------------------------------------------------------------------------------------------- -!> @brief subroutine for reading dataset of the type pReal with 5 dimensions +! creating a property list for data access properties + call h5pcreate_f(H5P_DATASET_ACCESS_F, aplist_id, hdferr) !-------------------------------------------------------------------------------------------------- -subroutine HDF5_read_pReal_5(dataset,loc_id,datasetName) + allocate(readSize(worldsize), source = 0_pInt) + readSize(worldrank+1) = localShape(5) + +#ifdef PETSc + if (present(parallel)) then; 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='HDF5_read_pInt5: h5pset_dxpl_mpio_f') + call MPI_allreduce(MPI_IN_PLACE,readSize,worldsize,MPI_LONG,MPI_SUM,PETSC_COMM_WORLD,ierr) ! get total output size over each process + if (ierr /= 0) call IO_error(894_pInt,ext_msg='HDF5_read_pInt5: MPI_allreduce') + endif; endif +#endif + + myStart = int([0,0,0,0,sum(readSize(1:worldrank))],HSIZE_T) + globalShape = [localShape(1:4),sum(readSize)] + +!-------------------------------------------------------------------------------------------------- +! create dataspace in memory (local shape) + call h5screate_simple_f(size(localShape), int(localShape,HSIZE_T), memspace_id, hdferr, & + int(localShape,HSIZE_T)) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_read_pInt{}: h5screate_simple_f/memspace_id') +!-------------------------------------------------------------------------------------------------- +! set I/O mode for read operations to collective + call h5pset_all_coll_metadata_ops_f(aplist_id, .true., hdferr) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_read_pInt5: h5pset_all_coll_metadata_ops_f') + +!-------------------------------------------------------------------------------------------------- +! open the dataset in the file + call h5dopen_f(loc_id,datasetName,dset_id,hdferr,dapl_id = aplist_id) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_read_pInt5: h5dopen_f') + +!-------------------------------------------------------------------------------------------------- +! get the space_id of dataset in the file + call h5dget_space_f(dset_id, filespace_id, hdferr) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_read_pInt5: h5dget_space_f') + +!-------------------------------------------------------------------------------------------------- +! select a hyperslab (the portion of the current process) in the file + call h5sselect_hyperslab_f(filespace_id, H5S_SELECT_SET_F, myStart, int(localShape,HSIZE_T), hdferr) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_read_pInt5: h5sselect_hyperslab_f') + +!-------------------------------------------------------------------------------------------------- +! read + call h5dread_f(dset_id, H5T_NATIVE_INTEGER,dataset,int(globalShape,HSIZE_T), 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') + +!-------------------------------------------------------------------------------------------------- +!close types, dataspaces + call h5pclose_f(plist_id, hdferr) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_read_pInt5: plist_id') + call h5dclose_f(dset_id, hdferr) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_read_pInt5: h5dclose_f') + call h5sclose_f(filespace_id, hdferr) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_read_pInt5: h5sclose_f/filespace_id') + call h5sclose_f(memspace_id, hdferr) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_read_pInt5: h5sclose_f/memspace_id') + +end subroutine HDF5_read_pInt5 + + +!-------------------------------------------------------------------------------------------------- +!> @brief subroutine for reading dataset of type pInt with 6 dimensions +!-------------------------------------------------------------------------------------------------- +subroutine HDF5_read_pInt6(loc_id,dataset,datasetName,parallel) + use numerics, only: & + worldrank, & + worldsize implicit none - real(pReal), intent(out), dimension(:,:,:,:,:) :: dataset + 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 - integer(pInt),dimension(:), allocatable :: myShape + logical, intent(in), optional :: parallel - integer :: hdferr - integer(HID_T) :: dset_id - myShape = shape(dataset) - call h5dopen_f(loc_id,datasetName,dset_id,hdferr) - if (hdferr /= 0) call IO_error(0_pInt,ext_msg='HDF5_read_pReal_shape5: h5dopen_f') - call h5dread_f(dset_id,H5T_NATIVE_DOUBLE,dataset,int(myShape,HSIZE_T),hdferr) - if (hdferr /= 0) call IO_error(0_pInt,ext_msg='HDF5_read_pReal_shape5: h5dread_f') - call h5dclose_f(dset_id,hdferr) - if (hdferr /= 0) call IO_error(0_pInt,ext_msg='HDF5_read_pReal_shape5: h5dclose_f') + integer(pInt), dimension(:), allocatable :: & + globalShape, & !< shape of the dataset (all processes) + localShape, & !< shape of the dataset (this process) + readSize !< contribution of all processes + integer :: ierr + integer(HDF5_ERR_TYPE) :: hdferr + integer(HID_T) :: dset_id, filespace_id, memspace_id, plist_id, aplist_id + integer(HSIZE_T), dimension(6) :: myStart -end subroutine HDF5_read_pReal_5 +!------------------------------------------------------------------------------------------------- +! determine shape of dataset + localShape = shape(dataset) + if (any(localShape(1:size(localShape)) == 0)) return + +!------------------------------------------------------------------------------------------------- +! creating a property list for transfer properties + call h5pcreate_f(H5P_DATASET_XFER_F, plist_id, hdferr) !-------------------------------------------------------------------------------------------------- -!> @brief subroutine for reading dataset of the type pReal with 6 dimensions +! creating a property list for data access properties + call h5pcreate_f(H5P_DATASET_ACCESS_F, aplist_id, hdferr) !-------------------------------------------------------------------------------------------------- -subroutine HDF5_read_pReal_6(dataset,loc_id,datasetName) + allocate(readSize(worldsize), source = 0_pInt) + readSize(worldrank+1) = localShape(6) + +#ifdef PETSc + if (present(parallel)) then; 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='HDF5_read_pInt6: h5pset_dxpl_mpio_f') + call MPI_allreduce(MPI_IN_PLACE,readSize,worldsize,MPI_LONG,MPI_SUM,PETSC_COMM_WORLD,ierr) ! get total output size over each process + if (ierr /= 0) call IO_error(894_pInt,ext_msg='HDF5_read_pInt6: MPI_allreduce') + endif; endif +#endif + + myStart = int([0,0,0,0,0,sum(readSize(1:worldrank))],HSIZE_T) + globalShape = [localShape(1:5),sum(readSize)] + +!-------------------------------------------------------------------------------------------------- +! create dataspace in memory (local shape) + call h5screate_simple_f(size(localShape), int(localShape,HSIZE_T), memspace_id, hdferr, & + int(localShape,HSIZE_T)) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_read_pInt{}: h5screate_simple_f/memspace_id') +!-------------------------------------------------------------------------------------------------- +! set I/O mode for read operations to collective + call h5pset_all_coll_metadata_ops_f(aplist_id, .true., hdferr) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_read_pInt6: h5pset_all_coll_metadata_ops_f') + +!-------------------------------------------------------------------------------------------------- +! open the dataset in the file + call h5dopen_f(loc_id,datasetName,dset_id,hdferr,dapl_id = aplist_id) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_read_pInt6: h5dopen_f') + +!-------------------------------------------------------------------------------------------------- +! get the space_id of dataset in the file + call h5dget_space_f(dset_id, filespace_id, hdferr) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_read_pInt6: h5dget_space_f') + +!-------------------------------------------------------------------------------------------------- +! select a hyperslab (the portion of the current process) in the file + call h5sselect_hyperslab_f(filespace_id, H5S_SELECT_SET_F, myStart, int(localShape,HSIZE_T), hdferr) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_read_pInt6: h5sselect_hyperslab_f') + +!-------------------------------------------------------------------------------------------------- +! read + call h5dread_f(dset_id, H5T_NATIVE_INTEGER,dataset,int(globalShape,HSIZE_T), 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') + +!-------------------------------------------------------------------------------------------------- +!close types, dataspaces + call h5pclose_f(plist_id, hdferr) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_read_pInt6: plist_id') + call h5dclose_f(dset_id, hdferr) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_read_pInt6: h5dclose_f') + call h5sclose_f(filespace_id, hdferr) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_read_pInt6: h5sclose_f/filespace_id') + call h5sclose_f(memspace_id, hdferr) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_read_pInt6: h5sclose_f/memspace_id') + +end subroutine HDF5_read_pInt6 + + +!-------------------------------------------------------------------------------------------------- +!> @brief subroutine for reading dataset of type pInt with 7 dimensions +!-------------------------------------------------------------------------------------------------- +subroutine HDF5_read_pInt7(loc_id,dataset,datasetName,parallel) + use numerics, only: & + worldrank, & + worldsize implicit none - real(pReal), intent(out), dimension(:,:,:,:,:,:) :: dataset + 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 - integer(pInt),dimension(:), allocatable :: myShape + logical, intent(in), optional :: parallel - integer :: hdferr - integer(HID_T) :: dset_id - myShape = shape(dataset) - call h5dopen_f(loc_id,datasetName,dset_id,hdferr) - if (hdferr /= 0) call IO_error(0_pInt,ext_msg='HDF5_read_pReal_shape6: h5dopen_f') - call h5dread_f(dset_id,H5T_NATIVE_DOUBLE,dataset,int(myShape,HSIZE_T),hdferr) - if (hdferr /= 0) call IO_error(0_pInt,ext_msg='HDF5_read_pReal_shape6: h5dread_f') - call h5dclose_f(dset_id,hdferr) - if (hdferr /= 0) call IO_error(0_pInt,ext_msg='HDF5_read_pReal_shape6: h5dclose_f') + integer(pInt), dimension(:), allocatable :: & + globalShape, & !< shape of the dataset (all processes) + localShape, & !< shape of the dataset (this process) + readSize !< contribution of all processes + integer :: ierr + integer(HDF5_ERR_TYPE) :: hdferr + integer(HID_T) :: dset_id, filespace_id, memspace_id, plist_id, aplist_id + integer(HSIZE_T), dimension(7) :: myStart -end subroutine HDF5_read_pReal_6 +!------------------------------------------------------------------------------------------------- +! determine shape of dataset + localShape = shape(dataset) + if (any(localShape(1:size(localShape)) == 0)) return + +!------------------------------------------------------------------------------------------------- +! creating a property list for transfer properties + call h5pcreate_f(H5P_DATASET_XFER_F, plist_id, hdferr) !-------------------------------------------------------------------------------------------------- -!> @brief subroutine for reading dataset of the type pReal with 7 dimensions +! creating a property list for data access properties + call h5pcreate_f(H5P_DATASET_ACCESS_F, aplist_id, hdferr) !-------------------------------------------------------------------------------------------------- -subroutine HDF5_read_pReal_7(dataset,loc_id,datasetName) + allocate(readSize(worldsize), source = 0_pInt) + readSize(worldrank+1) = localShape(7) + +#ifdef PETSc + if (present(parallel)) then; 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='HDF5_read_pInt7: h5pset_dxpl_mpio_f') + call MPI_allreduce(MPI_IN_PLACE,readSize,worldsize,MPI_LONG,MPI_SUM,PETSC_COMM_WORLD,ierr) ! get total output size over each process + if (ierr /= 0) call IO_error(894_pInt,ext_msg='HDF5_read_pInt7: MPI_allreduce') + endif; endif +#endif + + myStart = int([0,0,0,0,0,0,sum(readSize(1:worldrank))],HSIZE_T) + globalShape = [localShape(1:6),sum(readSize)] + +!-------------------------------------------------------------------------------------------------- +! create dataspace in memory (local shape) + call h5screate_simple_f(size(localShape), int(localShape,HSIZE_T), memspace_id, hdferr, & + int(localShape,HSIZE_T)) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_read_pInt{}: h5screate_simple_f/memspace_id') +!-------------------------------------------------------------------------------------------------- +! set I/O mode for read operations to collective + call h5pset_all_coll_metadata_ops_f(aplist_id, .true., hdferr) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_read_pInt7: h5pset_all_coll_metadata_ops_f') + +!-------------------------------------------------------------------------------------------------- +! open the dataset in the file + call h5dopen_f(loc_id,datasetName,dset_id,hdferr,dapl_id = aplist_id) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_read_pInt7: h5dopen_f') + +!-------------------------------------------------------------------------------------------------- +! get the space_id of dataset in the file + call h5dget_space_f(dset_id, filespace_id, hdferr) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_read_pInt7: h5dget_space_f') + +!-------------------------------------------------------------------------------------------------- +! select a hyperslab (the portion of the current process) in the file + call h5sselect_hyperslab_f(filespace_id, H5S_SELECT_SET_F, myStart, int(localShape,HSIZE_T), hdferr) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_read_pInt7: h5sselect_hyperslab_f') + +!-------------------------------------------------------------------------------------------------- +! read + call h5dread_f(dset_id, H5T_NATIVE_INTEGER,dataset,int(globalShape,HSIZE_T), 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') + +!-------------------------------------------------------------------------------------------------- +!close types, dataspaces + call h5pclose_f(plist_id, hdferr) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_read_pInt7: plist_id') + call h5dclose_f(dset_id, hdferr) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_read_pInt7: h5dclose_f') + call h5sclose_f(filespace_id, hdferr) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_read_pInt7: h5sclose_f/filespace_id') + call h5sclose_f(memspace_id, hdferr) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_read_pInt7: h5sclose_f/memspace_id') + +end subroutine HDF5_read_pInt7 + +!-------------------------------------------------------------------------------------------------- +!> @brief subroutine for writing dataset of type pReal with 1 dimensions +!-------------------------------------------------------------------------------------------------- +subroutine HDF5_write_pReal1(loc_id,dataset,datasetName,parallel) + use numerics, only: & + worldrank, & + worldsize implicit none - real(pReal), intent(out), dimension(:,:,:,:,:,:,:) :: dataset + 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 - integer(pInt),dimension(:), allocatable :: myShape + logical, intent(in), optional :: parallel - integer :: hdferr - integer(HID_T) :: dset_id - myShape = shape(dataset) - call h5dopen_f(loc_id,datasetName,dset_id,hdferr) - if (hdferr /= 0) call IO_error(0_pInt,ext_msg='HDF5_read_pReal_shape7: h5dopen_f') - call h5dread_f(dset_id,H5T_NATIVE_DOUBLE,dataset,int(myShape,HSIZE_T),hdferr) - if (hdferr /= 0) call IO_error(0_pInt,ext_msg='HDF5_read_pReal_shape7: h5dread_f') - call h5dclose_f(dset_id,hdferr) - if (hdferr /= 0) call IO_error(0_pInt,ext_msg='HDF5_read_pReal_shape7: h5dclose_f') + integer(pInt), dimension(:), allocatable :: & + globalShape, & !< shape of the dataset (all processes) + localShape, & !< shape of the dataset (this process) + outputSize !< contribution of all processes + integer :: ierr + integer(HDF5_ERR_TYPE) :: hdferr + integer(HID_T) :: dset_id, filespace_id, memspace_id, plist_id + integer(HSIZE_T), dimension(1) :: myStart -end subroutine HDF5_read_pReal_7 +!------------------------------------------------------------------------------------------------- +! determine shape of dataset + localShape = shape(dataset) + if (any(localShape(1:size(localShape)) == 0)) return + +!------------------------------------------------------------------------------------------------- +! creating a property list for transfer properties + call h5pcreate_f(H5P_DATASET_XFER_F, plist_id, hdferr) !-------------------------------------------------------------------------------------------------- -!> @brief subroutine for reading dataset of the type pInt with 1 dimension -!-------------------------------------------------------------------------------------------------- -subroutine HDF5_read_pInt_1(dataset,loc_id,datasetName) + allocate(outputSize(worldsize), source = 0_pInt) + outputSize(worldrank+1) = localShape(1) - implicit none - integer(pInt), intent(out), 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 - integer(pInt),dimension(:), allocatable :: myShape +#ifdef PETSc + if (present(parallel)) then; 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='HDF5_write_pReal1: h5pset_dxpl_mpio_f') + call MPI_allreduce(MPI_IN_PLACE,outputSize,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='HDF5_write_pReal1: MPI_allreduce') + endif; endif +#endif - integer :: hdferr - integer(HID_T) :: dset_id - myShape = shape(dataset) + myStart = int([sum(outputSize(1:worldrank))],HSIZE_T) + globalShape = [localShape(1:0),sum(outputSize)] - call h5dopen_f(loc_id,datasetName,dset_id,hdferr) - if (hdferr /= 0) call IO_error(0_pInt,ext_msg='HDF5_read_pInt__shape1: h5dopen_f') - call h5dread_f(dset_id,H5T_NATIVE_INTEGER,dataset,int(myShape,HSIZE_T),hdferr) - if (hdferr /= 0) call IO_error(0_pInt,ext_msg='HDF5_read_pInt__shape1: h5dread_f') - call h5dclose_f(dset_id,hdferr) - if (hdferr /= 0) call IO_error(0_pInt,ext_msg='HDF5_read_pInt__shape1: h5dclose_f') - -end subroutine HDF5_read_pInt_1 !-------------------------------------------------------------------------------------------------- -!> @brief subroutine for reading dataset of the type pInt with 2 dimensions -!-------------------------------------------------------------------------------------------------- -subroutine HDF5_read_pInt_2(dataset,loc_id,datasetName) - - implicit none - integer(pInt), intent(out), 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 - integer(pInt),dimension(:), allocatable :: myShape - - integer :: hdferr - integer(HID_T) :: dset_id - myShape = shape(dataset) - - call h5dopen_f(loc_id,datasetName,dset_id,hdferr) - if (hdferr /= 0) call IO_error(0_pInt,ext_msg='HDF5_read_pInt__shape2: h5dopen_f') - call h5dread_f(dset_id,H5T_NATIVE_INTEGER,dataset,int(myShape,HSIZE_T),hdferr) - if (hdferr /= 0) call IO_error(0_pInt,ext_msg='HDF5_read_pInt__shape2: h5dread_f') - call h5dclose_f(dset_id,hdferr) - if (hdferr /= 0) call IO_error(0_pInt,ext_msg='HDF5_read_pInt__shape2: h5dclose_f') - -end subroutine HDF5_read_pInt_2 +! create dataspace in memory (local shape) + call h5screate_simple_f(size(localShape), int(localShape,HSIZE_T), memspace_id, hdferr, & + int(localShape,HSIZE_T)) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_write_pReal1: h5dopen_f') !-------------------------------------------------------------------------------------------------- -!> @brief subroutine for reading dataset of the type pInt with 3 dimensions -!-------------------------------------------------------------------------------------------------- -subroutine HDF5_read_pInt_3(dataset,loc_id,datasetName) +! create dataspace in file (global shape) + call h5screate_simple_f(size(globalShape), int(globalShape,HSIZE_T), filespace_id, hdferr, & + int(globalShape,HSIZE_T)) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_write_pReal1: h5dget_space_f') - implicit none - integer(pInt), intent(out), 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 - integer(pInt),dimension(:), allocatable :: myShape - - integer :: hdferr - integer(HID_T) :: dset_id - myShape = shape(dataset) - - call h5dopen_f(loc_id,datasetName,dset_id,hdferr) - if (hdferr /= 0) call IO_error(0_pInt,ext_msg='HDF5_read_pInt__shape3: h5dopen_f') - call h5dread_f(dset_id,H5T_NATIVE_INTEGER,dataset,int(myShape,HSIZE_T),hdferr) - if (hdferr /= 0) call IO_error(0_pInt,ext_msg='HDF5_read_pInt__shape3: h5dread_f') - call h5dclose_f(dset_id,hdferr) - if (hdferr /= 0) call IO_error(0_pInt,ext_msg='HDF5_read_pInt__shape3: h5dclose_f') - -end subroutine HDF5_read_pInt_3 - -!-------------------------------------------------------------------------------------------------- -!> @brief subroutine for reading dataset of the type pInt with 4 dimensions -!-------------------------------------------------------------------------------------------------- -subroutine HDF5_read_pInt_4(dataset,loc_id,datasetName) - - implicit none - integer(pInt), intent(out), 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 - integer(pInt),dimension(:), allocatable :: myShape - - integer :: hdferr - integer(HID_T) :: dset_id - myShape = shape(dataset) - - call h5dopen_f(loc_id,datasetName,dset_id,hdferr) - if (hdferr /= 0) call IO_error(0_pInt,ext_msg='HDF5_read_pInt__shape4: h5dopen_f') - call h5dread_f(dset_id,H5T_NATIVE_INTEGER,dataset,int(myShape,HSIZE_T),hdferr) - if (hdferr /= 0) call IO_error(0_pInt,ext_msg='HDF5_read_pInt__shape4: h5dread_f') - call h5dclose_f(dset_id,hdferr) - if (hdferr /= 0) call IO_error(0_pInt,ext_msg='HDF5_read_pInt__shape4: h5dclose_f') - -end subroutine HDF5_read_pInt_4 - -!-------------------------------------------------------------------------------------------------- -!> @brief subroutine for reading dataset of the type pInt with 5 dimensions -!-------------------------------------------------------------------------------------------------- -subroutine HDF5_read_pInt_5(dataset,loc_id,datasetName) - - implicit none - integer(pInt), intent(out), 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 - integer(pInt),dimension(:), allocatable :: myShape - - integer :: hdferr - integer(HID_T) :: dset_id - myShape = shape(dataset) - - call h5dopen_f(loc_id,datasetName,dset_id,hdferr) - if (hdferr /= 0) call IO_error(0_pInt,ext_msg='HDF5_read_pInt__shape5: h5dopen_f') - call h5dread_f(dset_id,H5T_NATIVE_INTEGER,dataset,int(myShape,HSIZE_T),hdferr) - if (hdferr /= 0) call IO_error(0_pInt,ext_msg='HDF5_read_pInt__shape5: h5dread_f') - call h5dclose_f(dset_id,hdferr) - if (hdferr /= 0) call IO_error(0_pInt,ext_msg='HDF5_read_pInt__shape5: h5dclose_f') - -end subroutine HDF5_read_pInt_5 - -!-------------------------------------------------------------------------------------------------- -!> @brief subroutine for reading dataset of the type pInt with 6 dimensions -!-------------------------------------------------------------------------------------------------- -subroutine HDF5_read_pInt_6(dataset,loc_id,datasetName) - - implicit none - integer(pInt), intent(out), 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 - integer(pInt),dimension(:), allocatable :: myShape - - integer :: hdferr - integer(HID_T) :: dset_id - myShape = shape(dataset) - - call h5dopen_f(loc_id,datasetName,dset_id,hdferr) - if (hdferr /= 0) call IO_error(0_pInt,ext_msg='HDF5_read_pInt__shape6: h5dopen_f') - call h5dread_f(dset_id,H5T_NATIVE_INTEGER,dataset,int(myShape,HSIZE_T),hdferr) - if (hdferr /= 0) call IO_error(0_pInt,ext_msg='HDF5_read_pInt__shape6: h5dread_f') - call h5dclose_f(dset_id,hdferr) - if (hdferr /= 0) call IO_error(0_pInt,ext_msg='HDF5_read_pInt__shape6: h5dclose_f') - -end subroutine HDF5_read_pInt_6 - -!-------------------------------------------------------------------------------------------------- -!> @brief subroutine for reading dataset of the type pInt with 7 dimensions -!-------------------------------------------------------------------------------------------------- -subroutine HDF5_read_pInt_7(dataset,loc_id,datasetName) - - implicit none - integer(pInt), intent(out), 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 - integer(pInt),dimension(:), allocatable :: myShape - - integer :: hdferr - integer(HID_T) :: dset_id - myShape = shape(dataset) - - call h5dopen_f(loc_id,datasetName,dset_id,hdferr) - if (hdferr /= 0) call IO_error(0_pInt,ext_msg='HDF5_read_pInt__shape7: h5dopen_f') - call h5dread_f(dset_id,H5T_NATIVE_INTEGER,dataset,int(myShape,HSIZE_T),hdferr) - if (hdferr /= 0) call IO_error(0_pInt,ext_msg='HDF5_read_pInt__shape7: h5dread_f') - call h5dclose_f(dset_id,hdferr) - if (hdferr /= 0) call IO_error(0_pInt,ext_msg='HDF5_read_pInt__shape7: h5dclose_f') - -end subroutine HDF5_read_pInt_7 - -!-------------------------------------------------------------------------------------------------- -!> @brief subroutine for writing dataset of the type pReal with 1 dimensions -!-------------------------------------------------------------------------------------------------- -subroutine HDF5_write_pReal1(dataset,loc_id,datasetName) - - implicit none - real(pReal), intent(out), 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 - - integer(pInt), dimension(:), allocatable :: myShape ! @brief subroutine for writing dataset of the type pReal with 2 dimensions +!> @brief subroutine for writing dataset of type pReal with 2 dimensions !-------------------------------------------------------------------------------------------------- -subroutine HDF5_write_pReal2(dataset,loc_id,datasetName) +subroutine HDF5_write_pReal2(loc_id,dataset,datasetName,parallel) + use numerics, only: & + worldrank, & + worldsize implicit none - real(pReal), intent(out), dimension(:,:) :: dataset + 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(pInt), dimension(:), allocatable :: myShape ! @brief subroutine for writing dataset of the type pReal with 3 dimensions +!> @brief subroutine for writing dataset of type pReal with 3 dimensions !-------------------------------------------------------------------------------------------------- -subroutine HDF5_write_pReal3(dataset,loc_id,datasetName) +subroutine HDF5_write_pReal3(loc_id,dataset,datasetName,parallel) + use numerics, only: & + worldrank, & + worldsize implicit none - real(pReal), intent(out), dimension(:,:,:) :: dataset + 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(pInt), dimension(:), allocatable :: myShape ! @brief subroutine for writing dataset of the type pReal with 4 dimensions +!> @brief subroutine for writing dataset of type pReal with 4 dimensions !-------------------------------------------------------------------------------------------------- -subroutine HDF5_write_pReal4(dataset,loc_id,datasetName) +subroutine HDF5_write_pReal4(loc_id,dataset,datasetName,parallel) + use numerics, only: & + worldrank, & + worldsize implicit none - real(pReal), intent(out), dimension(:,:,:,:) :: dataset + 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(pInt), dimension(:), allocatable :: myShape ! @brief subroutine for writing dataset of the type pReal with 5 dimensions +!> @brief subroutine for writing dataset of type pReal with 5 dimensions !-------------------------------------------------------------------------------------------------- -subroutine HDF5_write_pReal5(dataset,loc_id,datasetName) +subroutine HDF5_write_pReal5(loc_id,dataset,datasetName,parallel) + use numerics, only: & + worldrank, & + worldsize implicit none - real(pReal), intent(out), dimension(:,:,:,:,:) :: dataset + 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(pInt), dimension(:), allocatable :: myShape ! @brief subroutine for writing dataset of the type pReal with 6 dimensions +!> @brief subroutine for writing dataset of type pReal with 6 dimensions !-------------------------------------------------------------------------------------------------- -subroutine HDF5_write_pReal6(dataset,loc_id,datasetName) +subroutine HDF5_write_pReal6(loc_id,dataset,datasetName,parallel) + use numerics, only: & + worldrank, & + worldsize implicit none - real(pReal), intent(out), dimension(:,:,:,:,:,:) :: dataset + 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(pInt), dimension(:), allocatable :: myShape ! @brief subroutine for writing dataset of the type pReal with 7 dimensions +!> @brief subroutine for writing dataset of type pReal with 7 dimensions !-------------------------------------------------------------------------------------------------- -subroutine HDF5_write_pReal7(dataset,loc_id,datasetName) +subroutine HDF5_write_pReal7(loc_id,dataset,datasetName,parallel) + use numerics, only: & + worldrank, & + worldsize implicit none - real(pReal), intent(out), dimension(:,:,:,:,:,:,:) :: dataset + 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(pInt), dimension(:), allocatable :: myShape ! @brief subroutine for writing dataset of the type pInt with 1 dimensions +!> @brief subroutine for writing dataset of type pInt with 1 dimensions !-------------------------------------------------------------------------------------------------- -subroutine HDF5_write_pInt1(dataset,loc_id,datasetName) +subroutine HDF5_write_pInt1(loc_id,dataset,datasetName,parallel) + use numerics, only: & + worldrank, & + worldsize implicit none - integer(pInt), intent(out), dimension(:) :: dataset + 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(pInt), dimension(:), allocatable :: myShape ! @brief subroutine for writing dataset of the type pInt with 2 dimensions +!> @brief subroutine for writing dataset of type pInt with 2 dimensions !-------------------------------------------------------------------------------------------------- -subroutine HDF5_write_pInt2(dataset,loc_id,datasetName) +subroutine HDF5_write_pInt2(loc_id,dataset,datasetName,parallel) + use numerics, only: & + worldrank, & + worldsize implicit none - integer(pInt), intent(out), dimension(:,:) :: dataset + 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(pInt), dimension(:), allocatable :: myShape ! @brief subroutine for writing dataset of the type pInt with 3 dimensions +!> @brief subroutine for writing dataset of type pInt with 3 dimensions !-------------------------------------------------------------------------------------------------- -subroutine HDF5_write_pInt3(dataset,loc_id,datasetName) +subroutine HDF5_write_pInt3(loc_id,dataset,datasetName,parallel) + use numerics, only: & + worldrank, & + worldsize implicit none - integer(pInt), intent(out), dimension(:,:,:) :: dataset + 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(pInt), dimension(:), allocatable :: myShape ! @brief subroutine for writing dataset of the type pInt with 4 dimensions +!> @brief subroutine for writing dataset of type pInt with 4 dimensions !-------------------------------------------------------------------------------------------------- -subroutine HDF5_write_pInt4(dataset,loc_id,datasetName) +subroutine HDF5_write_pInt4(loc_id,dataset,datasetName,parallel) + use numerics, only: & + worldrank, & + worldsize implicit none - integer(pInt), intent(out), dimension(:,:,:,:) :: dataset + 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(pInt), dimension(:), allocatable :: myShape ! @brief subroutine for writing dataset of the type pInt with 5 dimensions +!> @brief subroutine for writing dataset of type pInt with 5 dimensions !-------------------------------------------------------------------------------------------------- -subroutine HDF5_write_pInt5(dataset,loc_id,datasetName) +subroutine HDF5_write_pInt5(loc_id,dataset,datasetName,parallel) + use numerics, only: & + worldrank, & + worldsize implicit none - integer(pInt), intent(out), dimension(:,:,:,:,:) :: dataset + 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(pInt), dimension(:), allocatable :: myShape ! @brief subroutine for writing dataset of the type pInt with 6 dimensions +!> @brief subroutine for writing dataset of type pInt with 6 dimensions !-------------------------------------------------------------------------------------------------- -subroutine HDF5_write_pInt6(dataset,loc_id,datasetName) +subroutine HDF5_write_pInt6(loc_id,dataset,datasetName,parallel) + use numerics, only: & + worldrank, & + worldsize implicit none - integer(pInt), intent(out), dimension(:,:,:,:,:,:) :: dataset + 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(pInt), dimension(:), allocatable :: myShape ! @brief subroutine for writing dataset of the type pInt with 7 dimensions +!> @brief subroutine for writing dataset of type pInt with 7 dimensions !-------------------------------------------------------------------------------------------------- -subroutine HDF5_write_pInt7(dataset,loc_id,datasetName) +subroutine HDF5_write_pInt7(loc_id,dataset,datasetName,parallel) + use numerics, only: & + worldrank, & + worldsize implicit none - integer(pInt), intent(out), dimension(:,:,:,:,:,:,:) :: dataset + 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(pInt), dimension(:), allocatable :: myShape ! @brief adds a new scalar dataset to the given group location -!-------------------------------------------------------------------------------------------------- -subroutine HDF5_addScalarDataset(group,nnodes,label,SIunit) - use hdf5 - implicit none - integer(HID_T), intent(in) :: group - integer(pInt), intent(in) :: nnodes - character(len=*), intent(in) :: SIunit,label +end module HDF5_Utilities + + + + + - integer :: hdferr - integer(HID_T) :: space_id, dset_id -!-------------------------------------------------------------------------------------------------- -! 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='HDF5_addScalarDataset: h5screate_simple_f') -!-------------------------------------------------------------------------------------------------- -! create Dataset - call h5dcreate_f(group, trim(label),H5T_NATIVE_DOUBLE, space_id, dset_id, hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_addScalarDataset: h5dcreate_f') - call HDF5_addStringAttribute(dset_id,'unit',trim(SIunit)) -!-------------------------------------------------------------------------------------------------- -!close types, dataspaces - call h5dclose_f(dset_id, hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_addScalarDataset: h5dclose_f') - call h5sclose_f(space_id, hdferr) -end subroutine HDF5_addScalarDataset -!-------------------------------------------------------------------------------------------------- -!> @brief copies the current temp results to the actual results file -!-------------------------------------------------------------------------------------------------- -subroutine HDF5_forwardResults(time) - use hdf5 - use IO, only: & - IO_intOut - implicit none - integer :: hdferr - integer(HID_T) :: currentIncID - real(pReal), intent(in) :: time - character(len=1024) :: myName - currentInc = currentInc +1_pInt - write(6,*) 'forward results';flush(6) - write(myName,'(a,'//IO_intOut(currentInc)//')') 'inc',currentInc - currentIncID = HDF5_addGroup(myName) - call HDF5_setLink(myName,'current') -! call HDF5_flush(resultsFile) - call HDF5_closeGroup(currentIncID) -end subroutine HDF5_forwardResults -end module HDF5_Utilities \ No newline at end of file diff --git a/src/commercialFEM_fileList.f90 b/src/commercialFEM_fileList.f90 index 8d3e9c816..4feb52bed 100644 --- a/src/commercialFEM_fileList.f90 +++ b/src/commercialFEM_fileList.f90 @@ -4,12 +4,12 @@ !> @details List of files needed by MSC.Marc, Abaqus/Explicit, and Abaqus/Standard !-------------------------------------------------------------------------------------------------- #include "IO.f90" -#ifdef DAMASKHDF5 -#include "HDF5_utilities.f90" -#endif #include "numerics.f90" #include "debug.f90" #include "config.f90" +#ifdef DAMASKHDF5 +#include "HDF5_utilities.f90" +#endif #include "math.f90" #include "FEsolving.f90" #include "mesh.f90" diff --git a/src/config.f90 b/src/config.f90 index dcc14e015..b79442e62 100644 --- a/src/config.f90 +++ b/src/config.f90 @@ -1,4 +1,4 @@ -!-------------------------------------------------------------------------------------------------- +!------------------------------------------------------------------------------------------------- !> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH !> @brief Reads in the material configuration from file !> @details Reads the material configuration file, where solverJobName.materialConfig takes @@ -220,7 +220,7 @@ subroutine parseFile(sectionNames,part,line, & partPosition = [partPosition, i] ! needed when actually storing content do i = 1_pInt, size(partPosition) -1_pInt - sectionNames(i) = trim(adjustl(fileContent(partPosition(i)))) + sectionNames(i) = trim(adjustl(IO_getTag(fileContent(partPosition(i)),'[',']'))) do j = partPosition(i) + 1_pInt, partPosition(i+1) -1_pInt call part(i)%add(trim(adjustl(fileContent(j)))) enddo diff --git a/src/constitutive.f90 b/src/constitutive.f90 index b20d6e65d..7aba62c42 100644 --- a/src/constitutive.f90 +++ b/src/constitutive.f90 @@ -25,7 +25,8 @@ module constitutive constitutive_SandItsTangents, & constitutive_collectDotState, & constitutive_collectDeltaState, & - constitutive_postResults + constitutive_postResults, & + constitutive_results private :: & constitutive_hooke_SandItsTangents @@ -1126,4 +1127,43 @@ function constitutive_postResults(S6, Fi, FeArray, ipc, ip, el) end function constitutive_postResults + +!-------------------------------------------------------------------------------------------------- +!> @brief writes constitutive results to HDF5 output file +!-------------------------------------------------------------------------------------------------- +subroutine constitutive_results() + use material, only: & + PLASTICITY_ISOTROPIC_ID, & + PLASTICITY_PHENOPOWERLAW_ID, & + PLASTICITY_KINEHARDENING_ID, & + PLASTICITY_DISLOTWIN_ID, & + PLASTICITY_DISLOUCLA_ID, & + PLASTICITY_NONLOCAL_ID +#if defined(PETSc) || defined(DAMASKHDF5) + use results + use HDF5_utilities + use config, only: & + config_name_phase => phase_name ! anticipate logical name + + use material, only: & + phase_plasticityInstance, & + material_phase_plasticity_type => phase_plasticity + use plastic_phenopowerlaw, only: & + plastic_phenopowerlaw_results + + implicit none + integer(pInt) :: 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)))) + if (material_phase_plasticity_type(p) == PLASTICITY_PHENOPOWERLAW_ID) then + call plastic_phenopowerlaw_results(phase_plasticityInstance(p),'current/phase/'//trim(config_name_phase(p))) + endif + enddo + +#endif + + +end subroutine constitutive_results + end module constitutive diff --git a/src/plastic_phenopowerlaw.f90 b/src/plastic_phenopowerlaw.f90 index 82050086e..0fe63737e 100644 --- a/src/plastic_phenopowerlaw.f90 +++ b/src/plastic_phenopowerlaw.f90 @@ -92,13 +92,15 @@ module plastic_phenopowerlaw plastic_phenopowerlaw_init, & plastic_phenopowerlaw_LpAndItsTangent, & plastic_phenopowerlaw_dotState, & - plastic_phenopowerlaw_postResults + plastic_phenopowerlaw_postResults, & + plastic_phenopowerlaw_results private :: & kinetics_slip, & kinetics_twin contains + !-------------------------------------------------------------------------------------------------- !> @brief module initialization !> @details reads in material parameters, allocates arrays, and does sanity checks @@ -575,6 +577,35 @@ function plastic_phenopowerlaw_postResults(Mp,instance,of) result(postResults) 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 + + implicit none + integer(pInt), intent(in) :: instance + character(len=*) :: group + integer(pInt) :: o + + associate(prm => param(instance), stt => state(instance)) + outputsLoop: do o = 1_pInt,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 +#else + integer(pInt), intent(in) :: instance + character(len=*) :: group +#endif +end subroutine plastic_phenopowerlaw_results + + !-------------------------------------------------------------------------------------------------- !> @brief Shear rates on slip systems and their derivatives with respect to resolved stress !> @details Derivatives are calculated only optionally. diff --git a/src/quit.f90 b/src/quit.f90 index 4219830a5..ad61943e4 100644 --- a/src/quit.f90 +++ b/src/quit.f90 @@ -23,6 +23,8 @@ subroutine quit(stop_id) integer(pInt) :: error = 0_pInt PetscErrorCode :: ierr = 0 + call h5open_f(hdferr) + if (hdferr /= 0) write(6,'(a,i5)') ' Error in h5open_f',hdferr ! prevents error if not opened yet call h5close_f(hdferr) if (hdferr /= 0) write(6,'(a,i5)') ' Error in h5close_f',hdferr diff --git a/src/results.f90 b/src/results.f90 new file mode 100644 index 000000000..d38178993 --- /dev/null +++ b/src/results.f90 @@ -0,0 +1,974 @@ +!-------------------------------------------------------------------------------------------------- +!> @author Vitesh Shah, Max-Planck-Institut für Eisenforschung GmbH +!> @author Yi-Chin Yang, Max-Planck-Institut für Eisenforschung GmbH +!> @author Jennifer Nastola, Max-Planck-Institut für Eisenforschung GmbH +!> @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 +#endif + + implicit none + private + integer(HID_T), public, protected :: tempCoordinates, tempResults + integer(HID_T), private :: resultsFile, currentIncID, plist_id + + + public :: & + results_init, & + results_openJobFile, & + results_closeJobFile, & + results_addIncrement, & + results_addGroup, & + results_openGroup, & + results_writeVectorDataset, & + results_setLink, & + results_removeLink +contains + +subroutine results_init + use, intrinsic :: & + iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment) + use DAMASK_interface, only: & + getSolverJobName + implicit none + + write(6,'(/,a)') ' <<<+- results init -+>>>' +#include "compilation_info.f90" + + call HDF5_closeFile(HDF5_openFile(trim(getSolverJobName())//'.hdf5','w',.true.)) + +end subroutine results_init + + +!-------------------------------------------------------------------------------------------------- +!> @brief opens the results file to append data +!-------------------------------------------------------------------------------------------------- +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)) + +end subroutine results_openJobFile + + +!-------------------------------------------------------------------------------------------------- +!> @brief closes the results file +!-------------------------------------------------------------------------------------------------- +subroutine results_closeJobFile() + implicit none + + call HDF5_closeFile(resultsFile) + +end subroutine results_closeJobFile + + +!-------------------------------------------------------------------------------------------------- +!> @brief closes the results file +!-------------------------------------------------------------------------------------------------- +subroutine results_addIncrement(inc,time) + + implicit none + 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)))) + +end subroutine results_addIncrement + +!-------------------------------------------------------------------------------------------------- +!> @brief open a group from the results file +!-------------------------------------------------------------------------------------------------- +integer(HID_T) function results_openGroup(groupName) + + implicit none + character(len=*), intent(in) :: groupName + + results_openGroup = HDF5_openGroup(resultsFile,groupName) + +end function results_openGroup + + +!-------------------------------------------------------------------------------------------------- +!> @brief adds a new group to the results file +!-------------------------------------------------------------------------------------------------- +integer(HID_T) function results_addGroup(groupName) + + implicit none + character(len=*), intent(in) :: groupName + + results_addGroup = HDF5_addGroup(resultsFile,groupName) + +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 + + call HDF5_setLink(resultsFile,path,link) + +end subroutine results_setLink + + +!-------------------------------------------------------------------------------------------------- +!> @brief remove link to an object +!-------------------------------------------------------------------------------------------------- +subroutine results_removeLink(link) + use hdf5 + + implicit none + 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)//')') + +end subroutine results_removeLink + + +!-------------------------------------------------------------------------------------------------- +!> @brief stores a vector dataset in a group +!-------------------------------------------------------------------------------------------------- +subroutine results_writeVectorDataset(group,dataset,label,SIunit) + + implicit none + character(len=*), intent(in) :: SIunit,label,group + 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) + +end subroutine results_writeVectorDataset + + +!-------------------------------------------------------------------------------------------------- +!> @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 + + 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 + +!-------------------------------------------------------------------------------------------------- +! 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') + +!-------------------------------------------------------------------------------------------------- +! 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 +#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') +#endif + +!-------------------------------------------------------------------------------------------------- +! 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') + + 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') + +!-------------------------------------------------------------------------------------------------- +! 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 + +!-------------------------------------------------------------------------------------------------- +!> @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 + + 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 + + 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 + + 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 + +!-------------------------------------------------------------------------------------------------- +! 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(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 +#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') +#endif + +!-------------------------------------------------------------------------------------------------- +! 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') + + 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') + +!-------------------------------------------------------------------------------------------------- +!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 + +!-------------------------------------------------------------------------------------------------- +!> @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 + + +end module results diff --git a/src/spectral_mech_Basic.f90 b/src/spectral_mech_Basic.f90 index d6f353c91..003c9820d 100644 --- a/src/spectral_mech_Basic.f90 +++ b/src/spectral_mech_Basic.f90 @@ -80,6 +80,7 @@ subroutine basic_init #endif use IO, only: & IO_intOut, & + IO_error, & IO_read_realFile, & IO_timeStamp use debug, only: & @@ -173,7 +174,11 @@ subroutine basic_init call IO_read_realFile(777,'F_aimDot',trim(getSolverJobName()),size(F_aimDot)) read (777,rec=1) F_aimDot; close (777) F_aim = reshape(sum(sum(sum(F,dim=4),dim=3),dim=2) * wgt, [3,3]) ! average of F + call MPI_Allreduce(MPI_IN_PLACE,F_aim,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) + if(ierr /=0_pInt) call IO_error(894_pInt, ext_msg='F_aim') F_aim_lastInc = sum(sum(sum(F_lastInc,dim=5),dim=4),dim=3) * wgt ! average of F_lastInc + call MPI_Allreduce(MPI_IN_PLACE,F_aim_lastInc,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) + if(ierr /=0_pInt) call IO_error(894_pInt, ext_msg='F_aim_lastInc') elseif (restartInc == 0_pInt) then restart F_lastInc = spread(spread(spread(math_I3,3,grid(1)),4,grid(2)),5,grid3) ! initialize to identity F = reshape(F_lastInc,[9,grid(1),grid(2),grid3]) diff --git a/src/spectral_mech_Polarisation.f90 b/src/spectral_mech_Polarisation.f90 index 04f51cb35..b1da2a3f0 100644 --- a/src/spectral_mech_Polarisation.f90 +++ b/src/spectral_mech_Polarisation.f90 @@ -78,7 +78,6 @@ contains !-------------------------------------------------------------------------------------------------- !> @brief allocates all necessary fields and fills them with data, potentially from restart info -!> @todo use sourced allocation, e.g. allocate(Fdot,source = F_lastInc) !-------------------------------------------------------------------------------------------------- subroutine Polarisation_init #if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800 @@ -88,6 +87,7 @@ subroutine Polarisation_init #endif use IO, only: & IO_intOut, & + IO_error, & IO_read_realFile, & IO_timeStamp use debug, only: & @@ -191,7 +191,11 @@ subroutine Polarisation_init call IO_read_realFile(777,'F_aimDot',trim(getSolverJobName()),size(F_aimDot)) read (777,rec=1) F_aimDot; close (777) F_aim = reshape(sum(sum(sum(F,dim=4),dim=3),dim=2) * wgt, [3,3]) ! average of F + call MPI_Allreduce(MPI_IN_PLACE,F_aim,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) + if(ierr /=0_pInt) call IO_error(894_pInt, ext_msg='F_aim') F_aim_lastInc = sum(sum(sum(F_lastInc,dim=5),dim=4),dim=3) * wgt ! average of F_lastInc + call MPI_Allreduce(MPI_IN_PLACE,F_aim_lastInc,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) + if(ierr /=0_pInt) call IO_error(894_pInt, ext_msg='F_aim_lastInc') elseif (restartInc == 0_pInt) then restart F_lastInc = spread(spread(spread(math_I3,3,grid(1)),4,grid(2)),5,grid3) ! initialize to identity F = reshape(F_lastInc,[9,grid(1),grid(2),grid3])