From d585deee7e4f8fbdf631fafde3a1a4802b6630bf Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 25 Aug 2018 13:55:39 +0200 Subject: [PATCH 01/35] 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/35] 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/35] 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/35] 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/35] 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/35] 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/35] 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/35] 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/35] 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/35] 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/35] 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/35] 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/35] 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 From 3ef764da600842a925873b8ee71ddc103bc81fae Mon Sep 17 00:00:00 2001 From: Test User Date: Mon, 17 Sep 2018 01:53:58 +0200 Subject: [PATCH 14/35] [skip ci] updated version information after successful test of v2.0.2-558-gdfba3c9b --- VERSION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/VERSION b/VERSION index 119546fa7..429ed3632 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -v2.0.2-540-gce6e6679 +v2.0.2-558-gdfba3c9b From 245676728d7e3e6db0847fd1c28f6c5009fec22d Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Mon, 17 Sep 2018 23:28:00 +0200 Subject: [PATCH 15/35] not needed --- src/FEM_mech.f90 | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/FEM_mech.f90 b/src/FEM_mech.f90 index 14fbc6bbc..e10071eb8 100644 --- a/src/FEM_mech.f90 +++ b/src/FEM_mech.f90 @@ -5,7 +5,6 @@ !> @brief FEM PETSc solver !-------------------------------------------------------------------------------------------------- module FEM_mech -use PETSC #include #include #include @@ -15,7 +14,6 @@ use PETSC use PETScsnes use PETScDM use PETScDMplex - use PETSC use prec, only: & pInt, & pReal From 82cdc551a54979b49e5df231dc6b4d6d68190e69 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 18 Sep 2018 05:17:48 +0200 Subject: [PATCH 16/35] MPI Communicator only needed for PETSc 3.10. thx to Jaeyong for figuring this out --- src/FEM_mech.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/FEM_mech.f90 b/src/FEM_mech.f90 index e10071eb8..0fdaf6283 100644 --- a/src/FEM_mech.f90 +++ b/src/FEM_mech.f90 @@ -140,7 +140,7 @@ subroutine FEM_mech_init(fieldBC) nc = dimPlex call PetscQuadratureSetData(mechQuad,dimPlex,nc,nQuadrature,qPointsP,qWeightsP,ierr) CHKERRQ(ierr) - call PetscFECreateDefault(PETSC_COMM_SELF,mech_mesh,dimPlex,nc,PETSC_TRUE,prefix, & + call PetscFECreateDefault(mech_mesh,dimPlex,nc,PETSC_TRUE,prefix, & integrationOrder,mechFE,ierr); CHKERRQ(ierr) call PetscFESetQuadrature(mechFE,mechQuad,ierr); CHKERRQ(ierr) call PetscFEGetDimension(mechFE,nBasis,ierr); CHKERRQ(ierr) From 18db5f565214b183ae055888daf6adb3a067c608 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Mon, 27 Aug 2018 15:52:53 +0200 Subject: [PATCH 17/35] not needed --- src/prec.f90 | 10 ---------- 1 file changed, 10 deletions(-) diff --git a/src/prec.f90 b/src/prec.f90 index 857ec9559..0f942b3c1 100644 --- a/src/prec.f90 +++ b/src/prec.f90 @@ -87,16 +87,6 @@ module prec integer(pInt), pointer, dimension(:,:,:) :: p end type -#ifdef FEM - type, public :: tOutputData - integer(pInt) :: & - sizeIpCells = 0_pInt , & - sizeResults = 0_pInt - real(pReal), allocatable, dimension(:,:) :: & - output !< output data - end type -#endif - public :: & prec_init, & dEq, & From 20f0bee459f35f706ac21ae1470f5673576fbd24 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 19 Sep 2018 06:19:40 +0200 Subject: [PATCH 18/35] fallback dPdF not needed save a lot of memory --- src/crystallite.f90 | 5 ----- 1 file changed, 5 deletions(-) diff --git a/src/crystallite.f90 b/src/crystallite.f90 index 94ba805e8..af39fd572 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -71,8 +71,6 @@ module crystallite crystallite_dPdF, & !< current individual dPdF per grain (end of converged time step) crystallite_dPdF0, & !< individual dPdF per grain at start of FE inc crystallite_partioneddPdF0 !< individual dPdF per grain at start of homog inc - real(pReal), dimension(:,:,:,:,:,:,:), allocatable, private :: & - crystallite_fallbackdPdF !< dPdF fallback for non-converged grains (elastic prediction) logical, dimension(:,:,:), allocatable, public :: & crystallite_requested !< flag to request crystallite calculation logical, dimension(:,:,:), allocatable, public, protected :: & @@ -247,7 +245,6 @@ subroutine crystallite_init allocate(crystallite_dPdF(3,3,3,3,cMax,iMax,eMax), source=0.0_pReal) allocate(crystallite_dPdF0(3,3,3,3,cMax,iMax,eMax), source=0.0_pReal) allocate(crystallite_partioneddPdF0(3,3,3,3,cMax,iMax,eMax),source=0.0_pReal) - allocate(crystallite_fallbackdPdF(3,3,3,3,cMax,iMax,eMax), source=0.0_pReal) allocate(crystallite_dt(cMax,iMax,eMax), source=0.0_pReal) allocate(crystallite_subdt(cMax,iMax,eMax), source=0.0_pReal) allocate(crystallite_subFrac(cMax,iMax,eMax), source=0.0_pReal) @@ -454,7 +451,6 @@ subroutine crystallite_init ! enddo call crystallite_stressAndItsTangent(.true.) ! request elastic answers - crystallite_fallbackdPdF = crystallite_dPdF ! use initial elastic stiffness as fallback !-------------------------------------------------------------------------------------------------- ! debug output @@ -490,7 +486,6 @@ subroutine crystallite_init write(6,'(a35,1x,7(i8,1x))') 'crystallite_dPdF: ', shape(crystallite_dPdF) write(6,'(a35,1x,7(i8,1x))') 'crystallite_dPdF0: ', shape(crystallite_dPdF0) write(6,'(a35,1x,7(i8,1x))') 'crystallite_partioneddPdF0: ', shape(crystallite_partioneddPdF0) - write(6,'(a35,1x,7(i8,1x))') 'crystallite_fallbackdPdF: ', shape(crystallite_fallbackdPdF) write(6,'(a35,1x,7(i8,1x))') 'crystallite_orientation: ', shape(crystallite_orientation) write(6,'(a35,1x,7(i8,1x))') 'crystallite_orientation0: ', shape(crystallite_orientation0) write(6,'(a35,1x,7(i8,1x))') 'crystallite_rotation: ', shape(crystallite_rotation) From 431065e8928210b2568c5a176bf93177e2edcb74 Mon Sep 17 00:00:00 2001 From: Test User Date: Wed, 19 Sep 2018 10:22:02 +0200 Subject: [PATCH 19/35] [skip ci] updated version information after successful test of v2.0.2-560-g20f0bee4 --- VERSION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/VERSION b/VERSION index 429ed3632..cca202405 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -v2.0.2-558-gdfba3c9b +v2.0.2-560-g20f0bee4 From a8fb7d7adeaade3293da111c5eb88a049d127e22 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 19 Sep 2018 14:21:10 +0200 Subject: [PATCH 20/35] not needed but I'm under the impression that the compiler removes such things anyway --- src/crystallite.f90 | 9 --------- 1 file changed, 9 deletions(-) diff --git a/src/crystallite.f90 b/src/crystallite.f90 index af39fd572..d341f3653 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -57,7 +57,6 @@ module crystallite crystallite_Fe, & !< current "elastic" def grad (end of converged time step) crystallite_P !< 1st Piola-Kirchhoff stress per grain real(pReal), dimension(:,:,:,:,:), allocatable, private :: & - crystallite_subFe0,& !< "elastic" def grad at start of crystallite inc crystallite_invFp, & !< inverse of current plastic def grad (end of converged time step) crystallite_subFp0,& !< plastic def grad at start of crystallite inc crystallite_invFi, & !< inverse of current intermediate def grad (end of converged time step) @@ -233,7 +232,6 @@ subroutine crystallite_init allocate(crystallite_Fi(3,3,cMax,iMax,eMax), source=0.0_pReal) allocate(crystallite_invFi(3,3,cMax,iMax,eMax), source=0.0_pReal) allocate(crystallite_Fe(3,3,cMax,iMax,eMax), source=0.0_pReal) - allocate(crystallite_subFe0(3,3,cMax,iMax,eMax), source=0.0_pReal) allocate(crystallite_Lp0(3,3,cMax,iMax,eMax), source=0.0_pReal) allocate(crystallite_partionedLp0(3,3,cMax,iMax,eMax), source=0.0_pReal) allocate(crystallite_subLp0(3,3,cMax,iMax,eMax), source=0.0_pReal) @@ -473,7 +471,6 @@ subroutine crystallite_init write(6,'(a35,1x,7(i8,1x))') 'crystallite_partionedLi0: ', shape(crystallite_partionedLi0) write(6,'(a35,1x,7(i8,1x))') 'crystallite_subF: ', shape(crystallite_subF) write(6,'(a35,1x,7(i8,1x))') 'crystallite_subF0: ', shape(crystallite_subF0) - write(6,'(a35,1x,7(i8,1x))') 'crystallite_subFe0: ', shape(crystallite_subFe0) write(6,'(a35,1x,7(i8,1x))') 'crystallite_subFp0: ', shape(crystallite_subFp0) write(6,'(a35,1x,7(i8,1x))') 'crystallite_subFi0: ', shape(crystallite_subFi0) write(6,'(a35,1x,7(i8,1x))') 'crystallite_subLp0: ', shape(crystallite_subLp0) @@ -651,9 +648,6 @@ subroutine crystallite_stressAndItsTangent(updateJaco) crystallite_dPdF0(1:3,1:3,1:3,1:3,c,i,e) = crystallite_partioneddPdF0(1:3,1:3,1:3,1:3,c,i,e) ! ...stiffness crystallite_subF0(1:3,1:3,c,i,e) = crystallite_partionedF0(1:3,1:3,c,i,e) ! ...def grad crystallite_subTstar0_v(1:6,c,i,e) = crystallite_partionedTstar0_v(1:6,c,i,e) !...2nd PK stress - crystallite_subFe0(1:3,1:3,c,i,e) = math_mul33x33(math_mul33x33(crystallite_subF0(1:3,1:3,c,i,e), & - math_inv33(crystallite_subFp0(1:3,1:3,c,i,e))), & - math_inv33(crystallite_subFi0(1:3,1:3,c,i,e)))! only needed later on for stiffness calculation crystallite_subFrac(c,i,e) = 0.0_pReal crystallite_subStep(c,i,e) = 1.0_pReal/subStepSizeCryst crystallite_todo(c,i,e) = .true. @@ -915,9 +909,6 @@ subroutine crystallite_stressAndItsTangent(updateJaco) crystallite_subLi0(1:3,1:3,c,i,e) = crystallite_Li(1:3,1:3,c,i,e) ! ...intermediate velocity gradient crystallite_subFp0(1:3,1:3,c,i,e) = crystallite_Fp(1:3,1:3,c,i,e) ! ...plastic def grad crystallite_subFi0(1:3,1:3,c,i,e) = crystallite_Fi(1:3,1:3,c,i,e) ! ...intermediate def grad - crystallite_subFe0(1:3,1:3,c,i,e) = math_mul33x33(math_mul33x33(crystallite_subF (1:3,1:3,c,i,e), & - crystallite_invFp(1:3,1:3,c,i,e)), & - crystallite_invFi(1:3,1:3,c,i,e)) ! only needed later on for stiffness calculation !if abbrevation, make c and p private in omp plasticState (phaseAt(c,i,e))%subState0(:,phasememberAt(c,i,e)) = & plasticState (phaseAt(c,i,e))%state( :,phasememberAt(c,i,e)) From 0bf64645a1a4329a990e37695a732432192a78bd Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 19 Sep 2018 14:22:35 +0200 Subject: [PATCH 21/35] should be done by the plasticity laws (for the moment) --- src/crystallite.f90 | 13 ------------- 1 file changed, 13 deletions(-) diff --git a/src/crystallite.f90 b/src/crystallite.f90 index d341f3653..1c77037ad 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -435,19 +435,6 @@ subroutine crystallite_init enddo !$OMP END PARALLEL DO -!-------------------------------------------------------------------------------------------------- -! propagate dependent states to materialpoint and boundary value problem level -! do ph = 1_pInt,material_Nphase -! plasticState(ph)%partionedState0(plasticState(ph)%offsetDeltaState+plasticState(ph)%sizeDeltaState: & -! plasticState(ph)%sizeState,:) & -! = plasticState(ph)%state(plasticState(ph)%offsetDeltaState+plasticState(ph)%sizeDeltaState: & -! plasticState(ph)%sizeState,:) -! plasticState(ph)%state0 (plasticState(ph)%offsetDeltaState+plasticState(ph)%sizeDeltaState: & -! plasticState(ph)%sizeState,:) & -! = plasticState(ph)%state(plasticState(ph)%offsetDeltaState+plasticState(ph)%sizeDeltaState: & -! plasticState(ph)%sizeState,:) -! enddo - call crystallite_stressAndItsTangent(.true.) ! request elastic answers !-------------------------------------------------------------------------------------------------- From c313dc1675b235e137b6a342f9b3ae949cbe0751 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 19 Sep 2018 17:04:12 +0200 Subject: [PATCH 22/35] only read access --- src/crystallite.f90 | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/src/crystallite.f90 b/src/crystallite.f90 index 1c77037ad..3bad73b4d 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -38,6 +38,9 @@ module crystallite crystallite_orientation, & !< orientation as quaternion crystallite_orientation0, & !< initial orientation as quaternion crystallite_rotation !< grain rotation away from initial orientation as axis-angle (in degrees) in crystal reference frame + real(pReal), dimension(:,:,:,:,:), allocatable, public, protected :: & + crystallite_Fe, & !< current "elastic" def grad (end of converged time step) + crystallite_P !< 1st Piola-Kirchhoff stress per grain real(pReal), dimension(:,:,:,:,:), allocatable, public :: & crystallite_Fp, & !< current plastic def grad (end of converged time step) crystallite_Fp0, & !< plastic def grad at start of FE inc @@ -50,12 +53,10 @@ module crystallite crystallite_partionedF0, & !< def grad at start of homog inc crystallite_Lp, & !< current plastic velocitiy grad (end of converged time step) crystallite_Lp0, & !< plastic velocitiy grad at start of FE inc - crystallite_partionedLp0,& !< plastic velocity grad at start of homog inc + crystallite_partionedLp0, & !< plastic velocity grad at start of homog inc crystallite_Li, & !< current intermediate velocitiy grad (end of converged time step) crystallite_Li0, & !< intermediate velocitiy grad at start of FE inc - crystallite_partionedLi0,& !< intermediate velocity grad at start of homog inc - crystallite_Fe, & !< current "elastic" def grad (end of converged time step) - crystallite_P !< 1st Piola-Kirchhoff stress per grain + crystallite_partionedLi0 !< intermediate velocity grad at start of homog inc real(pReal), dimension(:,:,:,:,:), allocatable, private :: & crystallite_invFp, & !< inverse of current plastic def grad (end of converged time step) crystallite_subFp0,& !< plastic def grad at start of crystallite inc From b84476e681a7df20eb6a301801f6971844df61a3 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 19 Sep 2018 18:33:37 +0200 Subject: [PATCH 23/35] cleaning and debugging --- src/DAMASK_FEM.f90 | 22 +++------------------- src/FEM_mech.f90 | 7 ++----- src/FEM_utilities.f90 | 13 +++++-------- 3 files changed, 10 insertions(+), 32 deletions(-) diff --git a/src/DAMASK_FEM.f90 b/src/DAMASK_FEM.f90 index ee425585c..96204b621 100644 --- a/src/DAMASK_FEM.f90 +++ b/src/DAMASK_FEM.f90 @@ -42,9 +42,9 @@ program DAMASK_FEM restartWrite, & restartInc use numerics, only: & + worldrank, & maxCutBack, & - stagItMax, & - worldrank + stagItMax use mesh, only: & mesh_Nboundaries, & mesh_boundaries, & @@ -73,11 +73,7 @@ program DAMASK_FEM COMPONENT_SOLUTE_CVaH_ID, & COMPONENT_SOLUTE_CVaHPOT_ID, & COMPONENT_MGTWIN_PHI_ID, & - FIELD_MECH_label, & - FIELD_THERMAL_label, & - FIELD_DAMAGE_label, & - FIELD_SOLUTE_label, & - FIELD_MGTWIN_label + FIELD_MECH_label use FEM_mech implicit none @@ -405,18 +401,6 @@ program DAMASK_FEM case(FIELD_MECH_ID) write(6,'(2x,a)') 'Field '//trim(FIELD_MECH_label) - case(FIELD_THERMAL_ID) - write(6,'(2x,a)') 'Field '//trim(FIELD_THERMAL_label) - - case(FIELD_DAMAGE_ID) - write(6,'(2x,a)') 'Field '//trim(FIELD_DAMAGE_label) - - case(FIELD_MGTWIN_ID) - write(6,'(2x,a)') 'Field '//trim(FIELD_MGTWIN_label) - - case(FIELD_SOLUTE_ID) - write(6,'(2x,a)') 'Field '//trim(FIELD_SOLUTE_label) - end select do faceSet = 1_pInt, mesh_Nboundaries do component = 1_pInt, loadCases(currentLoadCase)%fieldBC(field)%nComponents diff --git a/src/FEM_mech.f90 b/src/FEM_mech.f90 index 0fdaf6283..09696b9ab 100644 --- a/src/FEM_mech.f90 +++ b/src/FEM_mech.f90 @@ -23,9 +23,6 @@ module FEM_mech tSolutionState, & tFieldBC, & tComponentBC - use numerics, only: & - worldrank, & - worldsize use mesh, only: & mesh_Nboundaries, & mesh_boundaries @@ -64,7 +61,7 @@ module FEM_mech FEM_mech_solution ,& FEM_mech_forward, & FEM_mech_destroy - + external :: PETScerrorf contains !-------------------------------------------------------------------------------------------------- @@ -80,7 +77,6 @@ subroutine FEM_mech_init(fieldBC) use mesh, only: & geomMesh use numerics, only: & - worldrank, & itmax, & integrationOrder use FEM_Zoo, only: & @@ -672,6 +668,7 @@ end subroutine FEM_mech_forward !-------------------------------------------------------------------------------------------------- subroutine FEM_mech_converged(snes_local,PETScIter,xnorm,snorm,fnorm,reason,dummy,ierr) use numerics, only: & + worldrank, & err_struct_tolAbs, & err_struct_tolRel use IO, only: & diff --git a/src/FEM_utilities.f90 b/src/FEM_utilities.f90 index 9feecbefd..117ffcafa 100644 --- a/src/FEM_utilities.f90 +++ b/src/FEM_utilities.f90 @@ -24,19 +24,16 @@ use PETScis real(pReal), public :: wgt !< weighting factor 1/Nelems real(pReal), public :: wgtDof !< weighting factor 1/Nelems real(pReal), public :: C_volAvg(3,3,3,3) - + !-------------------------------------------------------------------------------------------------- ! output data Vec, public :: coordinatesVec - !-------------------------------------------------------------------------------------------------- ! field labels information character(len=*), parameter, public :: & - FIELD_MECH_label = 'mechanical', & - FIELD_THERMAL_label = 'thermal', & - FIELD_DAMAGE_label = 'damage', & - FIELD_SOLUTE_label = 'solute', & - FIELD_MGTWIN_label = 'mgtwin' + FIELD_MECH_label = 'mechanical' + + integer(pInt), parameter :: structOrder = 2_pInt enum, bind(c) enumerator :: FIELD_UNDEFINED_ID, & @@ -207,7 +204,7 @@ subroutine utilities_init() call PetscOptionsInsertString(PETSC_NULL_OPTIONS,trim(petsc_defaultOptions),ierr) call PetscOptionsInsertString(PETSC_NULL_OPTIONS,trim(petsc_options),ierr) CHKERRQ(ierr) - !write(petsc_optionsPhysics,'(a,i0)') '-mechFE_petscspace_order ' , structOrder + write(petsc_optionsPhysics,'(a,i0)') '-mechFE_petscspace_order ' , structOrder call PetscOptionsInsertString(PETSC_NULL_OPTIONS,trim(petsc_optionsPhysics),ierr) CHKERRQ(ierr) From 11d4c28d885150d0f054f11a631e7dbabd23fc94 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 19 Sep 2018 19:45:57 +0200 Subject: [PATCH 24/35] flushes not needed + further cleaning --- src/crystallite.f90 | 57 +++++++++++---------------------------------- 1 file changed, 13 insertions(+), 44 deletions(-) diff --git a/src/crystallite.f90 b/src/crystallite.f90 index 3bad73b4d..f28c7cb00 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -422,11 +422,10 @@ subroutine crystallite_init call crystallite_orientations() crystallite_orientation0 = crystallite_orientation ! store initial orientations for calculation of grain rotations - !$OMP PARALLEL DO PRIVATE(myNcomponents) + !$OMP PARALLEL DO do e = FEsolving_execElem(1),FEsolving_execElem(2) - myNcomponents = homogenization_Ngrains(mesh_element(3,e)) do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) - do c = 1_pInt,myNcomponents + do c = 1_pInt,homogenization_Ngrains(mesh_element(3,e)) call constitutive_microstructure(crystallite_orientation, & ! pass orientation to constitutive module crystallite_Fe(1:3,1:3,c,i,e), & crystallite_Fp(1:3,1:3,c,i,e), & @@ -574,7 +573,6 @@ subroutine crystallite_stressAndItsTangent(updateJaco) neighboring_i, & o, & p, & - myNcomponents, & mySource ! local variables used for calculating analytic Jacobian real(pReal), dimension(3,3) :: temp_33 @@ -617,10 +615,9 @@ subroutine crystallite_stressAndItsTangent(updateJaco) ! initialize to starting condition crystallite_subStep = 0.0_pReal - !$OMP PARALLEL DO PRIVATE(myNcomponents) + !$OMP PARALLEL DO elementLooping1: do e = FEsolving_execElem(1),FEsolving_execElem(2) - myNcomponents = homogenization_Ngrains(mesh_element(3,e)) - do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e); do c = 1_pInt,myNcomponents + do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e); do c = 1_pInt,homogenization_Ngrains(mesh_element(3,e)) if (crystallite_requested(c,i,e)) then plasticState (phaseAt(c,i,e))%subState0( :,phasememberAt(c,i,e)) = & plasticState (phaseAt(c,i,e))%partionedState0(:,phasememberAt(c,i,e)) @@ -876,23 +873,20 @@ subroutine crystallite_stressAndItsTangent(updateJaco) endif timeSyncing1 - !$OMP PARALLEL DO PRIVATE(myNcomponents,formerSubStep) + !$OMP PARALLEL DO PRIVATE(formerSubStep) elementLooping3: do e = FEsolving_execElem(1),FEsolving_execElem(2) - 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,myNcomponents + do c = 1,homogenization_Ngrains(mesh_element(3,e)) ! --- wind forward --- if (crystallite_converged(c,i,e) .and. crystallite_clearToWindForward(i,e)) then formerSubStep = crystallite_subStep(c,i,e) crystallite_subFrac(c,i,e) = crystallite_subFrac(c,i,e) + crystallite_subStep(c,i,e) - !$OMP FLUSH(crystallite_subFrac) crystallite_subStep(c,i,e) = min(1.0_pReal - crystallite_subFrac(c,i,e), & stepIncreaseCryst * crystallite_subStep(c,i,e)) - !$OMP FLUSH(crystallite_subStep) + if (crystallite_subStep(c,i,e) > 0.0_pReal) then crystallite_subF0(1:3,1:3,c,i,e) = crystallite_subF(1:3,1:3,c,i,e) ! ...def grad - !$OMP FLUSH(crystallite_subF0) crystallite_subLp0(1:3,1:3,c,i,e) = crystallite_Lp(1:3,1:3,c,i,e) ! ...plastic velocity gradient crystallite_subLi0(1:3,1:3,c,i,e) = crystallite_Li(1:3,1:3,c,i,e) ! ...intermediate velocity gradient crystallite_subFp0(1:3,1:3,c,i,e) = crystallite_Fp(1:3,1:3,c,i,e) ! ...plastic def grad @@ -912,7 +906,6 @@ subroutine crystallite_stressAndItsTangent(updateJaco) else crystallite_todo(c,i,e) = .true. endif - !$OMP FLUSH(crystallite_todo) #ifdef DEBUG if (iand(debug_level(debug_crystallite),debug_levelBasic) /= 0_pInt & .and. ((e == debug_e .and. i == debug_i .and. c == debug_g) & @@ -923,7 +916,6 @@ subroutine crystallite_stressAndItsTangent(updateJaco) #endif else ! this crystallite just converged for the entire timestep crystallite_todo(c,i,e) = .false. ! so done here - !$OMP FLUSH(crystallite_todo) endif ! --- cutback --- @@ -934,15 +926,10 @@ subroutine crystallite_stressAndItsTangent(updateJaco) else crystallite_subStep(c,i,e) = subStepSizeCryst * crystallite_subStep(c,i,e) ! cut step in half and restore... endif - !$OMP FLUSH(crystallite_subStep) crystallite_Fp(1:3,1:3,c,i,e) = crystallite_subFp0(1:3,1:3,c,i,e) ! ...plastic def grad - !$OMP FLUSH(crystallite_Fp) crystallite_invFp(1:3,1:3,c,i,e) = math_inv33(crystallite_Fp(1:3,1:3,c,i,e)) - !$OMP FLUSH(crystallite_invFp) crystallite_Fi(1:3,1:3,c,i,e) = crystallite_subFi0(1:3,1:3,c,i,e) ! ...intermediate def grad - !$OMP FLUSH(crystallite_Fi) crystallite_invFi(1:3,1:3,c,i,e) = math_inv33(crystallite_Fi(1:3,1:3,c,i,e)) - !$OMP FLUSH(crystallite_invFi) crystallite_Lp(1:3,1:3,c,i,e) = crystallite_subLp0(1:3,1:3,c,i,e) ! ...plastic velocity grad crystallite_Li(1:3,1:3,c,i,e) = crystallite_subLi0(1:3,1:3,c,i,e) ! ...intermediate velocity grad plasticState (phaseAt(c,i,e))%state( :,phasememberAt(c,i,e)) = & @@ -955,7 +942,6 @@ subroutine crystallite_stressAndItsTangent(updateJaco) ! cant restore dotState here, since not yet calculated in first cutback after initialization crystallite_todo(c,i,e) = crystallite_subStep(c,i,e) > subStepMinCryst ! still on track or already done (beyond repair) - !$OMP FLUSH(crystallite_todo) #ifdef DEBUG if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt & .and. ((e == debug_e .and. i == debug_i .and. c == debug_g) & @@ -976,10 +962,8 @@ subroutine crystallite_stressAndItsTangent(updateJaco) if (crystallite_todo(c,i,e) .and. (crystallite_clearToWindForward(i,e) .or. crystallite_clearToCutback(i,e))) then crystallite_subF(1:3,1:3,c,i,e) = crystallite_subF0(1:3,1:3,c,i,e) & - + crystallite_subStep(c,i,e) & - * (crystallite_partionedF(1:3,1:3,c,i,e) & + + crystallite_subStep(c,i,e) * (crystallite_partionedF(1:3,1:3,c,i,e) & - crystallite_partionedF0(1:3,1:3,c,i,e)) - !$OMP FLUSH(crystallite_subF) crystallite_Fe(1:3,1:3,c,i,e) = math_mul33x33(math_mul33x33(crystallite_subF (1:3,1:3,c,i,e), & crystallite_invFp(1:3,1:3,c,i,e)), & crystallite_invFi(1:3,1:3,c,i,e)) @@ -997,9 +981,8 @@ subroutine crystallite_stressAndItsTangent(updateJaco) .and. crystallite_subStep <= subStepMinCryst)) then ! no way of rescuing a nonlocal ip that violated the lower time step limit, ... if (iand(debug_level(debug_crystallite),debug_levelExtensive) /= 0_pInt) then elementLooping4: do e = FEsolving_execElem(1),FEsolving_execElem(2) - myNcomponents = homogenization_Ngrains(mesh_element(3,e)) do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) - do c = 1,myNcomponents + do c = 1,homogenization_Ngrains(mesh_element(3,e)) if (.not. crystallite_localPlasticity(c,i,e) .and. .not. crystallite_todo(c,i,e) & .and. .not. crystallite_converged(c,i,e) .and. crystallite_subStep(c,i,e) <= subStepMinCryst) & write(6,'(a,i8,1x,i2,1x,i3)') '<< CRYST >> nonlocal violated minimum subStep at el ip ipc ',e,i,c @@ -1041,9 +1024,8 @@ subroutine crystallite_stressAndItsTangent(updateJaco) ! --+>> CHECK FOR NON-CONVERGED CRYSTALLITES <<+-- elementLooping5: do e = FEsolving_execElem(1),FEsolving_execElem(2) - 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,myNcomponents + do c = 1,homogenization_Ngrains(mesh_element(3,e)) if (.not. crystallite_converged(c,i,e)) then ! respond fully elastically (might be not required due to becoming terminally ill anyway) if(iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) & write(6,'(a,i8,1x,a,i8,a,1x,i2,1x,i3,/)') '<< CRYST >> no convergence: respond fully elastic at el (elFE) ip ipc ', & @@ -1080,11 +1062,10 @@ subroutine crystallite_stressAndItsTangent(updateJaco) computeJacobian: if(updateJaco) then !$OMP PARALLEL DO PRIVATE(dSdF,dSdFe,dSdFi,dLpdS,dLpdFi,dFpinvdF,dLidS,dLidFi,dFidS,& - !$OMP rhs_3333,lhs_3333,temp_99,temp_33,temp_3333,myNcomponents,error) + !$OMP rhs_3333,lhs_3333,temp_99,temp_33,temp_3333,error) elementLooping6: do e = FEsolving_execElem(1),FEsolving_execElem(2) - 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 + do c = 1_pInt,homogenization_Ngrains(mesh_element(3,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 @@ -1189,7 +1170,6 @@ subroutine crystallite_stressAndItsTangent(updateJaco) enddo elementLooping6 !$OMP END PARALLEL DO endif computeJacobian -!why not OMP? end subroutine crystallite_stressAndItsTangent @@ -2248,8 +2228,6 @@ subroutine crystallite_integrateStateAdaptiveEuler() + 0.5_pReal * sourceState(p)%p(mySource)%dotState(:,c) & * crystallite_subdt(g,i,e) ! contribution to absolute residuum in state enddo - !$OMP FLUSH(plasticStateResiduum) - !$OMP FLUSH(sourceStateResiduum) ! --- relative residui --- forall (s = 1_pInt:mySizePlasticDotState, abs(plasticState(p)%dotState(s,c)) > 0.0_pReal) & @@ -2261,11 +2239,8 @@ subroutine crystallite_integrateStateAdaptiveEuler() relSourceStateResiduum(s,mySource,g,i,e) = & sourceStateResiduum(s,mySource,g,i,e) / sourceState(p)%p(mySource)%dotState(s,c) enddo - !$OMP FLUSH(relPlasticStateResiduum) - !$OMP FLUSH(relSourceStateResiduum) #ifdef DEBUG - if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt & .and. ((e == debug_e .and. i == debug_i .and. g == debug_g)& .or. .not. iand(debug_level(debug_crystallite), debug_levelSelective) /= 0_pInt)) then @@ -2293,13 +2268,7 @@ subroutine crystallite_integrateStateAdaptiveEuler() abs(sourceStateResiduum(1:mySizeSourceDotState,mySource,g,i,e)) < & sourceState(p)%p(mySource)%aTolState(1:mySizeSourceDotState)) enddo - if (converged) then - crystallite_converged(g,i,e) = .true. ! ... converged per definitionem - if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) then - !$OMP CRITICAL (distributionState) - !$OMP END CRITICAL (distributionState) - endif - endif + if (converged) crystallite_converged(g,i,e) = .true. ! ... converged per definitionem endif enddo; enddo; enddo !$OMP ENDDO From e7e959af471ef75ee89d4f2dc409f37c1d80b282 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 19 Sep 2018 20:06:35 +0200 Subject: [PATCH 25/35] using explicit interface from PETScDT --- src/FEM_mech.f90 | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/FEM_mech.f90 b/src/FEM_mech.f90 index 09696b9ab..d77874515 100644 --- a/src/FEM_mech.f90 +++ b/src/FEM_mech.f90 @@ -14,6 +14,7 @@ module FEM_mech use PETScsnes use PETScDM use PETScDMplex + use PETScDT use prec, only: & pInt, & pReal @@ -100,8 +101,8 @@ subroutine FEM_mech_init(fieldBC) IS, pointer :: pBcComps(:), pBcPoints(:) PetscSection :: section PetscInt :: field, faceSet, topologDim, nNodalPoints - PetscReal, pointer :: qPointsP(:), qWeightsP(:), & - nodalPointsP(:), nodalWeightsP(:) + PetscReal, dimension(:) , pointer :: qPointsP, qWeightsP, & + nodalPointsP, nodalWeightsP PetscReal, allocatable, target :: nodalPoints(:), nodalWeights(:) PetscScalar, pointer :: px_scal(:) PetscScalar, allocatable, target :: x_scal(:) @@ -111,7 +112,8 @@ subroutine FEM_mech_init(fieldBC) PetscInt :: cellStart, cellEnd, cell, basis character(len=7) :: prefix = 'mechFE_' PetscErrorCode :: ierr - + PetscReal, allocatable, target, dimension(:) :: qWeights + write(6,'(/,a)') ' <<<+- FEM_mech init -+>>>' write(6,'(a15,a)') ' Current time: ',IO_timeStamp() #include "compilation_info.f90" @@ -123,8 +125,6 @@ subroutine FEM_mech_init(fieldBC) !-------------------------------------------------------------------------------------------------- ! Setup FEM mech discretization - allocate(qPoints(dimPlex*FEM_Zoo_nQuadrature(dimPlex,integrationOrder))) - allocate(qWeights(FEM_Zoo_nQuadrature(dimPlex,integrationOrder))) qPoints = FEM_Zoo_QuadraturePoints(dimPlex,integrationOrder)%p qWeights = FEM_Zoo_QuadratureWeights(dimPlex,integrationOrder)%p nQuadrature = FEM_Zoo_nQuadrature(dimPlex,integrationOrder) From 1bf1e13b46d143629fc0f05ddc66b89cbb447bbe Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 19 Sep 2018 20:23:28 +0200 Subject: [PATCH 26/35] some debug statements --- src/FEM_mech.f90 | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/src/FEM_mech.f90 b/src/FEM_mech.f90 index d77874515..d8d6102df 100644 --- a/src/FEM_mech.f90 +++ b/src/FEM_mech.f90 @@ -241,12 +241,17 @@ subroutine FEM_mech_init(fieldBC) call PetscFEGetDualSpace(mechFE,mechDualSpace,ierr); CHKERRQ(ierr) call DMPlexGetHeightStratum(mech_mesh,0,cellStart,cellEnd,ierr) CHKERRQ(ierr) + write(6,*) 'cellDof', cellDof;flush(6) + write(6,*) 'cell start and end-1',cellStart,cellEnd-1;flush(6) do cell = cellStart, cellEnd-1 !< loop over all elements + write(6,*) 'cell',cell;flush(6) x_scal = 0.0 call DMPlexComputeCellGeometryAffineFEM(mech_mesh,cell,pV0,pCellJ,pInvcellJ,detJ,ierr) CHKERRQ(ierr) cellJMat = reshape(pCellJ,shape=[dimPlex,dimPlex]) do basis = 0, nBasis-1 + write(6,*) 'nBasis-1',nBasis-1;flush(6) + write(6,*) 'basis',basis;flush(6) call PetscDualSpaceGetFunctional(mechDualSpace,basis,functional,ierr) CHKERRQ(ierr) call PetscQuadratureGetData(functional,dimPlex,nc,nNodalPoints,nodalPointsP,nodalWeightsP,ierr) From 191ad9df09c790c847cc8abb493efec4ac4d4842 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 19 Sep 2018 20:39:49 +0200 Subject: [PATCH 27/35] not used --- src/FEM_utilities.f90 | 31 +++++-------------------------- 1 file changed, 5 insertions(+), 26 deletions(-) diff --git a/src/FEM_utilities.f90 b/src/FEM_utilities.f90 index 117ffcafa..6242a4ebe 100644 --- a/src/FEM_utilities.f90 +++ b/src/FEM_utilities.f90 @@ -23,7 +23,6 @@ use PETScis ! grid related information information real(pReal), public :: wgt !< weighting factor 1/Nelems real(pReal), public :: wgtDof !< weighting factor 1/Nelems - real(pReal), public :: C_volAvg(3,3,3,3) !-------------------------------------------------------------------------------------------------- ! output data @@ -118,23 +117,12 @@ use PETScis utilities_indexActiveSet, & utilities_destroy, & FIELD_MECH_ID, & - FIELD_THERMAL_ID, & - FIELD_DAMAGE_ID, & - FIELD_SOLUTE_ID, & - FIELD_MGTWIN_ID, & COMPONENT_MECH_X_ID, & COMPONENT_MECH_Y_ID, & COMPONENT_MECH_Z_ID, & - COMPONENT_THERMAL_T_ID, & - COMPONENT_DAMAGE_PHI_ID, & - COMPONENT_SOLUTE_CV_ID, & - COMPONENT_SOLUTE_CVPOT_ID, & - COMPONENT_SOLUTE_CH_ID, & - COMPONENT_SOLUTE_CHPOT_ID, & - COMPONENT_SOLUTE_CVaH_ID, & - COMPONENT_SOLUTE_CVaHPOT_ID, & - COMPONENT_MGTWIN_PHI_ID + COMPONENT_THERMAL_T_ID + external :: PETScErrorF contains !-------------------------------------------------------------------------------------------------- @@ -173,14 +161,11 @@ subroutine utilities_init() implicit none - character(len=1024) :: petsc_optionsPhysics, grainStr + character(len=1024) :: petsc_optionsPhysics integer(pInt) :: dimPlex integer(pInt) :: headerID = 205_pInt - PetscInt, dimension(:), pointer :: points - PetscInt, allocatable :: nEntities(:), nOutputCells(:), nOutputNodes(:), mappingCells(:) - PetscInt :: cellStart, cellEnd, cell, ip, dim, ctr, qPt - PetscInt, allocatable :: connectivity(:,:) - Vec :: connectivityVec + PetscInt, allocatable :: nEntities(:), nOutputCells(:), nOutputNodes(:) + PetscInt :: dim PetscErrorCode :: ierr write(6,'(/,a)') ' <<<+- DAMASK_FEM_utilities init -+>>>' @@ -246,8 +231,6 @@ subroutine utilities_constitutiveResponse(timeinc,P_av,forwardData) use debug, only: & debug_reset, & debug_info - use numerics, only: & - worldrank use math, only: & math_transpose33, & math_rotate_forward33, & @@ -255,10 +238,8 @@ subroutine utilities_constitutiveResponse(timeinc,P_av,forwardData) use FEsolving, only: & restartWrite use homogenization, only: & - materialpoint_F0, & materialpoint_F, & materialpoint_P, & - materialpoint_dPdF, & materialpoint_stressAndItsTangent use mesh, only: & mesh_NcpElems @@ -311,9 +292,7 @@ subroutine utilities_constitutiveResponse(timeinc,P_av,forwardData) cutBack = .false. ! reset cutBack status P_av = sum(sum(materialpoint_P,dim=4),dim=3) * wgt ! average of P - C_volAvg = sum(sum(materialpoint_dPdF,dim=6),dim=5) * wgt call MPI_Allreduce(MPI_IN_PLACE,P_av,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) - call MPI_Allreduce(MPI_IN_PLACE,C_volAvg,81,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD, ierr) end subroutine utilities_constitutiveResponse From 4b4cd4e6c31587478ee731c53979ca34d9105aee Mon Sep 17 00:00:00 2001 From: Test User Date: Wed, 19 Sep 2018 21:08:45 +0200 Subject: [PATCH 28/35] [skip ci] updated version information after successful test of v2.0.2-565-g59043d58 --- VERSION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/VERSION b/VERSION index cca202405..d3cb228fe 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -v2.0.2-560-g20f0bee4 +v2.0.2-565-g59043d58 From 1623a33b486914b9c7ae2f2f3df4e8c99ebe694e Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 19 Sep 2018 21:45:12 +0200 Subject: [PATCH 29/35] cleaning (mainly OMP FLUSh) --- src/crystallite.f90 | 3 --- src/homogenization.f90 | 17 ++--------------- 2 files changed, 2 insertions(+), 18 deletions(-) diff --git a/src/crystallite.f90 b/src/crystallite.f90 index f28c7cb00..ee15cb67b 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -1855,9 +1855,6 @@ subroutine crystallite_integrateStateRKCK45() relSourceStateResiduum(s,mySource,g,i,e) = & sourceStateResiduum(s,mySource,g,i,e) / sourceState(p)%p(mySource)%state(s,cc) enddo - !$OMP FLUSH(relPlasticStateResiduum) - !$OMP FLUSH(relSourceStateResiduum) -! @Martin: do we need flushing? why..? crystallite_todo(g,i,e) = all(abs(relPlasticStateResiduum(1:mySizePlasticDotState,g,i,e)) < & rTol_crystalliteState .or. & abs(plasticStateResiduum(1:mySizePlasticDotState,g,i,e)) < & diff --git a/src/homogenization.f90 b/src/homogenization.f90 index 77d301400..a2b5d12e7 100644 --- a/src/homogenization.f90 +++ b/src/homogenization.f90 @@ -448,8 +448,6 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt) subStepSizeHomog, & stepIncreaseHomog, & nMPstate - use math, only: & - math_transpose33 use FEsolving, only: & FEsolving_execElem, & FEsolving_execIP, & @@ -524,9 +522,9 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt) write(6,'(/a,i5,1x,i2)') '<< HOMOG >> Material Point start at el ip ', debug_e, debug_i write(6,'(a,/,3(12x,3(f14.9,1x)/))') '<< HOMOG >> F0', & - math_transpose33(materialpoint_F0(1:3,1:3,debug_i,debug_e)) + transpose(materialpoint_F0(1:3,1:3,debug_i,debug_e)) write(6,'(a,/,3(12x,3(f14.9,1x)/))') '<< HOMOG >> F', & - math_transpose33(materialpoint_F(1:3,1:3,debug_i,debug_e)) + transpose(materialpoint_F(1:3,1:3,debug_i,debug_e)) !$OMP END CRITICAL (write2out) endif @@ -608,10 +606,8 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt) !--------------------------------------------------------------------------------------------------- ! calculate new subStep and new subFrac materialpoint_subFrac(i,e) = materialpoint_subFrac(i,e) + materialpoint_subStep(i,e) - !$OMP FLUSH(materialpoint_subFrac) materialpoint_subStep(i,e) = min(1.0_pReal-materialpoint_subFrac(i,e), & stepIncreaseHomog*materialpoint_subStep(i,e)) ! introduce flexibility for step increase/acceleration - !$OMP FLUSH(materialpoint_subStep) steppingNeeded: if (materialpoint_subStep(i,e) > subStepMinHomog) then @@ -671,7 +667,6 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt) hydrogenfluxState(mappingHomogenization(2,i,e))%subState0(:,mappingHomogenization(1,i,e)) = & hydrogenfluxState(mappingHomogenization(2,i,e))%State( :,mappingHomogenization(1,i,e))! ...internal hydrogen transport state materialpoint_subF0(1:3,1:3,i,e) = materialpoint_subF(1:3,1:3,i,e) ! ...def grad - !$OMP FLUSH(materialpoint_subF0) endif steppingNeeded else converged @@ -689,7 +684,6 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt) !$OMP END CRITICAL (setTerminallyIll) else ! cutback makes sense materialpoint_subStep(i,e) = subStepSizeHomog * materialpoint_subStep(i,e) ! crystallite had severe trouble, so do a significant cutback - !$OMP FLUSH(materialpoint_subStep) #ifdef DEBUG if (iand(debug_level(debug_homogenization), debug_levelExtensive) /= 0_pInt & @@ -810,13 +804,6 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt) materialpoint_doneAndHappy(1:2,i,e) = homogenization_updateState(i,e) materialpoint_converged(i,e) = all(materialpoint_doneAndHappy(1:2,i,e)) ! converged if done and happy endif - !$OMP FLUSH(materialpoint_converged) - if (materialpoint_converged(i,e)) then - if (iand(debug_level(debug_homogenization), debug_levelBasic) /= 0_pInt) then - !$OMP CRITICAL (distributionMPState) - !$OMP END CRITICAL (distributionMPState) - endif - endif endif enddo IpLooping3 enddo elementLooping3 From 6aa4dd842a1e5cfaadf6e455ffb2ac4482c76cf7 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 20 Sep 2018 06:09:02 +0200 Subject: [PATCH 30/35] define debug variables only if needed --- src/crystallite.f90 | 26 ++++++++++++++++++++++---- 1 file changed, 22 insertions(+), 4 deletions(-) diff --git a/src/crystallite.f90 b/src/crystallite.f90 index f28c7cb00..e1a444b2a 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -142,12 +142,14 @@ subroutine crystallite_init compiler_version, & compiler_options #endif +#ifdef DEBUG use debug, only: & debug_info, & debug_reset, & debug_level, & debug_crystallite, & debug_levelBasic +#endif use numerics, only: & numerics_integrator, & worldrank, & @@ -437,8 +439,7 @@ subroutine crystallite_init call crystallite_stressAndItsTangent(.true.) ! request elastic answers -!-------------------------------------------------------------------------------------------------- -! debug output +#ifdef DEBUG if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) then write(6,'(a35,1x,7(i8,1x))') 'crystallite_Fe: ', shape(crystallite_Fe) write(6,'(a35,1x,7(i8,1x))') 'crystallite_Fp: ', shape(crystallite_Fp) @@ -490,6 +491,7 @@ subroutine crystallite_init call debug_info call debug_reset +#endif end subroutine crystallite_init @@ -591,7 +593,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco) real(pReal), dimension(9,9):: temp_99 logical :: error - +#ifdef DEBUG if (iand(debug_level(debug_crystallite),debug_levelSelective) /= 0_pInt & .and. FEsolving_execElem(1) <= debug_e & .and. debug_e <= FEsolving_execElem(2)) then @@ -610,6 +612,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco) write(6,'(a,/,3(12x,3(f14.9,1x)/))') '<< CRYST >> Li0', & transpose(crystallite_partionedLi0(1:3,1:3,debug_g,debug_i,debug_e)) endif +#endif !-------------------------------------------------------------------------------------------------- ! initialize to starting condition @@ -741,8 +744,10 @@ subroutine crystallite_stressAndItsTangent(updateJaco) if (all(crystallite_localPlasticity .or. crystallite_converged)) then if (all(crystallite_localPlasticity .or. crystallite_subStep + crystallite_subFrac >= 1.0_pReal)) then crystallite_clearToWindForward = .true. ! final wind forward +#ifdef DEBUG if (iand(debug_level(debug_crystallite),debug_levelExtensive) /= 0_pInt) & write(6,'(a,i6)') '<< CRYST >> final wind forward' +#endif else !$OMP PARALLEL DO do e = FEsolving_execElem(1),FEsolving_execElem(2) @@ -751,8 +756,10 @@ subroutine crystallite_stressAndItsTangent(updateJaco) enddo enddo !$OMP END PARALLEL DO +#ifdef DEBUG if (iand(debug_level(debug_crystallite),debug_levelExtensive) /= 0_pInt) & write(6,'(a,i6)') '<< CRYST >> wind forward' +#endif endif else subFracIntermediate = maxval(crystallite_subFrac, mask=.not.crystallite_localPlasticity) @@ -838,8 +845,10 @@ subroutine crystallite_stressAndItsTangent(updateJaco) enddo enddo !$OMP END PARALLEL DO +#ifdef DEBUG if (iand(debug_level(debug_crystallite),debug_levelExtensive) /= 0_pInt) & write(6,'(a,i6)') '<< CRYST >> time synchronization: cutback' +#endif else !$OMP PARALLEL DO do e = FEsolving_execElem(1),FEsolving_execElem(2) @@ -848,8 +857,10 @@ subroutine crystallite_stressAndItsTangent(updateJaco) enddo enddo !$OMP END PARALLEL DO +#ifdef DEBUG if (iand(debug_level(debug_crystallite),debug_levelExtensive) /= 0_pInt) & write(6,'(a,i6)') '<< CRYST >> cutback' +#endif endif endif endif @@ -979,6 +990,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco) timeSyncing2: if(numerics_timeSyncing) then if (any(.not. crystallite_localPlasticity .and. .not. crystallite_todo .and. .not. crystallite_converged & .and. crystallite_subStep <= subStepMinCryst)) then ! no way of rescuing a nonlocal ip that violated the lower time step limit, ... +#ifdef DEBUG if (iand(debug_level(debug_crystallite),debug_levelExtensive) /= 0_pInt) then elementLooping4: do e = FEsolving_execElem(1),FEsolving_execElem(2) do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) @@ -990,6 +1002,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco) enddo enddo elementLooping4 endif +#endif where(.not. crystallite_localPlasticity) crystallite_todo = .false. ! ... so let all nonlocal ips die peacefully crystallite_subStep = 0.0_pReal @@ -997,6 +1010,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco) endif endif timeSyncing2 +#ifdef DEBUG if (iand(debug_level(debug_crystallite),debug_levelExtensive) /= 0_pInt) then write(6,'(/,a,f8.5)') '<< CRYST >> min(subStep) ',minval(crystallite_subStep) write(6,'(a,f8.5)') '<< CRYST >> max(subStep) ',maxval(crystallite_subStep) @@ -1009,9 +1023,9 @@ subroutine crystallite_stressAndItsTangent(updateJaco) flush(6) endif endif +#endif ! --- integrate --- requires fully defined state array (basic + dependent state) - if (any(crystallite_todo)) call integrateState() where(.not. crystallite_converged .and. crystallite_subStep > subStepMinCryst) & ! do not try non-converged & fully cutbacked any further crystallite_todo = .true. @@ -1027,9 +1041,11 @@ subroutine crystallite_stressAndItsTangent(updateJaco) do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed do c = 1,homogenization_Ngrains(mesh_element(3,e)) if (.not. crystallite_converged(c,i,e)) then ! respond fully elastically (might be not required due to becoming terminally ill anyway) +#ifdef DEBUG if(iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) & write(6,'(a,i8,1x,a,i8,a,1x,i2,1x,i3,/)') '<< CRYST >> no convergence: respond fully elastic at el (elFE) ip ipc ', & e,'(',mesh_element(1,e),')',i,c +#endif 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))) @@ -1037,6 +1053,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco) 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 +#ifdef DEBUG if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt & .and. ((e == debug_e .and. i == debug_i .and. c == debug_g) & .or. .not. iand(debug_level(debug_crystallite),debug_levelSelective) /= 0_pInt)) then @@ -1053,6 +1070,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco) transpose(crystallite_Li(1:3,1:3,c,i,e)) flush(6) endif +#endif enddo enddo enddo elementLooping5 From df0464c31bce615da018aca230ee4d702b8abd6f Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 20 Sep 2018 06:24:03 +0200 Subject: [PATCH 31/35] use (import) debug variables only when needed --- src/crystallite.f90 | 71 ++++++++++++++++++++++++++++----------------- 1 file changed, 44 insertions(+), 27 deletions(-) diff --git a/src/crystallite.f90 b/src/crystallite.f90 index e1a444b2a..a7463b9af 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -177,8 +177,7 @@ subroutine crystallite_init use config, only: & config_deallocate, & config_crystallite, & - crystallite_name, & - material_Nphase + crystallite_name use constitutive, only: & constitutive_initialFi, & constitutive_microstructure ! derived (shortcut) quantities of given state @@ -192,7 +191,6 @@ subroutine crystallite_init e, & !< counter in element loop o = 0_pInt, & !< counter in output loop r, & - ph, & !< counter in crystallite loop cMax, & !< maximum number of integration point components iMax, & !< maximum number of integration points eMax, & !< maximum number of elements @@ -508,6 +506,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco) subStepSizeCryst, & stepIncreaseCryst, & numerics_timeSyncing +#ifdef DEBUG use debug, only: & debug_level, & debug_crystallite, & @@ -517,6 +516,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco) debug_e, & debug_i, & debug_g +#endif use IO, only: & IO_warning, & IO_error @@ -656,10 +656,10 @@ subroutine crystallite_stressAndItsTangent(updateJaco) NiterationCrystallite = 0_pInt cutbackLooping: do while (any(crystallite_todo(:,startIP:endIP,FEsolving_execELem(1):FEsolving_execElem(2)))) - +#ifdef DEBUG if (iand(debug_level(debug_crystallite),debug_levelExtensive) /= 0_pInt) & write(6,'(a,i6)') '<< CRYST >> crystallite iteration ',NiterationCrystallite - +#endif timeSyncing1: if (any(.not. crystallite_localPlasticity) .and. numerics_timeSyncing) then ! Time synchronization can only be used for nonlocal calculations, and only there it makes sense. @@ -678,6 +678,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco) ! and its not clear how to fix this, so all nonlocals become terminally ill. if (any(crystallite_syncSubFrac .and. .not. crystallite_converged(1,:,:))) then +#ifdef DEBUG if (iand(debug_level(debug_crystallite),debug_levelExtensive) /= 0_pInt) then do e = FEsolving_execElem(1),FEsolving_execElem(2) do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) @@ -686,6 +687,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco) enddo enddo endif +#endif crystallite_syncSubFrac = .false. where(.not. crystallite_localPlasticity) crystallite_substep = 0.0_pReal @@ -700,8 +702,10 @@ subroutine crystallite_stressAndItsTangent(updateJaco) enddo enddo !$OMP END PARALLEL DO +#ifdef DEBUG if (iand(debug_level(debug_crystallite),debug_levelExtensive) /= 0_pInt) & write(6,'(a,i6)') '<< CRYST >> time synchronization: wind forward' +#endif endif elseif (any(crystallite_syncSubFracCompleted)) then @@ -717,8 +721,10 @@ subroutine crystallite_stressAndItsTangent(updateJaco) crystallite_clearToCutback(i,e) = crystallite_localPlasticity(1,i,e) .or. .not. crystallite_converged(1,i,e) enddo enddo +#ifdef DEBUG if (iand(debug_level(debug_crystallite),debug_levelExtensive) /= 0_pInt) & write(6,'(a,i6)') '<< CRYST >> time synchronization: done, proceed with cutback' +#endif else ! Normal calculation. @@ -1198,17 +1204,17 @@ end subroutine crystallite_stressAndItsTangent subroutine crystallite_integrateStateRK4() use, intrinsic :: & IEEE_arithmetic - use debug, only: & #ifdef DEBUG + use debug, only: & debug_e, & debug_i, & debug_g, & -#endif debug_level, & debug_crystallite, & debug_levelBasic, & debug_levelExtensive, & debug_levelSelective +#endif use FEsolving, only: & FEsolving_execElem, & FEsolving_execIP @@ -1488,17 +1494,17 @@ end subroutine crystallite_integrateStateRK4 subroutine crystallite_integrateStateRKCK45() use, intrinsic :: & IEEE_arithmetic - use debug, only: & #ifdef DEBUG + use debug, only: & debug_e, & debug_i, & debug_g, & -#endif debug_level, & debug_crystallite, & debug_levelBasic, & debug_levelExtensive, & debug_levelSelective +#endif use numerics, only: & rTol_crystalliteState use FEsolving, only: & @@ -1574,8 +1580,10 @@ subroutine crystallite_integrateStateRKCK45() singleRun ! flag indicating computation for single (g,i,e) triple eIter = FEsolving_execElem(1:2) +#ifdef DEBUG if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt) & write(6,'(a,1x,i1)') '<< CRYST >> Runge--Kutta step',1 +#endif ! --- LOOP ITERATOR FOR ELEMENT, GRAIN, IP --- do e = eIter(1),eIter(2) @@ -1970,9 +1978,10 @@ subroutine crystallite_integrateStateRKCK45() ! --- nonlocal convergence check --- - +#ifdef DEBUG if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt) & write(6,'(a,i8,a,i2,/)') '<< CRYST >> ', count(crystallite_converged(:,:,:)), ' grains converged' ! if not requesting Integration of just a single IP +#endif if ((.not. singleRun) .and. any(.not. crystallite_converged .and. .not. crystallite_localPlasticity)) & ! any non-local not yet converged (or broken)... crystallite_converged = crystallite_converged .and. crystallite_localPlasticity ! ...restart all non-local as not converged @@ -1985,17 +1994,17 @@ end subroutine crystallite_integrateStateRKCK45 subroutine crystallite_integrateStateAdaptiveEuler() use, intrinsic :: & IEEE_arithmetic - use debug, only: & #ifdef DEBUG + use debug, only: & debug_e, & debug_i, & debug_g, & -#endif debug_level, & debug_crystallite, & debug_levelBasic, & debug_levelExtensive, & debug_levelSelective +#endif use numerics, only: & rTol_crystalliteState use FEsolving, only: & @@ -2294,9 +2303,10 @@ subroutine crystallite_integrateStateAdaptiveEuler() ! --- NONLOCAL CONVERGENCE CHECK --- - +#ifdef DEBUG if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt) & write(6,'(a,i8,a,i2,/)') '<< CRYST >> ', count(crystallite_converged(:,:,:)), ' grains converged' +#endif if ((.not. singleRun) .and. any(.not. crystallite_converged .and. .not. crystallite_localPlasticity)) & ! any non-local not yet converged (or broken)... crystallite_converged = crystallite_converged .and. crystallite_localPlasticity ! ...restart all non-local as not converged @@ -2310,17 +2320,17 @@ end subroutine crystallite_integrateStateAdaptiveEuler subroutine crystallite_integrateStateEuler() use, intrinsic :: & IEEE_arithmetic - use debug, only: & #ifdef DEBUG + use debug, only: & debug_e, & debug_i, & debug_g, & -#endif debug_level, & debug_crystallite, & debug_levelBasic, & debug_levelExtensive, & debug_levelSelective +#endif use numerics, only: & numerics_timeSyncing use FEsolving, only: & @@ -2524,17 +2534,17 @@ end subroutine crystallite_integrateStateEuler subroutine crystallite_integrateStateFPI() use, intrinsic :: & IEEE_arithmetic - use debug, only: & #ifdef DEBUG + use debug, only: & debug_e, & debug_i, & debug_g, & -#endif debug_level,& debug_crystallite, & debug_levelBasic, & debug_levelExtensive, & debug_levelSelective +#endif use numerics, only: & nState, & rTol_crystalliteState @@ -2598,8 +2608,10 @@ subroutine crystallite_integrateStateFPI() singleRun = (eIter(1) == eIter(2) .and. iIter(1,eIter(1)) == iIter(2,eIter(2))) +#ifdef DEBUG if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt) & write(6,'(a,i8,a)') '<< CRYST >> ', count(crystallite_todo(:,:,:)),' grains todo at start of state integration' +#endif !-------------------------------------------------------------------------------------------------- ! initialize dotState @@ -2654,8 +2666,10 @@ subroutine crystallite_integrateStateFPI() NaN = NaN .or. any(IEEE_is_NaN(sourceState(p)%p(mySource)%dotState(:,c))) enddo if (NaN) then ! NaN occured in any dotState +#ifdef DEBUG if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt) & write(6,*) '<< CRYST >> dotstate ',plasticState(p)%dotState(:,c) +#endif if (.not. crystallite_localPlasticity(g,i,e)) then ! if broken is a non-local... !$OMP CRITICAL (checkTodo) crystallite_todo = crystallite_todo .and. crystallite_localPlasticity ! ...all non-locals done (and broken) @@ -2669,9 +2683,10 @@ subroutine crystallite_integrateStateFPI() !$OMP ENDDO ! --- UPDATE STATE --- +#ifdef DEBUG if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt) & write(6,'(a,i8,a)') '<< CRYST >> ', count(crystallite_todo(:,:,:)),' grains todo after preguess of state' - +#endif !$OMP DO PRIVATE(mySizePlasticDotState,mySizeSourceDotState,p,c) do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains @@ -2727,9 +2742,10 @@ subroutine crystallite_integrateStateFPI() ! --- STRESS INTEGRATION --- +#ifdef DEBUG if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt) & write(6,'(a,i8,a)') '<< CRYST >> ', count(crystallite_todo(:,:,:)),' grains todo before stress integration' - +#endif !$OMP DO do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains !$OMP FLUSH(crystallite_todo) @@ -2745,13 +2761,10 @@ subroutine crystallite_integrateStateFPI() enddo; enddo; enddo !$OMP ENDDO - !$OMP SINGLE - !$OMP CRITICAL (write2out) +#ifdef DEBUG if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt) & write(6,'(a,i8,a)') '<< CRYST >> ', count(crystallite_todo(:,:,:)),' grains todo after stress integration' - !$OMP END CRITICAL (write2out) - !$OMP END SINGLE - +#endif ! --- DOT STATE --- @@ -2940,10 +2953,11 @@ subroutine crystallite_integrateStateFPI() !$OMP ENDDO !$OMP END PARALLEL +#ifdef DEBUG if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt) & write(6,'(a,i8,a,i2)') '<< CRYST >> ', count(crystallite_converged(:,:,:)), & ' grains converged after state integration #', NiterationState - +#endif ! --- NON-LOCAL CONVERGENCE CHECK --- @@ -2952,12 +2966,15 @@ subroutine crystallite_integrateStateFPI() crystallite_converged = crystallite_converged .and. crystallite_localPlasticity ! ...restart all non-local as not converged endif +#ifdef DEBUG if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt) then write(6,'(a,i8,a)') '<< CRYST >> ', count(crystallite_converged(:,:,:)), & ' grains converged after non-local check' write(6,'(a,i8,a,i2,/)') '<< CRYST >> ', count(crystallite_todo(:,:,:)), & ' grains todo after state integration #', NiterationState endif +#endif + ! --- CHECK IF DONE WITH INTEGRATION --- doneWithIntegration = .true. @@ -3118,16 +3135,16 @@ logical function crystallite_integrateStress(& iJacoLpresiduum, & subStepSizeLp, & subStepSizeLi - use debug, only: debug_level, & #ifdef DEBUG + use debug, only: debug_level, & debug_e, & debug_i, & debug_g, & -#endif debug_crystallite, & debug_levelBasic, & debug_levelExtensive, & debug_levelSelective +#endif use constitutive, only: constitutive_LpAndItsTangents, & constitutive_LiAndItsTangents, & From 901355d2ae62409be8d3adfea0de93103c0a16e4 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 20 Sep 2018 06:27:53 +0200 Subject: [PATCH 32/35] don't use unnecessarily long names --- src/crystallite.f90 | 86 ++++++++++++++++++++++----------------------- 1 file changed, 43 insertions(+), 43 deletions(-) diff --git a/src/crystallite.f90 b/src/crystallite.f90 index a7463b9af..74f49920d 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -122,13 +122,13 @@ module crystallite crystallite_postResults private :: & integrateState, & - crystallite_integrateStateFPI, & - crystallite_integrateStateEuler, & - crystallite_integrateStateAdaptiveEuler, & - crystallite_integrateStateRK4, & - crystallite_integrateStateRKCK45, & - crystallite_integrateStress, & - crystallite_stateJump + integrateStateFPI, & + integrateStateEuler, & + integrateStateAdaptiveEuler, & + integrateStateRK4, & + integrateStateRKCK45, & + integrateStress, & + stateJump contains @@ -272,15 +272,15 @@ subroutine crystallite_init select case(numerics_integrator(1)) case(1_pInt) - integrateState => crystallite_integrateStateFPI + integrateState => integrateStateFPI case(2_pInt) - integrateState => crystallite_integrateStateEuler + integrateState => integrateStateEuler case(3_pInt) - integrateState => crystallite_integrateStateAdaptiveEuler + integrateState => integrateStateAdaptiveEuler case(4_pInt) - integrateState => crystallite_integrateStateRK4 + integrateState => integrateStateRK4 case(5_pInt) - integrateState => crystallite_integrateStateRKCK45 + integrateState => integrateStateRKCK45 end select @@ -1201,7 +1201,7 @@ end subroutine crystallite_stressAndItsTangent !-------------------------------------------------------------------------------------------------- !> @brief integrate stress, state with 4th order explicit Runge Kutta method !-------------------------------------------------------------------------------------------------- -subroutine crystallite_integrateStateRK4() +subroutine integrateStateRK4() use, intrinsic :: & IEEE_arithmetic #ifdef DEBUG @@ -1382,7 +1382,7 @@ subroutine crystallite_integrateStateRK4() do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains !$OMP FLUSH(crystallite_todo) if (crystallite_todo(g,i,e)) then - crystallite_todo(g,i,e) = crystallite_stateJump(g,i,e) + crystallite_todo(g,i,e) = stateJump(g,i,e) !$OMP FLUSH(crystallite_todo) if (.not. crystallite_todo(g,i,e) .and. .not. crystallite_localPlasticity(g,i,e)) then ! if broken non-local... !$OMP CRITICAL (checkTodo) @@ -1413,7 +1413,7 @@ subroutine crystallite_integrateStateRK4() do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains !$OMP FLUSH(crystallite_todo) if (crystallite_todo(g,i,e)) then - crystallite_todo(g,i,e) = crystallite_integrateStress(g,i,e,timeStepFraction(n)) ! fraction of original times step + crystallite_todo(g,i,e) = integrateStress(g,i,e,timeStepFraction(n)) ! fraction of original times step !$OMP FLUSH(crystallite_todo) if (.not. crystallite_todo(g,i,e) .and. .not. crystallite_localPlasticity(g,i,e)) then ! if broken non-local... !$OMP CRITICAL (checkTodo) @@ -1484,14 +1484,14 @@ subroutine crystallite_integrateStateRK4() endif endif -end subroutine crystallite_integrateStateRK4 +end subroutine integrateStateRK4 !-------------------------------------------------------------------------------------------------- !> @brief integrate stress, state with 5th order Runge-Kutta Cash-Karp method with !> adaptive step size (use 5th order solution to advance = "local extrapolation") !-------------------------------------------------------------------------------------------------- -subroutine crystallite_integrateStateRKCK45() +subroutine integrateStateRKCK45() use, intrinsic :: & IEEE_arithmetic #ifdef DEBUG @@ -1703,7 +1703,7 @@ subroutine crystallite_integrateStateRKCK45() do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains !$OMP FLUSH(crystallite_todo) if (crystallite_todo(g,i,e)) then - crystallite_todo(g,i,e) = crystallite_stateJump(g,i,e) + crystallite_todo(g,i,e) = stateJump(g,i,e) !$OMP FLUSH(crystallite_todo) if (.not. crystallite_todo(g,i,e) .and. .not. crystallite_localPlasticity(g,i,e)) then ! if broken non-local... !$OMP CRITICAL (checkTodo) @@ -1734,7 +1734,7 @@ subroutine crystallite_integrateStateRKCK45() do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains !$OMP FLUSH(crystallite_todo) if (crystallite_todo(g,i,e)) then - crystallite_todo(g,i,e) = crystallite_integrateStress(g,i,e,C(stage)) ! fraction of original time step + crystallite_todo(g,i,e) = integrateStress(g,i,e,C(stage)) ! fraction of original time step !$OMP FLUSH(crystallite_todo) if (.not. crystallite_todo(g,i,e) .and. .not. crystallite_localPlasticity(g,i,e)) then ! if broken non-local... !$OMP CRITICAL (checkTodo) @@ -1923,7 +1923,7 @@ subroutine crystallite_integrateStateRKCK45() do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains !$OMP FLUSH(crystallite_todo) if (crystallite_todo(g,i,e)) then - crystallite_todo(g,i,e) = crystallite_stateJump(g,i,e) + crystallite_todo(g,i,e) = stateJump(g,i,e) !$OMP FLUSH(crystallite_todo) if (.not. crystallite_todo(g,i,e) .and. .not. crystallite_localPlasticity(g,i,e)) then ! if broken non-local... !$OMP CRITICAL (checkTodo) @@ -1954,7 +1954,7 @@ subroutine crystallite_integrateStateRKCK45() do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains !$OMP FLUSH(crystallite_todo) if (crystallite_todo(g,i,e)) then - crystallite_todo(g,i,e) = crystallite_integrateStress(g,i,e) + crystallite_todo(g,i,e) = integrateStress(g,i,e) !$OMP FLUSH(crystallite_todo) if (.not. crystallite_todo(g,i,e) .and. .not. crystallite_localPlasticity(g,i,e)) then ! if broken non-local... !$OMP CRITICAL (checkTodo) @@ -1985,13 +1985,13 @@ subroutine crystallite_integrateStateRKCK45() if ((.not. singleRun) .and. any(.not. crystallite_converged .and. .not. crystallite_localPlasticity)) & ! any non-local not yet converged (or broken)... crystallite_converged = crystallite_converged .and. crystallite_localPlasticity ! ...restart all non-local as not converged -end subroutine crystallite_integrateStateRKCK45 +end subroutine integrateStateRKCK45 !-------------------------------------------------------------------------------------------------- !> @brief integrate stress, state with 1st order Euler method with adaptive step size !-------------------------------------------------------------------------------------------------- -subroutine crystallite_integrateStateAdaptiveEuler() +subroutine integrateStateAdaptiveEuler() use, intrinsic :: & IEEE_arithmetic #ifdef DEBUG @@ -2150,7 +2150,7 @@ subroutine crystallite_integrateStateAdaptiveEuler() do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains !$OMP FLUSH(crystallite_todo) if (crystallite_todo(g,i,e)) then - crystallite_todo(g,i,e) = crystallite_stateJump(g,i,e) + crystallite_todo(g,i,e) = stateJump(g,i,e) !$OMP FLUSH(crystallite_todo) if (.not. crystallite_todo(g,i,e) .and. .not. crystallite_localPlasticity(g,i,e)) then ! if broken non-local... !$OMP CRITICAL (checkTodo) @@ -2182,7 +2182,7 @@ subroutine crystallite_integrateStateAdaptiveEuler() do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains !$OMP FLUSH(crystallite_todo) if (crystallite_todo(g,i,e)) then - crystallite_todo(g,i,e) = crystallite_integrateStress(g,i,e) + crystallite_todo(g,i,e) = integrateStress(g,i,e) !$OMP FLUSH(crystallite_todo) if (.not. crystallite_todo(g,i,e) .and. .not. crystallite_localPlasticity(g,i,e)) then ! if broken non-local... !$OMP CRITICAL (checkTodo) @@ -2311,13 +2311,13 @@ subroutine crystallite_integrateStateAdaptiveEuler() crystallite_converged = crystallite_converged .and. crystallite_localPlasticity ! ...restart all non-local as not converged -end subroutine crystallite_integrateStateAdaptiveEuler +end subroutine integrateStateAdaptiveEuler !-------------------------------------------------------------------------------------------------- !> @brief integrate stress, and state with 1st order explicit Euler method !-------------------------------------------------------------------------------------------------- -subroutine crystallite_integrateStateEuler() +subroutine integrateStateEuler() use, intrinsic :: & IEEE_arithmetic #ifdef DEBUG @@ -2458,7 +2458,7 @@ eIter = FEsolving_execElem(1:2) do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains !$OMP FLUSH(crystallite_todo) if (crystallite_todo(g,i,e) .and. .not. crystallite_converged(g,i,e)) then - crystallite_todo(g,i,e) = crystallite_stateJump(g,i,e) + crystallite_todo(g,i,e) = stateJump(g,i,e) !$OMP FLUSH(crystallite_todo) if (.not. crystallite_todo(g,i,e) .and. .not. crystallite_localPlasticity(g,i,e) & ! if broken non-local... .and. .not. numerics_timeSyncing) then @@ -2492,7 +2492,7 @@ eIter = FEsolving_execElem(1:2) do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains !$OMP FLUSH(crystallite_todo) if (crystallite_todo(g,i,e) .and. .not. crystallite_converged(g,i,e)) then - crystallite_todo(g,i,e) = crystallite_integrateStress(g,i,e) + crystallite_todo(g,i,e) = integrateStress(g,i,e) !$OMP FLUSH(crystallite_todo) if (.not. crystallite_todo(g,i,e) .and. .not. crystallite_localPlasticity(g,i,e) & ! if broken non-local... .and. .not. numerics_timeSyncing) then @@ -2524,14 +2524,14 @@ eIter = FEsolving_execElem(1:2) crystallite_converged = crystallite_converged .and. crystallite_localPlasticity ! ...restart all non-local as not converged endif -end subroutine crystallite_integrateStateEuler +end subroutine integrateStateEuler !-------------------------------------------------------------------------------------------------- !> @brief integrate stress, state with adaptive 1st order explicit Euler method !> using Fixed Point Iteration to adapt the stepsize !-------------------------------------------------------------------------------------------------- -subroutine crystallite_integrateStateFPI() +subroutine integrateStateFPI() use, intrinsic :: & IEEE_arithmetic #ifdef DEBUG @@ -2750,7 +2750,7 @@ subroutine crystallite_integrateStateFPI() do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains !$OMP FLUSH(crystallite_todo) if (crystallite_todo(g,i,e) .and. .not. crystallite_converged(g,i,e)) then - crystallite_todo(g,i,e) = crystallite_integrateStress(g,i,e) + crystallite_todo(g,i,e) = integrateStress(g,i,e) !$OMP FLUSH(crystallite_todo) if (.not. crystallite_todo(g,i,e) .and. .not. crystallite_localPlasticity(g,i,e)) then ! broken non-local... !$OMP CRITICAL (checkTodo) @@ -2938,7 +2938,7 @@ subroutine crystallite_integrateStateFPI() do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains !$OMP FLUSH(crystallite_todo) if (crystallite_todo(g,i,e) .and. crystallite_converged(g,i,e)) then ! converged and still alive... - crystallite_todo(g,i,e) = crystallite_stateJump(g,i,e) + crystallite_todo(g,i,e) = stateJump(g,i,e) !$OMP FLUSH(crystallite_todo) if (.not. crystallite_todo(g,i,e)) then ! if state jump fails, then convergence is broken crystallite_converged(g,i,e) = .false. @@ -2988,14 +2988,14 @@ subroutine crystallite_integrateStateFPI() enddo elemLoop enddo crystalliteLooping -end subroutine crystallite_integrateStateFPI +end subroutine integrateStateFPI !-------------------------------------------------------------------------------------------------- !> @brief calculates a jump in the state according to the current state and the current stress !> returns true, if state jump was successfull or not needed. false indicates NaN in delta state !-------------------------------------------------------------------------------------------------- -logical function crystallite_stateJump(ipc,ip,el) +logical function stateJump(ipc,ip,el) use, intrinsic :: & IEEE_arithmetic use prec, only: & @@ -3045,7 +3045,7 @@ logical function crystallite_stateJump(ipc,ip,el) mySizePlasticDeltaState = plasticState(p)%sizeDeltaState if( any(IEEE_is_NaN(plasticState(p)%deltaState(1:mySizePlasticDeltaState,c)))) then ! NaN occured in deltaState - crystallite_stateJump = .false. + stateJump = .false. return endif @@ -3059,7 +3059,7 @@ logical function crystallite_stateJump(ipc,ip,el) myOffsetSourceDeltaState = sourceState(p)%p(mySource)%offsetDeltaState mySizeSourceDeltaState = sourceState(p)%p(mySource)%sizeDeltaState if (any(IEEE_is_NaN(sourceState(p)%p(mySource)%deltaState(1:mySizeSourceDeltaState,c)))) then ! NaN occured in deltaState - crystallite_stateJump = .false. + stateJump = .false. return endif sourceState(p)%p(mySource)%state(myOffsetSourceDeltaState + 1_pInt : & @@ -3082,9 +3082,9 @@ logical function crystallite_stateJump(ipc,ip,el) endif #endif - crystallite_stateJump = .true. + stateJump = .true. -end function crystallite_stateJump +end function stateJump !-------------------------------------------------------------------------------------------------- @@ -3118,7 +3118,7 @@ end function crystallite_push33ToRef !> @brief calculation of stress (P) with time integration based on a residuum in Lp and !> intermediate acceleration of the Newton-Raphson correction !-------------------------------------------------------------------------------------------------- -logical function crystallite_integrateStress(& +logical function integrateStress(& ipc,& ! grain number ip,& ! integration point number el,& ! element number @@ -3236,7 +3236,7 @@ logical function crystallite_integrateStress(& dgesv !* be pessimistic - crystallite_integrateStress = .false. + integrateStress = .false. #ifdef DEBUG if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt & .and. ((el == debug_e .and. ip == debug_i .and. ipc == debug_g) & @@ -3575,7 +3575,7 @@ logical function crystallite_integrateStress(& !* set return flag to true - crystallite_integrateStress = .true. + integrateStress = .true. #ifdef DEBUG if (iand(debug_level(debug_crystallite),debug_levelExtensive) /= 0_pInt & .and. ((el == debug_e .and. ip == debug_i .and. ipc == debug_g) & @@ -3590,7 +3590,7 @@ logical function crystallite_integrateStress(& endif #endif -end function crystallite_integrateStress +end function integrateStress !-------------------------------------------------------------------------------------------------- From fcff6b908a444c515a5c06c2cdfa96a9476b734f Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 20 Sep 2018 06:35:30 +0200 Subject: [PATCH 33/35] can be easily computed during post processing --- examples/ConfigFiles/Crystallite_All.config | 11 ++--- src/crystallite.f90 | 45 ++------------------- 2 files changed, 6 insertions(+), 50 deletions(-) diff --git a/examples/ConfigFiles/Crystallite_All.config b/examples/ConfigFiles/Crystallite_All.config index 761380fcd..d46c3e0e6 100644 --- a/examples/ConfigFiles/Crystallite_All.config +++ b/examples/ConfigFiles/Crystallite_All.config @@ -5,15 +5,10 @@ (output) orientation # quaternion (output) eulerangles # orientation as Bunge triple in degree (output) grainrotation # deviation from initial orientation as axis (1-3) and angle in degree (4) in crystal reference coordinates -(output) grainrotationx # deviation from initial orientation as angle in degrees around sample reference x axis -(output) grainrotationy # deviation from initial orientation as angle in degrees around sample reference y axis -(output) grainrotationz # deviation from initial orientation as angle in degrees around sample reference z axis -(output) f # deformation gradient tensor; synonyms: "defgrad" +(output) f # deformation gradient tensor (output) fe # elastic deformation gradient tensor (output) fp # plastic deformation gradient tensor -(output) e # total strain as Green-Lagrange tensor -(output) ee # elastic strain as Green-Lagrange tensor -(output) p # first Piola-Kichhoff stress tensor; synonyms: "firstpiola", "1stpiola" -(output) s # second Piola-Kichhoff stress tensor; synonyms: "tstar", "secondpiola", "2ndpiola" +(output) p # first Piola-Kichhoff stress tensor +(output) s # second Piola-Kichhoff stress tensor (output) lp # plastic velocity gradient tensor (output) elasmatrix # elastic stiffness matrix diff --git a/src/crystallite.f90 b/src/crystallite.f90 index 74f49920d..b079c7c72 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -90,9 +90,6 @@ module crystallite phase_ID, & texture_ID, & volume_ID, & - grainrotationx_ID, & - grainrotationy_ID, & - grainrotationz_ID, & orientation_ID, & grainrotation_ID, & eulerangles_ID, & @@ -102,8 +99,6 @@ module crystallite fi_ID, & lp_ID, & li_ID, & - e_ID, & - ee_ID, & p_ID, & s_ID, & elasmatrix_ID, & @@ -302,12 +297,6 @@ subroutine crystallite_init crystallite_outputID(o,c) = texture_ID case ('volume') outputName crystallite_outputID(o,c) = volume_ID - case ('grainrotationx') outputName - crystallite_outputID(o,c) = grainrotationx_ID - case ('grainrotationy') outputName - crystallite_outputID(o,c) = grainrotationy_ID - case ('grainrotationz') outputName - crystallite_outputID(o,c) = grainrotationx_ID case ('orientation') outputName crystallite_outputID(o,c) = orientation_ID case ('grainrotation') outputName @@ -326,10 +315,6 @@ subroutine crystallite_init crystallite_outputID(o,c) = lp_ID case ('li') outputName crystallite_outputID(o,c) = li_ID - case ('e') outputName - crystallite_outputID(o,c) = e_ID - case ('ee') outputName - crystallite_outputID(o,c) = ee_ID case ('p','firstpiola','1stpiola') outputName crystallite_outputID(o,c) = p_ID case ('s','tstar','secondpiola','2ndpiola') outputName @@ -350,13 +335,13 @@ subroutine crystallite_init do r = 1_pInt,size(config_crystallite) do o = 1_pInt,crystallite_Noutput(r) select case(crystallite_outputID(o,r)) - case(phase_ID,texture_ID,volume_ID,grainrotationx_ID,grainrotationy_ID,grainrotationz_ID) + case(phase_ID,texture_ID,volume_ID) mySize = 1_pInt case(orientation_ID,grainrotation_ID) mySize = 4_pInt case(eulerangles_ID) mySize = 3_pInt - case(defgrad_ID,fe_ID,fp_ID,fi_ID,lp_ID,li_ID,e_ID,ee_ID,p_ID,s_ID) + case(defgrad_ID,fe_ID,fp_ID,fi_ID,lp_ID,li_ID,p_ID,s_ID) mySize = 9_pInt case(elasmatrix_ID) mySize = 36_pInt @@ -3706,9 +3691,7 @@ function crystallite_postResults(ipc, ip, el) math_det33, & math_I3, & inDeg, & - math_Mandel6to33, & - math_qMul, & - math_qConj + math_Mandel6to33 use mesh, only: & mesh_element, & mesh_ipVolume, & @@ -3786,18 +3769,6 @@ function crystallite_postResults(ipc, ip, el) crystallite_postResults(c+1:c+mySize) = & math_qToEulerAxisAngle(crystallite_rotation(1:4,ipc,ip,el)) ! grain rotation away from initial orientation as axis-angle in sample reference coordinates crystallite_postResults(c+4) = inDeg * crystallite_postResults(c+4) ! angle in degree - case (grainrotationx_ID) - mySize = 1_pInt - rotation = math_qToEulerAxisAngle(crystallite_rotation(1:4,ipc,ip,el)) ! grain rotation away from initial orientation as axis-angle in sample reference coordinates - crystallite_postResults(c+1) = inDeg * rotation(1) * rotation(4) ! angle in degree - case (grainrotationy_ID) - mySize = 1_pInt - rotation = math_qToEulerAxisAngle(crystallite_rotation(1:4,ipc,ip,el)) ! grain rotation away from initial orientation as axis-angle in sample reference coordinates - crystallite_postResults(c+1) = inDeg * rotation(2) * rotation(4) ! angle in degree - case (grainrotationz_ID) - mySize = 1_pInt - rotation = math_qToEulerAxisAngle(crystallite_rotation(1:4,ipc,ip,el)) ! grain rotation away from initial orientation as axis-angle in sample reference coordinates - crystallite_postResults(c+1) = inDeg * rotation(3) * rotation(4) ! angle in degree ! remark: tensor output is of the form 11,12,13, 21,22,23, 31,32,33 ! thus row index i is slow, while column index j is fast. reminder: "row is slow" @@ -3806,20 +3777,10 @@ function crystallite_postResults(ipc, ip, el) mySize = 9_pInt crystallite_postResults(c+1:c+mySize) = & reshape(transpose(crystallite_partionedF(1:3,1:3,ipc,ip,el)),[mySize]) - case (e_ID) - mySize = 9_pInt - crystallite_postResults(c+1:c+mySize) = 0.5_pReal * reshape((math_mul33x33( & - transpose(crystallite_partionedF(1:3,1:3,ipc,ip,el)), & - crystallite_partionedF(1:3,1:3,ipc,ip,el)) - math_I3),[mySize]) case (fe_ID) mySize = 9_pInt crystallite_postResults(c+1:c+mySize) = & reshape(transpose(crystallite_Fe(1:3,1:3,ipc,ip,el)),[mySize]) - case (ee_ID) - Ee = 0.5_pReal *(math_mul33x33(transpose(crystallite_Fe(1:3,1:3,ipc,ip,el)), & - crystallite_Fe(1:3,1:3,ipc,ip,el)) - math_I3) - mySize = 9_pInt - crystallite_postResults(c+1:c+mySize) = reshape(Ee,[mySize]) case (fp_ID) mySize = 9_pInt crystallite_postResults(c+1:c+mySize) = & From 310ea62964a9dfeee83df3ba4cddfecb014239ea Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 20 Sep 2018 06:58:31 +0200 Subject: [PATCH 34/35] only print out the essential information --- .../Polycrystal/material.config | 2 - src/crystallite.f90 | 58 ++----------------- 2 files changed, 6 insertions(+), 54 deletions(-) diff --git a/examples/SpectralMethod/Polycrystal/material.config b/examples/SpectralMethod/Polycrystal/material.config index 5073f165e..eebb17473 100644 --- a/examples/SpectralMethod/Polycrystal/material.config +++ b/examples/SpectralMethod/Polycrystal/material.config @@ -18,8 +18,6 @@ mech none (output) f # deformation gradient tensor; synonyms: "defgrad" (output) fe # elastic deformation gradient tensor (output) fp # plastic deformation gradient tensor -(output) e # total strain as Green-Lagrange tensor -(output) ee # elastic strain as Green-Lagrange tensor (output) p # first Piola-Kichhoff stress tensor; synonyms: "firstpiola", "1stpiola" (output) lp # plastic velocity gradient tensor diff --git a/src/crystallite.f90 b/src/crystallite.f90 index b079c7c72..5f447664d 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -194,8 +194,6 @@ subroutine crystallite_init mySize character(len=65536), dimension(:), allocatable :: str - character(len=65536) :: & - tag = '' write(6,'(/,a)') ' <<<+- crystallite init -+>>>' write(6,'(a15,a)') ' Current time: ',IO_timeStamp() @@ -326,7 +324,7 @@ subroutine crystallite_init case ('neighboringelement') outputName crystallite_outputID(o,c) = neighboringelement_ID case default outputName - call IO_error(105_pInt,ext_msg=tag//' (Crystallite)') + call IO_error(105_pInt,ext_msg=trim(str(o))//' (Crystallite)') end select outputName enddo enddo @@ -424,51 +422,11 @@ subroutine crystallite_init #ifdef DEBUG if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) then - write(6,'(a35,1x,7(i8,1x))') 'crystallite_Fe: ', shape(crystallite_Fe) - write(6,'(a35,1x,7(i8,1x))') 'crystallite_Fp: ', shape(crystallite_Fp) - write(6,'(a35,1x,7(i8,1x))') 'crystallite_Fi: ', shape(crystallite_Fi) - write(6,'(a35,1x,7(i8,1x))') 'crystallite_Lp: ', shape(crystallite_Lp) - write(6,'(a35,1x,7(i8,1x))') 'crystallite_Li: ', shape(crystallite_Li) - write(6,'(a35,1x,7(i8,1x))') 'crystallite_F0: ', shape(crystallite_F0) - write(6,'(a35,1x,7(i8,1x))') 'crystallite_Fp0: ', shape(crystallite_Fp0) - write(6,'(a35,1x,7(i8,1x))') 'crystallite_Fi0: ', shape(crystallite_Fi0) - write(6,'(a35,1x,7(i8,1x))') 'crystallite_Lp0: ', shape(crystallite_Lp0) - write(6,'(a35,1x,7(i8,1x))') 'crystallite_Li0: ', shape(crystallite_Li0) - write(6,'(a35,1x,7(i8,1x))') 'crystallite_partionedF: ', shape(crystallite_partionedF) - write(6,'(a35,1x,7(i8,1x))') 'crystallite_partionedF0: ', shape(crystallite_partionedF0) - write(6,'(a35,1x,7(i8,1x))') 'crystallite_partionedFp0: ', shape(crystallite_partionedFp0) - write(6,'(a35,1x,7(i8,1x))') 'crystallite_partionedFi0: ', shape(crystallite_partionedFi0) - write(6,'(a35,1x,7(i8,1x))') 'crystallite_partionedLp0: ', shape(crystallite_partionedLp0) - write(6,'(a35,1x,7(i8,1x))') 'crystallite_partionedLi0: ', shape(crystallite_partionedLi0) - write(6,'(a35,1x,7(i8,1x))') 'crystallite_subF: ', shape(crystallite_subF) - write(6,'(a35,1x,7(i8,1x))') 'crystallite_subF0: ', shape(crystallite_subF0) - write(6,'(a35,1x,7(i8,1x))') 'crystallite_subFp0: ', shape(crystallite_subFp0) - write(6,'(a35,1x,7(i8,1x))') 'crystallite_subFi0: ', shape(crystallite_subFi0) - write(6,'(a35,1x,7(i8,1x))') 'crystallite_subLp0: ', shape(crystallite_subLp0) - write(6,'(a35,1x,7(i8,1x))') 'crystallite_subLi0: ', shape(crystallite_subLi0) - write(6,'(a35,1x,7(i8,1x))') 'crystallite_P: ', shape(crystallite_P) - write(6,'(a35,1x,7(i8,1x))') 'crystallite_Tstar_v: ', shape(crystallite_Tstar_v) - write(6,'(a35,1x,7(i8,1x))') 'crystallite_Tstar0_v: ', shape(crystallite_Tstar0_v) - write(6,'(a35,1x,7(i8,1x))') 'crystallite_partionedTstar0_v: ', shape(crystallite_partionedTstar0_v) - write(6,'(a35,1x,7(i8,1x))') 'crystallite_subTstar0_v: ', shape(crystallite_subTstar0_v) - write(6,'(a35,1x,7(i8,1x))') 'crystallite_dPdF: ', shape(crystallite_dPdF) - write(6,'(a35,1x,7(i8,1x))') 'crystallite_dPdF0: ', shape(crystallite_dPdF0) - write(6,'(a35,1x,7(i8,1x))') 'crystallite_partioneddPdF0: ', shape(crystallite_partioneddPdF0) - write(6,'(a35,1x,7(i8,1x))') 'crystallite_orientation: ', shape(crystallite_orientation) - write(6,'(a35,1x,7(i8,1x))') 'crystallite_orientation0: ', shape(crystallite_orientation0) - write(6,'(a35,1x,7(i8,1x))') 'crystallite_rotation: ', shape(crystallite_rotation) - write(6,'(a35,1x,7(i8,1x))') 'crystallite_disorientation: ', shape(crystallite_disorientation) - write(6,'(a35,1x,7(i8,1x))') 'crystallite_dt: ', shape(crystallite_dt) - write(6,'(a35,1x,7(i8,1x))') 'crystallite_subdt: ', shape(crystallite_subdt) - write(6,'(a35,1x,7(i8,1x))') 'crystallite_subFrac: ', shape(crystallite_subFrac) - write(6,'(a35,1x,7(i8,1x))') 'crystallite_subStep: ', shape(crystallite_subStep) - write(6,'(a35,1x,7(i8,1x))') 'crystallite_localPlasticity: ', shape(crystallite_localPlasticity) - write(6,'(a35,1x,7(i8,1x))') 'crystallite_requested: ', shape(crystallite_requested) - write(6,'(a35,1x,7(i8,1x))') 'crystallite_todo: ', shape(crystallite_todo) - write(6,'(a35,1x,7(i8,1x))') 'crystallite_converged: ', shape(crystallite_converged) - write(6,'(a35,1x,7(i8,1x))') 'crystallite_sizePostResults: ', shape(crystallite_sizePostResults) - write(6,'(a35,1x,7(i8,1x))') 'crystallite_sizePostResult: ', shape(crystallite_sizePostResult) - write(6,'(/,a35,1x,i10)') 'Number of nonlocal grains: ',count(.not. crystallite_localPlasticity) + write(6,'(a42,1x,i10)') ' # of elements: ', eMax + write(6,'(a42,1x,i10)') 'max # of integration points/element: ', iMax + write(6,'(a42,1x,i10)') 'max # of constituents/integration point: ', cMax + write(6,'(a42,1x,i10)') 'max # of neigbours/integration point: ', nMax + write(6,'(a42,1x,i10)') ' # of nonlocal constituents: ',count(.not. crystallite_localPlasticity) flush(6) endif @@ -3722,10 +3680,6 @@ function crystallite_postResults(ipc, ip, el) 1+plasticState(material_phase(ipc,ip,el))%sizePostResults + & sum(sourceState(material_phase(ipc,ip,el))%p(:)%sizePostResults)) :: & crystallite_postResults - real(pReal), dimension(3,3) :: & - Ee - real(pReal), dimension(4) :: & - rotation real(pReal) :: & detF integer(pInt) :: & From 3e7b80a3efb296f7d2fa0476f618daefdc3b80e2 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 20 Sep 2018 07:27:12 +0200 Subject: [PATCH 35/35] debug only available if compiled in debug mode --- src/homogenization.f90 | 15 +++++++-------- 1 file changed, 7 insertions(+), 8 deletions(-) diff --git a/src/homogenization.f90 b/src/homogenization.f90 index a2b5d12e7..82a97dc53 100644 --- a/src/homogenization.f90 +++ b/src/homogenization.f90 @@ -494,6 +494,7 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt) crystallite_converged, & crystallite_stressAndItsTangent, & crystallite_orientations +#ifdef DEBUG use debug, only: & debug_level, & debug_homogenization, & @@ -502,6 +503,7 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt) debug_levelSelective, & debug_e, & debug_i +#endif implicit none real(pReal), intent(in) :: dt !< time increment @@ -515,18 +517,16 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt) mySource, & myNgrains -!-------------------------------------------------------------------------------------------------- -! initialize to starting condition +#ifdef DEBUG if (iand(debug_level(debug_homogenization), debug_levelBasic) /= 0_pInt) then - !$OMP CRITICAL (write2out) write(6,'(/a,i5,1x,i2)') '<< HOMOG >> Material Point start at el ip ', debug_e, debug_i write(6,'(a,/,3(12x,3(f14.9,1x)/))') '<< HOMOG >> F0', & transpose(materialpoint_F0(1:3,1:3,debug_i,debug_e)) write(6,'(a,/,3(12x,3(f14.9,1x)/))') '<< HOMOG >> F', & transpose(materialpoint_F(1:3,1:3,debug_i,debug_e)) - !$OMP END CRITICAL (write2out) endif +#endif !-------------------------------------------------------------------------------------------------- ! initialize restoration points of ... @@ -746,8 +746,9 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt) if (materialpoint_subStep(i,e) > subStepMinHomog) then materialpoint_requested(i,e) = .true. - materialpoint_subF(1:3,1:3,i,e) = materialpoint_subF0(1:3,1:3,i,e) + & - materialpoint_subStep(i,e) * (materialpoint_F(1:3,1:3,i,e) - materialpoint_F0(1:3,1:3,i,e)) + materialpoint_subF(1:3,1:3,i,e) = materialpoint_subF0(1:3,1:3,i,e) & + + materialpoint_subStep(i,e) * (materialpoint_F(1:3,1:3,i,e) & + - materialpoint_F0(1:3,1:3,i,e)) materialpoint_subdt(i,e) = materialpoint_subStep(i,e) * dt materialpoint_doneAndHappy(1:2,i,e) = [.false.,.true.] endif @@ -825,9 +826,7 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt) enddo elementLooping4 !$OMP END PARALLEL DO else - !$OMP CRITICAL (write2out) write(6,'(/,a,/)') '<< HOMOG >> Material Point terminally ill' - !$OMP END CRITICAL (write2out) endif end subroutine materialpoint_stressAndItsTangent