From d585deee7e4f8fbdf631fafde3a1a4802b6630bf Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 25 Aug 2018 13:55:39 +0200 Subject: [PATCH 01/13] unified syntax _v for vector representation of tensor conflicts with 3333 suffix for 4-th order tensors. General idea: Mark symmetric second order tensors in vector notation with '6' and fourth order tensors in second order matrix notation with '99'. Append nothing for 'natural' representation, i.e. F and NOT F33, Lp_dTstar and NOT Lp_dTstar3333 --- PRIVATE | 2 +- src/constitutive.f90 | 91 ++++++++++++++++++++++---------------------- 2 files changed, 47 insertions(+), 46 deletions(-) diff --git a/PRIVATE b/PRIVATE index dfd67ea44..81fd7109f 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit dfd67ea44ba88ee1e0a33266a3986c64137908cf +Subproject commit 81fd7109fea8456b8eecaaef0eec041edcce7792 diff --git a/src/constitutive.f90 b/src/constitutive.f90 index 7833f70cf..9f827472c 100644 --- a/src/constitutive.f90 +++ b/src/constitutive.f90 @@ -345,6 +345,7 @@ end subroutine constitutive_init !-------------------------------------------------------------------------------------------------- !> @brief returns the homogenize elasticity matrix +!> ToDo: homogenizedC66 would be more consistent !-------------------------------------------------------------------------------------------------- function constitutive_homogenizedC(ipc,ip,el) use prec, only: & @@ -429,7 +430,7 @@ end subroutine constitutive_microstructure !-------------------------------------------------------------------------------------------------- !> @brief contains the constitutive equation for calculating the velocity gradient !-------------------------------------------------------------------------------------------------- -subroutine constitutive_LpAndItsTangent(Lp, dLp_dTstar3333, dLp_dFi3333, Tstar_v, Fi, ipc, ip, el) +subroutine constitutive_LpAndItsTangent(Lp, dLp_dTstar, dLp_dFi, Tstar6, Fi, ipc, ip, el) use prec, only: & pReal use math, only: & @@ -469,18 +470,18 @@ subroutine constitutive_LpAndItsTangent(Lp, dLp_dTstar3333, dLp_dFi3333, Tstar_v ip, & !< integration point el !< element real(pReal), intent(in), dimension(6) :: & - Tstar_v !< 2nd Piola-Kirchhoff stress + Tstar6 !< 2nd Piola-Kirchhoff stress real(pReal), intent(in), dimension(3,3) :: & Fi !< intermediate deformation gradient real(pReal), intent(out), dimension(3,3) :: & Lp !< plastic velocity gradient real(pReal), intent(out), dimension(3,3,3,3) :: & - dLp_dTstar3333, & !< derivative of Lp with respect to Tstar (4th-order tensor) - dLp_dFi3333 !< derivative of Lp with respect to Fi (4th-order tensor) + dLp_dTstar, & !< derivative of Lp with respect to Tstar (4th-order tensor) + dLp_dFi !< derivative of Lp with respect to Fi (4th-order tensor) real(pReal), dimension(6) :: & - Mstar_v !< Mandel stress work conjugate with Lp + Mstar6 !< Mandel stress work conjugate with Lp real(pReal), dimension(9,9) :: & - dLp_dMstar !< derivative of Lp with respect to Mstar (4th-order tensor) + dLp_dMstar99 !< derivative of Lp with respect to Mstar (4th-order tensor) real(pReal), dimension(3,3) :: & temp_33 integer(pInt) :: & @@ -492,39 +493,39 @@ subroutine constitutive_LpAndItsTangent(Lp, dLp_dTstar3333, dLp_dFi3333, Tstar_v ho = material_homog(ip,el) tme = thermalMapping(ho)%p(ip,el) - Mstar_v = math_Mandel33to6(math_mul33x33(math_mul33x33(transpose(Fi),Fi),math_Mandel6to33(Tstar_v))) + Mstar6 = math_Mandel33to6(math_mul33x33(math_mul33x33(transpose(Fi),Fi),math_Mandel6to33(Tstar6))) plasticityType: select case (phase_plasticity(material_phase(ipc,ip,el))) case (PLASTICITY_NONE_ID) plasticityType Lp = 0.0_pReal - dLp_dMstar = 0.0_pReal + dLp_dMstar99 = 0.0_pReal case (PLASTICITY_ISOTROPIC_ID) plasticityType - call plastic_isotropic_LpAndItsTangent (Lp,dLp_dMstar,Mstar_v,ipc,ip,el) + call plastic_isotropic_LpAndItsTangent (Lp,dLp_dMstar99,Mstar6,ipc,ip,el) case (PLASTICITY_PHENOPOWERLAW_ID) plasticityType - call plastic_phenopowerlaw_LpAndItsTangent (Lp,dLp_dMstar,Mstar_v,ipc,ip,el) + call plastic_phenopowerlaw_LpAndItsTangent (Lp,dLp_dMstar99,Mstar6,ipc,ip,el) case (PLASTICITY_KINEHARDENING_ID) plasticityType - call plastic_kinehardening_LpAndItsTangent (Lp,dLp_dMstar,Mstar_v,ipc,ip,el) + call plastic_kinehardening_LpAndItsTangent (Lp,dLp_dMstar99,Mstar6,ipc,ip,el) case (PLASTICITY_NONLOCAL_ID) plasticityType - call plastic_nonlocal_LpAndItsTangent (Lp,dLp_dMstar,Mstar_v, & + call plastic_nonlocal_LpAndItsTangent (Lp,dLp_dMstar99,Mstar6, & temperature(ho)%p(tme),ip,el) case (PLASTICITY_DISLOTWIN_ID) plasticityType - call plastic_dislotwin_LpAndItsTangent (Lp,dLp_dMstar,Mstar_v, & + call plastic_dislotwin_LpAndItsTangent (Lp,dLp_dMstar99,Mstar6, & temperature(ho)%p(tme),ipc,ip,el) case (PLASTICITY_DISLOUCLA_ID) plasticityType - call plastic_disloucla_LpAndItsTangent (Lp,dLp_dMstar,Mstar_v, & + call plastic_disloucla_LpAndItsTangent (Lp,dLp_dMstar99,Mstar6, & temperature(ho)%p(tme), ipc,ip,el) end select plasticityType - dLp_dTstar3333 = math_Plain99to3333(dLp_dMstar) - temp_33 = math_mul33x33(Fi,math_Mandel6to33(Tstar_v)) + dLp_dTstar = math_Plain99to3333(dLp_dMstar99) + temp_33 = math_mul33x33(Fi,math_Mandel6to33(Tstar6)) forall(i = 1_pInt:3_pInt, j = 1_pInt:3_pInt) & - dLp_dFi3333(i,j,1:3,1:3) = math_mul33x33(temp_33,transpose(dLp_dTstar3333(i,j,1:3,1:3))) + & - math_mul33x33(math_mul33x33(Fi,dLp_dTstar3333(i,j,1:3,1:3)),math_Mandel6to33(Tstar_v)) + dLp_dFi(i,j,1:3,1:3) = math_mul33x33(temp_33,transpose(dLp_dTstar(i,j,1:3,1:3))) + & + math_mul33x33(math_mul33x33(Fi,dLp_dTstar(i,j,1:3,1:3)),math_Mandel6to33(Tstar6)) temp_33 = math_mul33x33(transpose(Fi),Fi) forall(i = 1_pInt:3_pInt, j = 1_pInt:3_pInt) & - dLp_dTstar3333(i,j,1:3,1:3) = math_mul33x33(temp_33,dLp_dTstar3333(i,j,1:3,1:3)) + dLp_dTstar(i,j,1:3,1:3) = math_mul33x33(temp_33,dLp_dTstar(i,j,1:3,1:3)) end subroutine constitutive_LpAndItsTangent @@ -532,7 +533,7 @@ end subroutine constitutive_LpAndItsTangent !-------------------------------------------------------------------------------------------------- !> @brief contains the constitutive equation for calculating the velocity gradient !-------------------------------------------------------------------------------------------------- -subroutine constitutive_LiAndItsTangent(Li, dLi_dTstar3333, dLi_dFi3333, Tstar_v, Fi, ipc, ip, el) +subroutine constitutive_LiAndItsTangent(Li, dLi_dTstar3333, dLi_dFi3333, Tstar6, Fi, ipc, ip, el) use prec, only: & pReal use math, only: & @@ -570,7 +571,7 @@ subroutine constitutive_LiAndItsTangent(Li, dLi_dTstar3333, dLi_dFi3333, Tstar_v ip, & !< integration point el !< element real(pReal), intent(in), dimension(6) :: & - Tstar_v !< 2nd Piola-Kirchhoff stress + Tstar6 !< 2nd Piola-Kirchhoff stress real(pReal), intent(in), dimension(3,3) :: & Fi !< intermediate deformation gradient real(pReal), intent(out), dimension(3,3) :: & @@ -598,7 +599,7 @@ subroutine constitutive_LiAndItsTangent(Li, dLi_dTstar3333, dLi_dFi3333, Tstar_v plasticityType: select case (phase_plasticity(material_phase(ipc,ip,el))) case (PLASTICITY_isotropic_ID) plasticityType - call plastic_isotropic_LiAndItsTangent(my_Li, my_dLi_dTstar, Tstar_v, ipc, ip, el) + call plastic_isotropic_LiAndItsTangent(my_Li, my_dLi_dTstar, Tstar6, ipc, ip, el) case default plasticityType my_Li = 0.0_pReal my_dLi_dTstar = 0.0_pReal @@ -610,9 +611,9 @@ subroutine constitutive_LiAndItsTangent(Li, dLi_dTstar3333, dLi_dFi3333, Tstar_v KinematicsLoop: do k = 1_pInt, phase_Nkinematics(material_phase(ipc,ip,el)) kinematicsType: select case (phase_kinematics(k,material_phase(ipc,ip,el))) case (KINEMATICS_cleavage_opening_ID) kinematicsType - call kinematics_cleavage_opening_LiAndItsTangent(my_Li, my_dLi_dTstar, Tstar_v, ipc, ip, el) + call kinematics_cleavage_opening_LiAndItsTangent(my_Li, my_dLi_dTstar, Tstar6, ipc, ip, el) case (KINEMATICS_slipplane_opening_ID) kinematicsType - call kinematics_slipplane_opening_LiAndItsTangent(my_Li, my_dLi_dTstar, Tstar_v, ipc, ip, el) + call kinematics_slipplane_opening_LiAndItsTangent(my_Li, my_dLi_dTstar, Tstar6, ipc, ip, el) case (KINEMATICS_thermal_expansion_ID) kinematicsType call kinematics_thermal_expansion_LiAndItsTangent(my_Li, my_dLi_dTstar, ipc, ip, el) case (KINEMATICS_vacancy_strain_ID) kinematicsType @@ -798,7 +799,7 @@ end subroutine constitutive_hooke_TandItsTangent !-------------------------------------------------------------------------------------------------- !> @brief contains the constitutive equation for calculating the rate of change of microstructure !-------------------------------------------------------------------------------------------------- -subroutine constitutive_collectDotState(Tstar_v, FeArray, FpArray, subdt, subfracArray,ipc, ip, el) +subroutine constitutive_collectDotState(Tstar6, FeArray, FpArray, subdt, subfracArray,ipc, ip, el) use prec, only: & pReal, & pLongInt @@ -863,7 +864,7 @@ subroutine constitutive_collectDotState(Tstar_v, FeArray, FpArray, subdt, subfra FeArray, & !< elastic deformation gradient FpArray !< plastic deformation gradient real(pReal), intent(in), dimension(6) :: & - Tstar_v !< 2nd Piola Kirchhoff stress tensor (Mandel) + Tstar6 !< 2nd Piola Kirchhoff stress tensor (Mandel) integer(pInt) :: & ho, & !< homogenization tme, & !< thermal member position @@ -874,26 +875,26 @@ subroutine constitutive_collectDotState(Tstar_v, FeArray, FpArray, subdt, subfra plasticityType: select case (phase_plasticity(material_phase(ipc,ip,el))) case (PLASTICITY_ISOTROPIC_ID) plasticityType - call plastic_isotropic_dotState (Tstar_v,ipc,ip,el) + call plastic_isotropic_dotState (Tstar6,ipc,ip,el) case (PLASTICITY_PHENOPOWERLAW_ID) plasticityType - call plastic_phenopowerlaw_dotState(Tstar_v,ipc,ip,el) + call plastic_phenopowerlaw_dotState(Tstar6,ipc,ip,el) case (PLASTICITY_KINEHARDENING_ID) plasticityType - call plastic_kinehardening_dotState(Tstar_v,ipc,ip,el) + call plastic_kinehardening_dotState(Tstar6,ipc,ip,el) case (PLASTICITY_DISLOTWIN_ID) plasticityType - call plastic_dislotwin_dotState (Tstar_v,temperature(ho)%p(tme), & + call plastic_dislotwin_dotState (Tstar6,temperature(ho)%p(tme), & ipc,ip,el) case (PLASTICITY_DISLOUCLA_ID) plasticityType - call plastic_disloucla_dotState (Tstar_v,temperature(ho)%p(tme), & + call plastic_disloucla_dotState (Tstar6,temperature(ho)%p(tme), & ipc,ip,el) case (PLASTICITY_NONLOCAL_ID) plasticityType - call plastic_nonlocal_dotState (Tstar_v,FeArray,FpArray,temperature(ho)%p(tme), & + call plastic_nonlocal_dotState (Tstar6,FeArray,FpArray,temperature(ho)%p(tme), & subdt,subfracArray,ip,el) end select plasticityType SourceLoop: do s = 1_pInt, phase_Nsources(material_phase(ipc,ip,el)) sourceType: select case (phase_source(s,material_phase(ipc,ip,el))) case (SOURCE_damage_anisoBrittle_ID) sourceType - call source_damage_anisoBrittle_dotState (Tstar_v, ipc, ip, el) + call source_damage_anisoBrittle_dotState (Tstar6, ipc, ip, el) case (SOURCE_damage_isoDuctile_ID) sourceType call source_damage_isoDuctile_dotState ( ipc, ip, el) case (SOURCE_damage_anisoDuctile_ID) sourceType @@ -909,7 +910,7 @@ end subroutine constitutive_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 !-------------------------------------------------------------------------------------------------- -subroutine constitutive_collectDeltaState(Tstar_v, Fe, ipc, ip, el) +subroutine constitutive_collectDeltaState(Tstar6, Fe, ipc, ip, el) use prec, only: & pReal, & pLongInt @@ -944,7 +945,7 @@ subroutine constitutive_collectDeltaState(Tstar_v, Fe, ipc, ip, el) ip, & !< integration point el !< element real(pReal), intent(in), dimension(6) :: & - Tstar_v !< 2nd Piola-Kirchhoff stress + Tstar6 !< 2nd Piola-Kirchhoff stress real(pReal), intent(in), dimension(3,3) :: & Fe !< elastic deformation gradient integer(pInt) :: & @@ -952,9 +953,9 @@ subroutine constitutive_collectDeltaState(Tstar_v, Fe, ipc, ip, el) plasticityType: select case (phase_plasticity(material_phase(ipc,ip,el))) case (PLASTICITY_KINEHARDENING_ID) plasticityType - call plastic_kinehardening_deltaState(Tstar_v,ipc,ip,el) + call plastic_kinehardening_deltaState(Tstar6,ipc,ip,el) case (PLASTICITY_NONLOCAL_ID) plasticityType - call plastic_nonlocal_deltaState(Tstar_v,ip,el) + call plastic_nonlocal_deltaState(Tstar6,ip,el) end select plasticityType SourceLoop: do s = 1_pInt, phase_Nsources(material_phase(ipc,ip,el)) @@ -975,7 +976,7 @@ end subroutine constitutive_collectDeltaState !-------------------------------------------------------------------------------------------------- !> @brief returns array of constitutive results !-------------------------------------------------------------------------------------------------- -function constitutive_postResults(Tstar_v, FeArray, ipc, ip, el) +function constitutive_postResults(Tstar6, FeArray, ipc, ip, el) use prec, only: & pReal use mesh, only: & @@ -1035,7 +1036,7 @@ function constitutive_postResults(Tstar_v, FeArray, ipc, ip, el) real(pReal), intent(in), dimension(3,3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: & FeArray !< elastic deformation gradient real(pReal), intent(in), dimension(6) :: & - Tstar_v !< 2nd Piola Kirchhoff stress tensor (Mandel) + Tstar6 !< 2nd Piola Kirchhoff stress tensor (Mandel) integer(pInt) :: & startPos, endPos integer(pInt) :: & @@ -1053,22 +1054,22 @@ function constitutive_postResults(Tstar_v, FeArray, ipc, ip, el) plasticityType: select case (phase_plasticity(material_phase(ipc,ip,el))) case (PLASTICITY_ISOTROPIC_ID) plasticityType - constitutive_postResults(startPos:endPos) = plastic_isotropic_postResults(Tstar_v,ipc,ip,el) + constitutive_postResults(startPos:endPos) = plastic_isotropic_postResults(Tstar6,ipc,ip,el) case (PLASTICITY_PHENOPOWERLAW_ID) plasticityType constitutive_postResults(startPos:endPos) = & - plastic_phenopowerlaw_postResults(Tstar_v,ipc,ip,el) + plastic_phenopowerlaw_postResults(Tstar6,ipc,ip,el) case (PLASTICITY_KINEHARDENING_ID) plasticityType constitutive_postResults(startPos:endPos) = & - plastic_kinehardening_postResults(Tstar_v,ipc,ip,el) + plastic_kinehardening_postResults(Tstar6,ipc,ip,el) case (PLASTICITY_DISLOTWIN_ID) plasticityType constitutive_postResults(startPos:endPos) = & - plastic_dislotwin_postResults(Tstar_v,temperature(ho)%p(tme),ipc,ip,el) + plastic_dislotwin_postResults(Tstar6,temperature(ho)%p(tme),ipc,ip,el) case (PLASTICITY_DISLOUCLA_ID) plasticityType constitutive_postResults(startPos:endPos) = & - plastic_disloucla_postResults(Tstar_v,temperature(ho)%p(tme),ipc,ip,el) + plastic_disloucla_postResults(Tstar6,temperature(ho)%p(tme),ipc,ip,el) case (PLASTICITY_NONLOCAL_ID) plasticityType constitutive_postResults(startPos:endPos) = & - plastic_nonlocal_postResults (Tstar_v,FeArray,ip,el) + plastic_nonlocal_postResults (Tstar6,FeArray,ip,el) end select plasticityType SourceLoop: do s = 1_pInt, phase_Nsources(material_phase(ipc,ip,el)) From 953fff79ac323012c219049763ebaef226c554c7 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 25 Aug 2018 14:21:06 +0200 Subject: [PATCH 02/13] prepared LpAndItsTangent to remove superflous forward-backward conversion updating the individual plastic laws will be done in the respective branches --- src/constitutive.f90 | 47 ++++++++++++++++++++++++++++---------------- 1 file changed, 30 insertions(+), 17 deletions(-) diff --git a/src/constitutive.f90 b/src/constitutive.f90 index 9f827472c..786e1a302 100644 --- a/src/constitutive.f90 +++ b/src/constitutive.f90 @@ -470,19 +470,19 @@ subroutine constitutive_LpAndItsTangent(Lp, dLp_dTstar, dLp_dFi, Tstar6, Fi, ipc ip, & !< integration point el !< element real(pReal), intent(in), dimension(6) :: & - Tstar6 !< 2nd Piola-Kirchhoff stress + Tstar6 !< 2nd Piola-Kirchhoff stress (vector notation) real(pReal), intent(in), dimension(3,3) :: & Fi !< intermediate deformation gradient real(pReal), intent(out), dimension(3,3) :: & Lp !< plastic velocity gradient real(pReal), intent(out), dimension(3,3,3,3) :: & - dLp_dTstar, & !< derivative of Lp with respect to Tstar (4th-order tensor) - dLp_dFi !< derivative of Lp with respect to Fi (4th-order tensor) - real(pReal), dimension(6) :: & - Mstar6 !< Mandel stress work conjugate with Lp + dLp_dTstar, & !< derivative of Lp with respect to Tstar + dLp_dFi !< derivative of Lp with respect to Fi real(pReal), dimension(9,9) :: & - dLp_dMstar99 !< derivative of Lp with respect to Mstar (4th-order tensor) + dLp_dMstar99 !< derivative of Lp with respect to Mstar (matrix notation) real(pReal), dimension(3,3) :: & + Mstar, & !< Mandel stress work conjugate with Lp + Tstar, & !< 2nd Piola-Kirchhoff stress temp_33 integer(pInt) :: & ho, & !< homogenization @@ -493,34 +493,47 @@ subroutine constitutive_LpAndItsTangent(Lp, dLp_dTstar, dLp_dFi, Tstar6, Fi, ipc ho = material_homog(ip,el) tme = thermalMapping(ho)%p(ip,el) - Mstar6 = math_Mandel33to6(math_mul33x33(math_mul33x33(transpose(Fi),Fi),math_Mandel6to33(Tstar6))) + Tstar = math_Mandel6to33(Tstar6) + Mstar = math_mul33x33(math_mul33x33(transpose(Fi),Fi),Tstar)) plasticityType: select case (phase_plasticity(material_phase(ipc,ip,el))) case (PLASTICITY_NONE_ID) plasticityType Lp = 0.0_pReal - dLp_dMstar99 = 0.0_pReal + dLp_dTstar = 0.0_pReal + case (PLASTICITY_ISOTROPIC_ID) plasticityType - call plastic_isotropic_LpAndItsTangent (Lp,dLp_dMstar99,Mstar6,ipc,ip,el) + call plastic_isotropic_LpAndItsTangent (Lp,dLp_dMstar99, math_Mandel33to6(Mstar),ipc,ip,el) + dLp_dTstar = math_Plain99to3333(dLp_dMstar99) ! ToDo: We revert here the last statement in plastic_xx_LpAndItsTanget + case (PLASTICITY_PHENOPOWERLAW_ID) plasticityType - call plastic_phenopowerlaw_LpAndItsTangent (Lp,dLp_dMstar99,Mstar6,ipc,ip,el) + call plastic_phenopowerlaw_LpAndItsTangent (Lp,dLp_dMstar99, math_Mandel33to6(Mstar),ipc,ip,el) + dLp_dTstar = math_Plain99to3333(dLp_dMstar99) ! ToDo: We revert here the last statement in plastic_xx_LpAndItsTanget + case (PLASTICITY_KINEHARDENING_ID) plasticityType - call plastic_kinehardening_LpAndItsTangent (Lp,dLp_dMstar99,Mstar6,ipc,ip,el) + call plastic_kinehardening_LpAndItsTangent (Lp,dLp_dMstar99, math_Mandel33to6(Mstar),ipc,ip,el) + dLp_dTstar = math_Plain99to3333(dLp_dMstar99) ! ToDo: We revert here the last statement in plastic_xx_LpAndItsTanget + case (PLASTICITY_NONLOCAL_ID) plasticityType - call plastic_nonlocal_LpAndItsTangent (Lp,dLp_dMstar99,Mstar6, & + call plastic_nonlocal_LpAndItsTangent (Lp,dLp_dMstar99, math_Mandel33to6(Mstar), & temperature(ho)%p(tme),ip,el) + dLp_dTstar = math_Plain99to3333(dLp_dMstar99) ! ToDo: We revert here the last statement in plastic_xx_LpAndItsTanget + case (PLASTICITY_DISLOTWIN_ID) plasticityType - call plastic_dislotwin_LpAndItsTangent (Lp,dLp_dMstar99,Mstar6, & + call plastic_dislotwin_LpAndItsTangent (Lp,dLp_dMstar99, math_Mandel33to6(Mstar), & temperature(ho)%p(tme),ipc,ip,el) + dLp_dTstar = math_Plain99to3333(dLp_dMstar99) ! ToDo: We revert here the last statement in plastic_xx_LpAndItsTanget + case (PLASTICITY_DISLOUCLA_ID) plasticityType - call plastic_disloucla_LpAndItsTangent (Lp,dLp_dMstar99,Mstar6, & + call plastic_disloucla_LpAndItsTangent (Lp,dLp_dMstar99, math_Mandel33to6(Mstar), & temperature(ho)%p(tme), ipc,ip,el) + dLp_dTstar = math_Plain99to3333(dLp_dMstar99) ! ToDo: We revert here the last statement in plastic_xx_LpAndItsTanget + end select plasticityType - dLp_dTstar = math_Plain99to3333(dLp_dMstar99) - temp_33 = math_mul33x33(Fi,math_Mandel6to33(Tstar6)) + temp_33 = math_mul33x33(Fi,Tstar) forall(i = 1_pInt:3_pInt, j = 1_pInt:3_pInt) & dLp_dFi(i,j,1:3,1:3) = math_mul33x33(temp_33,transpose(dLp_dTstar(i,j,1:3,1:3))) + & - math_mul33x33(math_mul33x33(Fi,dLp_dTstar(i,j,1:3,1:3)),math_Mandel6to33(Tstar6)) + math_mul33x33(math_mul33x33(Fi,dLp_dTstar(i,j,1:3,1:3)),Tstar)) temp_33 = math_mul33x33(transpose(Fi),Fi) From 35688a6acf9030c9b1a51f7442c5bd4f09d235da Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 25 Aug 2018 14:29:03 +0200 Subject: [PATCH 03/13] temp33 not needed (compiler should be smart enough) --- src/constitutive.f90 | 10 +++------- 1 file changed, 3 insertions(+), 7 deletions(-) diff --git a/src/constitutive.f90 b/src/constitutive.f90 index 786e1a302..7a6abf91f 100644 --- a/src/constitutive.f90 +++ b/src/constitutive.f90 @@ -482,8 +482,7 @@ subroutine constitutive_LpAndItsTangent(Lp, dLp_dTstar, dLp_dFi, Tstar6, Fi, ipc dLp_dMstar99 !< derivative of Lp with respect to Mstar (matrix notation) real(pReal), dimension(3,3) :: & Mstar, & !< Mandel stress work conjugate with Lp - Tstar, & !< 2nd Piola-Kirchhoff stress - temp_33 + Tstar !< 2nd Piola-Kirchhoff stress integer(pInt) :: & ho, & !< homogenization tme !< thermal member position @@ -530,15 +529,12 @@ subroutine constitutive_LpAndItsTangent(Lp, dLp_dTstar, dLp_dFi, Tstar6, Fi, ipc end select plasticityType - temp_33 = math_mul33x33(Fi,Tstar) forall(i = 1_pInt:3_pInt, j = 1_pInt:3_pInt) & - dLp_dFi(i,j,1:3,1:3) = math_mul33x33(temp_33,transpose(dLp_dTstar(i,j,1:3,1:3))) + & + dLp_dFi(i,j,1:3,1:3) = math_mul33x33(math_mul33x33(Fi,Tstar),transpose(dLp_dTstar(i,j,1:3,1:3))) + & math_mul33x33(math_mul33x33(Fi,dLp_dTstar(i,j,1:3,1:3)),Tstar)) - temp_33 = math_mul33x33(transpose(Fi),Fi) - forall(i = 1_pInt:3_pInt, j = 1_pInt:3_pInt) & - dLp_dTstar(i,j,1:3,1:3) = math_mul33x33(temp_33,dLp_dTstar(i,j,1:3,1:3)) + dLp_dTstar(i,j,1:3,1:3) = math_mul33x33(math_mul33x33(transpose(Fi),Fi),dLp_dTstar(i,j,1:3,1:3)) end subroutine constitutive_LpAndItsTangent From 17d88184a7bded4dfd4a3858325eccbb3f5a06c4 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 25 Aug 2018 14:35:49 +0200 Subject: [PATCH 04/13] typos fixed --- src/constitutive.f90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/constitutive.f90 b/src/constitutive.f90 index 7a6abf91f..e001d6348 100644 --- a/src/constitutive.f90 +++ b/src/constitutive.f90 @@ -493,7 +493,7 @@ subroutine constitutive_LpAndItsTangent(Lp, dLp_dTstar, dLp_dFi, Tstar6, Fi, ipc tme = thermalMapping(ho)%p(ip,el) Tstar = math_Mandel6to33(Tstar6) - Mstar = math_mul33x33(math_mul33x33(transpose(Fi),Fi),Tstar)) + Mstar = math_mul33x33(math_mul33x33(transpose(Fi),Fi),Tstar) plasticityType: select case (phase_plasticity(material_phase(ipc,ip,el))) case (PLASTICITY_NONE_ID) plasticityType @@ -531,7 +531,7 @@ subroutine constitutive_LpAndItsTangent(Lp, dLp_dTstar, dLp_dFi, Tstar6, Fi, ipc forall(i = 1_pInt:3_pInt, j = 1_pInt:3_pInt) & dLp_dFi(i,j,1:3,1:3) = math_mul33x33(math_mul33x33(Fi,Tstar),transpose(dLp_dTstar(i,j,1:3,1:3))) + & - math_mul33x33(math_mul33x33(Fi,dLp_dTstar(i,j,1:3,1:3)),Tstar)) + math_mul33x33(math_mul33x33(Fi,dLp_dTstar(i,j,1:3,1:3)),Tstar) forall(i = 1_pInt:3_pInt, j = 1_pInt:3_pInt) & dLp_dTstar(i,j,1:3,1:3) = math_mul33x33(math_mul33x33(transpose(Fi),Fi),dLp_dTstar(i,j,1:3,1:3)) From e46605f0ef0e7a260b4e6b6ba35479b40e3efc1c Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 25 Aug 2018 14:36:21 +0200 Subject: [PATCH 05/13] forall is deprecated --- src/constitutive.f90 | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/src/constitutive.f90 b/src/constitutive.f90 index e001d6348..2f92e3f1e 100644 --- a/src/constitutive.f90 +++ b/src/constitutive.f90 @@ -529,12 +529,11 @@ subroutine constitutive_LpAndItsTangent(Lp, dLp_dTstar, dLp_dFi, Tstar6, Fi, ipc end select plasticityType - forall(i = 1_pInt:3_pInt, j = 1_pInt:3_pInt) & + do concurrent(i = 1_pInt:3_pInt, j = 1_pInt:3_pInt) dLp_dFi(i,j,1:3,1:3) = math_mul33x33(math_mul33x33(Fi,Tstar),transpose(dLp_dTstar(i,j,1:3,1:3))) + & math_mul33x33(math_mul33x33(Fi,dLp_dTstar(i,j,1:3,1:3)),Tstar) - - forall(i = 1_pInt:3_pInt, j = 1_pInt:3_pInt) & dLp_dTstar(i,j,1:3,1:3) = math_mul33x33(math_mul33x33(transpose(Fi),Fi),dLp_dTstar(i,j,1:3,1:3)) + enddo end subroutine constitutive_LpAndItsTangent From 7e002d1bfb483aabac321ca4d25cca7af5eeb424 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 25 Aug 2018 14:42:44 +0200 Subject: [PATCH 06/13] missing rename --- src/constitutive.f90 | 36 +++++++++++++++++++++++++----------- 1 file changed, 25 insertions(+), 11 deletions(-) diff --git a/src/constitutive.f90 b/src/constitutive.f90 index 2f92e3f1e..5b9662c02 100644 --- a/src/constitutive.f90 +++ b/src/constitutive.f90 @@ -496,6 +496,7 @@ subroutine constitutive_LpAndItsTangent(Lp, dLp_dTstar, dLp_dFi, Tstar6, Fi, ipc Mstar = math_mul33x33(math_mul33x33(transpose(Fi),Fi),Tstar) plasticityType: select case (phase_plasticity(material_phase(ipc,ip,el))) + case (PLASTICITY_NONE_ID) plasticityType Lp = 0.0_pReal dLp_dTstar = 0.0_pReal @@ -541,7 +542,7 @@ end subroutine constitutive_LpAndItsTangent !-------------------------------------------------------------------------------------------------- !> @brief contains the constitutive equation for calculating the velocity gradient !-------------------------------------------------------------------------------------------------- -subroutine constitutive_LiAndItsTangent(Li, dLi_dTstar3333, dLi_dFi3333, Tstar6, Fi, ipc, ip, el) +subroutine constitutive_LiAndItsTangent(Li, dLi_dTstar, dLi_dFi, Tstar6, Fi, ipc, ip, el) use prec, only: & pReal use math, only: & @@ -585,8 +586,8 @@ subroutine constitutive_LiAndItsTangent(Li, dLi_dTstar3333, dLi_dFi3333, Tstar6, real(pReal), intent(out), dimension(3,3) :: & Li !< intermediate velocity gradient real(pReal), intent(out), dimension(3,3,3,3) :: & - dLi_dTstar3333, & !< derivative of Li with respect to Tstar (4th-order tensor) - dLi_dFi3333 + dLi_dTstar, & !< derivative of Li with respect to Tstar + dLi_dFi real(pReal), dimension(3,3) :: & my_Li !< intermediate velocity gradient real(pReal), dimension(3,3,3,3) :: & @@ -602,8 +603,8 @@ subroutine constitutive_LiAndItsTangent(Li, dLi_dTstar3333, dLi_dFi3333, Tstar6, i, j Li = 0.0_pReal - dLi_dTstar3333 = 0.0_pReal - dLi_dFi3333 = 0.0_pReal + dLi_dTstar = 0.0_pReal + dLi_dFi = 0.0_pReal plasticityType: select case (phase_plasticity(material_phase(ipc,ip,el))) case (PLASTICITY_isotropic_ID) plasticityType @@ -614,7 +615,7 @@ subroutine constitutive_LiAndItsTangent(Li, dLi_dTstar3333, dLi_dFi3333, Tstar6, end select plasticityType Li = Li + my_Li - dLi_dTstar3333 = dLi_dTstar3333 + my_dLi_dTstar + dLi_dTstar = dLi_dTstar + my_dLi_dTstar KinematicsLoop: do k = 1_pInt, phase_Nkinematics(material_phase(ipc,ip,el)) kinematicsType: select case (phase_kinematics(k,material_phase(ipc,ip,el))) @@ -633,7 +634,7 @@ subroutine constitutive_LiAndItsTangent(Li, dLi_dTstar3333, dLi_dFi3333, Tstar6, my_dLi_dTstar = 0.0_pReal end select kinematicsType Li = Li + my_Li - dLi_dTstar3333 = dLi_dTstar3333 + my_dLi_dTstar + dLi_dTstar = dLi_dTstar + my_dLi_dTstar enddo KinematicsLoop FiInv = math_inv33(Fi) @@ -641,9 +642,9 @@ subroutine constitutive_LiAndItsTangent(Li, dLi_dTstar3333, dLi_dFi3333, Tstar6, Li = math_mul33x33(math_mul33x33(Fi,Li),FiInv)*detFi !< push forward to intermediate configuration temp_33 = math_mul33x33(FiInv,Li) forall(i = 1_pInt:3_pInt, j = 1_pInt:3_pInt) - dLi_dTstar3333(1:3,1:3,i,j) = math_mul33x33(math_mul33x33(Fi,dLi_dTstar3333(1:3,1:3,i,j)),FiInv)*detFi - dLi_dFi3333 (1:3,1:3,i,j) = dLi_dFi3333(1:3,1:3,i,j) + Li*FiInv(j,i) - dLi_dFi3333 (1:3,i,1:3,j) = dLi_dFi3333(1:3,i,1:3,j) + math_I3*temp_33(j,i) + Li*FiInv(j,i) + dLi_dTstar(1:3,1:3,i,j) = math_mul33x33(math_mul33x33(Fi,dLi_dTstar(1:3,1:3,i,j)),FiInv)*detFi + dLi_dFi (1:3,1:3,i,j) = dLi_dFi(1:3,1:3,i,j) + Li*FiInv(j,i) + dLi_dFi (1:3,i,1:3,j) = dLi_dFi(1:3,i,1:3,j) + math_I3*temp_33(j,i) + Li*FiInv(j,i) end forall end subroutine constitutive_LiAndItsTangent @@ -882,34 +883,47 @@ subroutine constitutive_collectDotState(Tstar6, FeArray, FpArray, subdt, subfrac tme = thermalMapping(ho)%p(ip,el) plasticityType: select case (phase_plasticity(material_phase(ipc,ip,el))) + case (PLASTICITY_ISOTROPIC_ID) plasticityType call plastic_isotropic_dotState (Tstar6,ipc,ip,el) + case (PLASTICITY_PHENOPOWERLAW_ID) plasticityType call plastic_phenopowerlaw_dotState(Tstar6,ipc,ip,el) + case (PLASTICITY_KINEHARDENING_ID) plasticityType call plastic_kinehardening_dotState(Tstar6,ipc,ip,el) + case (PLASTICITY_DISLOTWIN_ID) plasticityType call plastic_dislotwin_dotState (Tstar6,temperature(ho)%p(tme), & ipc,ip,el) + case (PLASTICITY_DISLOUCLA_ID) plasticityType call plastic_disloucla_dotState (Tstar6,temperature(ho)%p(tme), & ipc,ip,el) + case (PLASTICITY_NONLOCAL_ID) plasticityType call plastic_nonlocal_dotState (Tstar6,FeArray,FpArray,temperature(ho)%p(tme), & subdt,subfracArray,ip,el) end select plasticityType SourceLoop: do s = 1_pInt, phase_Nsources(material_phase(ipc,ip,el)) - sourceType: select case (phase_source(s,material_phase(ipc,ip,el))) + + sourceType: select case (phase_source(s,material_phase(ipc,ip,el))) + case (SOURCE_damage_anisoBrittle_ID) sourceType call source_damage_anisoBrittle_dotState (Tstar6, ipc, ip, el) + case (SOURCE_damage_isoDuctile_ID) sourceType call source_damage_isoDuctile_dotState ( ipc, ip, el) + case (SOURCE_damage_anisoDuctile_ID) sourceType call source_damage_anisoDuctile_dotState ( ipc, ip, el) + case (SOURCE_thermal_externalheat_ID) sourceType call source_thermal_externalheat_dotState( ipc, ip, el) + end select sourceType + enddo SourceLoop end subroutine constitutive_collectDotState From b8b5bac684aa0f8fe254a8f1b2cddca2d78d1b37 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 25 Aug 2018 15:59:34 +0200 Subject: [PATCH 07/13] dotState and deltaState parse Mstar instead of Tstar requires to parse in Fi --- src/constitutive.f90 | 70 +++++++++++++++++++++++++++++++------------- src/crystallite.f90 | 14 ++++++++- 2 files changed, 63 insertions(+), 21 deletions(-) diff --git a/src/constitutive.f90 b/src/constitutive.f90 index 5b9662c02..e43506c1d 100644 --- a/src/constitutive.f90 +++ b/src/constitutive.f90 @@ -580,7 +580,7 @@ subroutine constitutive_LiAndItsTangent(Li, dLi_dTstar, dLi_dFi, Tstar6, Fi, ipc ip, & !< integration point el !< element real(pReal), intent(in), dimension(6) :: & - Tstar6 !< 2nd Piola-Kirchhoff stress + Tstar6 !< 2nd Piola-Kirchhoff stress (vector notation) real(pReal), intent(in), dimension(3,3) :: & Fi !< intermediate deformation gradient real(pReal), intent(out), dimension(3,3) :: & @@ -740,11 +740,9 @@ subroutine constitutive_hooke_TandItsTangent(T, dT_dFe, dT_dFi, Fe, Fi, ipc, ip, use prec, only: & pReal use math, only : & - math_mul3x3, & math_mul33x33, & math_mul3333xx33, & math_Mandel66to3333, & - math_trace33, & math_I3 use material, only: & material_phase, & @@ -808,7 +806,7 @@ end subroutine constitutive_hooke_TandItsTangent !-------------------------------------------------------------------------------------------------- !> @brief contains the constitutive equation for calculating the rate of change of microstructure !-------------------------------------------------------------------------------------------------- -subroutine constitutive_collectDotState(Tstar6, FeArray, FpArray, subdt, subfracArray,ipc, ip, el) +subroutine constitutive_collectDotState(Tstar6, FeArray, Fi, FpArray, subdt, subfracArray,ipc, ip, el) use prec, only: & pReal, & pLongInt @@ -816,6 +814,10 @@ subroutine constitutive_collectDotState(Tstar6, FeArray, FpArray, subdt, subfrac debug_level, & debug_constitutive, & debug_levelBasic + use math, only: & + math_Mandel6to33, & + math_Mandel33to6, & + math_mul33x33 use mesh, only: & mesh_NcpElems, & mesh_maxNips @@ -872,8 +874,13 @@ subroutine constitutive_collectDotState(Tstar6, FeArray, FpArray, subdt, subfrac real(pReal), intent(in), dimension(3,3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: & FeArray, & !< elastic deformation gradient FpArray !< plastic deformation gradient + real(pReal), intent(in), dimension(3,3) :: & + Fi !< intermediate deformation gradient real(pReal), intent(in), dimension(6) :: & - Tstar6 !< 2nd Piola Kirchhoff stress tensor (Mandel) + Tstar6 !< 2nd Piola Kirchhoff stress (vector notation) + real(pReal), dimension(3,3) :: & + Tstar, & !< 2nd Piola Kirchhoff stress tensor + Mstar integer(pInt) :: & ho, & !< homogenization tme, & !< thermal member position @@ -882,27 +889,30 @@ subroutine constitutive_collectDotState(Tstar6, FeArray, FpArray, subdt, subfrac ho = material_homog( ip,el) tme = thermalMapping(ho)%p(ip,el) + Tstar = math_Mandel6to33(Tstar6) + Mstar = math_mul33x33(math_mul33x33(transpose(Fi),Fi),Tstar) + plasticityType: select case (phase_plasticity(material_phase(ipc,ip,el))) case (PLASTICITY_ISOTROPIC_ID) plasticityType - call plastic_isotropic_dotState (Tstar6,ipc,ip,el) + call plastic_isotropic_dotState (math_Mandel33to6(Mstar),ipc,ip,el) case (PLASTICITY_PHENOPOWERLAW_ID) plasticityType - call plastic_phenopowerlaw_dotState(Tstar6,ipc,ip,el) + call plastic_phenopowerlaw_dotState(math_Mandel33to6(Mstar),ipc,ip,el) case (PLASTICITY_KINEHARDENING_ID) plasticityType - call plastic_kinehardening_dotState(Tstar6,ipc,ip,el) + call plastic_kinehardening_dotState(math_Mandel33to6(Mstar),ipc,ip,el) case (PLASTICITY_DISLOTWIN_ID) plasticityType - call plastic_dislotwin_dotState (Tstar6,temperature(ho)%p(tme), & + call plastic_dislotwin_dotState (math_Mandel33to6(Mstar),temperature(ho)%p(tme), & ipc,ip,el) case (PLASTICITY_DISLOUCLA_ID) plasticityType - call plastic_disloucla_dotState (Tstar6,temperature(ho)%p(tme), & + call plastic_disloucla_dotState (math_Mandel33to6(Mstar),temperature(ho)%p(tme), & ipc,ip,el) case (PLASTICITY_NONLOCAL_ID) plasticityType - call plastic_nonlocal_dotState (Tstar6,FeArray,FpArray,temperature(ho)%p(tme), & + call plastic_nonlocal_dotState (math_Mandel33to6(Mstar),FeArray,FpArray,temperature(ho)%p(tme), & subdt,subfracArray,ip,el) end select plasticityType @@ -911,7 +921,7 @@ subroutine constitutive_collectDotState(Tstar6, FeArray, FpArray, subdt, subfrac sourceType: select case (phase_source(s,material_phase(ipc,ip,el))) case (SOURCE_damage_anisoBrittle_ID) sourceType - call source_damage_anisoBrittle_dotState (Tstar6, ipc, ip, el) + call source_damage_anisoBrittle_dotState (Tstar6, ipc, ip, el) !< correct stress? case (SOURCE_damage_isoDuctile_ID) sourceType call source_damage_isoDuctile_dotState ( ipc, ip, el) @@ -932,7 +942,7 @@ end subroutine constitutive_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 !-------------------------------------------------------------------------------------------------- -subroutine constitutive_collectDeltaState(Tstar6, Fe, ipc, ip, el) +subroutine constitutive_collectDeltaState(Tstar6, Fe, Fi, ipc, ip, el) use prec, only: & pReal, & pLongInt @@ -940,6 +950,10 @@ subroutine constitutive_collectDeltaState(Tstar6, Fe, ipc, ip, el) debug_level, & debug_constitutive, & debug_levelBasic + use math, only: & + math_Mandel6to33, & + math_Mandel33to6, & + math_mul33x33 use material, only: & phase_plasticity, & phase_source, & @@ -967,29 +981,45 @@ subroutine constitutive_collectDeltaState(Tstar6, Fe, ipc, ip, el) ip, & !< integration point el !< element real(pReal), intent(in), dimension(6) :: & - Tstar6 !< 2nd Piola-Kirchhoff stress + Tstar6 !< 2nd Piola Kirchhoff stress (vector notation) real(pReal), intent(in), dimension(3,3) :: & - Fe !< elastic deformation gradient + Fe, & !< elastic deformation gradient + Fi !< intermediate deformation gradient + real(pReal), dimension(3,3) :: & + Tstar, & !< 2nd Piola Kirchhoff stress tensor + Mstar integer(pInt) :: & s !< counter in source loop + Tstar = math_Mandel6to33(Tstar6) + Mstar = math_mul33x33(math_mul33x33(transpose(Fi),Fi),Tstar) + plasticityType: select case (phase_plasticity(material_phase(ipc,ip,el))) + case (PLASTICITY_KINEHARDENING_ID) plasticityType - call plastic_kinehardening_deltaState(Tstar6,ipc,ip,el) + call plastic_kinehardening_deltaState(math_Mandel33to6(Mstar),ipc,ip,el) + case (PLASTICITY_NONLOCAL_ID) plasticityType - call plastic_nonlocal_deltaState(Tstar6,ip,el) + call plastic_nonlocal_deltaState(math_Mandel33to6(Mstar),ip,el) + end select plasticityType - SourceLoop: do s = 1_pInt, phase_Nsources(material_phase(ipc,ip,el)) + sourceLoop: do s = 1_pInt, phase_Nsources(material_phase(ipc,ip,el)) + sourceType: select case (phase_source(s,material_phase(ipc,ip,el))) + case (SOURCE_damage_isoBrittle_ID) sourceType call source_damage_isoBrittle_deltaState (constitutive_homogenizedC(ipc,ip,el), Fe, & - ipc, ip, el) + ipc, ip, el) + case (SOURCE_vacancy_irradiation_ID) sourceType call source_vacancy_irradiation_deltaState(ipc, ip, el) + case (SOURCE_vacancy_thermalfluc_ID) sourceType call source_vacancy_thermalfluc_deltaState(ipc, ip, el) + end select sourceType + enddo SourceLoop end subroutine constitutive_collectDeltaState @@ -1058,7 +1088,7 @@ function constitutive_postResults(Tstar6, FeArray, ipc, ip, el) real(pReal), intent(in), dimension(3,3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: & FeArray !< elastic deformation gradient real(pReal), intent(in), dimension(6) :: & - Tstar6 !< 2nd Piola Kirchhoff stress tensor (Mandel) + Tstar6 !< 2nd Piola Kirchhoff stress (vector notation) integer(pInt) :: & startPos, endPos integer(pInt) :: & diff --git a/src/crystallite.f90 b/src/crystallite.f90 index 0ee71b5de..894744f28 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -1300,6 +1300,7 @@ subroutine crystallite_integrateStateRK4() if (crystallite_todo(g,i,e)) & call constitutive_collectDotState(crystallite_Tstar_v(1:6,g,i,e), & crystallite_Fe, & + crystallite_Fi(1:3,1:3,g,i,e), & crystallite_Fp, & crystallite_subdt(g,i,e), crystallite_subFrac, g,i,e) enddo; enddo; enddo @@ -1442,6 +1443,7 @@ subroutine crystallite_integrateStateRK4() if (crystallite_todo(g,i,e)) & call constitutive_collectDotState(crystallite_Tstar_v(1:6,g,i,e), & crystallite_Fe, & + crystallite_Fi(1:3,1:3,g,i,e), & crystallite_Fp, & timeStepFraction(n)*crystallite_subdt(g,i,e), & ! fraction of original timestep crystallite_subFrac, g,i,e) @@ -1609,6 +1611,7 @@ subroutine crystallite_integrateStateRKCK45() if (crystallite_todo(g,i,e)) & call constitutive_collectDotState(crystallite_Tstar_v(1:6,g,i,e), & crystallite_Fe, & + crystallite_Fi(1:3,1:3,g,i,e), & crystallite_Fp, & crystallite_subdt(g,i,e), crystallite_subFrac, g,i,e) enddo; enddo; enddo @@ -1761,6 +1764,7 @@ subroutine crystallite_integrateStateRKCK45() if (crystallite_todo(g,i,e)) & call constitutive_collectDotState(crystallite_Tstar_v(1:6,g,i,e), & crystallite_Fe, & + crystallite_Fi(1:3,1:3,g,i,e), & crystallite_Fp, & C(stage)*crystallite_subdt(g,i,e), & ! fraction of original timestep crystallite_subFrac, g,i,e) @@ -2088,6 +2092,7 @@ subroutine crystallite_integrateStateAdaptiveEuler() if (crystallite_todo(g,i,e)) & call constitutive_collectDotState(crystallite_Tstar_v(1:6,g,i,e), & crystallite_Fe, & + crystallite_Fi(1:3,1:3,g,i,e), & crystallite_Fp, & crystallite_subdt(g,i,e), crystallite_subFrac, g,i,e) enddo; enddo; enddo @@ -2209,6 +2214,7 @@ subroutine crystallite_integrateStateAdaptiveEuler() if (crystallite_todo(g,i,e)) & call constitutive_collectDotState(crystallite_Tstar_v(1:6,g,i,e), & crystallite_Fe, & + crystallite_Fi(1:3,1:3,g,i,e), & crystallite_Fp, & crystallite_subdt(g,i,e), crystallite_subFrac, g,i,e) enddo; enddo; enddo @@ -2417,6 +2423,7 @@ eIter = FEsolving_execElem(1:2) if (crystallite_todo(g,i,e) .and. .not. crystallite_converged(g,i,e)) & call constitutive_collectDotState(crystallite_Tstar_v(1:6,g,i,e), & crystallite_Fe, & + crystallite_Fi(1:3,1:3,g,i,e), & crystallite_Fp, & crystallite_subdt(g,i,e), crystallite_subFrac, g,i,e) enddo; enddo; enddo @@ -2677,6 +2684,7 @@ subroutine crystallite_integrateStateFPI() if (crystallite_todo(g,i,e)) & call constitutive_collectDotState(crystallite_Tstar_v(1:6,g,i,e), & crystallite_Fe, & + crystallite_Fi(1:3,1:3,g,i,e), & crystallite_Fp, & crystallite_subdt(g,i,e), crystallite_subFrac, g,i,e) enddo; enddo; enddo @@ -2799,6 +2807,7 @@ subroutine crystallite_integrateStateFPI() if (crystallite_todo(g,i,e) .and. .not. crystallite_converged(g,i,e)) & call constitutive_collectDotState(crystallite_Tstar_v(1:6,g,i,e), & crystallite_Fe, & + crystallite_Fi(1:3,1:3,g,i,e), & crystallite_Fp, & crystallite_subdt(g,i,e), crystallite_subFrac, g,i,e) enddo; enddo; enddo @@ -3057,7 +3066,10 @@ logical function crystallite_stateJump(ipc,ip,el) c = phasememberAt(ipc,ip,el) p = phaseAt(ipc,ip,el) - call constitutive_collectDeltaState(crystallite_Tstar_v(1:6,ipc,ip,el), crystallite_Fe(1:3,1:3,ipc,ip,el), ipc,ip,el) + call constitutive_collectDeltaState(crystallite_Tstar_v(1:6,ipc,ip,el), & + crystallite_Fe(1:3,1:3,ipc,ip,el), & + crystallite_Fi(1:3,1:3,ipc,ip,el), & + ipc,ip,el) myOffsetPlasticDeltaState = plasticState(p)%offsetDeltaState mySizePlasticDeltaState = plasticState(p)%sizeDeltaState From dd81fe806d575860a7090cb29a3e40e1213c71bd Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Mon, 27 Aug 2018 06:30:55 +0200 Subject: [PATCH 08/13] naming was not completely corrected --- src/constitutive.f90 | 32 ++++++++++++++------------------ 1 file changed, 14 insertions(+), 18 deletions(-) diff --git a/src/constitutive.f90 b/src/constitutive.f90 index e43506c1d..4275f0d8d 100644 --- a/src/constitutive.f90 +++ b/src/constitutive.f90 @@ -430,7 +430,7 @@ end subroutine constitutive_microstructure !-------------------------------------------------------------------------------------------------- !> @brief contains the constitutive equation for calculating the velocity gradient !-------------------------------------------------------------------------------------------------- -subroutine constitutive_LpAndItsTangent(Lp, dLp_dTstar, dLp_dFi, Tstar6, Fi, ipc, ip, el) +subroutine constitutive_LpAndItsTangent(Lp, dLp_dMstar, dLp_dFi, Tstar6, Fi, ipc, ip, el) use prec, only: & pReal use math, only: & @@ -476,7 +476,7 @@ subroutine constitutive_LpAndItsTangent(Lp, dLp_dTstar, dLp_dFi, Tstar6, Fi, ipc real(pReal), intent(out), dimension(3,3) :: & Lp !< plastic velocity gradient real(pReal), intent(out), dimension(3,3,3,3) :: & - dLp_dTstar, & !< derivative of Lp with respect to Tstar + dLp_dMstar, & !< derivative of Lp with respect to Mandel stress dLp_dFi !< derivative of Lp with respect to Fi real(pReal), dimension(9,9) :: & dLp_dMstar99 !< derivative of Lp with respect to Mstar (matrix notation) @@ -499,41 +499,41 @@ subroutine constitutive_LpAndItsTangent(Lp, dLp_dTstar, dLp_dFi, Tstar6, Fi, ipc case (PLASTICITY_NONE_ID) plasticityType Lp = 0.0_pReal - dLp_dTstar = 0.0_pReal + dLp_dMstar = 0.0_pReal case (PLASTICITY_ISOTROPIC_ID) plasticityType call plastic_isotropic_LpAndItsTangent (Lp,dLp_dMstar99, math_Mandel33to6(Mstar),ipc,ip,el) - dLp_dTstar = math_Plain99to3333(dLp_dMstar99) ! ToDo: We revert here the last statement in plastic_xx_LpAndItsTanget + dLp_dMstar = math_Plain99to3333(dLp_dMstar99) ! ToDo: We revert here the last statement in plastic_xx_LpAndItsTanget case (PLASTICITY_PHENOPOWERLAW_ID) plasticityType call plastic_phenopowerlaw_LpAndItsTangent (Lp,dLp_dMstar99, math_Mandel33to6(Mstar),ipc,ip,el) - dLp_dTstar = math_Plain99to3333(dLp_dMstar99) ! ToDo: We revert here the last statement in plastic_xx_LpAndItsTanget + dLp_dMstar = math_Plain99to3333(dLp_dMstar99) ! ToDo: We revert here the last statement in plastic_xx_LpAndItsTanget case (PLASTICITY_KINEHARDENING_ID) plasticityType call plastic_kinehardening_LpAndItsTangent (Lp,dLp_dMstar99, math_Mandel33to6(Mstar),ipc,ip,el) - dLp_dTstar = math_Plain99to3333(dLp_dMstar99) ! ToDo: We revert here the last statement in plastic_xx_LpAndItsTanget + dLp_dMstar = math_Plain99to3333(dLp_dMstar99) ! ToDo: We revert here the last statement in plastic_xx_LpAndItsTanget case (PLASTICITY_NONLOCAL_ID) plasticityType call plastic_nonlocal_LpAndItsTangent (Lp,dLp_dMstar99, math_Mandel33to6(Mstar), & temperature(ho)%p(tme),ip,el) - dLp_dTstar = math_Plain99to3333(dLp_dMstar99) ! ToDo: We revert here the last statement in plastic_xx_LpAndItsTanget + dLp_dMstar = math_Plain99to3333(dLp_dMstar99) ! ToDo: We revert here the last statement in plastic_xx_LpAndItsTanget case (PLASTICITY_DISLOTWIN_ID) plasticityType call plastic_dislotwin_LpAndItsTangent (Lp,dLp_dMstar99, math_Mandel33to6(Mstar), & temperature(ho)%p(tme),ipc,ip,el) - dLp_dTstar = math_Plain99to3333(dLp_dMstar99) ! ToDo: We revert here the last statement in plastic_xx_LpAndItsTanget + dLp_dMstar = math_Plain99to3333(dLp_dMstar99) ! ToDo: We revert here the last statement in plastic_xx_LpAndItsTanget case (PLASTICITY_DISLOUCLA_ID) plasticityType call plastic_disloucla_LpAndItsTangent (Lp,dLp_dMstar99, math_Mandel33to6(Mstar), & temperature(ho)%p(tme), ipc,ip,el) - dLp_dTstar = math_Plain99to3333(dLp_dMstar99) ! ToDo: We revert here the last statement in plastic_xx_LpAndItsTanget + dLp_dMstar = math_Plain99to3333(dLp_dMstar99) ! ToDo: We revert here the last statement in plastic_xx_LpAndItsTanget end select plasticityType do concurrent(i = 1_pInt:3_pInt, j = 1_pInt:3_pInt) - dLp_dFi(i,j,1:3,1:3) = math_mul33x33(math_mul33x33(Fi,Tstar),transpose(dLp_dTstar(i,j,1:3,1:3))) + & - math_mul33x33(math_mul33x33(Fi,dLp_dTstar(i,j,1:3,1:3)),Tstar) - dLp_dTstar(i,j,1:3,1:3) = math_mul33x33(math_mul33x33(transpose(Fi),Fi),dLp_dTstar(i,j,1:3,1:3)) + dLp_dFi(i,j,1:3,1:3) = math_mul33x33(math_mul33x33(Fi,Tstar),transpose(dLp_dMstar(i,j,1:3,1:3))) + & + math_mul33x33(math_mul33x33(Fi,dLp_dMstar(i,j,1:3,1:3)),Tstar) + dLp_dMstar(i,j,1:3,1:3) = math_mul33x33(math_mul33x33(transpose(Fi),Fi),dLp_dMstar(i,j,1:3,1:3)) enddo end subroutine constitutive_LpAndItsTangent @@ -879,7 +879,6 @@ subroutine constitutive_collectDotState(Tstar6, FeArray, Fi, FpArray, subdt, sub real(pReal), intent(in), dimension(6) :: & Tstar6 !< 2nd Piola Kirchhoff stress (vector notation) real(pReal), dimension(3,3) :: & - Tstar, & !< 2nd Piola Kirchhoff stress tensor Mstar integer(pInt) :: & ho, & !< homogenization @@ -889,8 +888,7 @@ subroutine constitutive_collectDotState(Tstar6, FeArray, Fi, FpArray, subdt, sub ho = material_homog( ip,el) tme = thermalMapping(ho)%p(ip,el) - Tstar = math_Mandel6to33(Tstar6) - Mstar = math_mul33x33(math_mul33x33(transpose(Fi),Fi),Tstar) + Mstar = math_mul33x33(math_mul33x33(transpose(Fi),Fi),math_Mandel6to33(Tstar6)) plasticityType: select case (phase_plasticity(material_phase(ipc,ip,el))) @@ -986,13 +984,11 @@ subroutine constitutive_collectDeltaState(Tstar6, Fe, Fi, ipc, ip, el) Fe, & !< elastic deformation gradient Fi !< intermediate deformation gradient real(pReal), dimension(3,3) :: & - Tstar, & !< 2nd Piola Kirchhoff stress tensor Mstar integer(pInt) :: & s !< counter in source loop - Tstar = math_Mandel6to33(Tstar6) - Mstar = math_mul33x33(math_mul33x33(transpose(Fi),Fi),Tstar) + Mstar = math_mul33x33(math_mul33x33(transpose(Fi),Fi),math_Mandel6to33(Tstar6)) plasticityType: select case (phase_plasticity(material_phase(ipc,ip,el))) From 4c14f988a30d2592b7af44be149f04de2957a9dd Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 28 Aug 2018 14:54:36 +0200 Subject: [PATCH 09/13] rename according to paper --- src/constitutive.f90 | 42 +++++++++++++++++++++--------------------- src/crystallite.f90 | 10 +++++----- 2 files changed, 26 insertions(+), 26 deletions(-) diff --git a/src/constitutive.f90 b/src/constitutive.f90 index 4275f0d8d..e781d77a7 100644 --- a/src/constitutive.f90 +++ b/src/constitutive.f90 @@ -22,13 +22,13 @@ module constitutive constitutive_LpAndItsTangent, & constitutive_LiAndItsTangent, & constitutive_initialFi, & - constitutive_TandItsTangent, & + constitutive_SandItsTangents, & constitutive_collectDotState, & constitutive_collectDeltaState, & constitutive_postResults private :: & - constitutive_hooke_TandItsTangent + constitutive_hooke_SandItsTangents contains @@ -705,10 +705,10 @@ end function constitutive_initialFi !-------------------------------------------------------------------------------------------------- !> @brief returns the 2nd Piola-Kirchhoff stress tensor and its tangent with respect to -!> the elastic deformation gradient depending on the selected elastic law (so far no case switch -!! because only Hooke is implemented +!> the elastic/intermediate deformation gradients depending on the selected elastic law +!! (so far no case switch because only Hooke is implemented) !-------------------------------------------------------------------------------------------------- -subroutine constitutive_TandItsTangent(T, dT_dFe, dT_dFi, Fe, Fi, ipc, ip, el) +subroutine constitutive_SandItsTangents(S, dS_dFe, dS_dFi, Fe, Fi, ipc, ip, el) use prec, only: & pReal @@ -721,22 +721,22 @@ subroutine constitutive_TandItsTangent(T, dT_dFe, dT_dFi, Fe, Fi, ipc, ip, el) Fe, & !< elastic deformation gradient Fi !< intermediate deformation gradient real(pReal), intent(out), dimension(3,3) :: & - T !< 2nd Piola-Kirchhoff stress tensor + S !< 2nd Piola-Kirchhoff stress tensor real(pReal), intent(out), dimension(3,3,3,3) :: & - dT_dFe, & !< derivative of 2nd P-K stress with respect to elastic deformation gradient - dT_dFi !< derivative of 2nd P-K stress with respect to intermediate deformation gradient + dS_dFe, & !< derivative of 2nd P-K stress with respect to elastic deformation gradient + dS_dFi !< derivative of 2nd P-K stress with respect to intermediate deformation gradient - call constitutive_hooke_TandItsTangent(T, dT_dFe, dT_dFi, Fe, Fi, ipc, ip, el) + call constitutive_hooke_SandItsTangents(S, dS_dFe, dS_dFi, Fe, Fi, ipc, ip, el) -end subroutine constitutive_TandItsTangent +end subroutine constitutive_SandItsTangents !-------------------------------------------------------------------------------------------------- !> @brief returns the 2nd Piola-Kirchhoff stress tensor and its tangent with respect to -!> the elastic deformation gradient using hookes law +!> the elastic and intermeidate deformation gradients using Hookes law !-------------------------------------------------------------------------------------------------- -subroutine constitutive_hooke_TandItsTangent(T, dT_dFe, dT_dFi, Fe, Fi, ipc, ip, el) +subroutine constitutive_hooke_SandItsTangents(S, dS_dFe, dS_dFi, Fe, Fi, ipc, ip, el) use prec, only: & pReal use math, only : & @@ -765,10 +765,10 @@ subroutine constitutive_hooke_TandItsTangent(T, dT_dFe, dT_dFi, Fe, Fi, ipc, ip, Fe, & !< elastic deformation gradient Fi !< intermediate deformation gradient real(pReal), intent(out), dimension(3,3) :: & - T !< 2nd Piola-Kirchhoff stress tensor in lattice configuration + S !< 2nd Piola-Kirchhoff stress tensor in lattice configuration real(pReal), intent(out), dimension(3,3,3,3) :: & - dT_dFe, & !< derivative of 2nd P-K stress with respect to elastic deformation gradient - dT_dFi !< derivative of 2nd P-K stress with respect to intermediate deformation gradient + dS_dFe, & !< derivative of 2nd P-K stress with respect to elastic deformation gradient + dS_dFi !< derivative of 2nd P-K stress with respect to intermediate deformation gradient real(pReal), dimension(3,3) :: E real(pReal), dimension(3,3,3,3) :: C integer(pInt) :: & @@ -791,16 +791,16 @@ subroutine constitutive_hooke_TandItsTangent(T, dT_dFe, dT_dFi, Fe, Fi, ipc, ip, enddo DegradationLoop E = 0.5_pReal*(math_mul33x33(transpose(Fe),Fe)-math_I3) !< Green-Lagrange strain in unloaded configuration - T = math_mul3333xx33(C,math_mul33x33(math_mul33x33(transpose(Fi),E),Fi)) !< 2PK stress in lattice configuration in work conjugate with GL strain pulled back to lattice configuration + S = math_mul3333xx33(C,math_mul33x33(math_mul33x33(transpose(Fi),E),Fi)) !< 2PK stress in lattice configuration in work conjugate with GL strain pulled back to lattice configuration - dT_dFe = 0.0_pReal + dS_dFe = 0.0_pReal forall (i=1_pInt:3_pInt, j=1_pInt:3_pInt) - dT_dFe(i,j,1:3,1:3) = & - math_mul33x33(Fe,math_mul33x33(math_mul33x33(Fi,C(i,j,1:3,1:3)),transpose(Fi))) !< dT_ij/dFe_kl = C_ijmn * Fi_lm * Fi_on * Fe_ko - dT_dFi(i,j,1:3,1:3) = 2.0_pReal*math_mul33x33(math_mul33x33(E,Fi),C(i,j,1:3,1:3)) !< dT_ij/dFi_kl = C_ijln * E_km * Fe_mn + dS_dFe(i,j,1:3,1:3) = & + math_mul33x33(Fe,math_mul33x33(math_mul33x33(Fi,C(i,j,1:3,1:3)),transpose(Fi))) !< dS_ij/dFe_kl = C_ijmn * Fi_lm * Fi_on * Fe_ko + dS_dFi(i,j,1:3,1:3) = 2.0_pReal*math_mul33x33(math_mul33x33(E,Fi),C(i,j,1:3,1:3)) !< dS_ij/dFi_kl = C_ijln * E_km * Fe_mn end forall -end subroutine constitutive_hooke_TandItsTangent +end subroutine constitutive_hooke_SandItsTangents !-------------------------------------------------------------------------------------------------- diff --git a/src/crystallite.f90 b/src/crystallite.f90 index 894744f28..139d40520 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -540,7 +540,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco) phase_Nsources, & phaseAt, phasememberAt use constitutive, only: & - constitutive_TandItsTangent, & + constitutive_SandItsTangents, & constitutive_LpAndItsTangent, & constitutive_LiAndItsTangent @@ -1065,7 +1065,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco) invFp = math_inv33(crystallite_partionedFp0(1:3,1:3,c,i,e)) Fe_guess = math_mul33x33(math_mul33x33(crystallite_partionedF(1:3,1:3,c,i,e), invFp), & math_inv33(crystallite_partionedFi0(1:3,1:3,c,i,e))) - call constitutive_TandItsTangent(Tstar,dSdFe,dSdFi,Fe_guess,crystallite_partionedFi0(1:3,1:3,c,i,e),c,i,e) + call constitutive_SandItsTangents(Tstar,dSdFe,dSdFi,Fe_guess,crystallite_partionedFi0(1:3,1:3,c,i,e),c,i,e) crystallite_P(1:3,1:3,c,i,e) = math_mul33x33(math_mul33x33(crystallite_partionedF(1:3,1:3,c,i,e), invFp), & math_mul33x33(Tstar,transpose(invFp))) endif @@ -1099,7 +1099,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco) myNcomponents = homogenization_Ngrains(mesh_element(3,e)) do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed do c = 1_pInt,myNcomponents - call constitutive_TandItsTangent(temp_33,dSdFe,dSdFi,crystallite_Fe(1:3,1:3,c,i,e), & + call constitutive_SandItsTangents(temp_33,dSdFe,dSdFi,crystallite_Fe(1:3,1:3,c,i,e), & crystallite_Fi(1:3,1:3,c,i,e),c,i,e) ! call constitutive law to calculate elastic stress tangent call constitutive_LiAndItsTangent(temp_33,dLidS,dLidFi,crystallite_Tstar_v(1:6,c,i,e), & @@ -3178,7 +3178,7 @@ logical function crystallite_integrateStress(& use constitutive, only: constitutive_LpAndItsTangent, & constitutive_LiAndItsTangent, & - constitutive_TandItsTangent + constitutive_SandItsTangents use math, only: math_mul33x33, & math_mul33xx33, & math_mul3333xx3333, & @@ -3369,7 +3369,7 @@ logical function crystallite_integrateStress(& B = math_I3 - dt*Lpguess Fe = math_mul33x33(math_mul33x33(A,B), invFi_new) ! current elastic deformation tensor - call constitutive_TandItsTangent(Tstar, dT_dFe3333, dT_dFi3333, & + call constitutive_SandItsTangents(Tstar, dT_dFe3333, dT_dFi3333, & Fe, Fi_new, ipc, ip, el) ! call constitutive law to calculate 2nd Piola-Kirchhoff stress and its derivative in unloaded configuration Tstar_v = math_Mandel33to6(Tstar) From 15e2d4a7cd61d484cda9a534ddff05dbf2fcc59b Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 28 Aug 2018 15:02:58 +0200 Subject: [PATCH 10/13] Tstar is S according to the DAMASK paper --- src/constitutive.f90 | 76 ++++++++++++++++++++++---------------------- 1 file changed, 38 insertions(+), 38 deletions(-) diff --git a/src/constitutive.f90 b/src/constitutive.f90 index e781d77a7..6f207f97c 100644 --- a/src/constitutive.f90 +++ b/src/constitutive.f90 @@ -430,7 +430,7 @@ end subroutine constitutive_microstructure !-------------------------------------------------------------------------------------------------- !> @brief contains the constitutive equation for calculating the velocity gradient !-------------------------------------------------------------------------------------------------- -subroutine constitutive_LpAndItsTangent(Lp, dLp_dMstar, dLp_dFi, Tstar6, Fi, ipc, ip, el) +subroutine constitutive_LpAndItsTangent(Lp, dLp_dMstar, dLp_dFi, S6, Fi, ipc, ip, el) use prec, only: & pReal use math, only: & @@ -470,7 +470,7 @@ subroutine constitutive_LpAndItsTangent(Lp, dLp_dMstar, dLp_dFi, Tstar6, Fi, ipc ip, & !< integration point el !< element real(pReal), intent(in), dimension(6) :: & - Tstar6 !< 2nd Piola-Kirchhoff stress (vector notation) + S6 !< 2nd Piola-Kirchhoff stress (vector notation) real(pReal), intent(in), dimension(3,3) :: & Fi !< intermediate deformation gradient real(pReal), intent(out), dimension(3,3) :: & @@ -482,7 +482,7 @@ subroutine constitutive_LpAndItsTangent(Lp, dLp_dMstar, dLp_dFi, Tstar6, Fi, ipc dLp_dMstar99 !< derivative of Lp with respect to Mstar (matrix notation) real(pReal), dimension(3,3) :: & Mstar, & !< Mandel stress work conjugate with Lp - Tstar !< 2nd Piola-Kirchhoff stress + S !< 2nd Piola-Kirchhoff stress integer(pInt) :: & ho, & !< homogenization tme !< thermal member position @@ -492,8 +492,8 @@ subroutine constitutive_LpAndItsTangent(Lp, dLp_dMstar, dLp_dFi, Tstar6, Fi, ipc ho = material_homog(ip,el) tme = thermalMapping(ho)%p(ip,el) - Tstar = math_Mandel6to33(Tstar6) - Mstar = math_mul33x33(math_mul33x33(transpose(Fi),Fi),Tstar) + S = math_Mandel6to33(S6) + Mstar = math_mul33x33(math_mul33x33(transpose(Fi),Fi),S) plasticityType: select case (phase_plasticity(material_phase(ipc,ip,el))) @@ -531,8 +531,8 @@ subroutine constitutive_LpAndItsTangent(Lp, dLp_dMstar, dLp_dFi, Tstar6, Fi, ipc end select plasticityType do concurrent(i = 1_pInt:3_pInt, j = 1_pInt:3_pInt) - dLp_dFi(i,j,1:3,1:3) = math_mul33x33(math_mul33x33(Fi,Tstar),transpose(dLp_dMstar(i,j,1:3,1:3))) + & - math_mul33x33(math_mul33x33(Fi,dLp_dMstar(i,j,1:3,1:3)),Tstar) + dLp_dFi(i,j,1:3,1:3) = math_mul33x33(math_mul33x33(Fi,S),transpose(dLp_dMstar(i,j,1:3,1:3))) + & + math_mul33x33(math_mul33x33(Fi,dLp_dMstar(i,j,1:3,1:3)),S) dLp_dMstar(i,j,1:3,1:3) = math_mul33x33(math_mul33x33(transpose(Fi),Fi),dLp_dMstar(i,j,1:3,1:3)) enddo @@ -542,7 +542,7 @@ end subroutine constitutive_LpAndItsTangent !-------------------------------------------------------------------------------------------------- !> @brief contains the constitutive equation for calculating the velocity gradient !-------------------------------------------------------------------------------------------------- -subroutine constitutive_LiAndItsTangent(Li, dLi_dTstar, dLi_dFi, Tstar6, Fi, ipc, ip, el) +subroutine constitutive_LiAndItsTangent(Li, dLi_dS, dLi_dFi, S6, Fi, ipc, ip, el) use prec, only: & pReal use math, only: & @@ -580,18 +580,18 @@ subroutine constitutive_LiAndItsTangent(Li, dLi_dTstar, dLi_dFi, Tstar6, Fi, ipc ip, & !< integration point el !< element real(pReal), intent(in), dimension(6) :: & - Tstar6 !< 2nd Piola-Kirchhoff stress (vector notation) + S6 !< 2nd Piola-Kirchhoff stress (vector notation) real(pReal), intent(in), dimension(3,3) :: & Fi !< intermediate deformation gradient real(pReal), intent(out), dimension(3,3) :: & Li !< intermediate velocity gradient real(pReal), intent(out), dimension(3,3,3,3) :: & - dLi_dTstar, & !< derivative of Li with respect to Tstar + dLi_dS, & !< derivative of Li with respect to S dLi_dFi real(pReal), dimension(3,3) :: & my_Li !< intermediate velocity gradient real(pReal), dimension(3,3,3,3) :: & - my_dLi_dTstar + my_dLi_dS real(pReal), dimension(3,3) :: & FiInv, & temp_33 @@ -603,38 +603,38 @@ subroutine constitutive_LiAndItsTangent(Li, dLi_dTstar, dLi_dFi, Tstar6, Fi, ipc i, j Li = 0.0_pReal - dLi_dTstar = 0.0_pReal + dLi_dS = 0.0_pReal dLi_dFi = 0.0_pReal plasticityType: select case (phase_plasticity(material_phase(ipc,ip,el))) case (PLASTICITY_isotropic_ID) plasticityType - call plastic_isotropic_LiAndItsTangent(my_Li, my_dLi_dTstar, Tstar6, ipc, ip, el) + call plastic_isotropic_LiAndItsTangent(my_Li, my_dLi_dS, S6, ipc, ip, el) case default plasticityType my_Li = 0.0_pReal - my_dLi_dTstar = 0.0_pReal + my_dLi_dS = 0.0_pReal end select plasticityType Li = Li + my_Li - dLi_dTstar = dLi_dTstar + my_dLi_dTstar + dLi_dS = dLi_dS + my_dLi_dS KinematicsLoop: do k = 1_pInt, phase_Nkinematics(material_phase(ipc,ip,el)) kinematicsType: select case (phase_kinematics(k,material_phase(ipc,ip,el))) case (KINEMATICS_cleavage_opening_ID) kinematicsType - call kinematics_cleavage_opening_LiAndItsTangent(my_Li, my_dLi_dTstar, Tstar6, ipc, ip, el) + call kinematics_cleavage_opening_LiAndItsTangent(my_Li, my_dLi_dS, S6, ipc, ip, el) case (KINEMATICS_slipplane_opening_ID) kinematicsType - call kinematics_slipplane_opening_LiAndItsTangent(my_Li, my_dLi_dTstar, Tstar6, ipc, ip, el) + call kinematics_slipplane_opening_LiAndItsTangent(my_Li, my_dLi_dS, S6, ipc, ip, el) case (KINEMATICS_thermal_expansion_ID) kinematicsType - call kinematics_thermal_expansion_LiAndItsTangent(my_Li, my_dLi_dTstar, ipc, ip, el) + call kinematics_thermal_expansion_LiAndItsTangent(my_Li, my_dLi_dS, ipc, ip, el) case (KINEMATICS_vacancy_strain_ID) kinematicsType - call kinematics_vacancy_strain_LiAndItsTangent(my_Li, my_dLi_dTstar, ipc, ip, el) + call kinematics_vacancy_strain_LiAndItsTangent(my_Li, my_dLi_dS, ipc, ip, el) case (KINEMATICS_hydrogen_strain_ID) kinematicsType - call kinematics_hydrogen_strain_LiAndItsTangent(my_Li, my_dLi_dTstar, ipc, ip, el) + call kinematics_hydrogen_strain_LiAndItsTangent(my_Li, my_dLi_dS, ipc, ip, el) case default kinematicsType my_Li = 0.0_pReal - my_dLi_dTstar = 0.0_pReal + my_dLi_dS = 0.0_pReal end select kinematicsType Li = Li + my_Li - dLi_dTstar = dLi_dTstar + my_dLi_dTstar + dLi_dS = dLi_dS + my_dLi_dS enddo KinematicsLoop FiInv = math_inv33(Fi) @@ -642,7 +642,7 @@ subroutine constitutive_LiAndItsTangent(Li, dLi_dTstar, dLi_dFi, Tstar6, Fi, ipc Li = math_mul33x33(math_mul33x33(Fi,Li),FiInv)*detFi !< push forward to intermediate configuration temp_33 = math_mul33x33(FiInv,Li) forall(i = 1_pInt:3_pInt, j = 1_pInt:3_pInt) - dLi_dTstar(1:3,1:3,i,j) = math_mul33x33(math_mul33x33(Fi,dLi_dTstar(1:3,1:3,i,j)),FiInv)*detFi + dLi_dS(1:3,1:3,i,j) = math_mul33x33(math_mul33x33(Fi,dLi_dS(1:3,1:3,i,j)),FiInv)*detFi dLi_dFi (1:3,1:3,i,j) = dLi_dFi(1:3,1:3,i,j) + Li*FiInv(j,i) dLi_dFi (1:3,i,1:3,j) = dLi_dFi(1:3,i,1:3,j) + math_I3*temp_33(j,i) + Li*FiInv(j,i) end forall @@ -806,7 +806,7 @@ end subroutine constitutive_hooke_SandItsTangents !-------------------------------------------------------------------------------------------------- !> @brief contains the constitutive equation for calculating the rate of change of microstructure !-------------------------------------------------------------------------------------------------- -subroutine constitutive_collectDotState(Tstar6, FeArray, Fi, FpArray, subdt, subfracArray,ipc, ip, el) +subroutine constitutive_collectDotState(S6, FeArray, Fi, FpArray, subdt, subfracArray,ipc, ip, el) use prec, only: & pReal, & pLongInt @@ -877,7 +877,7 @@ subroutine constitutive_collectDotState(Tstar6, FeArray, Fi, FpArray, subdt, sub real(pReal), intent(in), dimension(3,3) :: & Fi !< intermediate deformation gradient real(pReal), intent(in), dimension(6) :: & - Tstar6 !< 2nd Piola Kirchhoff stress (vector notation) + S6 !< 2nd Piola Kirchhoff stress (vector notation) real(pReal), dimension(3,3) :: & Mstar integer(pInt) :: & @@ -888,7 +888,7 @@ subroutine constitutive_collectDotState(Tstar6, FeArray, Fi, FpArray, subdt, sub ho = material_homog( ip,el) tme = thermalMapping(ho)%p(ip,el) - Mstar = math_mul33x33(math_mul33x33(transpose(Fi),Fi),math_Mandel6to33(Tstar6)) + Mstar = math_mul33x33(math_mul33x33(transpose(Fi),Fi),math_Mandel6to33(S6)) plasticityType: select case (phase_plasticity(material_phase(ipc,ip,el))) @@ -919,7 +919,7 @@ subroutine constitutive_collectDotState(Tstar6, FeArray, Fi, FpArray, subdt, sub sourceType: select case (phase_source(s,material_phase(ipc,ip,el))) case (SOURCE_damage_anisoBrittle_ID) sourceType - call source_damage_anisoBrittle_dotState (Tstar6, ipc, ip, el) !< correct stress? + call source_damage_anisoBrittle_dotState (S6, ipc, ip, el) !< correct stress? case (SOURCE_damage_isoDuctile_ID) sourceType call source_damage_isoDuctile_dotState ( ipc, ip, el) @@ -940,7 +940,7 @@ end subroutine constitutive_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 !-------------------------------------------------------------------------------------------------- -subroutine constitutive_collectDeltaState(Tstar6, Fe, Fi, ipc, ip, el) +subroutine constitutive_collectDeltaState(S6, Fe, Fi, ipc, ip, el) use prec, only: & pReal, & pLongInt @@ -979,7 +979,7 @@ subroutine constitutive_collectDeltaState(Tstar6, Fe, Fi, ipc, ip, el) ip, & !< integration point el !< element real(pReal), intent(in), dimension(6) :: & - Tstar6 !< 2nd Piola Kirchhoff stress (vector notation) + S6 !< 2nd Piola Kirchhoff stress (vector notation) real(pReal), intent(in), dimension(3,3) :: & Fe, & !< elastic deformation gradient Fi !< intermediate deformation gradient @@ -988,7 +988,7 @@ subroutine constitutive_collectDeltaState(Tstar6, Fe, Fi, ipc, ip, el) integer(pInt) :: & s !< counter in source loop - Mstar = math_mul33x33(math_mul33x33(transpose(Fi),Fi),math_Mandel6to33(Tstar6)) + Mstar = math_mul33x33(math_mul33x33(transpose(Fi),Fi),math_Mandel6to33(S6)) plasticityType: select case (phase_plasticity(material_phase(ipc,ip,el))) @@ -1024,7 +1024,7 @@ end subroutine constitutive_collectDeltaState !-------------------------------------------------------------------------------------------------- !> @brief returns array of constitutive results !-------------------------------------------------------------------------------------------------- -function constitutive_postResults(Tstar6, FeArray, ipc, ip, el) +function constitutive_postResults(S6, FeArray, ipc, ip, el) use prec, only: & pReal use mesh, only: & @@ -1084,7 +1084,7 @@ function constitutive_postResults(Tstar6, FeArray, ipc, ip, el) real(pReal), intent(in), dimension(3,3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: & FeArray !< elastic deformation gradient real(pReal), intent(in), dimension(6) :: & - Tstar6 !< 2nd Piola Kirchhoff stress (vector notation) + S6 !< 2nd Piola Kirchhoff stress (vector notation) integer(pInt) :: & startPos, endPos integer(pInt) :: & @@ -1102,22 +1102,22 @@ function constitutive_postResults(Tstar6, FeArray, ipc, ip, el) plasticityType: select case (phase_plasticity(material_phase(ipc,ip,el))) case (PLASTICITY_ISOTROPIC_ID) plasticityType - constitutive_postResults(startPos:endPos) = plastic_isotropic_postResults(Tstar6,ipc,ip,el) + constitutive_postResults(startPos:endPos) = plastic_isotropic_postResults(S6,ipc,ip,el) case (PLASTICITY_PHENOPOWERLAW_ID) plasticityType constitutive_postResults(startPos:endPos) = & - plastic_phenopowerlaw_postResults(Tstar6,ipc,ip,el) + plastic_phenopowerlaw_postResults(S6,ipc,ip,el) case (PLASTICITY_KINEHARDENING_ID) plasticityType constitutive_postResults(startPos:endPos) = & - plastic_kinehardening_postResults(Tstar6,ipc,ip,el) + plastic_kinehardening_postResults(S6,ipc,ip,el) case (PLASTICITY_DISLOTWIN_ID) plasticityType constitutive_postResults(startPos:endPos) = & - plastic_dislotwin_postResults(Tstar6,temperature(ho)%p(tme),ipc,ip,el) + plastic_dislotwin_postResults(S6,temperature(ho)%p(tme),ipc,ip,el) case (PLASTICITY_DISLOUCLA_ID) plasticityType constitutive_postResults(startPos:endPos) = & - plastic_disloucla_postResults(Tstar6,temperature(ho)%p(tme),ipc,ip,el) + plastic_disloucla_postResults(S6,temperature(ho)%p(tme),ipc,ip,el) case (PLASTICITY_NONLOCAL_ID) plasticityType constitutive_postResults(startPos:endPos) = & - plastic_nonlocal_postResults (Tstar6,FeArray,ip,el) + plastic_nonlocal_postResults (S6,FeArray,ip,el) end select plasticityType SourceLoop: do s = 1_pInt, phase_Nsources(material_phase(ipc,ip,el)) From fc3ce546673c9b67a07492b633ed47df528ffa2d Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 28 Aug 2018 15:07:39 +0200 Subject: [PATCH 11/13] return more than one tangent --- src/constitutive.f90 | 12 ++++++------ src/crystallite.f90 | 16 ++++++++-------- 2 files changed, 14 insertions(+), 14 deletions(-) diff --git a/src/constitutive.f90 b/src/constitutive.f90 index 6f207f97c..2321bed1d 100644 --- a/src/constitutive.f90 +++ b/src/constitutive.f90 @@ -19,8 +19,8 @@ module constitutive constitutive_init, & constitutive_homogenizedC, & constitutive_microstructure, & - constitutive_LpAndItsTangent, & - constitutive_LiAndItsTangent, & + constitutive_LpAndItsTangents, & + constitutive_LiAndItsTangents, & constitutive_initialFi, & constitutive_SandItsTangents, & constitutive_collectDotState, & @@ -430,7 +430,7 @@ end subroutine constitutive_microstructure !-------------------------------------------------------------------------------------------------- !> @brief contains the constitutive equation for calculating the velocity gradient !-------------------------------------------------------------------------------------------------- -subroutine constitutive_LpAndItsTangent(Lp, dLp_dMstar, dLp_dFi, S6, Fi, ipc, ip, el) +subroutine constitutive_LpAndItsTangents(Lp, dLp_dMstar, dLp_dFi, S6, Fi, ipc, ip, el) use prec, only: & pReal use math, only: & @@ -536,13 +536,13 @@ subroutine constitutive_LpAndItsTangent(Lp, dLp_dMstar, dLp_dFi, S6, Fi, ipc, ip dLp_dMstar(i,j,1:3,1:3) = math_mul33x33(math_mul33x33(transpose(Fi),Fi),dLp_dMstar(i,j,1:3,1:3)) enddo -end subroutine constitutive_LpAndItsTangent +end subroutine constitutive_LpAndItsTangents !-------------------------------------------------------------------------------------------------- !> @brief contains the constitutive equation for calculating the velocity gradient !-------------------------------------------------------------------------------------------------- -subroutine constitutive_LiAndItsTangent(Li, dLi_dS, dLi_dFi, S6, Fi, ipc, ip, el) +subroutine constitutive_LiAndItsTangents(Li, dLi_dS, dLi_dFi, S6, Fi, ipc, ip, el) use prec, only: & pReal use math, only: & @@ -647,7 +647,7 @@ subroutine constitutive_LiAndItsTangent(Li, dLi_dS, dLi_dFi, S6, Fi, ipc, ip, el dLi_dFi (1:3,i,1:3,j) = dLi_dFi(1:3,i,1:3,j) + math_I3*temp_33(j,i) + Li*FiInv(j,i) end forall -end subroutine constitutive_LiAndItsTangent +end subroutine constitutive_LiAndItsTangents !-------------------------------------------------------------------------------------------------- diff --git a/src/crystallite.f90 b/src/crystallite.f90 index 139d40520..575a9f014 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -541,8 +541,8 @@ subroutine crystallite_stressAndItsTangent(updateJaco) phaseAt, phasememberAt use constitutive, only: & constitutive_SandItsTangents, & - constitutive_LpAndItsTangent, & - constitutive_LiAndItsTangent + constitutive_LpAndItsTangents, & + constitutive_LiAndItsTangents implicit none logical, intent(in) :: & @@ -1102,7 +1102,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco) call constitutive_SandItsTangents(temp_33,dSdFe,dSdFi,crystallite_Fe(1:3,1:3,c,i,e), & crystallite_Fi(1:3,1:3,c,i,e),c,i,e) ! call constitutive law to calculate elastic stress tangent - call constitutive_LiAndItsTangent(temp_33,dLidS,dLidFi,crystallite_Tstar_v(1:6,c,i,e), & + call constitutive_LiAndItsTangents(temp_33,dLidS,dLidFi,crystallite_Tstar_v(1:6,c,i,e), & crystallite_Fi(1:3,1:3,c,i,e), & c,i,e) ! call constitutive law to calculate Li tangent in lattice configuration if (sum(abs(dLidS)) < tol_math_check) then @@ -1129,7 +1129,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco) dLidS = math_mul3333xx3333(dLidFi,dFidS) + dLidS endif - call constitutive_LpAndItsTangent(temp_33,dLpdS,dLpdFi,crystallite_Tstar_v(1:6,c,i,e), & + call constitutive_LpAndItsTangents(temp_33,dLpdS,dLpdFi,crystallite_Tstar_v(1:6,c,i,e), & crystallite_Fi(1:3,1:3,c,i,e),c,i,e) ! call constitutive law to calculate Lp tangent in lattice configuration dLpdS = math_mul3333xx3333(dLpdFi,dFidS) + dLpdS @@ -3176,8 +3176,8 @@ logical function crystallite_integrateStress(& debug_levelExtensive, & debug_levelSelective - use constitutive, only: constitutive_LpAndItsTangent, & - constitutive_LiAndItsTangent, & + use constitutive, only: constitutive_LpAndItsTangents, & + constitutive_LiAndItsTangents, & constitutive_SandItsTangents use math, only: math_mul33x33, & math_mul33xx33, & @@ -3386,7 +3386,7 @@ logical function crystallite_integrateStress(& write(6,'(a,/,6(e20.10,1x))') '<< CRYST >> Tstar', Tstar_v endif #endif - call constitutive_LpAndItsTangent(Lp_constitutive, dLp_dT3333, dLp_dFi3333, & + call constitutive_LpAndItsTangents(Lp_constitutive, dLp_dT3333, dLp_dFi3333, & Tstar_v, Fi_new, ipc, ip, el) #ifdef DEBUG @@ -3488,7 +3488,7 @@ logical function crystallite_integrateStress(& !* calculate intermediate velocity gradient and its tangent from constitutive law - call constitutive_LiAndItsTangent(Li_constitutive, dLi_dT3333, dLi_dFi3333, & + call constitutive_LiAndItsTangents(Li_constitutive, dLi_dT3333, dLi_dFi3333, & Tstar_v, Fi_new, ipc, ip, el) #ifdef DEBUG From b5e4ad247d454a83632237db350ec8681b9e17bf Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 28 Aug 2018 15:18:55 +0200 Subject: [PATCH 12/13] LpAndItsTangent naming convention following DAMASK paper --- src/constitutive.f90 | 56 +++++++++++++++++++++++--------------------- 1 file changed, 29 insertions(+), 27 deletions(-) diff --git a/src/constitutive.f90 b/src/constitutive.f90 index 2321bed1d..a16b1df30 100644 --- a/src/constitutive.f90 +++ b/src/constitutive.f90 @@ -430,7 +430,7 @@ end subroutine constitutive_microstructure !-------------------------------------------------------------------------------------------------- !> @brief contains the constitutive equation for calculating the velocity gradient !-------------------------------------------------------------------------------------------------- -subroutine constitutive_LpAndItsTangents(Lp, dLp_dMstar, dLp_dFi, S6, Fi, ipc, ip, el) +subroutine constitutive_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, S6, Fi, ipc, ip, el) use prec, only: & pReal use math, only: & @@ -476,12 +476,14 @@ subroutine constitutive_LpAndItsTangents(Lp, dLp_dMstar, dLp_dFi, S6, Fi, ipc, i real(pReal), intent(out), dimension(3,3) :: & Lp !< plastic velocity gradient real(pReal), intent(out), dimension(3,3,3,3) :: & - dLp_dMstar, & !< derivative of Lp with respect to Mandel stress + dLp_dS, & dLp_dFi !< derivative of Lp with respect to Fi + real(pReal), dimension(3,3,3,3) :: & + dLp_dMp !< derivative of Lp with respect to Mandel stress real(pReal), dimension(9,9) :: & - dLp_dMstar99 !< derivative of Lp with respect to Mstar (matrix notation) + dLp_dMp99 !< derivative of Lp with respect to Mstar (matrix notation) real(pReal), dimension(3,3) :: & - Mstar, & !< Mandel stress work conjugate with Lp + Mp, & !< Mandel stress work conjugate with Lp S !< 2nd Piola-Kirchhoff stress integer(pInt) :: & ho, & !< homogenization @@ -493,47 +495,47 @@ subroutine constitutive_LpAndItsTangents(Lp, dLp_dMstar, dLp_dFi, S6, Fi, ipc, i tme = thermalMapping(ho)%p(ip,el) S = math_Mandel6to33(S6) - Mstar = math_mul33x33(math_mul33x33(transpose(Fi),Fi),S) + Mp = math_mul33x33(math_mul33x33(transpose(Fi),Fi),S) plasticityType: select case (phase_plasticity(material_phase(ipc,ip,el))) case (PLASTICITY_NONE_ID) plasticityType Lp = 0.0_pReal - dLp_dMstar = 0.0_pReal + dLp_dMp = 0.0_pReal case (PLASTICITY_ISOTROPIC_ID) plasticityType - call plastic_isotropic_LpAndItsTangent (Lp,dLp_dMstar99, math_Mandel33to6(Mstar),ipc,ip,el) - dLp_dMstar = math_Plain99to3333(dLp_dMstar99) ! ToDo: We revert here the last statement in plastic_xx_LpAndItsTanget + call plastic_isotropic_LpAndItsTangent (Lp,dLp_dMp99, math_Mandel33to6(Mp),ipc,ip,el) + dLp_dMp = math_Plain99to3333(dLp_dMp99) ! ToDo: We revert here the last statement in plastic_xx_LpAndItsTanget case (PLASTICITY_PHENOPOWERLAW_ID) plasticityType - call plastic_phenopowerlaw_LpAndItsTangent (Lp,dLp_dMstar99, math_Mandel33to6(Mstar),ipc,ip,el) - dLp_dMstar = math_Plain99to3333(dLp_dMstar99) ! ToDo: We revert here the last statement in plastic_xx_LpAndItsTanget + call plastic_phenopowerlaw_LpAndItsTangent (Lp,dLp_dMp99, math_Mandel33to6(Mp),ipc,ip,el) + dLp_dMp = math_Plain99to3333(dLp_dMp99) ! ToDo: We revert here the last statement in plastic_xx_LpAndItsTanget case (PLASTICITY_KINEHARDENING_ID) plasticityType - call plastic_kinehardening_LpAndItsTangent (Lp,dLp_dMstar99, math_Mandel33to6(Mstar),ipc,ip,el) - dLp_dMstar = math_Plain99to3333(dLp_dMstar99) ! ToDo: We revert here the last statement in plastic_xx_LpAndItsTanget + call plastic_kinehardening_LpAndItsTangent (Lp,dLp_dMp99, math_Mandel33to6(Mp),ipc,ip,el) + dLp_dMp = math_Plain99to3333(dLp_dMp99) ! ToDo: We revert here the last statement in plastic_xx_LpAndItsTanget case (PLASTICITY_NONLOCAL_ID) plasticityType - call plastic_nonlocal_LpAndItsTangent (Lp,dLp_dMstar99, math_Mandel33to6(Mstar), & + call plastic_nonlocal_LpAndItsTangent (Lp,dLp_dMp99, math_Mandel33to6(Mp), & temperature(ho)%p(tme),ip,el) - dLp_dMstar = math_Plain99to3333(dLp_dMstar99) ! ToDo: We revert here the last statement in plastic_xx_LpAndItsTanget + dLp_dMp = math_Plain99to3333(dLp_dMp99) ! ToDo: We revert here the last statement in plastic_xx_LpAndItsTanget case (PLASTICITY_DISLOTWIN_ID) plasticityType - call plastic_dislotwin_LpAndItsTangent (Lp,dLp_dMstar99, math_Mandel33to6(Mstar), & + call plastic_dislotwin_LpAndItsTangent (Lp,dLp_dMp99, math_Mandel33to6(Mp), & temperature(ho)%p(tme),ipc,ip,el) - dLp_dMstar = math_Plain99to3333(dLp_dMstar99) ! ToDo: We revert here the last statement in plastic_xx_LpAndItsTanget + dLp_dMp = math_Plain99to3333(dLp_dMp99) ! ToDo: We revert here the last statement in plastic_xx_LpAndItsTanget case (PLASTICITY_DISLOUCLA_ID) plasticityType - call plastic_disloucla_LpAndItsTangent (Lp,dLp_dMstar99, math_Mandel33to6(Mstar), & + call plastic_disloucla_LpAndItsTangent (Lp,dLp_dMp99, math_Mandel33to6(Mp), & temperature(ho)%p(tme), ipc,ip,el) - dLp_dMstar = math_Plain99to3333(dLp_dMstar99) ! ToDo: We revert here the last statement in plastic_xx_LpAndItsTanget + dLp_dMp = math_Plain99to3333(dLp_dMp99) ! ToDo: We revert here the last statement in plastic_xx_LpAndItsTanget end select plasticityType do concurrent(i = 1_pInt:3_pInt, j = 1_pInt:3_pInt) - dLp_dFi(i,j,1:3,1:3) = math_mul33x33(math_mul33x33(Fi,S),transpose(dLp_dMstar(i,j,1:3,1:3))) + & - math_mul33x33(math_mul33x33(Fi,dLp_dMstar(i,j,1:3,1:3)),S) - dLp_dMstar(i,j,1:3,1:3) = math_mul33x33(math_mul33x33(transpose(Fi),Fi),dLp_dMstar(i,j,1:3,1:3)) + dLp_dFi(i,j,1:3,1:3) = math_mul33x33(math_mul33x33(Fi,S),transpose(dLp_dMp(i,j,1:3,1:3))) + & + math_mul33x33(math_mul33x33(Fi,dLp_dMp(i,j,1:3,1:3)),S) + dLp_dS(i,j,1:3,1:3) = math_mul33x33(math_mul33x33(transpose(Fi),Fi),dLp_dMp(i,j,1:3,1:3)) ! ToDo: @PS: why not: dLp_dMp:(FiT Fi) enddo end subroutine constitutive_LpAndItsTangents @@ -604,7 +606,7 @@ subroutine constitutive_LiAndItsTangents(Li, dLi_dS, dLi_dFi, S6, Fi, ipc, ip, e Li = 0.0_pReal dLi_dS = 0.0_pReal - dLi_dFi = 0.0_pReal + dLi_dFi = 0.0_pReal plasticityType: select case (phase_plasticity(material_phase(ipc,ip,el))) case (PLASTICITY_isotropic_ID) plasticityType @@ -641,11 +643,11 @@ subroutine constitutive_LiAndItsTangents(Li, dLi_dS, dLi_dFi, S6, Fi, ipc, ip, e detFi = math_det33(Fi) Li = math_mul33x33(math_mul33x33(Fi,Li),FiInv)*detFi !< push forward to intermediate configuration temp_33 = math_mul33x33(FiInv,Li) - forall(i = 1_pInt:3_pInt, j = 1_pInt:3_pInt) - dLi_dS(1:3,1:3,i,j) = math_mul33x33(math_mul33x33(Fi,dLi_dS(1:3,1:3,i,j)),FiInv)*detFi - dLi_dFi (1:3,1:3,i,j) = dLi_dFi(1:3,1:3,i,j) + Li*FiInv(j,i) - dLi_dFi (1:3,i,1:3,j) = dLi_dFi(1:3,i,1:3,j) + math_I3*temp_33(j,i) + Li*FiInv(j,i) - end forall + do concurrent(i = 1_pInt:3_pInt, j = 1_pInt:3_pInt) + dLi_dS(1:3,1:3,i,j) = math_mul33x33(math_mul33x33(Fi,dLi_dS(1:3,1:3,i,j)),FiInv)*detFi + dLi_dFi(1:3,1:3,i,j) = dLi_dFi(1:3,1:3,i,j) + Li*FiInv(j,i) + dLi_dFi(1:3,i,1:3,j) = dLi_dFi(1:3,i,1:3,j) + math_I3*temp_33(j,i) + Li*FiInv(j,i) + end do end subroutine constitutive_LiAndItsTangents From b884349e7b308e96a8d38e19df2ca65c83286f75 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 29 Aug 2018 13:16:37 +0200 Subject: [PATCH 13/13] only renaming 3333 not needed for dX_dY if X and Y are 3x3 tensors PK2 stress is S not T according to the DAMASK paper --- src/crystallite.f90 | 64 ++++++++++++++++++++++----------------------- 1 file changed, 32 insertions(+), 32 deletions(-) diff --git a/src/crystallite.f90 b/src/crystallite.f90 index 575a9f014..9eb6645de 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -3240,15 +3240,15 @@ logical function crystallite_integrateStress(& real(pReal), dimension(9,9) :: dRLp_dLp, & ! partial derivative of residuum (Jacobian for NEwton-Raphson scheme) dRLp_dLp2, & ! working copy of dRdLp dRLi_dLi ! partial derivative of residuumI (Jacobian for NEwton-Raphson scheme) - real(pReal), dimension(3,3,3,3):: dT_dFe3333, & ! partial derivative of 2nd Piola-Kirchhoff stress - dT_dFi3333, & - dFe_dLp3333, & ! partial derivative of elastic deformation gradient - dFe_dLi3333, & - dFi_dLi3333, & - dLp_dFi3333, & - dLi_dFi3333, & - dLp_dT3333, & - dLi_dT3333 + real(pReal), dimension(3,3,3,3):: dS_dFe, & ! partial derivative of 2nd Piola-Kirchhoff stress + dS_dFi, & + dFe_dLp, & ! partial derivative of elastic deformation gradient + dFe_dLi, & + dFi_dLi, & + dLp_dFi, & + dLi_dFi, & + dLp_dS, & + dLi_dS real(pReal) detInvFi, & ! determinant of InvFi steplengthLp, & steplengthLi, & @@ -3369,7 +3369,7 @@ logical function crystallite_integrateStress(& B = math_I3 - dt*Lpguess Fe = math_mul33x33(math_mul33x33(A,B), invFi_new) ! current elastic deformation tensor - call constitutive_SandItsTangents(Tstar, dT_dFe3333, dT_dFi3333, & + call constitutive_SandItsTangents(Tstar, dS_dFe, dS_dFi, & Fe, Fi_new, ipc, ip, el) ! call constitutive law to calculate 2nd Piola-Kirchhoff stress and its derivative in unloaded configuration Tstar_v = math_Mandel33to6(Tstar) @@ -3386,7 +3386,7 @@ logical function crystallite_integrateStress(& write(6,'(a,/,6(e20.10,1x))') '<< CRYST >> Tstar', Tstar_v endif #endif - call constitutive_LpAndItsTangents(Lp_constitutive, dLp_dT3333, dLp_dFi3333, & + call constitutive_LpAndItsTangents(Lp_constitutive, dLp_dS, dLp_dFi, & Tstar_v, Fi_new, ipc, ip, el) #ifdef DEBUG @@ -3437,18 +3437,18 @@ logical function crystallite_integrateStress(& !* calculate Jacobian for correction term if (mod(jacoCounterLp, iJacoLpresiduum) == 0_pInt) then - dFe_dLp3333 = 0.0_pReal + dFe_dLp = 0.0_pReal forall(o=1_pInt:3_pInt,p=1_pInt:3_pInt) & - dFe_dLp3333(o,1:3,p,1:3) = A(o,p)*transpose(invFi_new) ! dFe_dLp(i,j,k,l) = -dt * A(i,k) invFi(l,j) - dFe_dLp3333 = - dt * dFe_dLp3333 + dFe_dLp(o,1:3,p,1:3) = A(o,p)*transpose(invFi_new) ! dFe_dLp(i,j,k,l) = -dt * A(i,k) invFi(l,j) + dFe_dLp = - dt * dFe_dLp dRLp_dLp = math_identity2nd(9_pInt) & - - math_Plain3333to99(math_mul3333xx3333(math_mul3333xx3333(dLp_dT3333,dT_dFe3333),dFe_dLp3333)) + - math_Plain3333to99(math_mul3333xx3333(math_mul3333xx3333(dLp_dS,dS_dFe),dFe_dLp)) #ifdef DEBUG if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt & .and. ((el == debug_e .and. ip == debug_i .and. ipc == debug_g) & .or. .not. iand(debug_level(debug_crystallite), debug_levelSelective) /= 0_pInt)) then - write(6,'(a,/,9(12x,9(e12.4,1x)/))') '<< CRYST >> dLp_dT', math_Plain3333to99(dLp_dT3333) - write(6,'(a,1x,e20.10)') '<< CRYST >> dLp_dT norm', norm2(math_Plain3333to99(dLp_dT3333)) + write(6,'(a,/,9(12x,9(e12.4,1x)/))') '<< CRYST >> dLp_dS', math_Plain3333to99(dLp_dS) + write(6,'(a,1x,e20.10)') '<< CRYST >> dLp_dS norm', norm2(math_Plain3333to99(dLp_dS)) write(6,'(a,/,9(12x,9(e12.4,1x)/))') '<< CRYST >> dRLp_dLp', dRLp_dLp - math_identity2nd(9_pInt) write(6,'(a,1x,e20.10)') '<< CRYST >> dRLp_dLp norm', norm2(dRLp_dLp - math_identity2nd(9_pInt)) endif @@ -3466,9 +3466,9 @@ logical function crystallite_integrateStress(& .or. .not. iand(debug_level(debug_crystallite), debug_levelSelective) /= 0_pInt)) then write(6,*) write(6,'(a,/,9(12x,9(e15.3,1x)/))') '<< CRYST >> dR_dLp',transpose(dRLp_dLp) - write(6,'(a,/,9(12x,9(e15.3,1x)/))') '<< CRYST >> dFe_dLp',transpose(math_Plain3333to99(dFe_dLp3333)) - write(6,'(a,/,9(12x,9(e15.3,1x)/))') '<< CRYST >> dT_dFe_constitutive',transpose(math_Plain3333to99(dT_dFe3333)) - write(6,'(a,/,9(12x,9(e15.3,1x)/))') '<< CRYST >> dLp_dT_constitutive',transpose(math_Plain3333to99(dLp_dT3333)) + write(6,'(a,/,9(12x,9(e15.3,1x)/))') '<< CRYST >> dFe_dLp',transpose(math_Plain3333to99(dFe_dLp)) + write(6,'(a,/,9(12x,9(e15.3,1x)/))') '<< CRYST >> dS_dFe_constitutive',transpose(math_Plain3333to99(dS_dFe)) + write(6,'(a,/,9(12x,9(e15.3,1x)/))') '<< CRYST >> dLp_dS_constitutive',transpose(math_Plain3333to99(dLp_dS)) write(6,'(a,/,3(12x,3(e20.7,1x)/))') '<< CRYST >> A',transpose(A) write(6,'(a,/,3(12x,3(e20.7,1x)/))') '<< CRYST >> B',transpose(B) write(6,'(a,/,3(12x,3(e20.7,1x)/))') '<< CRYST >> Lp_constitutive',transpose(Lp_constitutive) @@ -3488,7 +3488,7 @@ logical function crystallite_integrateStress(& !* calculate intermediate velocity gradient and its tangent from constitutive law - call constitutive_LiAndItsTangents(Li_constitutive, dLi_dT3333, dLi_dFi3333, & + call constitutive_LiAndItsTangents(Li_constitutive, dLi_dS, dLi_dFi, & Tstar_v, Fi_new, ipc, ip, el) #ifdef DEBUG @@ -3523,19 +3523,19 @@ logical function crystallite_integrateStress(& if (mod(jacoCounterLi, iJacoLpresiduum) == 0_pInt) then temp_33 = math_mul33x33(math_mul33x33(A,B),invFi_current) - dFe_dLi3333 = 0.0_pReal - dFi_dLi3333 = 0.0_pReal + dFe_dLi = 0.0_pReal + dFi_dLi = 0.0_pReal forall(o=1_pInt:3_pInt,p=1_pInt:3_pInt) - dFe_dLi3333(1:3,o,1:3,p) = -dt*math_I3(o,p)*temp_33 ! dFe_dLp(i,j,k,l) = -dt * A(i,k) invFi(l,j) - dFi_dLi3333(1:3,o,1:3,p) = -dt*math_I3(o,p)*invFi_current + dFe_dLi(1:3,o,1:3,p) = -dt*math_I3(o,p)*temp_33 ! dFe_dLp(i,j,k,l) = -dt * A(i,k) invFi(l,j) + dFi_dLi(1:3,o,1:3,p) = -dt*math_I3(o,p)*invFi_current end forall forall(o=1_pInt:3_pInt,p=1_pInt:3_pInt) & - dFi_dLi3333(1:3,1:3,o,p) = math_mul33x33(math_mul33x33(Fi_new,dFi_dLi3333(1:3,1:3,o,p)),Fi_new) + dFi_dLi(1:3,1:3,o,p) = math_mul33x33(math_mul33x33(Fi_new,dFi_dLi(1:3,1:3,o,p)),Fi_new) dRLi_dLi = math_identity2nd(9_pInt) & - - math_Plain3333to99(math_mul3333xx3333(dLi_dT3333, math_mul3333xx3333(dT_dFe3333, dFe_dLi3333) + & - math_mul3333xx3333(dT_dFi3333, dFi_dLi3333))) & - - math_Plain3333to99(math_mul3333xx3333(dLi_dFi3333, dFi_dLi3333)) + - math_Plain3333to99(math_mul3333xx3333(dLi_dS, math_mul3333xx3333(dS_dFe, dFe_dLi) + & + math_mul3333xx3333(dS_dFi, dFi_dLi))) & + - math_Plain3333to99(math_mul3333xx3333(dLi_dFi, dFi_dLi)) work = math_plain33to9(residuumLi) call dgesv(9,1,dRLi_dLi,9,ipiv,work,9,ierr) ! solve dRLi/dLp * delta Li = -res for delta Li if (ierr /= 0_pInt) then @@ -3548,9 +3548,9 @@ logical function crystallite_integrateStress(& .or. .not. iand(debug_level(debug_crystallite), debug_levelSelective) /= 0_pInt)) then write(6,*) write(6,'(a,/,9(12x,9(e15.3,1x)/))') '<< CRYST >> dR_dLi',transpose(dRLi_dLi) - write(6,'(a,/,9(12x,9(e15.3,1x)/))') '<< CRYST >> dFe_dLi',transpose(math_Plain3333to99(dFe_dLi3333)) - write(6,'(a,/,9(12x,9(e15.3,1x)/))') '<< CRYST >> dT_dFi_constitutive',transpose(math_Plain3333to99(dT_dFi3333)) - write(6,'(a,/,9(12x,9(e15.3,1x)/))') '<< CRYST >> dLi_dT_constitutive',transpose(math_Plain3333to99(dLi_dT3333)) + write(6,'(a,/,9(12x,9(e15.3,1x)/))') '<< CRYST >> dFe_dLi',transpose(math_Plain3333to99(dFe_dLi)) + write(6,'(a,/,9(12x,9(e15.3,1x)/))') '<< CRYST >> dS_dFi_constitutive',transpose(math_Plain3333to99(dS_dFi)) + write(6,'(a,/,9(12x,9(e15.3,1x)/))') '<< CRYST >> dLi_dS_constitutive',transpose(math_Plain3333to99(dLi_dS)) write(6,'(a,/,3(12x,3(e20.7,1x)/))') '<< CRYST >> Li_constitutive',transpose(Li_constitutive) write(6,'(a,/,3(12x,3(e20.7,1x)/))') '<< CRYST >> Liguess',transpose(Liguess) endif