From 3aea8b39e9531fcf83c415efc24065eb5b88d672 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 12 Mar 2014 07:33:51 +0000 Subject: [PATCH] added some HDF5 functionality (needs to be activated with preprocessor makro) --- code/CPFEM.f90 | 14 +- code/IO.f90 | 277 +++++++++++++++++++++++++++++++++++++++- code/constitutive.f90 | 25 +++- code/homogenization.f90 | 2 +- code/material.f90 | 2 +- 5 files changed, 306 insertions(+), 14 deletions(-) diff --git a/code/CPFEM.f90 b/code/CPFEM.f90 index 243628e01..a636ce33e 100644 --- a/code/CPFEM.f90 +++ b/code/CPFEM.f90 @@ -544,8 +544,8 @@ subroutine CPFEM_general(mode, parallelExecution, ffn, ffn1, temperature, dt, el if (.not. parallelExecution) then FEsolving_execElem(1) = elCP FEsolving_execElem(2) = elCP - if (.not. microstructure_elemhomo(mesh_element(4,elCP)) .or. & ! calculate unless homogeneous - (microstructure_elemhomo(mesh_element(4,elCP)) .and. ip == 1_pInt)) then ! and then only first ip + if (.not. microstructure_elemhomo(mesh_element(4,elCP)) .or. & ! calculate unless homogeneous + (microstructure_elemhomo(mesh_element(4,elCP)) .and. ip == 1_pInt)) then ! and then only first ip FEsolving_execIP(1,elCP) = ip FEsolving_execIP(2,elCP) = ip if (iand(debug_level(debug_CPFEM), debug_levelExtensive) /= 0_pInt) then @@ -553,8 +553,8 @@ subroutine CPFEM_general(mode, parallelExecution, ffn, ffn1, temperature, dt, el write(6,'(a,i8,1x,i2)') '<< CPFEM >> calculation for elFE ip ',elCP,ip !$OMP END CRITICAL (write2out) endif - call materialpoint_stressAndItsTangent(updateJaco, dt) ! calculate stress and its tangent - call materialpoint_postResults() ! post results + call materialpoint_stressAndItsTangent(updateJaco, dt) ! calculate stress and its tangent + call materialpoint_postResults() endif !* parallel computation and calulation not yet done @@ -566,8 +566,8 @@ subroutine CPFEM_general(mode, parallelExecution, ffn, ffn1, temperature, dt, el ' to ',FEsolving_execElem(2) !$OMP END CRITICAL (write2out) endif - call materialpoint_stressAndItsTangent(updateJaco, dt) ! calculate stress and its tangent (parallel execution inside) - call materialpoint_postResults() ! post results + call materialpoint_stressAndItsTangent(updateJaco, dt) ! calculate stress and its tangent (parallel execution inside) + call materialpoint_postResults() CPFEM_calc_done = .true. endif @@ -579,7 +579,7 @@ subroutine CPFEM_general(mode, parallelExecution, ffn, ffn1, temperature, dt, el CPFEM_cs(1:6,ip,elCP) = rnd * CPFEM_odd_stress CPFEM_dcsde(1:6,1:6,ip,elCP) = CPFEM_odd_jacobian * math_identity2nd(6) else - if (microstructure_elemhomo(mesh_element(4,elCP)) .and. ip > 1_pInt) then ! me homogenous? --> copy from first ip + if (microstructure_elemhomo(mesh_element(4,elCP)) .and. ip > 1_pInt) then ! me homogenous? --> copy from first ip materialpoint_P(1:3,1:3,ip,elCP) = materialpoint_P(1:3,1:3,1,elCP) materialpoint_F(1:3,1:3,ip,elCP) = materialpoint_F(1:3,1:3,1,elCP) materialpoint_dPdF(1:3,1:3,1:3,1:3,ip,elCP) = materialpoint_dPdF(1:3,1:3,1:3,1:3,1,elCP) diff --git a/code/IO.f90 b/code/IO.f90 index 1c18b4312..c4c10268c 100644 --- a/code/IO.f90 +++ b/code/IO.f90 @@ -25,7 +25,11 @@ !> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH !> @brief input/output functions, partly depending on chosen solver !-------------------------------------------------------------------------------------------------- -module IO +module IO +#ifdef HDF + use hdf5, only: & + HID_T +#endif use prec, only: & pInt, & pReal @@ -34,6 +38,10 @@ module IO private character(len=5), parameter, public :: & IO_EOF = '#EOF#' !< end of file string +#ifdef HDF + integer(HID_T), public, protected :: resultsFile, tempResults +#endif + public :: & IO_init, & IO_read, & @@ -90,7 +98,10 @@ module IO private :: & abaqus_assembleInputFile #endif - +#ifdef HDF +public:: HDF5_mappingConstitutive, & + HDF5_closeJobFile +#endif contains @@ -105,6 +116,10 @@ subroutine IO_init write(6,'(a15,a)') ' Current time: ',IO_timeStamp() #include "compilation_info.f90" +#ifdef HDF + call HDF5_createJobFile +#endif + end subroutine IO_init @@ -1891,4 +1906,262 @@ recursive function abaqus_assembleInputFile(unit1,unit2) result(createSuccess) end function abaqus_assembleInputFile #endif +#ifdef HDF + + +!-------------------------------------------------------------------------------------------------- +!> @brief creates and initializes HDF5 output file +!-------------------------------------------------------------------------------------------------- +subroutine HDF5_createJobFile + use hdf5 + use DAMASK_interface, only: & + getSolverWorkingDirectoryName, & + getSolverJobName + + implicit none + integer :: hdferr + INTEGER(SIZE_T) :: typeSize + character(len=1024) :: path + +!-------------------------------------------------------------------------------------------------- +! initialize HDF5 library and check if integer and float type size match + call h5open_f(hdferr) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_createJobFile') + call h5tget_size_f(H5T_NATIVE_INTEGER,typeSize, hdferr) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_createJobFile: h5tget_size_f (int)') + if (int(pInt,SIZE_T)/=typeSize) call IO_error(0_pInt,ext_msg='pInt does not match H5T_NATIVE_INTEGER') + call h5tget_size_f(H5T_NATIVE_DOUBLE,typeSize, hdferr) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_createJobFile: h5tget_size_f (double)') + if (int(pReal,SIZE_T)/=typeSize) call IO_error(0_pInt,ext_msg='pReal does not match H5T_NATIVE_DOUBLE') + +!-------------------------------------------------------------------------------------------------- +! open file + path = trim(getSolverWorkingDirectoryName())//trim(getSolverJobName())//'.'//'DAMASKout' + call h5fcreate_f(path,H5F_ACC_TRUNC_F,resultsFile,hdferr) + if (hdferr < 0) call IO_error(100_pInt,ext_msg=path) + call HDF5_addStringAttribute(resultsFile,'createdBy','$Id$') + +!-------------------------------------------------------------------------------------------------- +! create mapping group in file + call HDF5_closeGroup(HDF5_addGroup("mapping")) + +end subroutine HDF5_createJobFile + + +!-------------------------------------------------------------------------------------------------- +!> @brief creates and initializes HDF5 output file +!-------------------------------------------------------------------------------------------------- +subroutine HDF5_closeJobFile() + use hdf5 + + implicit none + integer :: hdferr + call h5fclose_f(resultsFile,hdferr) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_closeJobFile') + +end subroutine HDF5_closeJobFile + + +!-------------------------------------------------------------------------------------------------- +!> @brief adds a new group to the results file +!-------------------------------------------------------------------------------------------------- +integer(HID_T) function HDF5_addGroup(path) + use hdf5 + + implicit none + character(len=*), intent(in) :: path + integer :: hdferr + + call h5gcreate_f(resultsFile, trim(path), HDF5_addGroup, hdferr) + if (hdferr < 0) call IO_error(1_pInt,ext_msg = 'HDF5_addGroup failed ('//trim(path)//' )') + +end function HDF5_addGroup + + +!-------------------------------------------------------------------------------------------------- +!> @brief adds a new group to the results file +!-------------------------------------------------------------------------------------------------- +integer(HID_T) function HDF5_openGroup(path) + use hdf5 + + implicit none + character(len=*), intent(in) :: path + integer :: hdferr + + call h5gopen_f(resultsFile, trim(path), HDF5_openGroup, hdferr) + if (hdferr < 0) call IO_error(1_pInt,ext_msg = 'HDF5_openGroup failed ('//trim(path)//' )') + +end function HDF5_openGroup + + +!-------------------------------------------------------------------------------------------------- +!> @brief closes a group +!-------------------------------------------------------------------------------------------------- +subroutine HDF5_closeGroup(ID) + use hdf5 + + implicit none + integer(HID_T), intent(in) :: ID + integer :: hdferr + + call h5gclose_f(ID, hdferr) + if (hdferr < 0) call IO_error(1_pInt,ext_msg = 'HDF5_closeGroup') + +end subroutine HDF5_closeGroup + + +!-------------------------------------------------------------------------------------------------- +!> @brief adds a new group to the results file +!-------------------------------------------------------------------------------------------------- +subroutine HDF5_addStringAttribute(entity,attrLabel,attrValue) + use hdf5 + + implicit none + integer(HID_T), intent(in) :: entity + character(len=*), intent(in) :: attrLabel, attrValue + integer :: hdferr + integer(HID_T) :: attr_id, space_id, type_id + + call h5screate_f(H5S_SCALAR_F,space_id,hdferr) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_addStringAttribute') + call h5tcopy_f(H5T_NATIVE_CHARACTER, type_id, hdferr) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_addStringAttribute') + 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') + 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') + 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') + call h5aclose_f(attr_id,hdferr) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_addStringAttribute') + call h5sclose_f(space_id,hdferr) + if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_addStringAttribute') + +end subroutine HDF5_addStringAttribute + + +!-------------------------------------------------------------------------------------------------- +!> @brief adds the unique mapping from spatial position and constituent ID to results +!-------------------------------------------------------------------------------------------------- +subroutine HDF5_mappingConstitutive(mapping) + use hdf5 + + implicit none + integer(pInt), intent(in), dimension(:,:,:) :: mapping + + integer :: hdf5err, NmatPoints,Nconstituents + integer(HID_T) :: mapping_id, dtype_id, dset_id, space_id,instance_id,position_id + + Nconstituents=size(mapping,1) + NmatPoints=size(mapping,2) + mapping_ID = HDF5_openGroup("mapping") + +!-------------------------------------------------------------------------------------------------- +! create dataspace + call h5screate_simple_f(2, int([Nconstituents,NmatPoints],HSIZE_T), space_id, hdf5err, & + int([Nconstituents,NmatPoints],HSIZE_T)) + if (hdf5err < 0) call IO_error(1_pInt,ext_msg='IO_addAttribute') + +!-------------------------------------------------------------------------------------------------- +! compound type + call h5tcreate_f(H5T_COMPOUND_F, 6_SIZE_T, dtype_id, hdf5err) + if (hdf5err < 0) call IO_error(1_pInt,ext_msg='IO_addAttribute: h5tcreate_f dtype_id') + call h5tinsert_f(dtype_id, "Constitutive Instance", 0_SIZE_T, H5T_STD_U16LE, hdf5err) + if (hdf5err < 0) call IO_error(1_pInt,ext_msg='IO_addAttribute: h5tinsert_f 0') + call h5tinsert_f(dtype_id, "Position in Instance Results", 2_SIZE_T, H5T_STD_U32LE, hdf5err) + if (hdf5err < 0) call IO_error(1_pInt,ext_msg='IO_addAttribute: h5tinsert_f 2') + +!-------------------------------------------------------------------------------------------------- +! create Dataset + call h5dcreate_f(mapping_id, "Constitutive", dtype_id, space_id, dset_id, hdf5err) + if (hdf5err < 0) call IO_error(1_pInt,ext_msg='IO_addAttribute') + +!-------------------------------------------------------------------------------------------------- +! Create memory types (one compound datatype for each member) + call h5tcreate_f(H5T_COMPOUND_F, int(pInt,SIZE_T), instance_id, hdf5err) + if (hdf5err < 0) call IO_error(1_pInt,ext_msg='IO_addAttribute: h5tcreate_f instance_id') + call h5tinsert_f(instance_id, "Constitutive Instance", 0_SIZE_T, H5T_NATIVE_INTEGER, hdf5err) + if (hdf5err < 0) call IO_error(1_pInt,ext_msg='IO_addAttribute: h5tinsert_f instance_id') + call h5tcreate_f(H5T_COMPOUND_F, int(pInt,SIZE_T), position_id, hdf5err) + if (hdf5err < 0) call IO_error(1_pInt,ext_msg='IO_addAttribute: h5tcreate_f position_id') + call h5tinsert_f(position_id, "Position in Instance Results", 0_SIZE_T, H5T_NATIVE_INTEGER, hdf5err) + if (hdf5err < 0) call IO_error(1_pInt,ext_msg='IO_addAttribute: h5tinsert_f position_id') +print*, mapping +print*, 'ddddddddd' +print*, mapping(1:Nconstituents,1:nmatpoints,1) +print*, mapping(1:Nconstituents,1:nmatpoints,2) + +!-------------------------------------------------------------------------------------------------- +! write data by fields in the datatype. Fields order is not important. + call h5dwrite_f(dset_id, position_id, mapping(1:Nconstituents,1:NmatPoints,1), & + int([Nconstituents, NmatPoints],HSIZE_T), hdf5err) + + call h5dwrite_f(dset_id, instance_id, mapping(1:Nconstituents,1:NmatPoints,2), & + int([Nconstituents, NmatPoints],HSIZE_T), hdf5err) + +!-------------------------------------------------------------------------------------------------- +!close types, dataspaces + call h5tclose_f(dtype_id, hdf5err) + call h5tclose_f(position_id, hdf5err) + call h5tclose_f(instance_id, hdf5err) + call h5dclose_f(dset_id, hdf5err) + call h5sclose_f(space_id, hdf5err) + call HDF5_closeGroup(mapping_ID) + +end subroutine HDF5_mappingConstitutive + + +subroutine HDF5_mappingHomogenization(Npoints) + use hdf5 + + implicit none + integer(pInt), intent(in) :: Npoints + integer(pInt) :: i + integer :: hdf5err + integer(HID_T) :: mapping_ID,dspace_id,dtype_id,dset_ID + INTEGER(HID_T) :: instance_id,position_id,elem_id,ip_id + + mapping_ID=HDF5_openGroup("mapping") + call h5screate_simple_f(1, [int(Npoints,HSIZE_T)], dspace_id, hdf5err) ! dataspace + +! compound type + CALL h5tcreate_f(H5T_COMPOUND_F, 11_SIZE_T, dtype_id, hdf5err) + CALL h5tinsert_f(dtype_id, "Homogenization Instance", 0_SIZE_T, H5T_STD_U16LE, hdf5err) + CALL h5tinsert_f(dtype_id, "Position in Instance", 2_SIZE_T, H5T_STD_U32LE, hdf5err) + CALL h5tinsert_f(dtype_id, "Element ID", 6_SIZE_T, H5T_STD_U32LE, hdf5err) + CALL h5tinsert_f(dtype_id, "Integration Point Number",10_SIZE_T, H5T_STD_U8LE, hdf5err) +! create Dataset + CALL h5dcreate_f(mapping_id, "Homogenization", dtype_id, dspace_id, dset_id, hdf5err) ! dataset + + ! Create memory types. We have to create a compound datatype + ! for each member we want to write. + ! + CALL h5tcreate_f(H5T_COMPOUND_F, int(pInt,SIZE_T), instance_id, hdf5err) + CALL h5tinsert_f(instance_id, "Homogenization Instance", 0_SIZE_T, H5T_NATIVE_INTEGER, hdf5err) + CALL h5tcreate_f(H5T_COMPOUND_F, int(pInt,SIZE_T), position_id, hdf5err) + CALL h5tinsert_f(position_id, "Position in Instance", 0_SIZE_T, H5T_NATIVE_INTEGER, hdf5err) + CALL h5tcreate_f(H5T_COMPOUND_F, int(pInt,SIZE_T), elem_id, hdf5err) + CALL h5tinsert_f(elem_id, "Element Number", 0_SIZE_T, H5T_NATIVE_INTEGER, hdf5err) + CALL h5tcreate_f(H5T_COMPOUND_F, int(pInt,SIZE_T), ip_id, hdf5err) + CALL h5tinsert_f(ip_id, "Integration Point Number",0_SIZE_T, H5T_NATIVE_INTEGER, hdf5err) + + ! Write data by fields in the datatype. Fields order is not important. + ! + CALL h5dwrite_f(dset_id, ip_id, spread(1_pInt,1,Npoints), [int(Npoints,HSIZE_T)], hdf5err) + CALL h5dwrite_f(dset_id, position_id, [(i,i=0_pInt,Npoints-1_pInt)], [int(Npoints,HSIZE_T)], hdf5err) + CALL h5dwrite_f(dset_id, instance_id, spread(1_pInt,1,Npoints), [int(Npoints,HSIZE_T)], hdf5err) + CALL h5dwrite_f(dset_id, elem_id, [(i,i=1_pInt,Npoints)], [int(Npoints,HSIZE_T)], hdf5err) + +!close + call h5tclose_f(ip_id, hdf5err) + call h5tclose_f(elem_id, hdf5err) + call h5tclose_f(position_id, hdf5err) + call h5tclose_f(instance_id, hdf5err) + call h5tclose_f(dtype_id, hdf5err) + call h5dclose_f(dset_id, hdf5err) + call h5sclose_f(dspace_id, hdf5err) + call HDF5_closeGroup(mapping_ID) + +end subroutine HDF5_mappingHomogenization +#endif end module IO diff --git a/code/constitutive.f90 b/code/constitutive.f90 index aa8928d32..9143c3b99 100644 --- a/code/constitutive.f90 +++ b/code/constitutive.f90 @@ -76,6 +76,14 @@ contains !> @brief allocates arrays pointing to array of the various constitutive modules !-------------------------------------------------------------------------------------------------- subroutine constitutive_init +#ifdef HDF + use hdf5, only: & + HID_T + use IO, only : & + HDF5_mappingConstitutive, & + HDF5_closeJobFile +#endif + use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment) use debug, only: & debug_level, & @@ -146,7 +154,12 @@ subroutine constitutive_init character(len=64), dimension(:,:), pointer :: thisOutput character(len=32) :: outputName !< name of output, intermediate fix until HDF5 output is ready logical :: knownPlasticity, nonlocalConstitutionPresent - +#ifdef HDF + integer(pInt), dimension(:,:,:), allocatable :: mapping + integer(pInt), dimension(:), allocatable :: InstancePosition + allocate(mapping(homogenization_maxngrains,mesh_ncpelems,2),source=0_pInt) + allocate(InstancePosition(material_nphase),source=0_pInt) +#endif nonlocalConstitutionPresent = .false. @@ -248,6 +261,10 @@ subroutine constitutive_init case default ! so far no output for elasticity end select instance = phase_plasticityInstance(material_phase(g,i,e)) +#ifdef HDF + InstancePosition(instance) = InstancePosition(instance)+1_pInt + mapping(g,e,1:2) = [instancePosition(instance),instance] +#endif select case(phase_plasticity(material_phase(g,i,e))) case (PLASTICITY_NONE_ID) allocate(constitutive_state0(g,i,e)%p(constitutive_none_sizeState(instance))) @@ -425,7 +442,10 @@ subroutine constitutive_init constitutive_state(g,i,e)%p = constitutive_state0(g,i,e)%p ! need to be defined for first call of constitutive_microstructure in crystallite_init endforall enddo - +#ifdef HDF + call HDF5_mappingConstitutive(mapping) + call HDF5_closeJobFile() +#endif !-------------------------------------------------------------------------------------------------- ! write out state size file call IO_write_jobIntFile(777,'sizeStateConst', size(constitutive_sizeState)) @@ -875,7 +895,6 @@ function constitutive_postResults(Tstar_v, Fe, temperature, ipc, ip, el) select case (phase_plasticity(material_phase(ipc,ip,el))) case (PLASTICITY_NONE_ID) - constitutive_postResults = 0.0_pReal case (PLASTICITY_TITANMOD_ID) constitutive_postResults = constitutive_titanmod_postResults(constitutive_state,ipc,ip,el) diff --git a/code/homogenization.f90 b/code/homogenization.f90 index 957570ff4..dd2f1450a 100644 --- a/code/homogenization.f90 +++ b/code/homogenization.f90 @@ -615,7 +615,7 @@ subroutine materialpoint_postResults grainLooping :do g = 1,myNgrains theSize = (1 + crystallite_sizePostResults(myCrystallite)) + (1 + constitutive_sizePostResults(g,i,e)) - materialpoint_results(thePos+1:thePos+theSize,i,e) = crystallite_postResults(g,i,e) ! tell crystallite results + materialpoint_results(thePos+1:thePos+theSize,i,e) = crystallite_postResults(g,i,e) ! tell crystallite results thePos = thePos + theSize enddo grainLooping enddo IpLooping diff --git a/code/material.f90 b/code/material.f90 index 6a332191e..ba3de0c0d 100644 --- a/code/material.f90 +++ b/code/material.f90 @@ -102,7 +102,7 @@ module material homogenization_typeInstance, & !< instance of particular type of each homogenization microstructure_crystallite !< crystallite setting ID of each microstructure - integer(pInt), dimension(:,:,:), allocatable, public:: & + integer(pInt), dimension(:,:,:), allocatable, public :: & material_phase !< phase (index) of each grain,IP,element integer(pInt), dimension(:,:,:), allocatable, public, protected :: & material_texture !< texture (index) of each grain,IP,element