diff --git a/PRIVATE b/PRIVATE index 45ef93dbf..591964dcf 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit 45ef93dbfa3e0e6fa830914b3632e188c308a099 +Subproject commit 591964dcf8521d95f6cccbfe840d462c430e63d9 diff --git a/src/CPFEM2.f90 b/src/CPFEM2.f90 index 5a500875d..e696858cf 100644 --- a/src/CPFEM2.f90 +++ b/src/CPFEM2.f90 @@ -74,9 +74,23 @@ 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 homogenization_restartRead(fileHandle) + call constitutive_restartRead(fileHandle) + + call HDF5_closeFile(fileHandle) + endif end subroutine CPFEM_init @@ -86,7 +100,19 @@ end subroutine CPFEM_init !-------------------------------------------------------------------------------------------------- subroutine CPFEM_restartWrite - call crystallite_restartWrite + integer(HID_T) :: fileHandle + character(len=pStringLen) :: fileName + + + 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 homogenization_restartWrite(fileHandle) + call constitutive_restartWrite(fileHandle) + + call HDF5_closeFile(fileHandle) end subroutine CPFEM_restartWrite diff --git a/src/constitutive.f90 b/src/constitutive.f90 index 696611549..7dac3e129 100644 --- a/src/constitutive.f90 +++ b/src/constitutive.f90 @@ -15,7 +15,6 @@ module constitutive use discretization use parallelization use HDF5_utilities - use DAMASK_interface use results implicit none @@ -42,42 +41,14 @@ 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 :: & - crystallite_F0, & !< def grad at start of FE inc - crystallite_Fe, & !< current "elastic" def grad (end of converged time step) - crystallite_subFp0,& !< plastic def grad at start of crystallite inc - crystallite_subFi0,& !< intermediate def grad at start of crystallite inc - crystallite_Lp0, & !< plastic velocitiy grad at start of FE inc - crystallite_partitionedLp0, & !< plastic velocity grad at start of homog inc - crystallite_S0, & !< 2nd Piola-Kirchhoff stress vector at start of FE inc - crystallite_partitionedS0 !< 2nd Piola-Kirchhoff stress vector at start of homog inc - real(pReal), dimension(:,:,:,:,:), allocatable, public :: & - crystallite_P, & !< 1st Piola-Kirchhoff stress per grain - crystallite_Lp, & !< current plastic velocitiy grad (end of converged time step) - crystallite_S, & !< current 2nd Piola-Kirchhoff stress vector (end of converged time step) - crystallite_partitionedF0, & !< def grad at start of homog inc - crystallite_F !< def grad to be reached at end of homog inc type :: tTensorContainer real(pReal), dimension(:,:,:), allocatable :: data end type - type(tTensorContainer), dimension(:), allocatable :: & - constitutive_mech_Fi, & - constitutive_mech_Fi0, & - constitutive_mech_partitionedFi0, & - constitutive_mech_Li, & - constitutive_mech_Li0, & - constitutive_mech_partitionedLi0, & - constitutive_mech_Fp, & - constitutive_mech_Fp0, & - constitutive_mech_partitionedFp0 - - type :: tNumerics integer :: & iJacoLpresiduum, & !< frequency of Jacobian update of residuum in Lp @@ -162,10 +133,6 @@ module constitutive end subroutine damage_results - module subroutine mech_restart_read(fileHandle) - integer(HID_T), intent(in) :: fileHandle - end subroutine mech_restart_read - module subroutine mech_initializeRestorationPoints(ph,me) integer, intent(in) :: ph, me end subroutine mech_initializeRestorationPoints @@ -185,6 +152,61 @@ 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 + + module subroutine mech_restartWrite(groupHandle,ph) + integer(HID_T), intent(in) :: groupHandle + integer, intent(in) :: ph + end subroutine mech_restartWrite + + module subroutine mech_restartRead(groupHandle,ph) + integer(HID_T), intent(in) :: groupHandle + integer, intent(in) :: ph + end subroutine mech_restartRead + + + module function constitutive_mech_getS(co,ip,el) result(S) + integer, intent(in) :: co, ip, el + real(pReal), dimension(3,3) :: S + end function constitutive_mech_getS + + module function constitutive_mech_getLp(co,ip,el) result(Lp) + integer, intent(in) :: co, ip, el + real(pReal), dimension(3,3) :: Lp + end function constitutive_mech_getLp + + module function constitutive_mech_getF(co,ip,el) result(F) + integer, intent(in) :: co, ip, el + real(pReal), dimension(3,3) :: F + end function constitutive_mech_getF + + module function constitutive_mech_getF_e(co,ip,el) result(F_e) + integer, intent(in) :: co, ip, el + real(pReal), dimension(3,3) :: F_e + end function constitutive_mech_getF_e + + module function constitutive_mech_getP(co,ip,el) result(P) + integer, intent(in) :: co, ip, el + real(pReal), dimension(3,3) :: P + end function constitutive_mech_getP + + module function constitutive_thermal_T(co,ip,el) result(T) + integer, intent(in) :: co, ip, el + real(pReal) :: T + end function constitutive_thermal_T + + module subroutine constitutive_mech_setF(F,co,ip,el) + real(pReal), dimension(3,3), intent(in) :: F + integer, intent(in) :: co, ip, el + end subroutine constitutive_mech_setF + ! == cleaned:end =================================================================================== module function crystallite_stress(dt,co,ip,el) result(converged_) @@ -238,15 +260,12 @@ module constitutive dPhiDot_dPhi end subroutine constitutive_damage_getRateAndItsTangents - module subroutine constitutive_thermal_getRateAndItsTangents(TDot, dTDot_dT, T, S, Lp, ip, el) + module subroutine constitutive_thermal_getRateAndItsTangents(TDot, dTDot_dT, T, ip, el) integer, intent(in) :: & ip, & !< integration point number el !< element number real(pReal), intent(in) :: & T - real(pReal), intent(in), dimension(:,:,:,:,:) :: & - S, & !< current 2nd Piola Kitchoff stress vector - Lp !< plastic velocity gradient real(pReal), intent(inout) :: & TDot, & dTDot_dT @@ -342,13 +361,11 @@ module constitutive end subroutine constitutive_plastic_LpAndItsTangents - module subroutine constitutive_plastic_dependentState(F, co, ip, el) + module subroutine constitutive_plastic_dependentState(co,ip,el) integer, intent(in) :: & co, & !< component-ID of integration point ip, & !< integration point el !< element - real(pReal), intent(in), dimension(3,3) :: & - F !< elastic deformation gradient end subroutine constitutive_plastic_dependentState @@ -390,12 +407,17 @@ module constitutive converged, & crystallite_init, & crystallite_stress, & - crystallite_stressTangent, & + constitutive_mech_dPdF, & crystallite_orientations, & crystallite_push33ToRef, & - crystallite_restartWrite, & + constitutive_restartWrite, & + constitutive_restartRead, & integrateSourceState, & - crystallite_restartRead, & + constitutive_mech_setF, & + constitutive_mech_getP, & + constitutive_mech_getLp, & + constitutive_mech_getF, & + constitutive_mech_getS, & constitutive_initializeRestorationPoints, & constitutive_windForward, & PLASTICITY_UNDEFINED_ID, & @@ -615,7 +637,7 @@ end subroutine constitutive_LiAndItsTangents !-------------------------------------------------------------------------------------------------- !> @brief contains the constitutive equation for calculating the rate of change of microstructure !-------------------------------------------------------------------------------------------------- -function constitutive_damage_collectDotState(S, co, ip, el,ph,of) result(broken) +function constitutive_damage_collectDotState(co,ip,el,ph,of) result(broken) integer, intent(in) :: & co, & !< component-ID of integration point @@ -623,8 +645,6 @@ function constitutive_damage_collectDotState(S, co, ip, el,ph,of) result(broken) el, & !< element ph, & of - real(pReal), intent(in), dimension(3,3) :: & - S !< 2nd Piola Kirchhoff stress (vector notation) integer :: & so !< counter in source loop logical :: broken @@ -637,7 +657,7 @@ function constitutive_damage_collectDotState(S, co, ip, el,ph,of) result(broken) sourceType: select case (phase_source(so,ph)) case (SOURCE_damage_anisoBrittle_ID) sourceType - call source_damage_anisoBrittle_dotState(S, co, ip, el) ! correct stress? + call source_damage_anisoBrittle_dotState(constitutive_mech_getS(co,ip,el), co, ip, el) ! correct stress? case (SOURCE_damage_isoDuctile_ID) sourceType call source_damage_isoDuctile_dotState(co, ip, el) @@ -748,7 +768,6 @@ subroutine constitutive_allocateState(state, & allocate(state%atol (sizeState), source=0.0_pReal) allocate(state%state0 (sizeState,Nconstituents), source=0.0_pReal) allocate(state%partitionedState0(sizeState,Nconstituents), source=0.0_pReal) - allocate(state%subState0 (sizeState,Nconstituents), source=0.0_pReal) allocate(state%state (sizeState,Nconstituents), source=0.0_pReal) allocate(state%dotState (sizeDotState,Nconstituents), source=0.0_pReal) @@ -792,17 +811,14 @@ end subroutine constitutive_restore !-------------------------------------------------------------------------------------------------- subroutine constitutive_forward - integer :: i, j + integer :: ph, so - crystallite_F0 = crystallite_F - crystallite_Lp0 = crystallite_Lp - crystallite_S0 = crystallite_S call constitutive_mech_forward() - do i = 1, size(sourceState) - do j = 1,phase_Nsources(i) - sourceState(i)%p(j)%state0 = sourceState(i)%p(j)%state + do ph = 1, size(sourceState) + do so = 1,phase_Nsources(ph) + sourceState(ph)%p(so)%state0 = sourceState(ph)%p(so)%state enddo; enddo end subroutine constitutive_forward @@ -838,12 +854,12 @@ end subroutine constitutive_results subroutine crystallite_init integer :: & - Nconstituents, & ph, & me, & co, & !< counter in integration point component loop ip, & !< counter in integration point loop el, & !< counter in element loop + so, & cMax, & !< maximum number of integration point components iMax, & !< maximum number of integration points eMax !< maximum number of elements @@ -866,19 +882,6 @@ subroutine crystallite_init iMax = discretization_nIPs eMax = discretization_Nelems - allocate(crystallite_F(3,3,cMax,iMax,eMax),source=0.0_pReal) - - allocate(crystallite_S0, & - crystallite_F0,crystallite_Lp0, & - crystallite_partitionedS0, & - crystallite_partitionedF0,& - crystallite_partitionedLp0, & - crystallite_S,crystallite_P, & - crystallite_Fe,crystallite_Lp, & - crystallite_subFp0,crystallite_subFi0, & - source = crystallite_F) - - allocate(crystallite_subdt(cMax,iMax,eMax),source=0.0_pReal) allocate(crystallite_orientation(cMax,iMax,eMax)) num_crystallite => config_numerics%get('crystallite',defaultVal=emptyDict) @@ -914,27 +917,10 @@ subroutine crystallite_init phases => config_material%get('phase') - allocate(constitutive_mech_Fi(phases%length)) - allocate(constitutive_mech_Fi0(phases%length)) - allocate(constitutive_mech_partitionedFi0(phases%length)) - allocate(constitutive_mech_Fp(phases%length)) - allocate(constitutive_mech_Fp0(phases%length)) - allocate(constitutive_mech_partitionedFp0(phases%length)) - allocate(constitutive_mech_Li(phases%length)) - allocate(constitutive_mech_Li0(phases%length)) - allocate(constitutive_mech_partitionedLi0(phases%length)) do ph = 1, phases%length - Nconstituents = count(material_phaseAt == ph) * discretization_nIPs - - allocate(constitutive_mech_Fi(ph)%data(3,3,Nconstituents)) - allocate(constitutive_mech_Fi0(ph)%data(3,3,Nconstituents)) - allocate(constitutive_mech_partitionedFi0(ph)%data(3,3,Nconstituents)) - allocate(constitutive_mech_Fp(ph)%data(3,3,Nconstituents)) - allocate(constitutive_mech_Fp0(ph)%data(3,3,Nconstituents)) - allocate(constitutive_mech_partitionedFp0(ph)%data(3,3,Nconstituents)) - allocate(constitutive_mech_Li(ph)%data(3,3,Nconstituents)) - allocate(constitutive_mech_Li0(ph)%data(3,3,Nconstituents)) - allocate(constitutive_mech_partitionedLi0(ph)%data(3,3,Nconstituents)) + do so = 1, phase_Nsources(ph) + allocate(sourceState(ph)%p(so)%subState0,source=sourceState(ph)%p(so)%state0) ! ToDo: hack + enddo enddo print'(a42,1x,i10)', ' # of elements: ', eMax @@ -942,34 +928,6 @@ subroutine crystallite_init print'(a42,1x,i10)', 'max # of constituents/integration point: ', cMax flush(IO_STDOUT) - !$OMP PARALLEL DO PRIVATE(ph,me) - do el = 1, size(material_phaseMemberAt,3); do ip = 1, size(material_phaseMemberAt,2) - do co = 1, homogenization_Nconstituents(material_homogenizationAt(el)) - - ph = material_phaseAt(co,el) - me = material_phaseMemberAt(co,ip,el) - constitutive_mech_Fp0(ph)%data(1:3,1:3,me) = material_orientation0(co,ip,el)%asMatrix() ! Fp reflects initial orientation (see 10.1016/j.actamat.2006.01.005) - constitutive_mech_Fp0(ph)%data(1:3,1:3,me) = constitutive_mech_Fp0(ph)%data(1:3,1:3,me) & - / math_det33(constitutive_mech_Fp0(ph)%data(1:3,1:3,me))**(1.0_pReal/3.0_pReal) - constitutive_mech_Fi0(ph)%data(1:3,1:3,me) = math_I3 - - crystallite_F0(1:3,1:3,co,ip,el) = math_I3 - - crystallite_Fe(1:3,1:3,co,ip,el) = math_inv33(matmul(constitutive_mech_Fi0(ph)%data(1:3,1:3,me), & - constitutive_mech_Fp0(ph)%data(1:3,1:3,me))) ! assuming that euler angles are given in internal strain free configuration - constitutive_mech_Fp(ph)%data(1:3,1:3,me) = constitutive_mech_Fp0(ph)%data(1:3,1:3,me) - constitutive_mech_Fi(ph)%data(1:3,1:3,me) = constitutive_mech_Fi0(ph)%data(1:3,1:3,me) - - constitutive_mech_partitionedFi0(ph)%data(1:3,1:3,me) = constitutive_mech_Fi0(ph)%data(1:3,1:3,me) - constitutive_mech_partitionedFp0(ph)%data(1:3,1:3,me) = constitutive_mech_Fp0(ph)%data(1:3,1:3,me) - - enddo - enddo; enddo - !$OMP END PARALLEL DO - - crystallite_partitionedF0 = crystallite_F0 - crystallite_F = crystallite_F0 - !$OMP PARALLEL DO PRIVATE(ph,me) do el = 1, size(material_phaseMemberAt,3) @@ -978,7 +936,7 @@ subroutine crystallite_init ph = material_phaseAt(co,el) me = material_phaseMemberAt(co,ip,el) call crystallite_orientations(co,ip,el) - call constitutive_plastic_dependentState(crystallite_partitionedF0(1:3,1:3,co,ip,el),co,ip,el) ! update dependent state variables to be consistent with basic states + call constitutive_plastic_dependentState(co,ip,el) ! update dependent state variables to be consistent with basic states enddo enddo enddo @@ -1003,15 +961,11 @@ subroutine constitutive_initializeRestorationPoints(ip,el) do co = 1,homogenization_Nconstituents(material_homogenizationAt(el)) ph = material_phaseAt(co,el) me = material_phaseMemberAt(co,ip,el) - crystallite_partitionedLp0(1:3,1:3,co,ip,el) = crystallite_Lp0(1:3,1:3,co,ip,el) - crystallite_partitionedF0(1:3,1:3,co,ip,el) = crystallite_F0(1:3,1:3,co,ip,el) - crystallite_partitionedS0(1:3,1:3,co,ip,el) = crystallite_S0(1:3,1:3,co,ip,el) call mech_initializeRestorationPoints(ph,me) do so = 1, phase_Nsources(material_phaseAt(co,el)) - sourceState(material_phaseAt(co,el))%p(so)%partitionedState0(:,material_phasememberAt(co,ip,el)) = & - sourceState(material_phaseAt(co,el))%p(so)%state0( :,material_phasememberAt(co,ip,el)) + sourceState(ph)%p(so)%partitionedState0(:,me) = sourceState(ph)%p(so)%state0(:,me) enddo enddo @@ -1035,9 +989,6 @@ subroutine constitutive_windForward(ip,el) do co = 1,homogenization_Nconstituents(material_homogenizationAt(el)) ph = material_phaseAt(co,el) me = material_phaseMemberAt(co,ip,el) - crystallite_partitionedF0 (1:3,1:3,co,ip,el) = crystallite_F (1:3,1:3,co,ip,el) - crystallite_partitionedLp0(1:3,1:3,co,ip,el) = crystallite_Lp(1:3,1:3,co,ip,el) - crystallite_partitionedS0 (1:3,1:3,co,ip,el) = crystallite_S (1:3,1:3,co,ip,el) call constitutive_mech_windForward(ph,me) do so = 1, phase_Nsources(material_phaseAt(co,el)) @@ -1048,136 +999,6 @@ subroutine constitutive_windForward(ip,el) end subroutine constitutive_windForward -!-------------------------------------------------------------------------------------------------- -!> @brief Calculate tangent (dPdF). -!-------------------------------------------------------------------------------------------------- -function crystallite_stressTangent(co,ip,el) result(dPdF) - - real(pReal), dimension(3,3,3,3) :: dPdF - integer, intent(in) :: & - co, & !< counter in constituent loop - ip, & !< counter in integration point loop - el !< counter in element loop - - 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, & - crystallite_Fe(1:3,1:3,co,ip,el), & - constitutive_mech_Fi(ph)%data(1:3,1:3,me),co,ip,el) - call constitutive_LiAndItsTangents(devNull,dLidS,dLidFi, & - crystallite_S (1:3,1:3,co,ip,el), & - 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(crystallite_subFp0(1:3,1:3,co,ip,el)) - invSubFi0 = math_inv33(crystallite_subFi0(1:3,1:3,co,ip,el)) - - 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) & - + crystallite_subdt(co,ip,el)*matmul(invSubFi0,dLidFi(1:3,1:3,o,p)) - 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)) - 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, & - crystallite_S (1:3,1:3,co,ip,el), & - 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(crystallite_F(1:3,1:3,co,ip,el),invSubFp0) - temp_33_3 = matmul(matmul(crystallite_F(1:3,1:3,co,ip,el),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 = crystallite_subdt(co,ip,el)*math_mul3333xx3333(dSdFe,temp_3333) & - + 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) = -crystallite_subdt(co,ip,el) & - * matmul(invSubFp0, matmul(temp_3333(1:3,1:3,p,o),invFi)) - enddo; enddo - -!-------------------------------------------------------------------------------------------------- -! assemble dPdF - temp_33_1 = matmul(crystallite_S(1:3,1:3,co,ip,el),transpose(invFp)) - temp_33_2 = matmul(crystallite_F(1:3,1:3,co,ip,el),invFp) - temp_33_3 = matmul(temp_33_2,crystallite_S(1:3,1:3,co,ip,el)) - - 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(crystallite_F(1:3,1:3,co,ip,el), & - 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 !-------------------------------------------------------------------------------------------------- @@ -1189,7 +1010,8 @@ subroutine crystallite_orientations(co,ip,el) el !< counter in element loop - call crystallite_orientation(co,ip,el)%fromMatrix(transpose(math_rotationalPart(crystallite_Fe(1:3,1:3,co,ip,el)))) + call crystallite_orientation(co,ip,el)%fromMatrix(transpose(math_rotationalPart(& + constitutive_mech_getF_e(co,ip,el)))) if (plasticState(material_phaseAt(1,el))%nonlocal) & call plastic_nonlocal_updateCompatibility(crystallite_orientation, & @@ -1205,17 +1027,17 @@ end subroutine crystallite_orientations function crystallite_push33ToRef(co,ip,el, tensor33) real(pReal), dimension(3,3), intent(in) :: tensor33 - real(pReal), dimension(3,3) :: T integer, intent(in):: & el, & ip, & co - real(pReal), dimension(3,3) :: crystallite_push33ToRef + real(pReal), dimension(3,3) :: T + + + T = matmul(material_orientation0(co,ip,el)%asMatrix(),transpose(math_inv33(constitutive_mech_getF(co,ip,el)))) ! ToDo: initial orientation correct? - T = matmul(material_orientation0(co,ip,el)%asMatrix(), & ! ToDo: initial orientation correct? - transpose(math_inv33(crystallite_F(1:3,1:3,co,ip,el)))) crystallite_push33ToRef = matmul(transpose(T),matmul(tensor33,T)) end function crystallite_push33ToRef @@ -1254,7 +1076,7 @@ function integrateSourceState(dt,co,ip,el) result(broken) converged_ = .true. broken = constitutive_thermal_collectDotState(ph,me) - broken = broken .or. constitutive_damage_collectDotState(crystallite_S(1:3,1:3,co,ip,el), co,ip,el,ph,me) + broken = broken .or. constitutive_damage_collectDotState(co,ip,el,ph,me) if(broken) return do so = 1, phase_Nsources(ph) @@ -1272,7 +1094,7 @@ function integrateSourceState(dt,co,ip,el) result(broken) enddo broken = constitutive_thermal_collectDotState(ph,me) - broken = broken .or. constitutive_damage_collectDotState(crystallite_S(1:3,1:3,co,ip,el), co,ip,el,ph,me) + broken = broken .or. constitutive_damage_collectDotState(co,ip,el,ph,me) if(broken) exit iteration do so = 1, phase_Nsources(ph) @@ -1292,7 +1114,7 @@ function integrateSourceState(dt,co,ip,el) result(broken) enddo if(converged_) then - broken = constitutive_damage_deltaState(crystallite_Fe(1:3,1:3,co,ip,el),co,ip,el,ph,me) + broken = constitutive_damage_deltaState(constitutive_mech_getF_e(co,ip,el),co,ip,el,ph,me) exit iteration endif @@ -1345,90 +1167,60 @@ 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 +! ToDo: Merge data into one file for MPI !-------------------------------------------------------------------------------------------------- -subroutine crystallite_restartWrite +subroutine constitutive_restartWrite(fileHandle) + integer(HID_T), intent(in) :: fileHandle + + integer(HID_T), dimension(2) :: groupHandle 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') + groupHandle(1) = HDF5_addGroup(fileHandle,'phase') - call HDF5_write(fileHandle,crystallite_F,'F') - call HDF5_write(fileHandle,crystallite_Lp, 'L_p') - call HDF5_write(fileHandle,crystallite_S, 'S') + do ph = 1, size(material_name_phase) + + groupHandle(2) = HDF5_addGroup(groupHandle(1),material_name_phase(ph)) + + call mech_restartWrite(groupHandle(2),ph) + + call HDF5_closeGroup(groupHandle(2)) - groupHandle = HDF5_addGroup(fileHandle,'phase') - do ph = 1,size(material_name_phase) - write(datasetName,'(i0,a)') ph,'_omega' - call HDF5_write(groupHandle,plasticState(ph)%state,datasetName) - write(datasetName,'(i0,a)') ph,'_F_i' - call HDF5_write(groupHandle,constitutive_mech_Fi(ph)%data,datasetName) - write(datasetName,'(i0,a)') ph,'_L_i' - call HDF5_write(groupHandle,constitutive_mech_Li(ph)%data,datasetName) - write(datasetName,'(i0,a)') ph,'_F_p' - call HDF5_write(groupHandle,constitutive_mech_Fp(ph)%data,datasetName) enddo - call HDF5_closeGroup(groupHandle) - groupHandle = HDF5_addGroup(fileHandle,'homogenization') - do ph = 1, size(material_name_homogenization) - write(datasetName,'(i0,a)') ph,'_omega' - call HDF5_write(groupHandle,homogState(ph)%state,datasetName) - enddo - call HDF5_closeGroup(groupHandle) + call HDF5_closeGroup(groupHandle(1)) - 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 +! ToDo: Merge data into one file for MPI !-------------------------------------------------------------------------------------------------- -subroutine crystallite_restartRead +subroutine constitutive_restartRead(fileHandle) + integer(HID_T), intent(in) :: fileHandle + + integer(HID_T), dimension(2) :: groupHandle integer :: ph - integer(HID_T) :: fileHandle, groupHandle - character(len=pStringLen) :: fileName, 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(1) = HDF5_openGroup(fileHandle,'phase') - call HDF5_read(fileHandle,crystallite_F0, 'F') - call HDF5_read(fileHandle,crystallite_Lp0,'L_p') - call HDF5_read(fileHandle,crystallite_S0, 'S') + do ph = 1, size(material_name_phase) + + groupHandle(2) = HDF5_openGroup(groupHandle(1),material_name_phase(ph)) + + call mech_restartRead(groupHandle(2),ph) + + call HDF5_closeGroup(groupHandle(2)) - groupHandle = HDF5_openGroup(fileHandle,'phase') - do ph = 1,size(material_name_phase) - write(datasetName,'(i0,a)') ph,'_omega' - call HDF5_read(groupHandle,plasticState(ph)%state0,datasetName) - write(datasetName,'(i0,a)') ph,'_F_i' - call HDF5_read(groupHandle,constitutive_mech_Fi0(ph)%data,datasetName) - write(datasetName,'(i0,a)') ph,'_L_i' - call HDF5_read(groupHandle,constitutive_mech_Li0(ph)%data,datasetName) - write(datasetName,'(i0,a)') ph,'_F_p' - call HDF5_read(groupHandle,constitutive_mech_Fp0(ph)%data,datasetName) enddo - call HDF5_closeGroup(groupHandle) - groupHandle = HDF5_openGroup(fileHandle,'homogenization') - do ph = 1,size(material_name_homogenization) - write(datasetName,'(i0,a)') ph,'_omega' - call HDF5_read(groupHandle,homogState(ph)%state0,datasetName) - enddo - call HDF5_closeGroup(groupHandle) + call HDF5_closeGroup(groupHandle(1)) - call HDF5_closeFile(fileHandle) - -end subroutine crystallite_restartRead +end subroutine constitutive_restartRead end module constitutive diff --git a/src/constitutive_mech.f90 b/src/constitutive_mech.f90 index c48c59ec9..44e777c80 100644 --- a/src/constitutive_mech.f90 +++ b/src/constitutive_mech.f90 @@ -15,6 +15,31 @@ submodule(constitutive) constitutive_mech integer(kind(SOURCE_undefined_ID)), dimension(:,:), allocatable :: & phase_stiffnessDegradation !< active stiffness degradation mechanisms of each phase + type(tTensorContainer), dimension(:), allocatable :: & + ! current value + constitutive_mech_Fe, & + constitutive_mech_Fi, & + constitutive_mech_Fp, & + constitutive_mech_F, & + constitutive_mech_Li, & + constitutive_mech_Lp, & + constitutive_mech_S, & + constitutive_mech_P, & + ! converged value at end of last solver increment + constitutive_mech_Fi0, & + constitutive_mech_Fp0, & + constitutive_mech_F0, & + constitutive_mech_Li0, & + constitutive_mech_Lp0, & + constitutive_mech_S0, & + ! converged value at end of last homogenization increment (RGC only) + constitutive_mech_partitionedFi0, & + constitutive_mech_partitionedFp0, & + constitutive_mech_partitionedF0, & + constitutive_mech_partitionedLi0, & + constitutive_mech_partitionedLp0, & + constitutive_mech_partitionedS0 + interface @@ -184,12 +209,9 @@ submodule(constitutive) constitutive_mech of end subroutine plastic_disloTungsten_dotState - module subroutine plastic_nonlocal_dotState(Mp, F, Temperature,timestep, & - instance,of,ip,el) + module subroutine plastic_nonlocal_dotState(Mp,Temperature,timestep,instance,of,ip,el) real(pReal), dimension(3,3), intent(in) :: & Mp !< MandelStress - real(pReal), dimension(3,3,homogenization_maxNconstituents,discretization_nIPs,discretization_Nelems), intent(in) :: & - F !< deformation gradient real(pReal), intent(in) :: & Temperature, & !< temperature timestep !< substepped crystallite time increment @@ -215,9 +237,7 @@ submodule(constitutive) constitutive_mech of end subroutine plastic_dislotungsten_dependentState - module subroutine plastic_nonlocal_dependentState(F, instance, of, ip, el) - real(pReal), dimension(3,3), intent(in) :: & - F !< deformation gradient + module subroutine plastic_nonlocal_dependentState(instance, of, ip, el) integer, intent(in) :: & instance, & of, & @@ -302,8 +322,13 @@ contains module subroutine mech_init integer :: & + el, & + ip, & + co, & ph, & - stiffDegradationCtr + me, & + stiffDegradationCtr, & + Nconstituents class(tNode), pointer :: & num_crystallite, & phases, & @@ -322,7 +347,51 @@ module subroutine mech_init allocate(phase_NstiffnessDegradations(phases%length),source=0) allocate(output_constituent(phases%length)) + allocate(constitutive_mech_Fe(phases%length)) + allocate(constitutive_mech_Fi(phases%length)) + allocate(constitutive_mech_Fi0(phases%length)) + allocate(constitutive_mech_partitionedFi0(phases%length)) + allocate(constitutive_mech_Fp(phases%length)) + allocate(constitutive_mech_Fp0(phases%length)) + allocate(constitutive_mech_partitionedFp0(phases%length)) + allocate(constitutive_mech_F(phases%length)) + allocate(constitutive_mech_F0(phases%length)) + allocate(constitutive_mech_partitionedF0(phases%length)) + allocate(constitutive_mech_Li(phases%length)) + allocate(constitutive_mech_Li0(phases%length)) + allocate(constitutive_mech_partitionedLi0(phases%length)) + allocate(constitutive_mech_partitionedLp0(phases%length)) + allocate(constitutive_mech_Lp0(phases%length)) + allocate(constitutive_mech_Lp(phases%length)) + allocate(constitutive_mech_S(phases%length)) + allocate(constitutive_mech_P(phases%length)) + allocate(constitutive_mech_S0(phases%length)) + allocate(constitutive_mech_partitionedS0(phases%length)) + do ph = 1, phases%length + Nconstituents = count(material_phaseAt == ph) * discretization_nIPs + + allocate(constitutive_mech_Fi(ph)%data(3,3,Nconstituents)) + allocate(constitutive_mech_Fe(ph)%data(3,3,Nconstituents)) + allocate(constitutive_mech_Fi0(ph)%data(3,3,Nconstituents)) + allocate(constitutive_mech_partitionedFi0(ph)%data(3,3,Nconstituents)) + allocate(constitutive_mech_Fp(ph)%data(3,3,Nconstituents)) + allocate(constitutive_mech_Fp0(ph)%data(3,3,Nconstituents)) + allocate(constitutive_mech_partitionedFp0(ph)%data(3,3,Nconstituents)) + allocate(constitutive_mech_Li(ph)%data(3,3,Nconstituents)) + allocate(constitutive_mech_Li0(ph)%data(3,3,Nconstituents)) + allocate(constitutive_mech_partitionedLi0(ph)%data(3,3,Nconstituents)) + allocate(constitutive_mech_partitionedLp0(ph)%data(3,3,Nconstituents)) + allocate(constitutive_mech_Lp0(ph)%data(3,3,Nconstituents)) + allocate(constitutive_mech_Lp(ph)%data(3,3,Nconstituents)) + allocate(constitutive_mech_S(ph)%data(3,3,Nconstituents),source=0.0_pReal) + allocate(constitutive_mech_P(ph)%data(3,3,Nconstituents),source=0.0_pReal) + allocate(constitutive_mech_S0(ph)%data(3,3,Nconstituents),source=0.0_pReal) + allocate(constitutive_mech_partitionedS0(ph)%data(3,3,Nconstituents),source=0.0_pReal) + allocate(constitutive_mech_F(ph)%data(3,3,Nconstituents)) + allocate(constitutive_mech_F0(ph)%data(3,3,Nconstituents)) + allocate(constitutive_mech_partitionedF0(ph)%data(3,3,Nconstituents)) + phase => phases%get(ph) mech => phase%get('mechanics') #if defined(__GFORTRAN__) @@ -355,6 +424,33 @@ module subroutine mech_init enddo endif + !$OMP PARALLEL DO PRIVATE(ph,me) + do el = 1, size(material_phaseMemberAt,3); do ip = 1, size(material_phaseMemberAt,2) + do co = 1, homogenization_Nconstituents(material_homogenizationAt(el)) + + ph = material_phaseAt(co,el) + me = material_phaseMemberAt(co,ip,el) + + constitutive_mech_Fp0(ph)%data(1:3,1:3,me) = material_orientation0(co,ip,el)%asMatrix() ! Fp reflects initial orientation (see 10.1016/j.actamat.2006.01.005) + constitutive_mech_Fp0(ph)%data(1:3,1:3,me) = constitutive_mech_Fp0(ph)%data(1:3,1:3,me) & + / math_det33(constitutive_mech_Fp0(ph)%data(1:3,1:3,me))**(1.0_pReal/3.0_pReal) + constitutive_mech_Fi0(ph)%data(1:3,1:3,me) = math_I3 + constitutive_mech_F0(ph)%data(1:3,1:3,me) = math_I3 + + constitutive_mech_Fe(ph)%data(1:3,1:3,me) = math_inv33(matmul(constitutive_mech_Fi0(ph)%data(1:3,1:3,me), & + constitutive_mech_Fp0(ph)%data(1:3,1:3,me))) ! assuming that euler angles are given in internal strain free configuration + constitutive_mech_Fp(ph)%data(1:3,1:3,me) = constitutive_mech_Fp0(ph)%data(1:3,1:3,me) + constitutive_mech_Fi(ph)%data(1:3,1:3,me) = constitutive_mech_Fi0(ph)%data(1:3,1:3,me) + constitutive_mech_F(ph)%data(1:3,1:3,me) = constitutive_mech_F0(ph)%data(1:3,1:3,me) + + constitutive_mech_partitionedFi0(ph)%data(1:3,1:3,me) = constitutive_mech_Fi0(ph)%data(1:3,1:3,me) + constitutive_mech_partitionedFp0(ph)%data(1:3,1:3,me) = constitutive_mech_Fp0(ph)%data(1:3,1:3,me) + constitutive_mech_partitionedF0(ph)%data(1:3,1:3,me) = constitutive_mech_F0(ph)%data(1:3,1:3,me) + + enddo + enddo; enddo + !$OMP END PARALLEL DO + ! initialize plasticity allocate(plasticState(phases%length)) @@ -480,32 +576,35 @@ end subroutine constitutive_hooke_SandItsTangents !-------------------------------------------------------------------------------------------------- !> @brief calls microstructure function of the different plasticity constitutive models !-------------------------------------------------------------------------------------------------- -module subroutine constitutive_plastic_dependentState(F, co, ip, el) +module subroutine constitutive_plastic_dependentState(co, ip, el) integer, intent(in) :: & co, & !< component-ID of integration point ip, & !< integration point el !< element - real(pReal), intent(in), dimension(3,3) :: & - F !< deformation gradient integer :: & ho, & !< homogenization tme, & !< thermal member position - instance, of + instance, me + ho = material_homogenizationAt(el) tme = material_homogenizationMemberAt(ip,el) - of = material_phasememberAt(co,ip,el) + me = material_phasememberAt(co,ip,el) instance = phase_plasticityInstance(material_phaseAt(co,el)) plasticityType: select case (phase_plasticity(material_phaseAt(co,el))) + case (PLASTICITY_DISLOTWIN_ID) plasticityType - call plastic_dislotwin_dependentState(temperature(ho)%p(tme),instance,of) + call plastic_dislotwin_dependentState(temperature(ho)%p(tme),instance,me) + case (PLASTICITY_DISLOTUNGSTEN_ID) plasticityType - call plastic_dislotungsten_dependentState(instance,of) + call plastic_dislotungsten_dependentState(instance,me) + case (PLASTICITY_NONLOCAL_ID) plasticityType - call plastic_nonlocal_dependentState (F,instance,of,ip,el) + call plastic_nonlocal_dependentState(instance,me,ip,el) + end select plasticityType end subroutine constitutive_plastic_dependentState @@ -539,13 +638,13 @@ module subroutine constitutive_plastic_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, & ho, & !< homogenization tme !< thermal member position integer :: & - i, j, instance, of + i, j, instance, me ho = material_homogenizationAt(el) tme = material_homogenizationMemberAt(ip,el) Mp = matmul(matmul(transpose(Fi),Fi),S) - of = material_phasememberAt(co,ip,el) + me = material_phasememberAt(co,ip,el) instance = phase_plasticityInstance(material_phaseAt(co,el)) plasticityType: select case (phase_plasticity(material_phaseAt(co,el))) @@ -555,22 +654,22 @@ module subroutine constitutive_plastic_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, & dLp_dMp = 0.0_pReal case (PLASTICITY_ISOTROPIC_ID) plasticityType - call plastic_isotropic_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,of) + call plastic_isotropic_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,me) case (PLASTICITY_PHENOPOWERLAW_ID) plasticityType - call plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,of) + call plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,me) case (PLASTICITY_KINEHARDENING_ID) plasticityType - call plastic_kinehardening_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,of) + call plastic_kinehardening_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,me) case (PLASTICITY_NONLOCAL_ID) plasticityType - call plastic_nonlocal_LpAndItsTangent(Lp,dLp_dMp,Mp, temperature(ho)%p(tme),instance,of,ip,el) + call plastic_nonlocal_LpAndItsTangent(Lp,dLp_dMp,Mp, temperature(ho)%p(tme),instance,me,ip,el) case (PLASTICITY_DISLOTWIN_ID) plasticityType - call plastic_dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,temperature(ho)%p(tme),instance,of) + call plastic_dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,temperature(ho)%p(tme),instance,me) case (PLASTICITY_DISLOTUNGSTEN_ID) plasticityType - call plastic_dislotungsten_LpAndItsTangent(Lp,dLp_dMp,Mp,temperature(ho)%p(tme),instance,of) + call plastic_dislotungsten_LpAndItsTangent(Lp,dLp_dMp,Mp,temperature(ho)%p(tme),instance,me) end select plasticityType @@ -586,7 +685,7 @@ end subroutine constitutive_plastic_LpAndItsTangents !-------------------------------------------------------------------------------------------------- !> @brief contains the constitutive equation for calculating the rate of change of microstructure !-------------------------------------------------------------------------------------------------- -function mech_collectDotState(subdt, co, ip, el,ph,of) result(broken) +function mech_collectDotState(subdt,co,ip,el,ph,of) result(broken) integer, intent(in) :: & co, & !< component-ID of integration point @@ -601,15 +700,15 @@ function mech_collectDotState(subdt, co, ip, el,ph,of) result(broken) integer :: & ho, & !< homogenization tme, & !< thermal member position - i, & !< counter in source loop instance logical :: broken + ho = material_homogenizationAt(el) tme = material_homogenizationMemberAt(ip,el) instance = phase_plasticityInstance(ph) Mp = matmul(matmul(transpose(constitutive_mech_Fi(ph)%data(1:3,1:3,of)),& - constitutive_mech_Fi(ph)%data(1:3,1:3,of)),crystallite_S(1:3,1:3,co,ip,el)) + constitutive_mech_Fi(ph)%data(1:3,1:3,of)),constitutive_mech_S(ph)%data(1:3,1:3,of)) plasticityType: select case (phase_plasticity(ph)) @@ -629,8 +728,7 @@ function mech_collectDotState(subdt, co, ip, el,ph,of) result(broken) call plastic_disloTungsten_dotState(Mp,temperature(ho)%p(tme),instance,of) case (PLASTICITY_NONLOCAL_ID) plasticityType - call plastic_nonlocal_dotState(Mp,crystallite_partitionedF0,temperature(ho)%p(tme),subdt, & - instance,of,ip,el) + call plastic_nonlocal_dotState(Mp,temperature(ho)%p(tme),subdt,instance,of,ip,el) end select plasticityType broken = any(IEEE_is_NaN(plasticState(ph)%dotState(:,of))) @@ -642,7 +740,7 @@ end function mech_collectDotState !> @brief for constitutive models having an instantaneous change of state !> will return false if delta state is not needed/supported by the constitutive model !-------------------------------------------------------------------------------------------------- -function constitutive_deltaState(S, Fi, co, ip, el, ph, of) result(broken) +function constitutive_deltaState(co, ip, el, ph, of) result(broken) integer, intent(in) :: & co, & !< component-ID of integration point @@ -650,19 +748,19 @@ function constitutive_deltaState(S, Fi, co, ip, el, ph, of) result(broken) el, & !< element ph, & of - real(pReal), intent(in), dimension(3,3) :: & - S, & !< 2nd Piola Kirchhoff stress - Fi !< intermediate deformation gradient + logical :: & + broken + real(pReal), dimension(3,3) :: & Mp integer :: & instance, & myOffset, & mySize - logical :: & - broken - Mp = matmul(matmul(transpose(Fi),Fi),S) + + Mp = matmul(matmul(transpose(constitutive_mech_Fi(ph)%data(1:3,1:3,of)),& + constitutive_mech_Fi(ph)%data(1:3,1:3,of)),constitutive_mech_S(ph)%data(1:3,1:3,of)) instance = phase_plasticityInstance(ph) plasticityType: select case (phase_plasticity(ph)) @@ -728,18 +826,14 @@ module subroutine mech_results(group,ph) end subroutine mech_results - module subroutine mech_restart_read(fileHandle) - integer(HID_T), intent(in) :: fileHandle - end subroutine mech_restart_read - !-------------------------------------------------------------------------------------------------- !> @brief calculation of stress (P) with time integration based on a residuum in Lp and !> intermediate acceleration of the Newton-Raphson correction !-------------------------------------------------------------------------------------------------- -function integrateStress(F,Delta_t,co,ip,el) result(broken) +function integrateStress(F,subFp0,subFi0,Delta_t,co,ip,el) result(broken) - real(pReal), dimension(3,3), intent(in) :: F + real(pReal), dimension(3,3), intent(in) :: F,subFp0,subFi0 real(pReal), intent(in) :: Delta_t integer, intent(in):: el, & ! element index ip, & ! integration point index @@ -798,19 +892,20 @@ function integrateStress(F,Delta_t,co,ip,el) result(broken) jacoCounterLi ! counters to check for Jacobian update logical :: error,broken - broken = .true. - call constitutive_plastic_dependentState(crystallite_F(1:3,1:3,co,ip,el),co,ip,el) + broken = .true. ph = material_phaseAt(co,el) me = material_phaseMemberAt(co,ip,el) - Lpguess = crystallite_Lp(1:3,1:3,co,ip,el) ! take as first guess + call constitutive_plastic_dependentState(co,ip,el) + + Lpguess = constitutive_mech_Lp(ph)%data(1:3,1:3,me) ! take as first guess Liguess = constitutive_mech_Li(ph)%data(1:3,1:3,me) ! take as first guess - call math_invert33(invFp_current,devNull,error,crystallite_subFp0(1:3,1:3,co,ip,el)) + call math_invert33(invFp_current,devNull,error,subFp0) if (error) return ! error - call math_invert33(invFi_current,devNull,error,crystallite_subFi0(1:3,1:3,co,ip,el)) + call math_invert33(invFi_current,devNull,error,subFi0) if (error) return ! error A = matmul(F,invFp_current) ! intermediate tensor needed later to calculate dFe_dLp @@ -935,13 +1030,13 @@ function integrateStress(F,Delta_t,co,ip,el) result(broken) call math_invert33(Fp_new,devNull,error,invFp_new) if (error) return ! error - crystallite_P (1:3,1:3,co,ip,el) = matmul(matmul(F,invFp_new),matmul(S,transpose(invFp_new))) - crystallite_S (1:3,1:3,co,ip,el) = S - crystallite_Lp (1:3,1:3,co,ip,el) = Lpguess + constitutive_mech_P(ph)%data(1:3,1:3,me) = matmul(matmul(F,invFp_new),matmul(S,transpose(invFp_new))) + constitutive_mech_S(ph)%data(1:3,1:3,me) = S + constitutive_mech_Lp(ph)%data(1:3,1:3,me) = Lpguess constitutive_mech_Li(ph)%data(1:3,1:3,me) = Liguess - constitutive_mech_Fp(ph)%data(1:3,1:3,me) = Fp_new / math_det33(Fp_new)**(1.0_pReal/3.0_pReal) ! regularize + constitutive_mech_Fp(ph)%data(1:3,1:3,me) = Fp_new / math_det33(Fp_new)**(1.0_pReal/3.0_pReal) ! regularize constitutive_mech_Fi(ph)%data(1:3,1:3,me) = Fi_new - crystallite_Fe (1:3,1:3,co,ip,el) = matmul(matmul(F,invFp_new),invFi_new) + constitutive_mech_Fe(ph)%data(1:3,1:3,me) = matmul(matmul(F,invFp_new),invFi_new) broken = .false. end function integrateStress @@ -951,9 +1046,10 @@ end function integrateStress !> @brief integrate stress, state with adaptive 1st order explicit Euler method !> using Fixed Point Iteration to adapt the stepsize !-------------------------------------------------------------------------------------------------- -function integrateStateFPI(F_0,F,Delta_t,co,ip,el) result(broken) +function integrateStateFPI(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el) result(broken) - real(pReal), intent(in),dimension(3,3) :: F_0,F + real(pReal), intent(in),dimension(3,3) :: F_0,F,subFp0,subFi0 + real(pReal), intent(in),dimension(:) :: subState0 real(pReal), intent(in) :: Delta_t integer, intent(in) :: & el, & !< element index in element loop @@ -982,7 +1078,7 @@ function integrateStateFPI(F_0,F,Delta_t,co,ip,el) result(broken) if(broken) return sizeDotState = plasticState(ph)%sizeDotState - plasticState(ph)%state(1:sizeDotState,me) = plasticState(ph)%subState0(1:sizeDotState,me) & + plasticState(ph)%state(1:sizeDotState,me) = subState0 & + plasticState(ph)%dotState (1:sizeDotState,me) * Delta_t dotState(1:sizeDotState,2) = 0.0_pReal @@ -991,7 +1087,7 @@ function integrateStateFPI(F_0,F,Delta_t,co,ip,el) result(broken) if(nIterationState > 1) dotState(1:sizeDotState,2) = dotState(1:sizeDotState,1) dotState(1:sizeDotState,1) = plasticState(ph)%dotState(:,me) - broken = integrateStress(F,Delta_t,co,ip,el) + broken = integrateStress(F,subFp0,subFi0,Delta_t,co,ip,el) if(broken) exit iteration broken = mech_collectDotState(Delta_t, co,ip,el,ph,me) @@ -1002,13 +1098,12 @@ function integrateStateFPI(F_0,F,Delta_t,co,ip,el) result(broken) plasticState(ph)%dotState(:,me) = plasticState(ph)%dotState(:,me) * zeta & + dotState(1:sizeDotState,1) * (1.0_pReal - zeta) r(1:sizeDotState) = plasticState(ph)%state (1:sizeDotState,me) & - - plasticState(ph)%subState0(1:sizeDotState,me) & + - subState0 & - plasticState(ph)%dotState (1:sizeDotState,me) * Delta_t plasticState(ph)%state(1:sizeDotState,me) = plasticState(ph)%state(1:sizeDotState,me) & - r(1:sizeDotState) if (converged(r(1:sizeDotState),plasticState(ph)%state(1:sizeDotState,me),plasticState(ph)%atol(1:sizeDotState))) then - broken = constitutive_deltaState(crystallite_S(1:3,1:3,co,ip,el), & - constitutive_mech_Fi(ph)%data(1:3,1:3,me),co,ip,el,ph,me) + broken = constitutive_deltaState(co,ip,el,ph,me) exit iteration endif @@ -1043,9 +1138,10 @@ end function integrateStateFPI !-------------------------------------------------------------------------------------------------- !> @brief integrate state with 1st order explicit Euler method !-------------------------------------------------------------------------------------------------- -function integrateStateEuler(F_0,F,Delta_t,co,ip,el) result(broken) +function integrateStateEuler(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el) result(broken) - real(pReal), intent(in),dimension(3,3) :: F_0,F + real(pReal), intent(in),dimension(3,3) :: F_0,F,subFp0,subFi0 + real(pReal), intent(in),dimension(:) :: subState0 real(pReal), intent(in) :: Delta_t integer, intent(in) :: & el, & !< element index in element loop @@ -1067,14 +1163,13 @@ function integrateStateEuler(F_0,F,Delta_t,co,ip,el) result(broken) if(broken) return sizeDotState = plasticState(ph)%sizeDotState - plasticState(ph)%state(1:sizeDotState,me) = plasticState(ph)%subState0(1:sizeDotState,me) & - + plasticState(ph)%dotState (1:sizeDotState,me) * Delta_t + plasticState(ph)%state(1:sizeDotState,me) = subState0 & + + plasticState(ph)%dotState(1:sizeDotState,me) * Delta_t - broken = constitutive_deltaState(crystallite_S(1:3,1:3,co,ip,el), & - constitutive_mech_Fi(ph)%data(1:3,1:3,me),co,ip,el,ph,me) + broken = constitutive_deltaState(co,ip,el,ph,me) if(broken) return - broken = integrateStress(F,Delta_t,co,ip,el) + broken = integrateStress(F,subFp0,subFi0,Delta_t,co,ip,el) end function integrateStateEuler @@ -1082,9 +1177,10 @@ end function integrateStateEuler !-------------------------------------------------------------------------------------------------- !> @brief integrate stress, state with 1st order Euler method with adaptive step size !-------------------------------------------------------------------------------------------------- -function integrateStateAdaptiveEuler(F_0,F,Delta_t,co,ip,el) result(broken) +function integrateStateAdaptiveEuler(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el) result(broken) - real(pReal), intent(in),dimension(3,3) :: F_0,F + real(pReal), intent(in),dimension(3,3) :: F_0,F,subFp0,subFi0 + real(pReal), intent(in),dimension(:) :: subState0 real(pReal), intent(in) :: Delta_t integer, intent(in) :: & el, & !< element index in element loop @@ -1109,14 +1205,13 @@ function integrateStateAdaptiveEuler(F_0,F,Delta_t,co,ip,el) result(broken) sizeDotState = plasticState(ph)%sizeDotState residuum_plastic(1:sizeDotState) = - plasticState(ph)%dotstate(1:sizeDotState,me) * 0.5_pReal * Delta_t - plasticState(ph)%state(1:sizeDotState,me) = plasticState(ph)%subState0(1:sizeDotState,me) & + plasticState(ph)%state(1:sizeDotState,me) = subState0 & + plasticState(ph)%dotstate(1:sizeDotState,me) * Delta_t - broken = constitutive_deltaState(crystallite_S(1:3,1:3,co,ip,el), & - constitutive_mech_Fi(ph)%data(1:3,1:3,me),co,ip,el,ph,me) + broken = constitutive_deltaState(co,ip,el,ph,me) if(broken) return - broken = integrateStress(F,Delta_t,co,ip,el) + broken = integrateStress(F,subFp0,subFi0,Delta_t,co,ip,el) if(broken) return broken = mech_collectDotState(Delta_t, co,ip,el,ph,me) @@ -1132,9 +1227,10 @@ end function integrateStateAdaptiveEuler !--------------------------------------------------------------------------------------------------- !> @brief Integrate state (including stress integration) with the classic Runge Kutta method !--------------------------------------------------------------------------------------------------- -function integrateStateRK4(F_0,F,Delta_t,co,ip,el) result(broken) +function integrateStateRK4(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el) result(broken) - real(pReal), intent(in),dimension(3,3) :: F_0,F + real(pReal), intent(in),dimension(3,3) :: F_0,F,subFp0,subFi0 + real(pReal), intent(in),dimension(:) :: subState0 real(pReal), intent(in) :: Delta_t integer, intent(in) :: co,ip,el logical :: broken @@ -1151,7 +1247,7 @@ function integrateStateRK4(F_0,F,Delta_t,co,ip,el) result(broken) B = [1.0_pReal/6.0_pReal, 1.0_pReal/3.0_pReal, 1.0_pReal/3.0_pReal, 1.0_pReal/6.0_pReal] - broken = integrateStateRK(F_0,F,Delta_t,co,ip,el,A,B,C) + broken = integrateStateRK(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el,A,B,C) end function integrateStateRK4 @@ -1159,9 +1255,10 @@ end function integrateStateRK4 !--------------------------------------------------------------------------------------------------- !> @brief Integrate state (including stress integration) with the Cash-Carp method !--------------------------------------------------------------------------------------------------- -function integrateStateRKCK45(F_0,F,Delta_t,co,ip,el) result(broken) +function integrateStateRKCK45(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el) result(broken) - real(pReal), intent(in),dimension(3,3) :: F_0,F + real(pReal), intent(in),dimension(3,3) :: F_0,F,subFp0,subFi0 + real(pReal), intent(in),dimension(:) :: subState0 real(pReal), intent(in) :: Delta_t integer, intent(in) :: co,ip,el logical :: broken @@ -1185,7 +1282,7 @@ function integrateStateRKCK45(F_0,F,Delta_t,co,ip,el) result(broken) 13525.0_pReal/55296.0_pReal, 277.0_pReal/14336.0_pReal, 1._pReal/4._pReal] - broken = integrateStateRK(F_0,F,Delta_t,co,ip,el,A,B,C,DB) + broken = integrateStateRK(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el,A,B,C,DB) end function integrateStateRKCK45 @@ -1194,9 +1291,10 @@ end function integrateStateRKCK45 !> @brief Integrate state (including stress integration) with an explicit Runge-Kutta method or an !! embedded explicit Runge-Kutta method !-------------------------------------------------------------------------------------------------- -function integrateStateRK(F_0,F,Delta_t,co,ip,el,A,B,C,DB) result(broken) +function integrateStateRK(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el,A,B,C,DB) result(broken) - real(pReal), intent(in),dimension(3,3) :: F_0,F + real(pReal), intent(in),dimension(3,3) :: F_0,F,subFp0,subFi0 + real(pReal), intent(in),dimension(:) :: subState0 real(pReal), intent(in) :: Delta_t real(pReal), dimension(:,:), intent(in) :: A real(pReal), dimension(:), intent(in) :: B, C @@ -1234,10 +1332,10 @@ function integrateStateRK(F_0,F,Delta_t,co,ip,el,A,B,C,DB) result(broken) + A(n,stage) * plastic_RKdotState(1:sizeDotState,n) enddo - plasticState(ph)%state(1:sizeDotState,me) = plasticState(ph)%subState0(1:sizeDotState,me) & + plasticState(ph)%state(1:sizeDotState,me) = subState0 & + plasticState(ph)%dotState (1:sizeDotState,me) * Delta_t - broken = integrateStress(F_0 + (F - F_0) * Delta_t * C(stage),Delta_t * C(stage),co,ip,el) + broken = integrateStress(F_0 + (F - F_0) * Delta_t * C(stage),subFp0,subFi0,Delta_t * C(stage),co,ip,el) if(broken) exit broken = mech_collectDotState(Delta_t*C(stage),co,ip,el,ph,me) @@ -1249,7 +1347,7 @@ function integrateStateRK(F_0,F,Delta_t,co,ip,el,A,B,C,DB) result(broken) plastic_RKdotState(1:sizeDotState,size(B)) = plasticState (ph)%dotState(:,me) plasticState(ph)%dotState(:,me) = matmul(plastic_RKdotState(1:sizeDotState,1:size(B)),B) - plasticState(ph)%state(1:sizeDotState,me) = plasticState(ph)%subState0(1:sizeDotState,me) & + plasticState(ph)%state(1:sizeDotState,me) = subState0 & + plasticState(ph)%dotState (1:sizeDotState,me) * Delta_t if(present(DB)) & @@ -1259,11 +1357,10 @@ function integrateStateRK(F_0,F,Delta_t,co,ip,el,A,B,C,DB) result(broken) if(broken) return - broken = constitutive_deltaState(crystallite_S(1:3,1:3,co,ip,el), & - constitutive_mech_Fi(ph)%data(1:3,1:3,me),co,ip,el,ph,me) + broken = constitutive_deltaState(co,ip,el,ph,me) if(broken) return - broken = integrateStress(F,Delta_t,co,ip,el) + broken = integrateStress(F,subFp0,subFi0,Delta_t,co,ip,el) end function integrateStateRK @@ -1288,33 +1385,28 @@ subroutine crystallite_results(group,ph) select case (output_constituent(ph)%label(ou)) case('F') - selected_tensors = select_tensors(crystallite_F,ph) - call results_writeDataset(group//'/mechanics/',selected_tensors,output_constituent(ph)%label(ou),& + call results_writeDataset(group//'/mechanics/',constitutive_mech_F(ph)%data,'F',& 'deformation gradient','1') case('F_e') - selected_tensors = select_tensors(crystallite_Fe,ph) - call results_writeDataset(group//'/mechanics/',selected_tensors,output_constituent(ph)%label(ou),& + call results_writeDataset(group//'/mechanics/',constitutive_mech_Fe(ph)%data,'F_e',& 'elastic deformation gradient','1') case('F_p') - call results_writeDataset(group//'/mechanics/',constitutive_mech_Fp(ph)%data,output_constituent(ph)%label(ou),& + call results_writeDataset(group//'/mechanics/',constitutive_mech_Fp(ph)%data,'F_p', & 'plastic deformation gradient','1') case('F_i') - call results_writeDataset(group//'/mechanics/',constitutive_mech_Fi(ph)%data,output_constituent(ph)%label(ou),& + call results_writeDataset(group//'/mechanics/',constitutive_mech_Fi(ph)%data,'F_i', & 'inelastic deformation gradient','1') case('L_p') - selected_tensors = select_tensors(crystallite_Lp,ph) - call results_writeDataset(group//'/mechanics/',selected_tensors,output_constituent(ph)%label(ou),& + call results_writeDataset(group//'/mechanics/',constitutive_mech_Lp(ph)%data,'L_p', & 'plastic velocity gradient','1/s') case('L_i') - call results_writeDataset(group//'/mechanics/',constitutive_mech_Li(ph)%data,output_constituent(ph)%label(ou),& + call results_writeDataset(group//'/mechanics/',constitutive_mech_Li(ph)%data,'L_i', & 'inelastic velocity gradient','1/s') case('P') - selected_tensors = select_tensors(crystallite_P,ph) - call results_writeDataset(group//'/mechanics/',selected_tensors,output_constituent(ph)%label(ou),& + call results_writeDataset(group//'/mechanics/',constitutive_mech_P(ph)%data,'P', & 'First Piola-Kirchhoff stress','Pa') case('S') - selected_tensors = select_tensors(crystallite_S,ph) - call results_writeDataset(group//'/mechanics/',selected_tensors,output_constituent(ph)%label(ou),& + call results_writeDataset(group//'/mechanics/',constitutive_mech_S(ph)%data,'S', & 'Second Piola-Kirchhoff stress','Pa') case('O') select case(lattice_structure(ph)) @@ -1341,33 +1433,6 @@ subroutine crystallite_results(group,ph) contains - !------------------------------------------------------------------------------------------------ - !> @brief select tensors for output - !------------------------------------------------------------------------------------------------ - function select_tensors(dataset,ph) - - integer, intent(in) :: ph - real(pReal), dimension(:,:,:,:,:), intent(in) :: dataset - real(pReal), allocatable, dimension(:,:,:) :: select_tensors - integer :: el,ip,co,j - - allocate(select_tensors(3,3,count(material_phaseAt==ph)*discretization_nIPs)) - - j=0 - do el = 1, size(material_phaseAt,2) - do ip = 1, discretization_nIPs - do co = 1, size(material_phaseAt,1) !ToDo: this needs to be changed for varying Ngrains - if (material_phaseAt(co,el) == ph) then - j = j + 1 - select_tensors(1:3,1:3,j) = dataset(1:3,1:3,co,ip,el) - endif - enddo - enddo - enddo - - end function select_tensors - - !-------------------------------------------------------------------------------------------------- !> @brief select rotations for output !-------------------------------------------------------------------------------------------------- @@ -1407,7 +1472,11 @@ module subroutine mech_initializeRestorationPoints(ph,me) constitutive_mech_partitionedFi0(ph)%data(1:3,1:3,me) = constitutive_mech_Fi0(ph)%data(1:3,1:3,me) constitutive_mech_partitionedFp0(ph)%data(1:3,1:3,me) = constitutive_mech_Fp0(ph)%data(1:3,1:3,me) + constitutive_mech_partitionedF0(ph)%data(1:3,1:3,me) = constitutive_mech_F0(ph)%data(1:3,1:3,me) constitutive_mech_partitionedLi0(ph)%data(1:3,1:3,me) = constitutive_mech_Li0(ph)%data(1:3,1:3,me) + constitutive_mech_partitionedLp0(ph)%data(1:3,1:3,me) = constitutive_mech_Lp0(ph)%data(1:3,1:3,me) + constitutive_mech_partitionedS0(ph)%data(1:3,1:3,me) = constitutive_mech_S0(ph)%data(1:3,1:3,me) + plasticState(ph)%partitionedState0(:,me) = plasticState(ph)%state0(:,me) end subroutine mech_initializeRestorationPoints @@ -1423,7 +1492,10 @@ module subroutine constitutive_mech_windForward(ph,me) constitutive_mech_partitionedFp0(ph)%data(1:3,1:3,me) = constitutive_mech_Fp(ph)%data(1:3,1:3,me) constitutive_mech_partitionedFi0(ph)%data(1:3,1:3,me) = constitutive_mech_Fi(ph)%data(1:3,1:3,me) + constitutive_mech_partitionedF0(ph)%data(1:3,1:3,me) = constitutive_mech_F(ph)%data(1:3,1:3,me) constitutive_mech_partitionedLi0(ph)%data(1:3,1:3,me) = constitutive_mech_Li(ph)%data(1:3,1:3,me) + constitutive_mech_partitionedLp0(ph)%data(1:3,1:3,me) = constitutive_mech_Lp(ph)%data(1:3,1:3,me) + constitutive_mech_partitionedS0(ph)%data(1:3,1:3,me) = constitutive_mech_S(ph)%data(1:3,1:3,me) plasticState(ph)%partitionedState0(:,me) = plasticState(ph)%state(:,me) @@ -1440,10 +1512,13 @@ module subroutine constitutive_mech_forward() do ph = 1, size(plasticState) - plasticState(ph)%state0 = plasticState(ph)%state constitutive_mech_Fi0(ph) = constitutive_mech_Fi(ph) constitutive_mech_Fp0(ph) = constitutive_mech_Fp(ph) + constitutive_mech_F0(ph) = constitutive_mech_F(ph) constitutive_mech_Li0(ph) = constitutive_mech_Li(ph) + constitutive_mech_Lp0(ph) = constitutive_mech_Lp(ph) + constitutive_mech_S0(ph) = constitutive_mech_S(ph) + plasticState(ph)%state0 = plasticState(ph)%state enddo end subroutine constitutive_mech_forward @@ -1487,41 +1562,42 @@ module function crystallite_stress(dt,co,ip,el) result(converged_) real(pReal) :: & formerSubStep integer :: & - NiterationCrystallite, & ! number of iterations in crystallite loop - so, ph, me + so, ph, me, sizeDotState logical :: todo real(pReal) :: subFrac,subStep real(pReal), dimension(3,3) :: & - subLp0, & !< plastic velocity grad at start of crystallite inc - subLi0, & !< intermediate velocity grad at start of crystallite inc + subFp0, & + subFi0, & + subLp0, & + subLi0, & subF0, & subF + real(pReal), dimension(:), allocatable :: subState0 ph = material_phaseAt(co,el) me = material_phaseMemberAt(co,ip,el) - subLi0 = constitutive_mech_partitionedLi0(ph)%data(1:3,1:3,me) - subLp0 = crystallite_partitionedLp0(1:3,1:3,co,ip,el) + sizeDotState = plasticState(ph)%sizeDotState + + subLi0 = constitutive_mech_partitionedLi0(ph)%data(1:3,1:3,me) + subLp0 = constitutive_mech_partitionedLp0(ph)%data(1:3,1:3,me) + subState0 = plasticState(ph)%partitionedState0(:,me) + - plasticState(ph)%subState0(:,me) = plasticState(ph)%partitionedState0(:,me) do so = 1, phase_Nsources(ph) sourceState(ph)%p(so)%subState0(:,me) = sourceState(ph)%p(so)%partitionedState0(:,me) enddo - crystallite_subFp0(1:3,1:3,co,ip,el) = constitutive_mech_partitionedFp0(ph)%data(1:3,1:3,me) - crystallite_subFi0(1:3,1:3,co,ip,el) = constitutive_mech_partitionedFi0(ph)%data(1:3,1:3,me) - subF0 = crystallite_partitionedF0(1:3,1:3,co,ip,el) + subFp0 = constitutive_mech_partitionedFp0(ph)%data(1:3,1:3,me) + subFi0 = constitutive_mech_partitionedFi0(ph)%data(1:3,1:3,me) + subF0 = constitutive_mech_partitionedF0(ph)%data(1:3,1:3,me) subFrac = 0.0_pReal subStep = 1.0_pReal/num%subStepSizeCryst todo = .true. converged_ = .false. ! pretend failed step of 1/subStepSizeCryst todo = .true. - NiterationCrystallite = 0 cutbackLooping: do while (todo) - NiterationCrystallite = NiterationCrystallite + 1 -!-------------------------------------------------------------------------------------------------- -! wind forward if (converged_) then formerSubStep = subStep subFrac = subFrac + subStep @@ -1531,11 +1607,11 @@ module function crystallite_stress(dt,co,ip,el) result(converged_) if (todo) then subF0 = subF - subLp0 = crystallite_Lp (1:3,1:3,co,ip,el) + subLp0 = constitutive_mech_Lp(ph)%data(1:3,1:3,me) subLi0 = constitutive_mech_Li(ph)%data(1:3,1:3,me) - crystallite_subFp0(1:3,1:3,co,ip,el) = constitutive_mech_Fp(ph)%data(1:3,1:3,me) - crystallite_subFi0(1:3,1:3,co,ip,el) = constitutive_mech_Fi(ph)%data(1:3,1:3,me) - plasticState(ph)%subState0(:,me) = plasticState(ph)%state(:,me) + subFp0 = constitutive_mech_Fp(ph)%data(1:3,1:3,me) + subFi0 = constitutive_mech_Fi(ph)%data(1:3,1:3,me) + subState0 = plasticState(ph)%state(:,me) do so = 1, phase_Nsources(ph) sourceState(ph)%p(so)%subState0(:,me) = sourceState(ph)%p(so)%state(:,me) enddo @@ -1544,14 +1620,14 @@ module function crystallite_stress(dt,co,ip,el) result(converged_) ! cut back (reduced time and restore) else subStep = num%subStepSizeCryst * subStep - constitutive_mech_Fp(ph)%data(1:3,1:3,me) = crystallite_subFp0(1:3,1:3,co,ip,el) - constitutive_mech_Fi(ph)%data(1:3,1:3,me) = crystallite_subFi0(1:3,1:3,co,ip,el) - crystallite_S (1:3,1:3,co,ip,el) = crystallite_S0 (1:3,1:3,co,ip,el) - if (subStep < 1.0_pReal) then ! actual (not initial) cutback - crystallite_Lp (1:3,1:3,co,ip,el) = subLp0 - constitutive_mech_Li(ph)%data(1:3,1:3,me) = subLi0 + constitutive_mech_Fp(ph)%data(1:3,1:3,me) = subFp0 + constitutive_mech_Fi(ph)%data(1:3,1:3,me) = subFi0 + constitutive_mech_S(ph)%data(1:3,1:3,me) = constitutive_mech_S0(ph)%data(1:3,1:3,me) ! why no subS0 ? is S0 of any use? + if (subStep < 1.0_pReal) then ! actual (not initial) cutback + constitutive_mech_Lp(ph)%data(1:3,1:3,me) = subLp0 + constitutive_mech_Li(ph)%data(1:3,1:3,me) = subLi0 endif - plasticState(ph)%state(:,me) = plasticState(ph)%subState0(:,me) + plasticState(ph)%state(:,me) = subState0 do so = 1, phase_Nsources(ph) sourceState(ph)%p(so)%state(:,me) = sourceState(ph)%p(so)%subState0(:,me) enddo @@ -1563,11 +1639,10 @@ module function crystallite_stress(dt,co,ip,el) result(converged_) ! prepare for integration if (todo) then subF = subF0 & - + subStep * (crystallite_F(1:3,1:3,co,ip,el) - crystallite_partitionedF0(1:3,1:3,co,ip,el)) - crystallite_Fe(1:3,1:3,co,ip,el) = matmul(subF,math_inv33(matmul(constitutive_mech_Fi(ph)%data(1:3,1:3,me), & - constitutive_mech_Fp(ph)%data(1:3,1:3,me)))) - crystallite_subdt(co,ip,el) = subStep * dt - converged_ = .not. integrateState(subF0,subF,subStep * dt,co,ip,el) + + subStep * (constitutive_mech_F(ph)%data(1:3,1:3,me) - constitutive_mech_partitionedF0(ph)%data(1:3,1:3,me)) + constitutive_mech_Fe(ph)%data(1:3,1:3,me) = matmul(subF,math_inv33(matmul(constitutive_mech_Fi(ph)%data(1:3,1:3,me), & + constitutive_mech_Fp(ph)%data(1:3,1:3,me)))) + converged_ = .not. integrateState(subF0,subF,subFp0,subFi0,subState0(1:sizeDotState),subStep * dt,co,ip,el) converged_ = converged_ .and. .not. integrateSourceState(subStep * dt,co,ip,el) endif @@ -1586,20 +1661,22 @@ module subroutine mech_restore(ip,el,includeL) el !< element number logical, intent(in) :: & includeL !< protect agains fake cutback + integer :: & - co, p, m !< constituent number + co, ph, me + do co = 1,homogenization_Nconstituents(material_homogenizationAt(el)) - p = material_phaseAt(co,el) - m = material_phaseMemberAt(co,ip,el) + ph = material_phaseAt(co,el) + me = material_phaseMemberAt(co,ip,el) if (includeL) then - crystallite_Lp(1:3,1:3,co,ip,el) = crystallite_partitionedLp0(1:3,1:3,co,ip,el) - constitutive_mech_Li(p)%data(1:3,1:3,m) = constitutive_mech_partitionedLi0(p)%data(1:3,1:3,m) + constitutive_mech_Lp(ph)%data(1:3,1:3,me) = constitutive_mech_partitionedLp0(ph)%data(1:3,1:3,me) + constitutive_mech_Li(ph)%data(1:3,1:3,me) = constitutive_mech_partitionedLi0(ph)%data(1:3,1:3,me) endif ! maybe protecting everything from overwriting makes more sense - constitutive_mech_Fp(p)%data(1:3,1:3,m) = constitutive_mech_partitionedFp0(p)%data(1:3,1:3,m) - constitutive_mech_Fi(p)%data(1:3,1:3,m) = constitutive_mech_partitionedFi0(p)%data(1:3,1:3,m) - crystallite_S (1:3,1:3,co,ip,el) = crystallite_partitionedS0 (1:3,1:3,co,ip,el) + constitutive_mech_Fp(ph)%data(1:3,1:3,me) = constitutive_mech_partitionedFp0(ph)%data(1:3,1:3,me) + constitutive_mech_Fi(ph)%data(1:3,1:3,me) = constitutive_mech_partitionedFi0(ph)%data(1:3,1:3,me) + constitutive_mech_S(ph)%data(1:3,1:3,me) = constitutive_mech_partitionedS0(ph)%data(1:3,1:3,me) plasticState (material_phaseAt(co,el))%state( :,material_phasememberAt(co,ip,el)) = & plasticState (material_phaseAt(co,el))%partitionedState0(:,material_phasememberAt(co,ip,el)) @@ -1607,5 +1684,239 @@ 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 + + +module subroutine mech_restartWrite(groupHandle,ph) + + integer(HID_T), intent(in) :: groupHandle + integer, intent(in) :: ph + + + call HDF5_write(groupHandle,plasticState(ph)%state,'omega') + call HDF5_write(groupHandle,constitutive_mech_Fi(ph)%data,'F_i') + call HDF5_write(groupHandle,constitutive_mech_Li(ph)%data,'L_i') + call HDF5_write(groupHandle,constitutive_mech_Lp(ph)%data,'L_p') + call HDF5_write(groupHandle,constitutive_mech_Fp(ph)%data,'F_p') + call HDF5_write(groupHandle,constitutive_mech_S(ph)%data,'S') + call HDF5_write(groupHandle,constitutive_mech_F(ph)%data,'F') + +end subroutine mech_restartWrite + + +module subroutine mech_restartRead(groupHandle,ph) + + integer(HID_T), intent(in) :: groupHandle + integer, intent(in) :: ph + + + call HDF5_read(groupHandle,plasticState(ph)%state0,'omega') + call HDF5_read(groupHandle,constitutive_mech_Fi0(ph)%data,'F_i') + call HDF5_read(groupHandle,constitutive_mech_Li0(ph)%data,'L_i') + call HDF5_read(groupHandle,constitutive_mech_Lp0(ph)%data,'L_p') + call HDF5_read(groupHandle,constitutive_mech_Fp0(ph)%data,'F_p') + call HDF5_read(groupHandle,constitutive_mech_S0(ph)%data,'S') + call HDF5_read(groupHandle,constitutive_mech_F0(ph)%data,'F') + +end subroutine mech_restartRead + + +! getter for non-mech (e.g. thermal) +module function constitutive_mech_getS(co,ip,el) result(S) + + integer, intent(in) :: co, ip, el + real(pReal), dimension(3,3) :: S + + + S = constitutive_mech_S(material_phaseAt(co,el))%data(1:3,1:3,material_phaseMemberAt(co,ip,el)) + +end function constitutive_mech_getS + + +! getter for non-mech (e.g. thermal) +module function constitutive_mech_getLp(co,ip,el) result(Lp) + + integer, intent(in) :: co, ip, el + real(pReal), dimension(3,3) :: Lp + + + Lp = constitutive_mech_Lp(material_phaseAt(co,el))%data(1:3,1:3,material_phaseMemberAt(co,ip,el)) + +end function constitutive_mech_getLp + + +! getter for non-mech (e.g. thermal) +module function constitutive_mech_getF(co,ip,el) result(F) + + integer, intent(in) :: co, ip, el + real(pReal), dimension(3,3) :: F + + + F = constitutive_mech_F(material_phaseAt(co,el))%data(1:3,1:3,material_phaseMemberAt(co,ip,el)) + +end function constitutive_mech_getF + + +! getter for non-mech (e.g. thermal) +module function constitutive_mech_getF_e(co,ip,el) result(F_e) + + integer, intent(in) :: co, ip, el + real(pReal), dimension(3,3) :: F_e + + + F_e = constitutive_mech_Fe(material_phaseAt(co,el))%data(1:3,1:3,material_phaseMemberAt(co,ip,el)) + +end function constitutive_mech_getF_e + + + +! getter for non-mech (e.g. thermal) +module function constitutive_mech_getP(co,ip,el) result(P) + + integer, intent(in) :: co, ip, el + real(pReal), dimension(3,3) :: P + + + P = constitutive_mech_P(material_phaseAt(co,el))%data(1:3,1:3,material_phaseMemberAt(co,ip,el)) + +end function constitutive_mech_getP + + +! setter for homogenization +module subroutine constitutive_mech_setF(F,co,ip,el) + + real(pReal), dimension(3,3), intent(in) :: F + integer, intent(in) :: co, ip, el + + + constitutive_mech_F(material_phaseAt(co,el))%data(1:3,1:3,material_phaseMemberAt(co,ip,el)) = F + +end subroutine constitutive_mech_setF + end submodule constitutive_mech diff --git a/src/constitutive_plastic_nonlocal.f90 b/src/constitutive_plastic_nonlocal.f90 index 0d7875291..2244eb7ad 100644 --- a/src/constitutive_plastic_nonlocal.f90 +++ b/src/constitutive_plastic_nonlocal.f90 @@ -552,10 +552,8 @@ end function plastic_nonlocal_init !-------------------------------------------------------------------------------------------------- !> @brief calculates quantities characterizing the microstructure !-------------------------------------------------------------------------------------------------- -module subroutine plastic_nonlocal_dependentState(F, instance, of, ip, el) +module subroutine plastic_nonlocal_dependentState(instance, of, ip, el) - real(pReal), dimension(3,3), intent(in) :: & - F integer, intent(in) :: & instance, & of, & @@ -647,7 +645,7 @@ module subroutine plastic_nonlocal_dependentState(F, instance, of, ip, el) ph = material_phaseAt(1,el) me = material_phaseMemberAt(1,ip,el) invFp = math_inv33(constitutive_mech_Fp(ph)%data(1:3,1:3,me)) - invFe = matmul(constitutive_mech_Fp(ph)%data(1:3,1:3,me),math_inv33(F)) + invFe = math_inv33(constitutive_mech_Fe(ph)%data(1:3,1:3,me)) rho_edg_delta = rho0(:,mob_edg_pos) - rho0(:,mob_edg_neg) rho_scr_delta = rho0(:,mob_scr_pos) - rho0(:,mob_scr_neg) @@ -976,13 +974,11 @@ end subroutine plastic_nonlocal_deltaState !--------------------------------------------------------------------------------------------------- !> @brief calculates the rate of change of microstructure !--------------------------------------------------------------------------------------------------- -module subroutine plastic_nonlocal_dotState(Mp, F, Temperature,timestep, & +module subroutine plastic_nonlocal_dotState(Mp, Temperature,timestep, & instance,of,ip,el) real(pReal), dimension(3,3), intent(in) :: & Mp !< MandelStress - real(pReal), dimension(3,3,homogenization_maxNconstituents,discretization_nIPs,discretization_Nelems), intent(in) :: & - F !< Deformation gradient real(pReal), intent(in) :: & Temperature, & !< temperature timestep !< substepped crystallite time increment @@ -1149,7 +1145,7 @@ module subroutine plastic_nonlocal_dotState(Mp, F, Temperature,timestep, & - rhoDip(s,1) / timestep - rhoDotAthermalAnnihilation(s,9) & - rhoDotSingle2DipoleGlide(s,9)) ! make sure that we do not annihilate more dipoles than we have - rhoDot = rhoDotFlux(F,timestep, instance,of,ip,el) & + rhoDot = rhoDotFlux(timestep, instance,of,ip,el) & + rhoDotMultiplication & + rhoDotSingle2DipoleGlide & + rhoDotAthermalAnnihilation & @@ -1178,10 +1174,8 @@ end subroutine plastic_nonlocal_dotState !--------------------------------------------------------------------------------------------------- !> @brief calculates the rate of change of microstructure !--------------------------------------------------------------------------------------------------- -function rhoDotFlux(F,timestep, instance,of,ip,el) +function rhoDotFlux(timestep,instance,of,ip,el) - real(pReal), dimension(3,3,homogenization_maxNconstituents,discretization_nIPs,discretization_Nelems), intent(in) :: & - F !< Deformation gradient real(pReal), intent(in) :: & timestep !< substepped crystallite time increment integer, intent(in) :: & @@ -1293,7 +1287,7 @@ function rhoDotFlux(F,timestep, instance,of,ip,el) m(1:3,:,3) = -prm%slip_transverse m(1:3,:,4) = prm%slip_transverse - my_F = F(1:3,1:3,1,ip,el) + my_F = constitutive_mech_F(ph)%data(1:3,1:3,of) my_Fe = matmul(my_F, math_inv33(constitutive_mech_Fp(ph)%data(1:3,1:3,of))) neighbors: do n = 1,nIPneighbors @@ -1311,7 +1305,7 @@ function rhoDotFlux(F,timestep, instance,of,ip,el) if (neighbor_n > 0) then ! if neighbor exists, average deformation gradient neighbor_instance = phase_plasticityInstance(material_phaseAt(1,neighbor_el)) - neighbor_F = F(1:3,1:3,1,neighbor_ip,neighbor_el) + neighbor_F = constitutive_mech_F(np)%data(1:3,1:3,no) neighbor_Fe = matmul(neighbor_F, math_inv33(constitutive_mech_Fp(np)%data(1:3,1:3,no))) Favg = 0.5_pReal * (my_F + neighbor_F) else ! if no neighbor, take my value as average diff --git a/src/constitutive_thermal.f90 b/src/constitutive_thermal.f90 index a7d5d3259..1a05b983f 100644 --- a/src/constitutive_thermal.f90 +++ b/src/constitutive_thermal.f90 @@ -9,7 +9,7 @@ submodule(constitutive) constitutive_thermal integer, intent(in) :: source_length logical, dimension(:,:), allocatable :: mySources end function source_thermal_dissipation_init - + module function source_thermal_externalheat_init(source_length) result(mySources) integer, intent(in) :: source_length logical, dimension(:,:), allocatable :: mySources @@ -55,8 +55,8 @@ module subroutine thermal_init if(maxval(phase_Nsources) /= 0) then where(source_thermal_dissipation_init (maxval(phase_Nsources))) phase_source = SOURCE_thermal_dissipation_ID where(source_thermal_externalheat_init(maxval(phase_Nsources))) phase_source = SOURCE_thermal_externalheat_ID - endif - + endif + !-------------------------------------------------------------------------------------------------- !initialize kinematic mechanisms if(maxval(phase_Nkinematics) /= 0) where(kinematics_thermal_expansion_init(maxval(phase_Nkinematics))) & @@ -68,15 +68,13 @@ end subroutine thermal_init !---------------------------------------------------------------------------------------------- !< @brief calculates thermal dissipation rate !---------------------------------------------------------------------------------------------- -module subroutine constitutive_thermal_getRateAndItsTangents(TDot, dTDot_dT, T, S, Lp, ip, el) +module subroutine constitutive_thermal_getRateAndItsTangents(TDot, dTDot_dT, T, ip, el) + integer, intent(in) :: & ip, & !< integration point number el !< element number real(pReal), intent(in) :: & - T - real(pReal), intent(in), dimension(:,:,:,:,:) :: & - S, & !< current 2nd Piola Kirchhoff stress - Lp !< plastic velocity gradient + T !< plastic velocity gradient real(pReal), intent(inout) :: & TDot, & dTDot_dT @@ -84,6 +82,7 @@ module subroutine constitutive_thermal_getRateAndItsTangents(TDot, dTDot_dT, T, real(pReal) :: & my_Tdot, & my_dTdot_dT + real(pReal), dimension(3,3) :: Lp, S integer :: & phase, & homog, & @@ -101,10 +100,10 @@ module subroutine constitutive_thermal_getRateAndItsTangents(TDot, dTDot_dT, T, do source = 1, phase_Nsources(phase) select case(phase_source(source,phase)) case (SOURCE_thermal_dissipation_ID) + Lp = constitutive_mech_getLp(grain,ip,el) + S = constitutive_mech_getS(grain,ip,el) call source_thermal_dissipation_getRateAndItsTangent(my_Tdot, my_dTdot_dT, & - S(1:3,1:3,grain,ip,el), & - Lp(1:3,1:3,grain,ip,el), & - phase) + S, Lp, phase) case (SOURCE_thermal_externalheat_ID) call source_thermal_externalheat_getRateAndItsTangent(my_Tdot, my_dTdot_dT, & @@ -122,4 +121,22 @@ module subroutine constitutive_thermal_getRateAndItsTangents(TDot, dTDot_dT, T, end subroutine constitutive_thermal_getRateAndItsTangents + + +! getter for non-thermal (e.g. mech) +module function constitutive_thermal_T(co,ip,el) result(T) + + integer, intent(in) :: co, ip, el + real(pReal) :: T + + integer :: ho, tme + + ho = material_homogenizationAt(el) + tme = material_homogenizationMemberAt(ip,el) + + T = temperature(ho)%p(tme) + +end function constitutive_thermal_T + + end submodule constitutive_thermal diff --git a/src/homogenization.f90 b/src/homogenization.f90 index 52553b57b..686bb9885 100644 --- a/src/homogenization.f90 +++ b/src/homogenization.f90 @@ -16,6 +16,7 @@ module homogenization use thermal_conduction use damage_none use damage_nonlocal + use HDF5_utilities use results implicit none @@ -63,7 +64,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 @@ -91,7 +93,9 @@ module homogenization homogenization_init, & materialpoint_stressAndItsTangent, & homogenization_forward, & - homogenization_results + homogenization_results, & + homogenization_restartRead, & + homogenization_restartWrite contains @@ -168,10 +172,8 @@ subroutine materialpoint_stressAndItsTangent(dt,FEsolving_execIP,FEsolving_execE converged = .false. ! pretend failed step ... subStep = 1.0_pReal/num%subStepSizeHomog ! ... larger then the requested calculation - if (homogState(ho)%sizeState > 0) & - homogState(ho)%subState0(:,me) = homogState(ho)%State0(:,me) - if (damageState(ho)%sizeState > 0) & - damageState(ho)%subState0(:,me) = damageState(ho)%State0(:,me) + if (homogState(ho)%sizeState > 0) homogState(ho)%subState0(:,me) = homogState(ho)%State0(:,me) + if (damageState(ho)%sizeState > 0) damageState(ho)%subState0(:,me) = damageState(ho)%State0(:,me) cutBackLooping: do while (.not. terminallyIll .and. subStep > num%subStepMinHomog) @@ -184,10 +186,8 @@ subroutine materialpoint_stressAndItsTangent(dt,FEsolving_execIP,FEsolving_execE ! wind forward grain starting point call constitutive_windForward(ip,el) - if(homogState(ho)%sizeState > 0) & - homogState(ho)%subState0(:,me) = homogState(ho)%State(:,me) - if(damageState(ho)%sizeState > 0) & - damageState(ho)%subState0(:,me) = damageState(ho)%State(:,me) + if(homogState(ho)%sizeState > 0) homogState(ho)%subState0(:,me) = homogState(ho)%State(:,me) + if(damageState(ho)%sizeState > 0) damageState(ho)%subState0(:,me) = damageState(ho)%State(:,me) endif steppingNeeded elseif ( (myNgrains == 1 .and. subStep <= 1.0 ) .or. & ! single grain already tried internal subStepping in crystallite @@ -201,10 +201,8 @@ subroutine materialpoint_stressAndItsTangent(dt,FEsolving_execIP,FEsolving_execE call constitutive_restore(ip,el,subStep < 1.0_pReal) - if(homogState(ho)%sizeState > 0) & - homogState(ho)%State(:,me) = homogState(ho)%subState0(:,me) - if(damageState(ho)%sizeState > 0) & - damageState(ho)%State(:,me) = damageState(ho)%subState0(:,me) + if(homogState(ho)%sizeState > 0) homogState(ho)%State(:,me) = homogState(ho)%subState0(:,me) + if(damageState(ho)%sizeState > 0) damageState(ho)%State(:,me) = damageState(ho)%subState0(:,me) endif if (subStep > num%subStepMinHomog) doneAndHappy = [.false.,.true.] @@ -257,7 +255,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 @@ -320,4 +318,59 @@ subroutine homogenization_forward end subroutine homogenization_forward + +!-------------------------------------------------------------------------------------------------- +!-------------------------------------------------------------------------------------------------- +subroutine homogenization_restartWrite(fileHandle) + + integer(HID_T), intent(in) :: fileHandle + + integer(HID_T), dimension(2) :: groupHandle + integer :: ho + + + groupHandle(1) = HDF5_addGroup(fileHandle,'homogenization') + + do ho = 1, size(material_name_homogenization) + + groupHandle(2) = HDF5_addGroup(groupHandle(1),material_name_homogenization(ho)) + + call HDF5_read(groupHandle(2),homogState(ho)%state,'omega') ! ToDo: should be done by mech + + call HDF5_closeGroup(groupHandle(2)) + + enddo + + call HDF5_closeGroup(groupHandle(1)) + +end subroutine homogenization_restartWrite + + +!-------------------------------------------------------------------------------------------------- +!-------------------------------------------------------------------------------------------------- +subroutine homogenization_restartRead(fileHandle) + + integer(HID_T), intent(in) :: fileHandle + + integer(HID_T), dimension(2) :: groupHandle + integer :: ho + + + groupHandle(1) = HDF5_openGroup(fileHandle,'homogenization') + + do ho = 1, size(material_name_homogenization) + + groupHandle(2) = HDF5_openGroup(groupHandle(1),material_name_homogenization(ho)) + + call HDF5_write(groupHandle(2),homogState(ho)%state,'omega') ! ToDo: should be done by mech + + call HDF5_closeGroup(groupHandle(2)) + + enddo + + call HDF5_closeGroup(groupHandle(1)) + +end subroutine homogenization_restartRead + + end module homogenization diff --git a/src/homogenization_mech.f90 b/src/homogenization_mech.f90 index 641e960fd..783d08dd1 100644 --- a/src/homogenization_mech.f90 +++ b/src/homogenization_mech.f90 @@ -52,12 +52,11 @@ submodule(homogenization) homogenization_mech end subroutine mech_RGC_averageStressAndItsTangent - module function mech_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el) result(doneAndHappy) + module function mech_RGC_updateState(P,F,avgF,dt,dPdF,ip,el) result(doneAndHappy) logical, dimension(2) :: doneAndHappy real(pReal), dimension(:,:,:), intent(in) :: & P,& !< partitioned stresses - F,& !< partitioned deformation gradients - F0 !< partitioned initial deformation gradients + F !< partitioned deformation gradients real(pReal), dimension(:,:,:,:,:), intent(in) :: dPdF !< partitioned stiffnesses real(pReal), dimension(3,3), intent(in) :: avgF !< average F real(pReal), intent(in) :: dt !< time increment @@ -113,67 +112,73 @@ module subroutine mech_partition(subF,ip,el) ip, & !< integration point el !< element number + integer :: co + real(pReal) :: F(3,3,homogenization_Nconstituents(material_homogenizationAt(el))) + + chosenHomogenization: select case(homogenization_type(material_homogenizationAt(el))) case (HOMOGENIZATION_NONE_ID) chosenHomogenization - crystallite_F(1:3,1:3,1,ip,el) = subF + F(1:3,1:3,1) = subF case (HOMOGENIZATION_ISOSTRAIN_ID) chosenHomogenization - call mech_isostrain_partitionDeformation(& - crystallite_F(1:3,1:3,1:homogenization_Nconstituents(material_homogenizationAt(el)),ip,el), & - subF) + call mech_isostrain_partitionDeformation(F,subF) case (HOMOGENIZATION_RGC_ID) chosenHomogenization - call mech_RGC_partitionDeformation(& - crystallite_F(1:3,1:3,1:homogenization_Nconstituents(material_homogenizationAt(el)),ip,el), & - subF,& - ip, & - el) + call mech_RGC_partitionDeformation(F,subF,ip,el) end select chosenHomogenization + do co = 1,homogenization_Nconstituents(material_homogenizationAt(el)) + call constitutive_mech_setF(F(1:3,1:3,co),co,ip,el) + enddo + + 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))) + real(pReal) :: Ps(3,3,homogenization_Nconstituents(material_homogenizationAt(el))) ce = (el-1)* discretization_nIPs + ip chosenHomogenization: select case(homogenization_type(material_homogenizationAt(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_P(1:3,1:3,ce) = constitutive_mech_getP(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(co,ip,el) + dPdFs(:,:,:,:,co) = constitutive_mech_dPdF(dt,co,ip,el) + Ps(:,:,co) = constitutive_mech_getP(co,ip,el) enddo call mech_isostrain_averageStressAndItsTangent(& homogenization_P(1:3,1:3,ce), & homogenization_dPdF(1:3,1:3,1:3,1:3,ce),& - crystallite_P(1:3,1:3,1:homogenization_Nconstituents(material_homogenizationAt(el)),ip,el), & - dPdFs, & + Ps,dPdFs, & homogenization_typeInstance(material_homogenizationAt(el))) case (HOMOGENIZATION_RGC_ID) chosenHomogenization do co = 1, homogenization_Nconstituents(material_homogenizationAt(el)) - dPdFs(:,:,:,:,co) = crystallite_stressTangent(co,ip,el) + dPdFs(:,:,:,:,co) = constitutive_mech_dPdF(dt,co,ip,el) + Ps(:,:,co) = constitutive_mech_getP(co,ip,el) enddo call mech_RGC_averageStressAndItsTangent(& homogenization_P(1:3,1:3,ce), & homogenization_dPdF(1:3,1:3,1:3,1:3,ce),& - crystallite_P(1:3,1:3,1:homogenization_Nconstituents(material_homogenizationAt(el)),ip,el), & - dPdFs, & + Ps,dPdFs, & homogenization_typeInstance(material_homogenizationAt(el))) end select chosenHomogenization @@ -198,21 +203,17 @@ module function mech_updateState(subdt,subF,ip,el) result(doneAndHappy) integer :: co real(pReal) :: dPdFs(3,3,3,3,homogenization_Nconstituents(material_homogenizationAt(el))) + real(pReal) :: Fs(3,3,homogenization_Nconstituents(material_homogenizationAt(el))) + real(pReal) :: Ps(3,3,homogenization_Nconstituents(material_homogenizationAt(el))) 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) = constitutive_mech_dPdF(subdt,co,ip,el) + Fs(:,:,co) = constitutive_mech_getF(co,ip,el) + Ps(:,:,co) = constitutive_mech_getP(co,ip,el) enddo - doneAndHappy = & - mech_RGC_updateState(crystallite_P(1:3,1:3,1:homogenization_Nconstituents(material_homogenizationAt(el)),ip,el), & - crystallite_F(1:3,1:3,1:homogenization_Nconstituents(material_homogenizationAt(el)),ip,el), & - crystallite_partitionedF0(1:3,1:3,1:homogenization_Nconstituents(material_homogenizationAt(el)),ip,el),& - subF,& - subdt, & - dPdFs, & - ip, & - el) + doneAndHappy = mech_RGC_updateState(Ps,Fs,subF,subdt,dPdFs,ip,el) else doneAndHappy = .true. endif diff --git a/src/homogenization_mech_RGC.f90 b/src/homogenization_mech_RGC.f90 index 04ec73845..580bb0268 100644 --- a/src/homogenization_mech_RGC.f90 +++ b/src/homogenization_mech_RGC.f90 @@ -24,9 +24,6 @@ submodule(homogenization:homogenization_mech) homogenization_mech_RGC end type tParameters type :: tRGCstate - real(pReal), pointer, dimension(:) :: & - work, & - penaltyEnergy real(pReal), pointer, dimension(:,:) :: & relaxationVector end type tRGCstate @@ -170,8 +167,7 @@ module subroutine mech_RGC_init(num_homogMech) nIntFaceTot = 3*( (prm%N_constituents(1)-1)*prm%N_constituents(2)*prm%N_constituents(3) & + prm%N_constituents(1)*(prm%N_constituents(2)-1)*prm%N_constituents(3) & + prm%N_constituents(1)*prm%N_constituents(2)*(prm%N_constituents(3)-1)) - sizeState = nIntFaceTot & - + size(['avg constitutive work ','average penalty energy']) + sizeState = nIntFaceTot homogState(h)%sizeState = sizeState allocate(homogState(h)%state0 (sizeState,Nmaterialpoints), source=0.0_pReal) @@ -180,8 +176,6 @@ module subroutine mech_RGC_init(num_homogMech) stt%relaxationVector => homogState(h)%state(1:nIntFaceTot,:) st0%relaxationVector => homogState(h)%state0(1:nIntFaceTot,:) - stt%work => homogState(h)%state(nIntFaceTot+1,:) - stt%penaltyEnergy => homogState(h)%state(nIntFaceTot+2,:) allocate(dst%volumeDiscrepancy( Nmaterialpoints), source=0.0_pReal) allocate(dst%relaxationRate_avg( Nmaterialpoints), source=0.0_pReal) @@ -243,12 +237,11 @@ end subroutine mech_RGC_partitionDeformation !> @brief update the internal state of the homogenization scheme and tell whether "done" and ! "happy" with result !-------------------------------------------------------------------------------------------------- -module function mech_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el) result(doneAndHappy) +module function mech_RGC_updateState(P,F,avgF,dt,dPdF,ip,el) result(doneAndHappy) logical, dimension(2) :: doneAndHappy real(pReal), dimension(:,:,:), intent(in) :: & P,& !< partitioned stresses - F,& !< partitioned deformation gradients - F0 !< partitioned initial deformation gradients + F !< partitioned deformation gradients real(pReal), dimension(:,:,:,:,:), intent(in) :: dPdF !< partitioned stiffnesses real(pReal), dimension(3,3), intent(in) :: avgF !< average F real(pReal), intent(in) :: dt !< time increment @@ -287,8 +280,8 @@ module function mech_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el) result(doneAndHa !-------------------------------------------------------------------------------------------------- ! allocate the size of the global relaxation arrays/jacobian matrices depending on the size of the cluster - allocate(resid(3*nIntFaceTot), source=0.0_pReal) - allocate(tract(nIntFaceTot,3), source=0.0_pReal) + allocate(resid(3*nIntFaceTot), source=0.0_pReal) + allocate(tract(nIntFaceTot,3), source=0.0_pReal) relax = stt%relaxationVector(:,of) drelax = stt%relaxationVector(:,of) - st0%relaxationVector(:,of) @@ -346,17 +339,6 @@ module function mech_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el) result(doneAndHa if (residMax < num%rtol*stresMax .or. residMax < num%atol) then doneAndHappy = .true. -!-------------------------------------------------------------------------------------------------- -! compute/update the state for postResult, i.e., all energy densities computed by time-integration - do iGrain = 1,product(prm%N_constituents) - do i = 1,3;do j = 1,3 - stt%work(of) = stt%work(of) & - + P(i,j,iGrain)*(F(i,j,iGrain) - F0(i,j,iGrain))/real(nGrain,pReal) - stt%penaltyEnergy(of) = stt%penaltyEnergy(of) & - + R(i,j,iGrain)*(F(i,j,iGrain) - F0(i,j,iGrain))/real(nGrain,pReal) - enddo; enddo - enddo - dst%mismatch(1:3,of) = sum(NN,2)/real(nGrain,pReal) dst%relaxationRate_avg(of) = sum(abs(drelax))/dt/real(3*nIntFaceTot,pReal) dst%relaxationRate_max(of) = maxval(abs(drelax))/dt @@ -523,7 +505,6 @@ module function mech_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el) result(doneAndHa integer, dimension (3) :: iGrain3,iGNghb3,nGDim real(pReal), dimension (3,3) :: gDef,nDef real(pReal), dimension (3) :: nVect,surfCorr - real(pReal), dimension (2) :: Gmoduli integer :: iGrain,iGNghb,iFace,i,j,k,l real(pReal) :: muGrain,muGNghb,nDefNorm real(pReal), parameter :: & @@ -755,15 +736,9 @@ module subroutine mech_RGC_results(instance,group) associate(stt => state(instance), dst => dependentState(instance), prm => param(instance)) outputsLoop: do o = 1,size(prm%output) select case(trim(prm%output(o))) - case('W') - call results_writeDataset(group,stt%work,trim(prm%output(o)), & - 'work density','J/m³') case('M') call results_writeDataset(group,dst%mismatch,trim(prm%output(o)), & 'average mismatch tensor','1') - case('R') - call results_writeDataset(group,stt%penaltyEnergy,trim(prm%output(o)), & - 'mismatch penalty density','J/m³') case('Delta_V') call results_writeDataset(group,dst%volumeDiscrepancy,trim(prm%output(o)), & 'volume discrepancy','m³') diff --git a/src/thermal_conduction.f90 b/src/thermal_conduction.f90 index d30e50677..f98d36d3b 100644 --- a/src/thermal_conduction.f90 +++ b/src/thermal_conduction.f90 @@ -94,7 +94,7 @@ subroutine thermal_conduction_getSourceAndItsTangent(Tdot, dTdot_dT, T, ip, el) dTdot_dT = 0.0_pReal homog = material_homogenizationAt(el) - call constitutive_thermal_getRateAndItsTangents(TDot, dTDot_dT, T, crystallite_S,crystallite_Lp ,ip, el) + call constitutive_thermal_getRateAndItsTangents(TDot, dTDot_dT, T, ip, el) Tdot = Tdot/real(homogenization_Nconstituents(homog),pReal) dTdot_dT = dTdot_dT/real(homogenization_Nconstituents(homog),pReal) @@ -112,14 +112,16 @@ function thermal_conduction_getConductivity(ip,el) el !< element number real(pReal), dimension(3,3) :: & thermal_conduction_getConductivity + integer :: & - grain + co thermal_conduction_getConductivity = 0.0_pReal - do grain = 1, homogenization_Nconstituents(material_homogenizationAt(el)) + + do co = 1, homogenization_Nconstituents(material_homogenizationAt(el)) thermal_conduction_getConductivity = thermal_conduction_getConductivity + & - crystallite_push33ToRef(grain,ip,el,lattice_K(:,:,material_phaseAt(grain,el))) + crystallite_push33ToRef(co,ip,el,lattice_K(:,:,material_phaseAt(co,el))) enddo thermal_conduction_getConductivity = thermal_conduction_getConductivity & @@ -138,14 +140,16 @@ function thermal_conduction_getSpecificHeat(ip,el) el !< element number real(pReal) :: & thermal_conduction_getSpecificHeat + integer :: & - grain + co + thermal_conduction_getSpecificHeat = 0.0_pReal - do grain = 1, homogenization_Nconstituents(material_homogenizationAt(el)) + do co = 1, homogenization_Nconstituents(material_homogenizationAt(el)) thermal_conduction_getSpecificHeat = thermal_conduction_getSpecificHeat & - + lattice_c_p(material_phaseAt(grain,el)) + + lattice_c_p(material_phaseAt(co,el)) enddo thermal_conduction_getSpecificHeat = thermal_conduction_getSpecificHeat & @@ -164,15 +168,16 @@ function thermal_conduction_getMassDensity(ip,el) el !< element number real(pReal) :: & thermal_conduction_getMassDensity + integer :: & - grain + co + thermal_conduction_getMassDensity = 0.0_pReal - - do grain = 1, homogenization_Nconstituents(material_homogenizationAt(el)) + do co = 1, homogenization_Nconstituents(material_homogenizationAt(el)) thermal_conduction_getMassDensity = thermal_conduction_getMassDensity & - + lattice_rho(material_phaseAt(grain,el)) + + lattice_rho(material_phaseAt(co,el)) enddo thermal_conduction_getMassDensity = thermal_conduction_getMassDensity &