diff --git a/src/constitutive.f90 b/src/constitutive.f90 index 808870059..fe09b3662 100644 --- a/src/constitutive.f90 +++ b/src/constitutive.f90 @@ -42,8 +42,7 @@ module constitutive KINEMATICS_SLIPPLANE_OPENING_ID, & KINEMATICS_THERMAL_EXPANSION_ID end enum - real(pReal), dimension(:,:,:), allocatable :: & - crystallite_subdt !< substepped time increment of each grain + type(rotation), dimension(:,:,:), allocatable :: & crystallite_orientation !< current orientation real(pReal), dimension(:,:,:,:,:), allocatable :: & @@ -874,7 +873,6 @@ subroutine crystallite_init crystallite_Fe,crystallite_Lp, & source = crystallite_F) - allocate(crystallite_subdt(cMax,iMax,eMax),source=0.0_pReal) allocate(crystallite_orientation(cMax,iMax,eMax)) @@ -1103,11 +1101,11 @@ function crystallite_stressTangent(co,ip,el) result(dPdF) 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) & - + crystallite_subdt(co,ip,el)*matmul(invSubFi0,dLidFi(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) & - - crystallite_subdt(co,ip,el)*matmul(invSubFi0,dLidS(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 @@ -1136,7 +1134,7 @@ function crystallite_stressTangent(co,ip,el) result(dPdF) 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 = crystallite_subdt(co,ip,el)*math_mul3333xx3333(dSdFe,temp_3333) & + 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)) @@ -1152,8 +1150,7 @@ function crystallite_stressTangent(co,ip,el) result(dPdF) ! calculate dFpinvdF temp_3333 = math_mul3333xx3333(dLpdS,dSdF) do o=1,3; do p=1,3 - dFpinvdF(1:3,1:3,p,o) = -crystallite_subdt(co,ip,el) & - * matmul(invSubFp0, matmul(temp_3333(1:3,1:3,p,o),invFi)) + dFpinvdF(1:3,1:3,p,o) = - matmul(invSubFp0, matmul(temp_3333(1:3,1:3,p,o),invFi)) * dt enddo; enddo !-------------------------------------------------------------------------------------------------- diff --git a/src/constitutive_mech.f90 b/src/constitutive_mech.f90 index 11ced6f40..7c819c480 100644 --- a/src/constitutive_mech.f90 +++ b/src/constitutive_mech.f90 @@ -1526,7 +1526,6 @@ module function crystallite_stress(dt,co,ip,el) result(converged_) todo = .true. converged_ = .false. ! pretend failed step of 1/subStepSizeCryst - crystallite_subdt(co,ip,el) = dt todo = .true. NiterationCrystallite = 0 cutbackLooping: do while (todo) diff --git a/src/homogenization.f90 b/src/homogenization.f90 index 52553b57b..d61fa57e8 100644 --- a/src/homogenization.f90 +++ b/src/homogenization.f90 @@ -63,7 +63,8 @@ module homogenization el !< element number end subroutine mech_partition - module subroutine mech_homogenize(ip,el) + module subroutine mech_homogenize(dt,ip,el) + real(pReal), intent(in) :: dt integer, intent(in) :: & ip, & !< integration point el !< element number @@ -257,7 +258,7 @@ subroutine materialpoint_stressAndItsTangent(dt,FEsolving_execIP,FEsolving_execE do co = 1, myNgrains call crystallite_orientations(co,ip,el) enddo - call mech_homogenize(ip,el) + call mech_homogenize(dt,ip,el) enddo IpLooping3 enddo elementLooping3 !$OMP END PARALLEL DO diff --git a/src/homogenization_mech.f90 b/src/homogenization_mech.f90 index 641e960fd..e3e9cfb3e 100644 --- a/src/homogenization_mech.f90 +++ b/src/homogenization_mech.f90 @@ -138,11 +138,13 @@ end subroutine mech_partition !-------------------------------------------------------------------------------------------------- !> @brief Average P and dPdF from the individual constituents. !-------------------------------------------------------------------------------------------------- -module subroutine mech_homogenize(ip,el) +module subroutine mech_homogenize(dt,ip,el) + real(pReal), intent(in) :: dt integer, intent(in) :: & ip, & !< integration point el !< element number + integer :: co,ce real(pReal) :: dPdFs(3,3,3,3,homogenization_Nconstituents(material_homogenizationAt(el))) @@ -152,11 +154,11 @@ module subroutine mech_homogenize(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(1,ip,el) + homogenization_dPdF(1:3,1:3,1:3,1:3,ce) = crystallite_stressTangent(dt,1,ip,el) case (HOMOGENIZATION_ISOSTRAIN_ID) chosenHomogenization do co = 1, homogenization_Nconstituents(material_homogenizationAt(el)) - dPdFs(:,:,:,:,co) = crystallite_stressTangent(co,ip,el) + dPdFs(:,:,:,:,co) = crystallite_stressTangent(dt,co,ip,el) enddo call mech_isostrain_averageStressAndItsTangent(& homogenization_P(1:3,1:3,ce), & @@ -167,7 +169,7 @@ module subroutine mech_homogenize(ip,el) case (HOMOGENIZATION_RGC_ID) chosenHomogenization do co = 1, homogenization_Nconstituents(material_homogenizationAt(el)) - dPdFs(:,:,:,:,co) = crystallite_stressTangent(co,ip,el) + dPdFs(:,:,:,:,co) = crystallite_stressTangent(dt,co,ip,el) enddo call mech_RGC_averageStressAndItsTangent(& homogenization_P(1:3,1:3,ce), & @@ -202,7 +204,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(co,ip,el) + dPdFs(:,:,:,:,co) = crystallite_stressTangent(subdt,co,ip,el) enddo doneAndHappy = & mech_RGC_updateState(crystallite_P(1:3,1:3,1:homogenization_Nconstituents(material_homogenizationAt(el)),ip,el), &