diff --git a/src/commercialFEM_fileList.f90 b/src/commercialFEM_fileList.f90 index d145d3965..621a091d8 100644 --- a/src/commercialFEM_fileList.f90 +++ b/src/commercialFEM_fileList.f90 @@ -32,15 +32,15 @@ #include "phase_mechanics_eigendeformation.f90" #include "phase_mechanics_eigendeformation_cleavageopening.f90" #include "phase_mechanics_eigendeformation_slipplaneopening.f90" +#include "phase_thermal.f90" +#include "phase_thermal_dissipation.f90" +#include "phase_thermal_externalheat.f90" +#include "phase_mechanics_eigendeformation_thermalexpansion.f90" #include "phase_damage.f90" #include "phase_damage_isobrittle.f90" #include "phase_damage_isoductile.f90" #include "phase_damage_anisobrittle.f90" #include "phase_damage_anisoductile.f90" -#include "phase_thermal.f90" -#include "phase_mechanics_eigendeformation_thermalexpansion.f90" -#include "phase_thermal_dissipation.f90" -#include "phase_thermal_externalheat.f90" #include "damage_none.f90" #include "damage_nonlocal.f90" #include "homogenization.f90" diff --git a/src/phase_mechanics_plastic.f90 b/src/phase_mechanics_plastic.f90 index 035b03171..f6c010485 100644 --- a/src/phase_mechanics_plastic.f90 +++ b/src/phase_mechanics_plastic.f90 @@ -2,7 +2,7 @@ submodule(phase:mechanics) plastic interface - module subroutine isotropic_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,me) + module subroutine isotropic_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,me) real(pReal), dimension(3,3), intent(out) :: & Lp real(pReal), dimension(3,3,3,3), intent(out) :: & @@ -10,11 +10,11 @@ submodule(phase:mechanics) plastic real(pReal), dimension(3,3), intent(in) :: & Mp integer, intent(in) :: & - instance, & + ph, & me end subroutine isotropic_LpAndItsTangent - pure module subroutine phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,me) + pure module subroutine phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,me) real(pReal), dimension(3,3), intent(out) :: & Lp real(pReal), dimension(3,3,3,3), intent(out) :: & @@ -22,11 +22,11 @@ submodule(phase:mechanics) plastic real(pReal), dimension(3,3), intent(in) :: & Mp integer, intent(in) :: & - instance, & + ph, & me end subroutine phenopowerlaw_LpAndItsTangent - pure module subroutine kinehardening_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,me) + pure module subroutine kinehardening_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,me) real(pReal), dimension(3,3), intent(out) :: & Lp real(pReal), dimension(3,3,3,3), intent(out) :: & @@ -34,11 +34,11 @@ submodule(phase:mechanics) plastic real(pReal), dimension(3,3), intent(in) :: & Mp integer, intent(in) :: & - instance, & + ph, & me end subroutine kinehardening_LpAndItsTangent - module subroutine dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,instance,me) + module subroutine dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,ph,me) real(pReal), dimension(3,3), intent(out) :: & Lp real(pReal), dimension(3,3,3,3), intent(out) :: & @@ -49,11 +49,11 @@ submodule(phase:mechanics) plastic real(pReal), intent(in) :: & T integer, intent(in) :: & - instance, & + ph, & me end subroutine dislotwin_LpAndItsTangent - pure module subroutine dislotungsten_LpAndItsTangent(Lp,dLp_dMp,Mp,T,instance,me) + pure module subroutine dislotungsten_LpAndItsTangent(Lp,dLp_dMp,Mp,T,ph,me) real(pReal), dimension(3,3), intent(out) :: & Lp real(pReal), dimension(3,3,3,3), intent(out) :: & @@ -64,12 +64,12 @@ submodule(phase:mechanics) plastic real(pReal), intent(in) :: & T integer, intent(in) :: & - instance, & + ph, & me end subroutine dislotungsten_LpAndItsTangent module subroutine nonlocal_LpAndItsTangent(Lp,dLp_dMp, & - Mp,Temperature,instance,me,ip,el) + Mp,Temperature,ph,me,ip,el) real(pReal), dimension(3,3), intent(out) :: & Lp real(pReal), dimension(3,3,3,3), intent(out) :: & @@ -80,71 +80,71 @@ submodule(phase:mechanics) plastic real(pReal), intent(in) :: & Temperature integer, intent(in) :: & - instance, & + ph, & me, & ip, & !< current integration point el !< current element number end subroutine nonlocal_LpAndItsTangent - module subroutine isotropic_dotState(Mp,instance,me) + module subroutine isotropic_dotState(Mp,ph,me) real(pReal), dimension(3,3), intent(in) :: & Mp !< Mandel stress integer, intent(in) :: & - instance, & + ph, & me end subroutine isotropic_dotState - module subroutine phenopowerlaw_dotState(Mp,instance,me) + module subroutine phenopowerlaw_dotState(Mp,ph,me) real(pReal), dimension(3,3), intent(in) :: & Mp !< Mandel stress integer, intent(in) :: & - instance, & + ph, & me end subroutine phenopowerlaw_dotState - module subroutine plastic_kinehardening_dotState(Mp,instance,me) + module subroutine plastic_kinehardening_dotState(Mp,ph,me) real(pReal), dimension(3,3), intent(in) :: & Mp !< Mandel stress integer, intent(in) :: & - instance, & + ph, & me end subroutine plastic_kinehardening_dotState - module subroutine dislotwin_dotState(Mp,T,instance,me) + module subroutine dislotwin_dotState(Mp,T,ph,me) real(pReal), dimension(3,3), intent(in) :: & Mp !< Mandel stress real(pReal), intent(in) :: & T integer, intent(in) :: & - instance, & + ph, & me end subroutine dislotwin_dotState - module subroutine dislotungsten_dotState(Mp,T,instance,me) + module subroutine dislotungsten_dotState(Mp,T,ph,me) real(pReal), dimension(3,3), intent(in) :: & Mp !< Mandel stress real(pReal), intent(in) :: & T integer, intent(in) :: & - instance, & + ph, & me end subroutine dislotungsten_dotState - module subroutine nonlocal_dotState(Mp,Temperature,timestep,instance,me,ip,el) + module subroutine nonlocal_dotState(Mp,Temperature,timestep,ph,me,ip,el) real(pReal), dimension(3,3), intent(in) :: & Mp !< MandelStress real(pReal), intent(in) :: & Temperature, & !< temperature timestep !< substepped crystallite time increment integer, intent(in) :: & - instance, & + ph, & me, & ip, & !< current integration point el !< current element number end subroutine nonlocal_dotState - module subroutine dislotwin_dependentState(T,instance,me) + module subroutine dislotwin_dependentState(T,instance,me) integer, intent(in) :: & instance, & me @@ -213,13 +213,12 @@ module subroutine plastic_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, & real(pReal), dimension(3,3) :: & Mp !< Mandel stress work conjugate with Lp integer :: & - i, j, instance, me, ph + i, j, me, ph Mp = matmul(matmul(transpose(Fi),Fi),S) me = material_phasememberAt(co,ip,el) ph = material_phaseAt(co,el) - instance = phase_plasticityInstance(ph) plasticityType: select case (phase_plasticity(material_phaseAt(co,el))) @@ -228,22 +227,22 @@ module subroutine plastic_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, & dLp_dMp = 0.0_pReal case (PLASTICITY_ISOTROPIC_ID) plasticityType - call isotropic_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,me) + call isotropic_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,me) case (PLASTICITY_PHENOPOWERLAW_ID) plasticityType - call phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,me) + call phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,me) case (PLASTICITY_KINEHARDENING_ID) plasticityType - call kinehardening_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,me) + call kinehardening_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,me) case (PLASTICITY_NONLOCAL_ID) plasticityType - call nonlocal_LpAndItsTangent(Lp,dLp_dMp,Mp, thermal_T(ph,me),instance,me,ip,el) + call nonlocal_LpAndItsTangent(Lp,dLp_dMp,Mp, thermal_T(ph,me),ph,me,ip,el) case (PLASTICITY_DISLOTWIN_ID) plasticityType - call dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp, thermal_T(ph,me),instance,me) + call dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp, thermal_T(ph,me),ph,me) case (PLASTICITY_DISLOTUNGSTEN_ID) plasticityType - call dislotungsten_LpAndItsTangent(Lp,dLp_dMp,Mp, thermal_T(ph,me),instance,me) + call dislotungsten_LpAndItsTangent(Lp,dLp_dMp,Mp, thermal_T(ph,me),ph,me) end select plasticityType @@ -271,35 +270,31 @@ module function plastic_dotState(subdt,co,ip,el,ph,me) result(broken) subdt !< timestep real(pReal), dimension(3,3) :: & Mp - integer :: & - instance logical :: broken - instance = phase_plasticityInstance(ph) - Mp = matmul(matmul(transpose(constitutive_mech_Fi(ph)%data(1:3,1:3,me)),& constitutive_mech_Fi(ph)%data(1:3,1:3,me)),constitutive_mech_S(ph)%data(1:3,1:3,me)) plasticityType: select case (phase_plasticity(ph)) case (PLASTICITY_ISOTROPIC_ID) plasticityType - call isotropic_dotState(Mp,instance,me) + call isotropic_dotState(Mp,ph,me) case (PLASTICITY_PHENOPOWERLAW_ID) plasticityType - call phenopowerlaw_dotState(Mp,instance,me) + call phenopowerlaw_dotState(Mp,ph,me) case (PLASTICITY_KINEHARDENING_ID) plasticityType - call plastic_kinehardening_dotState(Mp,instance,me) + call plastic_kinehardening_dotState(Mp,ph,me) case (PLASTICITY_DISLOTWIN_ID) plasticityType - call dislotwin_dotState(Mp,thermal_T(ph,me),instance,me) + call dislotwin_dotState(Mp,thermal_T(ph,me),ph,me) case (PLASTICITY_DISLOTUNGSTEN_ID) plasticityType - call dislotungsten_dotState(Mp,thermal_T(ph,me),instance,me) + call dislotungsten_dotState(Mp,thermal_T(ph,me),ph,me) case (PLASTICITY_NONLOCAL_ID) plasticityType - call nonlocal_dotState(Mp,thermal_T(ph,me),subdt,instance,me,ip,el) + call nonlocal_dotState(Mp,thermal_T(ph,me),subdt,ph,me,ip,el) end select plasticityType broken = any(IEEE_is_NaN(plasticState(ph)%dotState(:,me))) diff --git a/src/phase_mechanics_plastic_dislotungsten.f90 b/src/phase_mechanics_plastic_dislotungsten.f90 index c57669edb..ee0dfa9da 100644 --- a/src/phase_mechanics_plastic_dislotungsten.f90 +++ b/src/phase_mechanics_plastic_dislotungsten.f90 @@ -273,7 +273,7 @@ end function plastic_dislotungsten_init !> @brief Calculate plastic velocity gradient and its tangent. !-------------------------------------------------------------------------------------------------- pure module subroutine dislotungsten_LpAndItsTangent(Lp,dLp_dMp, & - Mp,T,instance,me) + Mp,T,ph,me) real(pReal), dimension(3,3), intent(out) :: & Lp !< plastic velocity gradient real(pReal), dimension(3,3,3,3), intent(out) :: & @@ -284,21 +284,21 @@ pure module subroutine dislotungsten_LpAndItsTangent(Lp,dLp_dMp, & real(pReal), intent(in) :: & T !< temperature integer, intent(in) :: & - instance, & + ph, & me integer :: & i,k,l,m,n - real(pReal), dimension(param(instance)%sum_N_sl) :: & + real(pReal), dimension(param(phase_plasticityInstance(ph))%sum_N_sl) :: & dot_gamma_pos,dot_gamma_neg, & ddot_gamma_dtau_pos,ddot_gamma_dtau_neg Lp = 0.0_pReal dLp_dMp = 0.0_pReal - associate(prm => param(instance)) + associate(prm => param(phase_plasticityInstance(ph))) - call kinetics(Mp,T,instance,me,dot_gamma_pos,dot_gamma_neg,ddot_gamma_dtau_pos,ddot_gamma_dtau_neg) + call kinetics(Mp,T,phase_plasticityInstance(ph),me,dot_gamma_pos,dot_gamma_neg,ddot_gamma_dtau_pos,ddot_gamma_dtau_neg) do i = 1, prm%sum_N_sl Lp = Lp + (dot_gamma_pos(i)+dot_gamma_neg(i))*prm%P_sl(1:3,1:3,i) forall (k=1:3,l=1:3,m=1:3,n=1:3) & @@ -315,19 +315,19 @@ end subroutine dislotungsten_LpAndItsTangent !-------------------------------------------------------------------------------------------------- !> @brief Calculate the rate of change of microstructure. !-------------------------------------------------------------------------------------------------- -module subroutine dislotungsten_dotState(Mp,T,instance,me) +module subroutine dislotungsten_dotState(Mp,T,ph,me) real(pReal), dimension(3,3), intent(in) :: & Mp !< Mandel stress real(pReal), intent(in) :: & T !< temperature integer, intent(in) :: & - instance, & + ph, & me real(pReal) :: & VacancyDiffusion - real(pReal), dimension(param(instance)%sum_N_sl) :: & + real(pReal), dimension(param(phase_plasticityInstance(ph))%sum_N_sl) :: & gdot_pos, gdot_neg,& tau_pos,& tau_neg, & @@ -336,9 +336,10 @@ module subroutine dislotungsten_dotState(Mp,T,instance,me) dot_rho_dip_climb, & dip_distance - associate(prm => param(instance), stt => state(instance),dot => dotState(instance), dst => dependentState(instance)) + associate(prm => param(phase_plasticityInstance(ph)), stt => state(phase_plasticityInstance(ph)),& + dot => dotState(phase_plasticityInstance(ph)), dst => dependentState(phase_plasticityInstance(ph))) - call kinetics(Mp,T,instance,me,& + call kinetics(Mp,T,phase_plasticityInstance(ph),me,& gdot_pos,gdot_neg, & tau_pos_out = tau_pos,tau_neg_out = tau_neg) diff --git a/src/phase_mechanics_plastic_dislotwin.f90 b/src/phase_mechanics_plastic_dislotwin.f90 index 4889b7698..1ad1957ec 100644 --- a/src/phase_mechanics_plastic_dislotwin.f90 +++ b/src/phase_mechanics_plastic_dislotwin.f90 @@ -521,12 +521,12 @@ end function plastic_dislotwin_homogenizedC !-------------------------------------------------------------------------------------------------- !> @brief Calculate plastic velocity gradient and its tangent. !-------------------------------------------------------------------------------------------------- -module subroutine dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,instance,me) +module subroutine dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,ph,me) real(pReal), dimension(3,3), intent(out) :: Lp real(pReal), dimension(3,3,3,3), intent(out) :: dLp_dMp real(pReal), dimension(3,3), intent(in) :: Mp - integer, intent(in) :: instance,me + integer, intent(in) :: ph,me real(pReal), intent(in) :: T integer :: i,k,l,m,n @@ -535,11 +535,11 @@ module subroutine dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,instance,me) BoltzmannRatio, & ddot_gamma_dtau, & tau - real(pReal), dimension(param(instance)%sum_N_sl) :: & + real(pReal), dimension(param(phase_plasticityInstance(ph))%sum_N_sl) :: & dot_gamma_sl,ddot_gamma_dtau_slip - real(pReal), dimension(param(instance)%sum_N_tw) :: & + real(pReal), dimension(param(phase_plasticityInstance(ph))%sum_N_tw) :: & dot_gamma_twin,ddot_gamma_dtau_twin - real(pReal), dimension(param(instance)%sum_N_tr) :: & + real(pReal), dimension(param(phase_plasticityInstance(ph))%sum_N_tr) :: & dot_gamma_tr,ddot_gamma_dtau_trans real(pReal):: dot_gamma_sb real(pReal), dimension(3,3) :: eigVectors, P_sb @@ -564,7 +564,7 @@ module subroutine dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,instance,me) 0, 1, 1 & ],pReal),[ 3,6]) - associate(prm => param(instance), stt => state(instance)) + associate(prm => param(phase_plasticityInstance(ph)), stt => state(phase_plasticityInstance(ph))) f_unrotated = 1.0_pReal & - sum(stt%f_tw(1:prm%sum_N_tw,me)) & @@ -573,7 +573,7 @@ module subroutine dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,instance,me) Lp = 0.0_pReal dLp_dMp = 0.0_pReal - call kinetics_slip(Mp,T,instance,me,dot_gamma_sl,ddot_gamma_dtau_slip) + call kinetics_slip(Mp,T,phase_plasticityInstance(ph),me,dot_gamma_sl,ddot_gamma_dtau_slip) slipContribution: do i = 1, prm%sum_N_sl Lp = Lp + dot_gamma_sl(i)*prm%P_sl(1:3,1:3,i) forall (k=1:3,l=1:3,m=1:3,n=1:3) & @@ -581,7 +581,7 @@ module subroutine dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,instance,me) + ddot_gamma_dtau_slip(i) * prm%P_sl(k,l,i) * prm%P_sl(m,n,i) enddo slipContribution - call kinetics_twin(Mp,T,dot_gamma_sl,instance,me,dot_gamma_twin,ddot_gamma_dtau_twin) + call kinetics_twin(Mp,T,dot_gamma_sl,phase_plasticityInstance(ph),me,dot_gamma_twin,ddot_gamma_dtau_twin) twinContibution: do i = 1, prm%sum_N_tw Lp = Lp + dot_gamma_twin(i)*prm%P_tw(1:3,1:3,i) forall (k=1:3,l=1:3,m=1:3,n=1:3) & @@ -589,7 +589,7 @@ module subroutine dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,instance,me) + ddot_gamma_dtau_twin(i)* prm%P_tw(k,l,i)*prm%P_tw(m,n,i) enddo twinContibution - call kinetics_trans(Mp,T,dot_gamma_sl,instance,me,dot_gamma_tr,ddot_gamma_dtau_trans) + call kinetics_trans(Mp,T,dot_gamma_sl,phase_plasticityInstance(ph),me,dot_gamma_tr,ddot_gamma_dtau_trans) transContibution: do i = 1, prm%sum_N_tr Lp = Lp + dot_gamma_tr(i)*prm%P_tr(1:3,1:3,i) forall (k=1:3,l=1:3,m=1:3,n=1:3) & @@ -634,14 +634,14 @@ end subroutine dislotwin_LpAndItsTangent !-------------------------------------------------------------------------------------------------- !> @brief Calculate the rate of change of microstructure. !-------------------------------------------------------------------------------------------------- -module subroutine dislotwin_dotState(Mp,T,instance,me) +module subroutine dislotwin_dotState(Mp,T,ph,me) real(pReal), dimension(3,3), intent(in):: & Mp !< Mandel stress real(pReal), intent(in) :: & T !< temperature at integration point integer, intent(in) :: & - instance, & + ph, & me integer :: i @@ -653,24 +653,24 @@ module subroutine dislotwin_dotState(Mp,T,instance,me) tau, & sigma_cl, & !< climb stress b_d !< ratio of Burgers vector to stacking fault width - real(pReal), dimension(param(instance)%sum_N_sl) :: & + real(pReal), dimension(param(phase_plasticityInstance(ph))%sum_N_sl) :: & dot_rho_dip_formation, & dot_rho_dip_climb, & rho_dip_distance_min, & dot_gamma_sl - real(pReal), dimension(param(instance)%sum_N_tw) :: & + real(pReal), dimension(param(phase_plasticityInstance(ph))%sum_N_tw) :: & dot_gamma_twin - real(pReal), dimension(param(instance)%sum_N_tr) :: & + real(pReal), dimension(param(phase_plasticityInstance(ph))%sum_N_tr) :: & dot_gamma_tr - associate(prm => param(instance), stt => state(instance), & - dot => dotState(instance), dst => dependentState(instance)) + associate(prm => param(phase_plasticityInstance(ph)), stt => state(phase_plasticityInstance(ph)), & + dot => dotState(phase_plasticityInstance(ph)), dst => dependentState(phase_plasticityInstance(ph))) f_unrotated = 1.0_pReal & - sum(stt%f_tw(1:prm%sum_N_tw,me)) & - sum(stt%f_tr(1:prm%sum_N_tr,me)) - call kinetics_slip(Mp,T,instance,me,dot_gamma_sl) + call kinetics_slip(Mp,T,phase_plasticityInstance(ph),me,dot_gamma_sl) dot%gamma_sl(:,me) = abs(dot_gamma_sl) rho_dip_distance_min = prm%D_a*prm%b_sl @@ -721,10 +721,10 @@ module subroutine dislotwin_dotState(Mp,T,instance,me) - 2.0_pReal*rho_dip_distance_min/prm%b_sl * stt%rho_dip(:,me)*abs(dot_gamma_sl) & - dot_rho_dip_climb - call kinetics_twin(Mp,T,dot_gamma_sl,instance,me,dot_gamma_twin) + call kinetics_twin(Mp,T,dot_gamma_sl,phase_plasticityInstance(ph),me,dot_gamma_twin) dot%f_tw(:,me) = f_unrotated*dot_gamma_twin/prm%gamma_char - call kinetics_trans(Mp,T,dot_gamma_sl,instance,me,dot_gamma_tr) + call kinetics_trans(Mp,T,dot_gamma_sl,phase_plasticityInstance(ph),me,dot_gamma_tr) dot%f_tr(:,me) = f_unrotated*dot_gamma_tr end associate diff --git a/src/phase_mechanics_plastic_isotropic.f90 b/src/phase_mechanics_plastic_isotropic.f90 index 080a56cfd..e7b72fcd2 100644 --- a/src/phase_mechanics_plastic_isotropic.f90 +++ b/src/phase_mechanics_plastic_isotropic.f90 @@ -168,7 +168,7 @@ end function plastic_isotropic_init !-------------------------------------------------------------------------------------------------- !> @brief Calculate plastic velocity gradient and its tangent. !-------------------------------------------------------------------------------------------------- -module subroutine isotropic_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,me) +module subroutine isotropic_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,me) real(pReal), dimension(3,3), intent(out) :: & Lp !< plastic velocity gradient @@ -178,7 +178,7 @@ module subroutine isotropic_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,me) real(pReal), dimension(3,3), intent(in) :: & Mp !< Mandel stress integer, intent(in) :: & - instance, & + ph, & me real(pReal), dimension(3,3) :: & @@ -190,7 +190,7 @@ module subroutine isotropic_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,me) integer :: & k, l, m, n - associate(prm => param(instance), stt => state(instance)) + associate(prm => param(phase_plasticityInstance(ph)), stt => state(phase_plasticityInstance(ph))) Mp_dev = math_deviatoric33(Mp) squarenorm_Mp_dev = math_tensordot(Mp_dev,Mp_dev) @@ -262,12 +262,12 @@ module subroutine plastic_isotropic_LiAndItsTangent(Li,dLi_dMi,Mi,instance,me) !-------------------------------------------------------------------------------------------------- !> @brief Calculate the rate of change of microstructure. !-------------------------------------------------------------------------------------------------- -module subroutine isotropic_dotState(Mp,instance,me) +module subroutine isotropic_dotState(Mp,ph,me) real(pReal), dimension(3,3), intent(in) :: & Mp !< Mandel stress integer, intent(in) :: & - instance, & + ph, & me real(pReal) :: & @@ -275,7 +275,8 @@ module subroutine isotropic_dotState(Mp,instance,me) xi_inf_star, & !< saturation xi norm_Mp !< norm of the (deviatoric) Mandel stress - associate(prm => param(instance), stt => state(instance), dot => dotState(instance)) + associate(prm => param(phase_plasticityInstance(ph)), stt => state(phase_plasticityInstance(ph)), & + dot => dotState(phase_plasticityInstance(ph))) if (prm%dilatation) then norm_Mp = sqrt(math_tensordot(Mp,Mp)) diff --git a/src/phase_mechanics_plastic_kinehardening.f90 b/src/phase_mechanics_plastic_kinehardening.f90 index de80a646d..034e5fd32 100644 --- a/src/phase_mechanics_plastic_kinehardening.f90 +++ b/src/phase_mechanics_plastic_kinehardening.f90 @@ -240,7 +240,7 @@ end function plastic_kinehardening_init !-------------------------------------------------------------------------------------------------- !> @brief Calculate plastic velocity gradient and its tangent. !-------------------------------------------------------------------------------------------------- -pure module subroutine kinehardening_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,me) +pure module subroutine kinehardening_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,me) real(pReal), dimension(3,3), intent(out) :: & Lp !< plastic velocity gradient @@ -250,21 +250,21 @@ pure module subroutine kinehardening_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,me) real(pReal), dimension(3,3), intent(in) :: & Mp !< Mandel stress integer, intent(in) :: & - instance, & + ph, & me integer :: & i,k,l,m,n - real(pReal), dimension(param(instance)%sum_N_sl) :: & + real(pReal), dimension(param(phase_plasticityInstance(ph))%sum_N_sl) :: & gdot_pos,gdot_neg, & dgdot_dtau_pos,dgdot_dtau_neg Lp = 0.0_pReal dLp_dMp = 0.0_pReal - associate(prm => param(instance)) + associate(prm => param(phase_plasticityInstance(ph))) - call kinetics(Mp,instance,me,gdot_pos,gdot_neg,dgdot_dtau_pos,dgdot_dtau_neg) + call kinetics(Mp,phase_plasticityInstance(ph),me,gdot_pos,gdot_neg,dgdot_dtau_pos,dgdot_dtau_neg) do i = 1, prm%sum_N_sl Lp = Lp + (gdot_pos(i)+gdot_neg(i))*prm%P(1:3,1:3,i) @@ -282,23 +282,24 @@ end subroutine kinehardening_LpAndItsTangent !-------------------------------------------------------------------------------------------------- !> @brief Calculate the rate of change of microstructure. !-------------------------------------------------------------------------------------------------- -module subroutine plastic_kinehardening_dotState(Mp,instance,me) +module subroutine plastic_kinehardening_dotState(Mp,ph,me) real(pReal), dimension(3,3), intent(in) :: & Mp !< Mandel stress integer, intent(in) :: & - instance, & + ph, & me real(pReal) :: & sumGamma - real(pReal), dimension(param(instance)%sum_N_sl) :: & + real(pReal), dimension(param(phase_plasticityInstance(ph))%sum_N_sl) :: & gdot_pos,gdot_neg - associate(prm => param(instance), stt => state(instance), dot => dotState(instance)) + associate(prm => param(phase_plasticityInstance(ph)), stt => state(phase_plasticityInstance(ph)),& + dot => dotState(phase_plasticityInstance(ph))) - call kinetics(Mp,instance,me,gdot_pos,gdot_neg) + call kinetics(Mp,phase_plasticityInstance(ph),me,gdot_pos,gdot_neg) dot%accshear(:,me) = abs(gdot_pos+gdot_neg) sumGamma = sum(stt%accshear(:,me)) diff --git a/src/phase_mechanics_plastic_nonlocal.f90 b/src/phase_mechanics_plastic_nonlocal.f90 index 884bc6bc0..d0ed01cf0 100644 --- a/src/phase_mechanics_plastic_nonlocal.f90 +++ b/src/phase_mechanics_plastic_nonlocal.f90 @@ -758,13 +758,13 @@ end subroutine nonlocal_dependentState !> @brief calculates plastic velocity gradient and its tangent !-------------------------------------------------------------------------------------------------- module subroutine nonlocal_LpAndItsTangent(Lp,dLp_dMp, & - Mp,Temperature,instance,me,ip,el) + Mp,Temperature,ph,me,ip,el) real(pReal), dimension(3,3), intent(out) :: & Lp !< plastic velocity gradient real(pReal), dimension(3,3,3,3), intent(out) :: & dLp_dMp integer, intent(in) :: & - instance, & + ph, & me, & ip, & !< current integration point el !< current element number @@ -782,24 +782,25 @@ module subroutine nonlocal_LpAndItsTangent(Lp,dLp_dMp, & l, & t, & !< dislocation type s !< index of my current slip system - real(pReal), dimension(param(instance)%sum_N_sl,8) :: & + real(pReal), dimension(param(phase_plasticityInstance(ph))%sum_N_sl,8) :: & rhoSgl !< single dislocation densities (including blocked) - real(pReal), dimension(param(instance)%sum_N_sl,10) :: & + real(pReal), dimension(param(phase_plasticityInstance(ph))%sum_N_sl,10) :: & rho - real(pReal), dimension(param(instance)%sum_N_sl,4) :: & + real(pReal), dimension(param(phase_plasticityInstance(ph))%sum_N_sl,4) :: & v, & !< velocity tauNS, & !< resolved shear stress including non Schmid and backstress terms dv_dtau, & !< velocity derivative with respect to the shear stress dv_dtauNS !< velocity derivative with respect to the shear stress - real(pReal), dimension(param(instance)%sum_N_sl) :: & + real(pReal), dimension(param(phase_plasticityInstance(ph))%sum_N_sl) :: & tau, & !< resolved shear stress including backstress terms gdotTotal !< shear rate - associate(prm => param(instance),dst=>microstructure(instance),stt=>state(instance)) + associate(prm => param(phase_plasticityInstance(ph)),dst=>microstructure(phase_plasticityInstance(ph)),& + stt=>state(phase_plasticityInstance(ph))) ns = prm%sum_N_sl !*** shortcut to state variables - rho = getRho(instance,me,ip,el) + rho = getRho(phase_plasticityInstance(ph),me,ip,el) rhoSgl = rho(:,sgl) do s = 1,ns @@ -819,7 +820,7 @@ module subroutine nonlocal_LpAndItsTangent(Lp,dLp_dMp, & ! edges call kinetics(v(:,1), dv_dtau(:,1), dv_dtauNS(:,1), & - tau, tauNS(:,1), dst%tau_pass(:,me),1,Temperature, instance) + tau, tauNS(:,1), dst%tau_pass(:,me),1,Temperature, phase_plasticityInstance(ph)) v(:,2) = v(:,1) dv_dtau(:,2) = dv_dtau(:,1) dv_dtauNS(:,2) = dv_dtauNS(:,1) @@ -832,7 +833,7 @@ module subroutine nonlocal_LpAndItsTangent(Lp,dLp_dMp, & else do t = 3,4 call kinetics(v(:,t), dv_dtau(:,t), dv_dtauNS(:,t), & - tau, tauNS(:,t), dst%tau_pass(:,me),2,Temperature, instance) + tau, tauNS(:,t), dst%tau_pass(:,me),2,Temperature, phase_plasticityInstance(ph)) enddo endif @@ -973,7 +974,7 @@ end subroutine plastic_nonlocal_deltaState !> @brief calculates the rate of change of microstructure !--------------------------------------------------------------------------------------------------- module subroutine nonlocal_dotState(Mp, Temperature,timestep, & - instance,me,ip,el) + ph,me,ip,el) real(pReal), dimension(3,3), intent(in) :: & Mp !< MandelStress @@ -981,7 +982,7 @@ module subroutine nonlocal_dotState(Mp, Temperature,timestep, & Temperature, & !< temperature timestep !< substepped crystallite time increment integer, intent(in) :: & - instance, & + ph, & me, & ip, & !< current integration point el !< current element number @@ -992,7 +993,7 @@ module subroutine nonlocal_dotState(Mp, Temperature,timestep, & c, & !< character of dislocation t, & !< type of dislocation s !< index of my current slip system - real(pReal), dimension(param(instance)%sum_N_sl,10) :: & + real(pReal), dimension(param(phase_plasticityInstance(ph))%sum_N_sl,10) :: & rho, & rho0, & !< dislocation density at beginning of time step rhoDot, & !< density evolution @@ -1000,45 +1001,44 @@ module subroutine nonlocal_dotState(Mp, Temperature,timestep, & rhoDotSingle2DipoleGlide, & !< density evolution by dipole formation (by glide) rhoDotAthermalAnnihilation, & !< density evolution by athermal annihilation rhoDotThermalAnnihilation !< density evolution by thermal annihilation - real(pReal), dimension(param(instance)%sum_N_sl,8) :: & + real(pReal), dimension(param(phase_plasticityInstance(ph))%sum_N_sl,8) :: & rhoSgl, & !< current single dislocation densities (positive/negative screw and edge without dipoles) my_rhoSgl0 !< single dislocation densities of central ip (positive/negative screw and edge without dipoles) - real(pReal), dimension(param(instance)%sum_N_sl,4) :: & + real(pReal), dimension(param(phase_plasticityInstance(ph))%sum_N_sl,4) :: & v, & !< current dislocation glide velocity v0, & gdot !< shear rates - real(pReal), dimension(param(instance)%sum_N_sl) :: & + real(pReal), dimension(param(phase_plasticityInstance(ph))%sum_N_sl) :: & tau, & !< current resolved shear stress vClimb !< climb velocity of edge dipoles - real(pReal), dimension(param(instance)%sum_N_sl,2) :: & + real(pReal), dimension(param(phase_plasticityInstance(ph))%sum_N_sl,2) :: & rhoDip, & !< current dipole dislocation densities (screw and edge dipoles) dLower, & !< minimum stable dipole distance for edges and screws dUpper !< current maximum stable dipole distance for edges and screws real(pReal) :: & selfDiffusion !< self diffusion - ph = material_phaseAt(1,el) if (timestep <= 0.0_pReal) then plasticState(ph)%dotState = 0.0_pReal return endif - associate(prm => param(instance), & - dst => microstructure(instance), & - dot => dotState(instance), & - stt => state(instance)) + associate(prm => param(phase_plasticityInstance(ph)), & + dst => microstructure(phase_plasticityInstance(ph)), & + dot => dotState(phase_plasticityInstance(ph)), & + stt => state(phase_plasticityInstance(ph))) ns = prm%sum_N_sl tau = 0.0_pReal gdot = 0.0_pReal - rho = getRho(instance,me,ip,el) + rho = getRho(phase_plasticityInstance(ph),me,ip,el) rhoSgl = rho(:,sgl) rhoDip = rho(:,dip) - rho0 = getRho0(instance,me,ip,el) + rho0 = getRho0(phase_plasticityInstance(ph),me,ip,el) my_rhoSgl0 = rho0(:,sgl) - forall (s = 1:ns, t = 1:4) v(s,t) = plasticState(ph)%state(iV(s,t,instance),me) + forall (s = 1:ns, t = 1:4) v(s,t) = plasticState(ph)%state(iV(s,t,phase_plasticityInstance(ph)),me) gdot = rhoSgl(:,1:4) * v * spread(prm%b_sl,2,4) #ifdef DEBUG @@ -1087,7 +1087,7 @@ module subroutine nonlocal_dotState(Mp, Temperature,timestep, & * sqrt(stt%rho_forest(:,me)) / prm%i_sl / prm%b_sl, 2, 4) endif isBCC - forall (s = 1:ns, t = 1:4) v0(s,t) = plasticState(ph)%state0(iV(s,t,instance),me) + forall (s = 1:ns, t = 1:4) v0(s,t) = plasticState(ph)%state0(iV(s,t,phase_plasticityInstance(ph)),me) !**************************************************************************** @@ -1143,7 +1143,7 @@ module subroutine nonlocal_dotState(Mp, Temperature,timestep, & - rhoDip(s,1) / timestep - rhoDotAthermalAnnihilation(s,9) & - rhoDotSingle2DipoleGlide(s,9)) ! make sure that we do not annihilate more dipoles than we have - rhoDot = rhoDotFlux(timestep, instance,me,ip,el) & + rhoDot = rhoDotFlux(timestep, phase_plasticityInstance(ph),me,ip,el) & + rhoDotMultiplication & + rhoDotSingle2DipoleGlide & + rhoDotAthermalAnnihilation & diff --git a/src/phase_mechanics_plastic_phenopowerlaw.f90 b/src/phase_mechanics_plastic_phenopowerlaw.f90 index a9538aeec..1d91a97fb 100644 --- a/src/phase_mechanics_plastic_phenopowerlaw.f90 +++ b/src/phase_mechanics_plastic_phenopowerlaw.f90 @@ -285,7 +285,7 @@ end function plastic_phenopowerlaw_init !> @details asummes that deformation by dislocation glide affects twinned and untwinned volume ! equally (Taylor assumption). Twinning happens only in untwinned volume !-------------------------------------------------------------------------------------------------- -pure module subroutine phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,me) +pure module subroutine phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,me) real(pReal), dimension(3,3), intent(out) :: & Lp !< plastic velocity gradient @@ -295,23 +295,23 @@ pure module subroutine phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,me) real(pReal), dimension(3,3), intent(in) :: & Mp !< Mandel stress integer, intent(in) :: & - instance, & + ph, & me integer :: & i,k,l,m,n - real(pReal), dimension(param(instance)%sum_N_sl) :: & + real(pReal), dimension(param(phase_plasticityInstance(ph))%sum_N_sl) :: & gdot_slip_pos,gdot_slip_neg, & dgdot_dtauslip_pos,dgdot_dtauslip_neg - real(pReal), dimension(param(instance)%sum_N_tw) :: & + real(pReal), dimension(param(phase_plasticityInstance(ph))%sum_N_tw) :: & gdot_twin,dgdot_dtautwin Lp = 0.0_pReal dLp_dMp = 0.0_pReal - associate(prm => param(instance)) + associate(prm => param(phase_plasticityInstance(ph))) - call kinetics_slip(Mp,instance,me,gdot_slip_pos,gdot_slip_neg,dgdot_dtauslip_pos,dgdot_dtauslip_neg) + call kinetics_slip(Mp,phase_plasticityInstance(ph),me,gdot_slip_pos,gdot_slip_neg,dgdot_dtauslip_pos,dgdot_dtauslip_neg) slipSystems: do i = 1, prm%sum_N_sl Lp = Lp + (gdot_slip_pos(i)+gdot_slip_neg(i))*prm%P_sl(1:3,1:3,i) forall (k=1:3,l=1:3,m=1:3,n=1:3) & @@ -320,7 +320,7 @@ pure module subroutine phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,me) + dgdot_dtauslip_neg(i) * prm%P_sl(k,l,i) * prm%nonSchmid_neg(m,n,i) enddo slipSystems - call kinetics_twin(Mp,instance,me,gdot_twin,dgdot_dtautwin) + call kinetics_twin(Mp,phase_plasticityInstance(ph),me,gdot_twin,dgdot_dtautwin) twinSystems: do i = 1, prm%sum_N_tw Lp = Lp + gdot_twin(i)*prm%P_tw(1:3,1:3,i) forall (k=1:3,l=1:3,m=1:3,n=1:3) & @@ -336,23 +336,24 @@ end subroutine phenopowerlaw_LpAndItsTangent !-------------------------------------------------------------------------------------------------- !> @brief Calculate the rate of change of microstructure. !-------------------------------------------------------------------------------------------------- -module subroutine phenopowerlaw_dotState(Mp,instance,me) +module subroutine phenopowerlaw_dotState(Mp,ph,me) real(pReal), dimension(3,3), intent(in) :: & Mp !< Mandel stress integer, intent(in) :: & - instance, & + ph, & me real(pReal) :: & c_SlipSlip,c_TwinSlip,c_TwinTwin, & xi_slip_sat_offset,& sumGamma,sumF - real(pReal), dimension(param(instance)%sum_N_sl) :: & + real(pReal), dimension(param(phase_plasticityInstance(ph))%sum_N_sl) :: & left_SlipSlip,right_SlipSlip, & gdot_slip_pos,gdot_slip_neg - associate(prm => param(instance), stt => state(instance), dot => dotState(instance)) + associate(prm => param(phase_plasticityInstance(ph)), stt => state(phase_plasticityInstance(ph)), & + dot => dotState(phase_plasticityInstance(ph))) sumGamma = sum(stt%gamma_slip(:,me)) sumF = sum(stt%gamma_twin(:,me)/prm%gamma_tw_char) @@ -372,9 +373,9 @@ module subroutine phenopowerlaw_dotState(Mp,instance,me) !-------------------------------------------------------------------------------------------------- ! shear rates - call kinetics_slip(Mp,instance,me,gdot_slip_pos,gdot_slip_neg) + call kinetics_slip(Mp,phase_plasticityInstance(ph),me,gdot_slip_pos,gdot_slip_neg) dot%gamma_slip(:,me) = abs(gdot_slip_pos+gdot_slip_neg) - call kinetics_twin(Mp,instance,me,dot%gamma_twin(:,me)) + call kinetics_twin(Mp,phase_plasticityInstance(ph),me,dot%gamma_twin(:,me)) !-------------------------------------------------------------------------------------------------- ! hardening