diff --git a/src/CPFEM2.f90 b/src/CPFEM2.f90 index 5a500875d..b1e03659b 100644 --- a/src/CPFEM2.f90 +++ b/src/CPFEM2.f90 @@ -74,9 +74,22 @@ end subroutine CPFEM_initAll !-------------------------------------------------------------------------------------------------- subroutine CPFEM_init + integer(HID_T) :: fileHandle + character(len=pStringLen) :: fileName + + print'(/,a)', ' <<<+- CPFEM init -+>>>'; flush(IO_STDOUT) - if (interface_restartInc > 0) call crystallite_restartRead + + if (interface_restartInc > 0) then + print'(/,a,i0,a)', ' reading restart information of increment from file'; flush(IO_STDOUT) + write(fileName,'(a,i0,a)') trim(getSolverJobName())//'_',worldrank,'.hdf5' + fileHandle = HDF5_openFile(fileName) + + call constitutive_restartRead(fileHandle) + + call HDF5_closeFile(fileHandle) + endif end subroutine CPFEM_init @@ -85,8 +98,19 @@ end subroutine CPFEM_init !> @brief Write restart information. !-------------------------------------------------------------------------------------------------- subroutine CPFEM_restartWrite + + integer(HID_T) :: fileHandle + character(len=pStringLen) :: fileName + - call crystallite_restartWrite + print*, ' writing field and constitutive data required for restart to file';flush(IO_STDOUT) + + write(fileName,'(a,i0,a)') trim(getSolverJobName())//'_',worldrank,'.hdf5' + fileHandle = HDF5_openFile(fileName,'a') + + call constitutive_restartWrite(fileHandle) + + call HDF5_closeFile(fileHandle) end subroutine CPFEM_restartWrite diff --git a/src/constitutive.f90 b/src/constitutive.f90 index 261e9e304..67e8b33c8 100644 --- a/src/constitutive.f90 +++ b/src/constitutive.f90 @@ -184,6 +184,15 @@ module constitutive includeL end subroutine mech_restore + module function constitutive_mech_dPdF(dt,co,ip,el) result(dPdF) + real(pReal), intent(in) :: dt + integer, intent(in) :: & + co, & !< counter in constituent loop + ip, & !< counter in integration point loop + el !< counter in element loop + real(pReal), dimension(3,3,3,3) :: dPdF + end function constitutive_mech_dPdF + ! == cleaned:end =================================================================================== module function crystallite_stress(dt,co,ip,el) result(converged_) @@ -384,16 +393,16 @@ module constitutive converged, & crystallite_init, & crystallite_stress, & - crystallite_stressTangent, & + constitutive_mech_dPdF, & crystallite_orientations, & crystallite_push33ToRef, & - crystallite_restartWrite, & + constitutive_restartWrite, & + constitutive_restartRead, & integrateSourceState, & constitutive_mech_setF, & constitutive_mech_getLp, & constitutive_mech_getF, & constitutive_mech_getS, & - crystallite_restartRead, & constitutive_initializeRestorationPoints, & constitutive_windForward, & PLASTICITY_UNDEFINED_ID, & @@ -975,134 +984,6 @@ subroutine constitutive_windForward(ip,el) end subroutine constitutive_windForward -!-------------------------------------------------------------------------------------------------- -!> @brief Calculate tangent (dPdF). -!-------------------------------------------------------------------------------------------------- -function crystallite_stressTangent(dt,co,ip,el) result(dPdF) - - real(pReal), intent(in) :: dt - integer, intent(in) :: & - co, & !< counter in constituent loop - ip, & !< counter in integration point loop - el !< counter in element loop - real(pReal), dimension(3,3,3,3) :: dPdF - - integer :: & - o, & - p, ph, me - real(pReal), dimension(3,3) :: devNull, & - invSubFp0,invSubFi0,invFp,invFi, & - temp_33_1, temp_33_2, temp_33_3 - real(pReal), dimension(3,3,3,3) :: dSdFe, & - dSdF, & - dSdFi, & - dLidS, & ! tangent in lattice configuration - dLidFi, & - dLpdS, & - dLpdFi, & - dFidS, & - dFpinvdF, & - rhs_3333, & - lhs_3333, & - temp_3333 - real(pReal), dimension(9,9):: temp_99 - logical :: error - - - ph = material_phaseAt(co,el) - me = material_phaseMemberAt(co,ip,el) - - call constitutive_hooke_SandItsTangents(devNull,dSdFe,dSdFi, & - constitutive_mech_Fe(ph)%data(1:3,1:3,me), & - constitutive_mech_Fi(ph)%data(1:3,1:3,me),co,ip,el) - call constitutive_LiAndItsTangents(devNull,dLidS,dLidFi, & - constitutive_mech_S(ph)%data(1:3,1:3,me), & - constitutive_mech_Fi(ph)%data(1:3,1:3,me), & - co,ip,el) - - invFp = math_inv33(constitutive_mech_Fp(ph)%data(1:3,1:3,me)) - invFi = math_inv33(constitutive_mech_Fi(ph)%data(1:3,1:3,me)) - invSubFp0 = math_inv33(constitutive_mech_partitionedFp0(ph)%data(1:3,1:3,me)) - invSubFi0 = math_inv33(constitutive_mech_partitionedFi0(ph)%data(1:3,1:3,me)) - - if (sum(abs(dLidS)) < tol_math_check) then - dFidS = 0.0_pReal - else - lhs_3333 = 0.0_pReal; rhs_3333 = 0.0_pReal - do o=1,3; do p=1,3 - lhs_3333(1:3,1:3,o,p) = lhs_3333(1:3,1:3,o,p) & - + matmul(invSubFi0,dLidFi(1:3,1:3,o,p)) * dt - lhs_3333(1:3,o,1:3,p) = lhs_3333(1:3,o,1:3,p) & - + invFi*invFi(p,o) - rhs_3333(1:3,1:3,o,p) = rhs_3333(1:3,1:3,o,p) & - - matmul(invSubFi0,dLidS(1:3,1:3,o,p)) * dt - enddo; enddo - call math_invert(temp_99,error,math_3333to99(lhs_3333)) - if (error) then - call IO_warning(warning_ID=600,el=el,ip=ip,g=co, & - ext_msg='inversion error in analytic tangent calculation') - dFidS = 0.0_pReal - else - dFidS = math_mul3333xx3333(math_99to3333(temp_99),rhs_3333) - endif - dLidS = math_mul3333xx3333(dLidFi,dFidS) + dLidS - endif - - call constitutive_plastic_LpAndItsTangents(devNull,dLpdS,dLpdFi, & - constitutive_mech_S(ph)%data(1:3,1:3,me), & - constitutive_mech_Fi(ph)%data(1:3,1:3,me),co,ip,el) - dLpdS = math_mul3333xx3333(dLpdFi,dFidS) + dLpdS - -!-------------------------------------------------------------------------------------------------- -! calculate dSdF - temp_33_1 = transpose(matmul(invFp,invFi)) - temp_33_2 = matmul(constitutive_mech_F(ph)%data(1:3,1:3,me),invSubFp0) - temp_33_3 = matmul(matmul(constitutive_mech_F(ph)%data(1:3,1:3,me),invFp), invSubFi0) - - do o=1,3; do p=1,3 - rhs_3333(p,o,1:3,1:3) = matmul(dSdFe(p,o,1:3,1:3),temp_33_1) - temp_3333(1:3,1:3,p,o) = matmul(matmul(temp_33_2,dLpdS(1:3,1:3,p,o)), invFi) & - + matmul(temp_33_3,dLidS(1:3,1:3,p,o)) - enddo; enddo - lhs_3333 = math_mul3333xx3333(dSdFe,temp_3333) * dt & - + math_mul3333xx3333(dSdFi,dFidS) - - call math_invert(temp_99,error,math_eye(9)+math_3333to99(lhs_3333)) - if (error) then - call IO_warning(warning_ID=600,el=el,ip=ip,g=co, & - ext_msg='inversion error in analytic tangent calculation') - dSdF = rhs_3333 - else - dSdF = math_mul3333xx3333(math_99to3333(temp_99),rhs_3333) - endif - -!-------------------------------------------------------------------------------------------------- -! calculate dFpinvdF - temp_3333 = math_mul3333xx3333(dLpdS,dSdF) - do o=1,3; do p=1,3 - dFpinvdF(1:3,1:3,p,o) = - matmul(invSubFp0, matmul(temp_3333(1:3,1:3,p,o),invFi)) * dt - enddo; enddo - -!-------------------------------------------------------------------------------------------------- -! assemble dPdF - temp_33_1 = matmul(constitutive_mech_S(ph)%data(1:3,1:3,me),transpose(invFp)) - temp_33_2 = matmul(constitutive_mech_F(ph)%data(1:3,1:3,me),invFp) - temp_33_3 = matmul(temp_33_2,constitutive_mech_S(ph)%data(1:3,1:3,me)) - - dPdF = 0.0_pReal - do p=1,3 - dPdF(p,1:3,p,1:3) = transpose(matmul(invFp,temp_33_1)) - enddo - do o=1,3; do p=1,3 - dPdF(1:3,1:3,p,o) = dPdF(1:3,1:3,p,o) & - + matmul(matmul(constitutive_mech_F(ph)%data(1:3,1:3,me),dFpinvdF(1:3,1:3,p,o)),temp_33_1) & - + matmul(matmul(temp_33_2,dSdF(1:3,1:3,p,o)),transpose(invFp)) & - + matmul(temp_33_3,transpose(dFpinvdF(1:3,1:3,p,o))) - enddo; enddo - -end function crystallite_stressTangent - - !-------------------------------------------------------------------------------------------------- !> @brief calculates orientations !-------------------------------------------------------------------------------------------------- @@ -1273,16 +1154,12 @@ end function converged !> @brief Write current restart information (Field and constitutive data) to file. ! ToDo: Merge data into one file for MPI, move state to constitutive and homogenization, respectively !-------------------------------------------------------------------------------------------------- -subroutine crystallite_restartWrite +subroutine constitutive_restartWrite(fileHandle) + integer(HID_T), intent(in) :: fileHandle integer :: ph - integer(HID_T) :: fileHandle, groupHandle - character(len=pStringLen) :: fileName, datasetName - - print*, ' writing field and constitutive data required for restart to file';flush(IO_STDOUT) - - write(fileName,'(a,i0,a)') trim(getSolverJobName())//'_',worldrank,'.hdf5' - fileHandle = HDF5_openFile(fileName,'a') + integer(HID_T) :: groupHandle + character(len=pStringLen) :: datasetName groupHandle = HDF5_addGroup(fileHandle,'phase') do ph = 1,size(material_name_phase) @@ -1310,25 +1187,21 @@ subroutine crystallite_restartWrite enddo call HDF5_closeGroup(groupHandle) - call HDF5_closeFile(fileHandle) -end subroutine crystallite_restartWrite +end subroutine constitutive_restartWrite !-------------------------------------------------------------------------------------------------- !> @brief Read data for restart ! ToDo: Merge data into one file for MPI, move state to constitutive and homogenization, respectively !-------------------------------------------------------------------------------------------------- -subroutine crystallite_restartRead +subroutine constitutive_restartRead(fileHandle) + integer(HID_T), intent(in) :: fileHandle integer :: ph - integer(HID_T) :: fileHandle, groupHandle - character(len=pStringLen) :: fileName, datasetName + integer(HID_T) :: groupHandle + character(len=pStringLen) ::datasetName - print'(/,a,i0,a)', ' reading restart information of increment from file' - - write(fileName,'(a,i0,a)') trim(getSolverJobName())//'_',worldrank,'.hdf5' - fileHandle = HDF5_openFile(fileName) groupHandle = HDF5_openGroup(fileHandle,'phase') do ph = 1,size(material_name_phase) @@ -1356,9 +1229,7 @@ subroutine crystallite_restartRead enddo call HDF5_closeGroup(groupHandle) - call HDF5_closeFile(fileHandle) - -end subroutine crystallite_restartRead +end subroutine constitutive_restartRead ! getter for non-mech (e.g. thermal) diff --git a/src/constitutive_mech.f90 b/src/constitutive_mech.f90 index fa9a5eda6..6392fb0ee 100644 --- a/src/constitutive_mech.f90 +++ b/src/constitutive_mech.f90 @@ -1688,5 +1688,132 @@ module subroutine mech_restore(ip,el,includeL) end subroutine mech_restore +!-------------------------------------------------------------------------------------------------- +!> @brief Calculate tangent (dPdF). +!-------------------------------------------------------------------------------------------------- +module function constitutive_mech_dPdF(dt,co,ip,el) result(dPdF) + + real(pReal), intent(in) :: dt + integer, intent(in) :: & + co, & !< counter in constituent loop + ip, & !< counter in integration point loop + el !< counter in element loop + real(pReal), dimension(3,3,3,3) :: dPdF + + integer :: & + o, & + p, ph, me + real(pReal), dimension(3,3) :: devNull, & + invSubFp0,invSubFi0,invFp,invFi, & + temp_33_1, temp_33_2, temp_33_3 + real(pReal), dimension(3,3,3,3) :: dSdFe, & + dSdF, & + dSdFi, & + dLidS, & ! tangent in lattice configuration + dLidFi, & + dLpdS, & + dLpdFi, & + dFidS, & + dFpinvdF, & + rhs_3333, & + lhs_3333, & + temp_3333 + real(pReal), dimension(9,9):: temp_99 + logical :: error + + + ph = material_phaseAt(co,el) + me = material_phaseMemberAt(co,ip,el) + + call constitutive_hooke_SandItsTangents(devNull,dSdFe,dSdFi, & + constitutive_mech_Fe(ph)%data(1:3,1:3,me), & + constitutive_mech_Fi(ph)%data(1:3,1:3,me),co,ip,el) + call constitutive_LiAndItsTangents(devNull,dLidS,dLidFi, & + constitutive_mech_S(ph)%data(1:3,1:3,me), & + constitutive_mech_Fi(ph)%data(1:3,1:3,me), & + co,ip,el) + + invFp = math_inv33(constitutive_mech_Fp(ph)%data(1:3,1:3,me)) + invFi = math_inv33(constitutive_mech_Fi(ph)%data(1:3,1:3,me)) + invSubFp0 = math_inv33(constitutive_mech_partitionedFp0(ph)%data(1:3,1:3,me)) + invSubFi0 = math_inv33(constitutive_mech_partitionedFi0(ph)%data(1:3,1:3,me)) + + if (sum(abs(dLidS)) < tol_math_check) then + dFidS = 0.0_pReal + else + lhs_3333 = 0.0_pReal; rhs_3333 = 0.0_pReal + do o=1,3; do p=1,3 + lhs_3333(1:3,1:3,o,p) = lhs_3333(1:3,1:3,o,p) & + + matmul(invSubFi0,dLidFi(1:3,1:3,o,p)) * dt + lhs_3333(1:3,o,1:3,p) = lhs_3333(1:3,o,1:3,p) & + + invFi*invFi(p,o) + rhs_3333(1:3,1:3,o,p) = rhs_3333(1:3,1:3,o,p) & + - matmul(invSubFi0,dLidS(1:3,1:3,o,p)) * dt + enddo; enddo + call math_invert(temp_99,error,math_3333to99(lhs_3333)) + if (error) then + call IO_warning(warning_ID=600,el=el,ip=ip,g=co, & + ext_msg='inversion error in analytic tangent calculation') + dFidS = 0.0_pReal + else + dFidS = math_mul3333xx3333(math_99to3333(temp_99),rhs_3333) + endif + dLidS = math_mul3333xx3333(dLidFi,dFidS) + dLidS + endif + + call constitutive_plastic_LpAndItsTangents(devNull,dLpdS,dLpdFi, & + constitutive_mech_S(ph)%data(1:3,1:3,me), & + constitutive_mech_Fi(ph)%data(1:3,1:3,me),co,ip,el) + dLpdS = math_mul3333xx3333(dLpdFi,dFidS) + dLpdS + +!-------------------------------------------------------------------------------------------------- +! calculate dSdF + temp_33_1 = transpose(matmul(invFp,invFi)) + temp_33_2 = matmul(constitutive_mech_F(ph)%data(1:3,1:3,me),invSubFp0) + temp_33_3 = matmul(matmul(constitutive_mech_F(ph)%data(1:3,1:3,me),invFp), invSubFi0) + + do o=1,3; do p=1,3 + rhs_3333(p,o,1:3,1:3) = matmul(dSdFe(p,o,1:3,1:3),temp_33_1) + temp_3333(1:3,1:3,p,o) = matmul(matmul(temp_33_2,dLpdS(1:3,1:3,p,o)), invFi) & + + matmul(temp_33_3,dLidS(1:3,1:3,p,o)) + enddo; enddo + lhs_3333 = math_mul3333xx3333(dSdFe,temp_3333) * dt & + + math_mul3333xx3333(dSdFi,dFidS) + + call math_invert(temp_99,error,math_eye(9)+math_3333to99(lhs_3333)) + if (error) then + call IO_warning(warning_ID=600,el=el,ip=ip,g=co, & + ext_msg='inversion error in analytic tangent calculation') + dSdF = rhs_3333 + else + dSdF = math_mul3333xx3333(math_99to3333(temp_99),rhs_3333) + endif + +!-------------------------------------------------------------------------------------------------- +! calculate dFpinvdF + temp_3333 = math_mul3333xx3333(dLpdS,dSdF) + do o=1,3; do p=1,3 + dFpinvdF(1:3,1:3,p,o) = - matmul(invSubFp0, matmul(temp_3333(1:3,1:3,p,o),invFi)) * dt + enddo; enddo + +!-------------------------------------------------------------------------------------------------- +! assemble dPdF + temp_33_1 = matmul(constitutive_mech_S(ph)%data(1:3,1:3,me),transpose(invFp)) + temp_33_2 = matmul(constitutive_mech_F(ph)%data(1:3,1:3,me),invFp) + temp_33_3 = matmul(temp_33_2,constitutive_mech_S(ph)%data(1:3,1:3,me)) + + dPdF = 0.0_pReal + do p=1,3 + dPdF(p,1:3,p,1:3) = transpose(matmul(invFp,temp_33_1)) + enddo + do o=1,3; do p=1,3 + dPdF(1:3,1:3,p,o) = dPdF(1:3,1:3,p,o) & + + matmul(matmul(constitutive_mech_F(ph)%data(1:3,1:3,me),dFpinvdF(1:3,1:3,p,o)),temp_33_1) & + + matmul(matmul(temp_33_2,dSdF(1:3,1:3,p,o)),transpose(invFp)) & + + matmul(temp_33_3,transpose(dFpinvdF(1:3,1:3,p,o))) + enddo; enddo + +end function constitutive_mech_dPdF + end submodule constitutive_mech diff --git a/src/homogenization_mech.f90 b/src/homogenization_mech.f90 index 4a9e1856f..1d0942f3e 100644 --- a/src/homogenization_mech.f90 +++ b/src/homogenization_mech.f90 @@ -156,11 +156,11 @@ module subroutine mech_homogenize(dt,ip,el) case (HOMOGENIZATION_NONE_ID) chosenHomogenization homogenization_P(1:3,1:3,ce) = crystallite_P(1:3,1:3,1,ip,el) - homogenization_dPdF(1:3,1:3,1:3,1:3,ce) = crystallite_stressTangent(dt,1,ip,el) + homogenization_dPdF(1:3,1:3,1:3,1:3,ce) = constitutive_mech_dPdF(dt,1,ip,el) case (HOMOGENIZATION_ISOSTRAIN_ID) chosenHomogenization do co = 1, homogenization_Nconstituents(material_homogenizationAt(el)) - dPdFs(:,:,:,:,co) = crystallite_stressTangent(dt,co,ip,el) + dPdFs(:,:,:,:,co) = constitutive_mech_dPdF(dt,co,ip,el) enddo call mech_isostrain_averageStressAndItsTangent(& homogenization_P(1:3,1:3,ce), & @@ -171,7 +171,7 @@ module subroutine mech_homogenize(dt,ip,el) case (HOMOGENIZATION_RGC_ID) chosenHomogenization do co = 1, homogenization_Nconstituents(material_homogenizationAt(el)) - dPdFs(:,:,:,:,co) = crystallite_stressTangent(dt,co,ip,el) + dPdFs(:,:,:,:,co) = constitutive_mech_dPdF(dt,co,ip,el) enddo call mech_RGC_averageStressAndItsTangent(& homogenization_P(1:3,1:3,ce), & @@ -207,7 +207,7 @@ module function mech_updateState(subdt,subF,ip,el) result(doneAndHappy) if (homogenization_type(material_homogenizationAt(el)) == HOMOGENIZATION_RGC_ID) then do co = 1, homogenization_Nconstituents(material_homogenizationAt(el)) - dPdFs(:,:,:,:,co) = crystallite_stressTangent(subdt,co,ip,el) + dPdFs(:,:,:,:,co) = constitutive_mech_dPdF(subdt,co,ip,el) Fs(:,:,co) = constitutive_mech_getF(co,ip,el) enddo doneAndHappy = &