diff --git a/CMakeLists.txt b/CMakeLists.txt index 198528b59..533b14e85 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -4,7 +4,7 @@ include (FindPkgConfig REQUIRED) # Dummy project to determine compiler names and version project (Prerequisites LANGUAGES) set(ENV{PKG_CONFIG_PATH} "$ENV{PETSC_DIR}/$ENV{PETSC_ARCH}/lib/pkgconfig") -pkg_check_modules (PETSC REQUIRED PETSc>=3.12.0 PETSc<3.15.0) +pkg_check_modules (PETSC REQUIRED PETSc>=3.12.0 PETSc<3.16.0) pkg_get_variable (CMAKE_Fortran_COMPILER PETSc fcompiler) pkg_get_variable (CMAKE_C_COMPILER PETSc ccompiler) diff --git a/src/DAMASK_interface.f90 b/src/DAMASK_interface.f90 index b7cb1d80f..e01c050fa 100644 --- a/src/DAMASK_interface.f90 +++ b/src/DAMASK_interface.f90 @@ -51,9 +51,7 @@ subroutine DAMASK_interface_init #include #if PETSC_VERSION_MAJOR!=3 || PETSC_VERSION_MINORPETSC_MINOR_MAX -=================================================================================================== --- WRONG PETSc VERSION --- WRONG PETSc VERSION --- WRONG PETSc VERSION --- WRONG PETSc VERSION -- -=================================================================================================== +-- UNSUPPORTED PETSc VERSION --- UNSUPPORTED PETSc VERSION --- UNSUPPORTED PETSc VERSION --- #endif character(len=pPathLen*3+pStringLen) :: & diff --git a/src/phase.f90 b/src/phase.f90 index b65bf334f..93743d64a 100644 --- a/src/phase.f90 +++ b/src/phase.f90 @@ -161,13 +161,13 @@ module phase real(pReal), dimension(3,3) :: P end function phase_P - module function thermal_T(ph,me) result(T) - integer, intent(in) :: ph,me + module function thermal_T(ph,en) result(T) + integer, intent(in) :: ph,en real(pReal) :: T end function thermal_T - module function thermal_dot_T(ph,me) result(dot_T) - integer, intent(in) :: ph,me + module function thermal_dot_T(ph,en) result(dot_T) + integer, intent(in) :: ph,en real(pReal) :: dot_T end function thermal_dot_T @@ -216,10 +216,10 @@ module phase ! == cleaned:end =================================================================================== - module function thermal_stress(Delta_t,ph,me) result(converged_) + module function thermal_stress(Delta_t,ph,en) result(converged_) real(pReal), intent(in) :: Delta_t - integer, intent(in) :: ph, me + integer, intent(in) :: ph, en logical :: converged_ end function thermal_stress @@ -253,8 +253,8 @@ module phase f end function phase_f_phi - module function phase_f_T(ph,me) result(f) - integer, intent(in) :: ph, me + module function phase_f_T(ph,en) result(f) + integer, intent(in) :: ph, en real(pReal) :: f end function phase_f_T diff --git a/src/phase_mechanical_plastic.f90 b/src/phase_mechanical_plastic.f90 index 3ec275795..33dfd681f 100644 --- a/src/phase_mechanical_plastic.f90 +++ b/src/phase_mechanical_plastic.f90 @@ -37,7 +37,7 @@ submodule(phase:mechanical) plastic myPlasticity end function plastic_nonlocal_init - module subroutine isotropic_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,me) + module subroutine isotropic_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,en) real(pReal), dimension(3,3), intent(out) :: & Lp real(pReal), dimension(3,3,3,3), intent(out) :: & @@ -46,10 +46,10 @@ submodule(phase:mechanical) plastic Mp integer, intent(in) :: & ph, & - me + en end subroutine isotropic_LpAndItsTangent - pure module subroutine phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,me) + pure module subroutine phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,en) real(pReal), dimension(3,3), intent(out) :: & Lp real(pReal), dimension(3,3,3,3), intent(out) :: & @@ -58,10 +58,10 @@ submodule(phase:mechanical) plastic Mp integer, intent(in) :: & ph, & - me + en end subroutine phenopowerlaw_LpAndItsTangent - pure module subroutine kinehardening_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,me) + pure module subroutine kinehardening_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,en) real(pReal), dimension(3,3), intent(out) :: & Lp real(pReal), dimension(3,3,3,3), intent(out) :: & @@ -70,10 +70,10 @@ submodule(phase:mechanical) plastic Mp integer, intent(in) :: & ph, & - me + en end subroutine kinehardening_LpAndItsTangent - module subroutine dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,ph,me) + module subroutine dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,ph,en) real(pReal), dimension(3,3), intent(out) :: & Lp real(pReal), dimension(3,3,3,3), intent(out) :: & @@ -85,10 +85,10 @@ submodule(phase:mechanical) plastic T integer, intent(in) :: & ph, & - me + en end subroutine dislotwin_LpAndItsTangent - pure module subroutine dislotungsten_LpAndItsTangent(Lp,dLp_dMp,Mp,T,ph,me) + pure module subroutine dislotungsten_LpAndItsTangent(Lp,dLp_dMp,Mp,T,ph,en) real(pReal), dimension(3,3), intent(out) :: & Lp real(pReal), dimension(3,3,3,3), intent(out) :: & @@ -100,11 +100,11 @@ submodule(phase:mechanical) plastic T integer, intent(in) :: & ph, & - me + en end subroutine dislotungsten_LpAndItsTangent module subroutine nonlocal_LpAndItsTangent(Lp,dLp_dMp, & - Mp,Temperature,ph,me) + Mp,Temperature,ph,en) real(pReal), dimension(3,3), intent(out) :: & Lp real(pReal), dimension(3,3,3,3), intent(out) :: & @@ -116,55 +116,55 @@ submodule(phase:mechanical) plastic Temperature integer, intent(in) :: & ph, & - me + en end subroutine nonlocal_LpAndItsTangent - module subroutine isotropic_dotState(Mp,ph,me) + module subroutine isotropic_dotState(Mp,ph,en) real(pReal), dimension(3,3), intent(in) :: & Mp !< Mandel stress integer, intent(in) :: & ph, & - me + en end subroutine isotropic_dotState - module subroutine phenopowerlaw_dotState(Mp,ph,me) + module subroutine phenopowerlaw_dotState(Mp,ph,en) real(pReal), dimension(3,3), intent(in) :: & Mp !< Mandel stress integer, intent(in) :: & ph, & - me + en end subroutine phenopowerlaw_dotState - module subroutine plastic_kinehardening_dotState(Mp,ph,me) + module subroutine plastic_kinehardening_dotState(Mp,ph,en) real(pReal), dimension(3,3), intent(in) :: & Mp !< Mandel stress integer, intent(in) :: & ph, & - me + en end subroutine plastic_kinehardening_dotState - module subroutine dislotwin_dotState(Mp,T,ph,me) + module subroutine dislotwin_dotState(Mp,T,ph,en) real(pReal), dimension(3,3), intent(in) :: & Mp !< Mandel stress real(pReal), intent(in) :: & T integer, intent(in) :: & ph, & - me + en end subroutine dislotwin_dotState - module subroutine dislotungsten_dotState(Mp,T,ph,me) + module subroutine dislotungsten_dotState(Mp,T,ph,en) real(pReal), dimension(3,3), intent(in) :: & Mp !< Mandel stress real(pReal), intent(in) :: & T integer, intent(in) :: & ph, & - me + en end subroutine dislotungsten_dotState - module subroutine nonlocal_dotState(Mp,Temperature,timestep,ph,me,ip,el) + module subroutine nonlocal_dotState(Mp,Temperature,timestep,ph,en,ip,el) real(pReal), dimension(3,3), intent(in) :: & Mp !< MandelStress real(pReal), intent(in) :: & @@ -172,47 +172,47 @@ submodule(phase:mechanical) plastic timestep !< substepped crystallite time increment integer, intent(in) :: & ph, & - me, & + en, & ip, & !< current integration point el !< current element number end subroutine nonlocal_dotState - module subroutine dislotwin_dependentState(T,ph,me) + module subroutine dislotwin_dependentState(T,ph,en) integer, intent(in) :: & ph, & - me + en real(pReal), intent(in) :: & T end subroutine dislotwin_dependentState - module subroutine dislotungsten_dependentState(ph,me) + module subroutine dislotungsten_dependentState(ph,en) integer, intent(in) :: & ph, & - me + en end subroutine dislotungsten_dependentState - module subroutine nonlocal_dependentState(ph, me, ip, el) + module subroutine nonlocal_dependentState(ph, en, ip, el) integer, intent(in) :: & ph, & - me, & + en, & ip, & !< current integration point el !< current element number end subroutine nonlocal_dependentState - module subroutine plastic_kinehardening_deltaState(Mp,ph,me) + module subroutine plastic_kinehardening_deltaState(Mp,ph,en) real(pReal), dimension(3,3), intent(in) :: & Mp !< Mandel stress integer, intent(in) :: & ph, & - me + en end subroutine plastic_kinehardening_deltaState - module subroutine plastic_nonlocal_deltaState(Mp,ph,me) + module subroutine plastic_nonlocal_deltaState(Mp,ph,en) real(pReal), dimension(3,3), intent(in) :: & Mp integer, intent(in) :: & ph, & - me + en end subroutine plastic_nonlocal_deltaState end interface @@ -357,22 +357,22 @@ module subroutine plastic_dependentState(co, ip, el) integer :: & ph, & - me + en ph = material_phaseAt(co,el) - me = material_phasememberAt(co,ip,el) + en = material_phasememberAt(co,ip,el) plasticType: select case (phase_plasticity(material_phaseAt(co,el))) case (PLASTICITY_DISLOTWIN_ID) plasticType - call dislotwin_dependentState(thermal_T(ph,me),ph,me) + call dislotwin_dependentState(thermal_T(ph,en),ph,en) case (PLASTICITY_DISLOTUNGSTEN_ID) plasticType - call dislotungsten_dependentState(ph,me) + call dislotungsten_dependentState(ph,en) case (PLASTICITY_NONLOCAL_ID) plasticType - call nonlocal_dependentState(ph,me,ip,el) + call nonlocal_dependentState(ph,en,ip,el) end select plasticType diff --git a/src/phase_mechanical_plastic_dislotungsten.f90 b/src/phase_mechanical_plastic_dislotungsten.f90 index 2e7130756..bb53684e1 100644 --- a/src/phase_mechanical_plastic_dislotungsten.f90 +++ b/src/phase_mechanical_plastic_dislotungsten.f90 @@ -17,7 +17,7 @@ submodule(phase:plastic) dislotungsten D_0 = 1.0_pReal, & !< prefactor for self-diffusion coefficient Q_cl = 1.0_pReal !< activation energy for dislocation climb real(pReal), allocatable, dimension(:) :: & - b_sl, & !< magnitude me Burgers vector [m] + b_sl, & !< magnitude of Burgers vector [m] D_a, & i_sl, & !< Adj. parameter for distance between 2 forest dislocations f_at, & !< factor to calculate atomic volume @@ -271,7 +271,7 @@ end function plastic_dislotungsten_init !> @brief Calculate plastic velocity gradient and its tangent. !-------------------------------------------------------------------------------------------------- pure module subroutine dislotungsten_LpAndItsTangent(Lp,dLp_dMp, & - Mp,T,ph,me) + Mp,T,ph,en) real(pReal), dimension(3,3), intent(out) :: & Lp !< plastic velocity gradient real(pReal), dimension(3,3,3,3), intent(out) :: & @@ -283,7 +283,7 @@ pure module subroutine dislotungsten_LpAndItsTangent(Lp,dLp_dMp, & T !< temperature integer, intent(in) :: & ph, & - me + en integer :: & i,k,l,m,n @@ -296,7 +296,7 @@ pure module subroutine dislotungsten_LpAndItsTangent(Lp,dLp_dMp, & associate(prm => param(ph)) - call kinetics(Mp,T,ph,me,dot_gamma_pos,dot_gamma_neg,ddot_gamma_dtau_pos,ddot_gamma_dtau_neg) + call kinetics(Mp,T,ph,en,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) & @@ -313,7 +313,7 @@ end subroutine dislotungsten_LpAndItsTangent !-------------------------------------------------------------------------------------------------- !> @brief Calculate the rate of change of microstructure. !-------------------------------------------------------------------------------------------------- -module subroutine dislotungsten_dotState(Mp,T,ph,me) +module subroutine dislotungsten_dotState(Mp,T,ph,en) real(pReal), dimension(3,3), intent(in) :: & Mp !< Mandel stress @@ -321,7 +321,7 @@ module subroutine dislotungsten_dotState(Mp,T,ph,me) T !< temperature integer, intent(in) :: & ph, & - me + en real(pReal) :: & VacancyDiffusion @@ -337,11 +337,11 @@ module subroutine dislotungsten_dotState(Mp,T,ph,me) associate(prm => param(ph), stt => state(ph),& dot => dotState(ph), dst => dependentState(ph)) - call kinetics(Mp,T,ph,me,& + call kinetics(Mp,T,ph,en,& gdot_pos,gdot_neg, & tau_pos_out = tau_pos,tau_neg_out = tau_neg) - dot%gamma_sl(:,me) = (gdot_pos+gdot_neg) ! ToDo: needs to be abs + dot%gamma_sl(:,en) = (gdot_pos+gdot_neg) ! ToDo: needs to be abs VacancyDiffusion = prm%D_0*exp(-prm%Q_cl/(kB*T)) where(dEq0(tau_pos)) ! ToDo: use avg of pos and neg @@ -350,20 +350,20 @@ module subroutine dislotungsten_dotState(Mp,T,ph,me) else where dip_distance = math_clip(3.0_pReal*prm%mu*prm%b_sl/(16.0_pReal*PI*abs(tau_pos)), & prm%D_a, & ! lower limit - dst%Lambda_sl(:,me)) ! upper limit - dot_rho_dip_formation = merge(2.0_pReal*dip_distance* stt%rho_mob(:,me)*abs(dot%gamma_sl(:,me))/prm%b_sl, & ! ToDo: ignore region of spontaneous annihilation + dst%Lambda_sl(:,en)) ! upper limit + dot_rho_dip_formation = merge(2.0_pReal*dip_distance* stt%rho_mob(:,en)*abs(dot%gamma_sl(:,en))/prm%b_sl, & ! ToDo: ignore region of spontaneous annihilation 0.0_pReal, & prm%dipoleformation) v_cl = (3.0_pReal*prm%mu*VacancyDiffusion*prm%f_at/(2.0_pReal*pi*kB*T)) & * (1.0_pReal/(dip_distance+prm%D_a)) - dot_rho_dip_climb = (4.0_pReal*v_cl*stt%rho_dip(:,me))/(dip_distance-prm%D_a) ! ToDo: Discuss with Franz: Stress dependency? + dot_rho_dip_climb = (4.0_pReal*v_cl*stt%rho_dip(:,en))/(dip_distance-prm%D_a) ! ToDo: Discuss with Franz: Stress dependency? end where - dot%rho_mob(:,me) = abs(dot%gamma_sl(:,me))/(prm%b_sl*dst%Lambda_sl(:,me)) & ! multiplication + dot%rho_mob(:,en) = abs(dot%gamma_sl(:,en))/(prm%b_sl*dst%Lambda_sl(:,en)) & ! multiplication - dot_rho_dip_formation & - - (2.0_pReal*prm%D_a)/prm%b_sl*stt%rho_mob(:,me)*abs(dot%gamma_sl(:,me)) ! Spontaneous annihilation of 2 single edge dislocations - dot%rho_dip(:,me) = dot_rho_dip_formation & - - (2.0_pReal*prm%D_a)/prm%b_sl*stt%rho_dip(:,me)*abs(dot%gamma_sl(:,me)) & ! Spontaneous annihilation of a single edge dislocation with a dipole constituent + - (2.0_pReal*prm%D_a)/prm%b_sl*stt%rho_mob(:,en)*abs(dot%gamma_sl(:,en)) ! Spontaneous annihilation of 2 single edge dislocations + dot%rho_dip(:,en) = dot_rho_dip_formation & + - (2.0_pReal*prm%D_a)/prm%b_sl*stt%rho_dip(:,en)*abs(dot%gamma_sl(:,en)) & ! Spontaneous annihilation of a single edge dislocation with a dipole constituent - dot_rho_dip_climb end associate @@ -374,22 +374,22 @@ end subroutine dislotungsten_dotState !-------------------------------------------------------------------------------------------------- !> @brief Calculate derived quantities from state. !-------------------------------------------------------------------------------------------------- -module subroutine dislotungsten_dependentState(ph,me) +module subroutine dislotungsten_dependentState(ph,en) integer, intent(in) :: & ph, & - me + en real(pReal), dimension(param(ph)%sum_N_sl) :: & dislocationSpacing associate(prm => param(ph), stt => state(ph),dst => dependentState(ph)) - dislocationSpacing = sqrt(matmul(prm%forestProjection,stt%rho_mob(:,me)+stt%rho_dip(:,me))) - dst%threshold_stress(:,me) = prm%mu*prm%b_sl & - * sqrt(matmul(prm%h_sl_sl,stt%rho_mob(:,me)+stt%rho_dip(:,me))) + dislocationSpacing = sqrt(matmul(prm%forestProjection,stt%rho_mob(:,en)+stt%rho_dip(:,en))) + dst%threshold_stress(:,en) = prm%mu*prm%b_sl & + * sqrt(matmul(prm%h_sl_sl,stt%rho_mob(:,en)+stt%rho_dip(:,en))) - dst%Lambda_sl(:,me) = prm%D/(1.0_pReal+prm%D*dislocationSpacing/prm%i_sl) + dst%Lambda_sl(:,en) = prm%D/(1.0_pReal+prm%D*dislocationSpacing/prm%i_sl) end associate @@ -438,7 +438,7 @@ end subroutine plastic_dislotungsten_results ! NOTE: Against the common convention, the result (i.e. intent(out)) variables are the last to ! have the optional arguments at the end !-------------------------------------------------------------------------------------------------- -pure subroutine kinetics(Mp,T,ph,me, & +pure subroutine kinetics(Mp,T,ph,en, & dot_gamma_pos,dot_gamma_neg,ddot_gamma_dtau_pos,ddot_gamma_dtau_neg,tau_pos_out,tau_neg_out) real(pReal), dimension(3,3), intent(in) :: & @@ -447,7 +447,7 @@ pure subroutine kinetics(Mp,T,ph,me, & T !< temperature integer, intent(in) :: & ph, & - me + en real(pReal), intent(out), dimension(param(ph)%sum_N_sl) :: & dot_gamma_pos, & @@ -478,11 +478,11 @@ pure subroutine kinetics(Mp,T,ph,me, & if (present(tau_neg_out)) tau_neg_out = tau_neg associate(BoltzmannRatio => prm%Q_s/(kB*T), & - dot_gamma_0 => stt%rho_mob(:,me)*prm%b_sl*prm%v_0, & - effectiveLength => dst%Lambda_sl(:,me) - prm%w) + dot_gamma_0 => stt%rho_mob(:,en)*prm%b_sl*prm%v_0, & + effectiveLength => dst%Lambda_sl(:,en) - prm%w) - significantPositiveTau: where(abs(tau_pos)-dst%threshold_stress(:,me) > tol_math_check) - StressRatio = (abs(tau_pos)-dst%threshold_stress(:,me))/prm%tau_Peierls + significantPositiveTau: where(abs(tau_pos)-dst%threshold_stress(:,en) > tol_math_check) + StressRatio = (abs(tau_pos)-dst%threshold_stress(:,en))/prm%tau_Peierls StressRatio_p = StressRatio** prm%p StressRatio_pminus1 = StressRatio**(prm%p-1.0_pReal) needsGoodName = exp(-BoltzmannRatio*(1-StressRatio_p) ** prm%q) @@ -498,7 +498,7 @@ pure subroutine kinetics(Mp,T,ph,me, & end where significantPositiveTau if (present(ddot_gamma_dtau_pos)) then - significantPositiveTau2: where(abs(tau_pos)-dst%threshold_stress(:,me) > tol_math_check) + significantPositiveTau2: where(abs(tau_pos)-dst%threshold_stress(:,en) > tol_math_check) dtn = -1.0_pReal * t_n * BoltzmannRatio * prm%p * prm%q * (1.0_pReal-StressRatio_p)**(prm%q - 1.0_pReal) & * (StressRatio)**(prm%p - 1.0_pReal) / prm%tau_Peierls dtk = -1.0_pReal * t_k / tau_pos @@ -511,8 +511,8 @@ pure subroutine kinetics(Mp,T,ph,me, & end where significantPositiveTau2 endif - significantNegativeTau: where(abs(tau_neg)-dst%threshold_stress(:,me) > tol_math_check) - StressRatio = (abs(tau_neg)-dst%threshold_stress(:,me))/prm%tau_Peierls + significantNegativeTau: where(abs(tau_neg)-dst%threshold_stress(:,en) > tol_math_check) + StressRatio = (abs(tau_neg)-dst%threshold_stress(:,en))/prm%tau_Peierls StressRatio_p = StressRatio** prm%p StressRatio_pminus1 = StressRatio**(prm%p-1.0_pReal) needsGoodName = exp(-BoltzmannRatio*(1-StressRatio_p) ** prm%q) @@ -528,7 +528,7 @@ pure subroutine kinetics(Mp,T,ph,me, & end where significantNegativeTau if (present(ddot_gamma_dtau_neg)) then - significantNegativeTau2: where(abs(tau_neg)-dst%threshold_stress(:,me) > tol_math_check) + significantNegativeTau2: where(abs(tau_neg)-dst%threshold_stress(:,en) > tol_math_check) dtn = -1.0_pReal * t_n * BoltzmannRatio * prm%p * prm%q * (1.0_pReal-StressRatio_p)**(prm%q - 1.0_pReal) & * (StressRatio)**(prm%p - 1.0_pReal) / prm%tau_Peierls dtk = -1.0_pReal * t_k / tau_neg diff --git a/src/phase_mechanical_plastic_dislotwin.f90 b/src/phase_mechanical_plastic_dislotwin.f90 index b04f28f0e..151250f14 100644 --- a/src/phase_mechanical_plastic_dislotwin.f90 +++ b/src/phase_mechanical_plastic_dislotwin.f90 @@ -519,12 +519,12 @@ end function plastic_dislotwin_homogenizedC !-------------------------------------------------------------------------------------------------- !> @brief Calculate plastic velocity gradient and its tangent. !-------------------------------------------------------------------------------------------------- -module subroutine dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,ph,me) +module subroutine dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,ph,en) 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) :: ph,me + integer, intent(in) :: ph,en real(pReal), intent(in) :: T integer :: i,k,l,m,n @@ -565,13 +565,13 @@ module subroutine dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,ph,me) associate(prm => param(ph), stt => state(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)) + - sum(stt%f_tw(1:prm%sum_N_tw,en)) & + - sum(stt%f_tr(1:prm%sum_N_tr,en)) Lp = 0.0_pReal dLp_dMp = 0.0_pReal - call kinetics_slip(Mp,T,ph,me,dot_gamma_sl,ddot_gamma_dtau_slip) + call kinetics_slip(Mp,T,ph,en,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) & @@ -579,7 +579,7 @@ module subroutine dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,ph,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,ph,me,dot_gamma_tw,ddot_gamma_dtau_tw) + call kinetics_twin(Mp,T,dot_gamma_sl,ph,en,dot_gamma_tw,ddot_gamma_dtau_tw) twinContibution: do i = 1, prm%sum_N_tw Lp = Lp + dot_gamma_tw(i)*prm%P_tw(1:3,1:3,i) forall (k=1:3,l=1:3,m=1:3,n=1:3) & @@ -587,7 +587,7 @@ module subroutine dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,ph,me) + ddot_gamma_dtau_tw(i)* prm%P_tw(k,l,i)*prm%P_tw(m,n,i) enddo twinContibution - call kinetics_trans(Mp,T,dot_gamma_sl,ph,me,dot_gamma_tr,ddot_gamma_dtau_tr) + call kinetics_trans(Mp,T,dot_gamma_sl,ph,en,dot_gamma_tr,ddot_gamma_dtau_tr) 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) & @@ -632,7 +632,7 @@ end subroutine dislotwin_LpAndItsTangent !-------------------------------------------------------------------------------------------------- !> @brief Calculate the rate of change of microstructure. !-------------------------------------------------------------------------------------------------- -module subroutine dislotwin_dotState(Mp,T,ph,me) +module subroutine dislotwin_dotState(Mp,T,ph,en) real(pReal), dimension(3,3), intent(in):: & Mp !< Mandel stress @@ -640,7 +640,7 @@ module subroutine dislotwin_dotState(Mp,T,ph,me) T !< temperature at integration point integer, intent(in) :: & ph, & - me + en integer :: i real(pReal) :: & @@ -664,11 +664,11 @@ module subroutine dislotwin_dotState(Mp,T,ph,me) dot => dotState(ph), dst => dependentState(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)) + - sum(stt%f_tw(1:prm%sum_N_tw,en)) & + - sum(stt%f_tr(1:prm%sum_N_tr,en)) - call kinetics_slip(Mp,T,ph,me,dot_gamma_sl) - dot%gamma_sl(:,me) = abs(dot_gamma_sl) + call kinetics_slip(Mp,T,ph,en,dot_gamma_sl) + dot%gamma_sl(:,en) = abs(dot_gamma_sl) rho_dip_distance_min = prm%D_a*prm%b_sl @@ -680,11 +680,11 @@ module subroutine dislotwin_dotState(Mp,T,ph,me) dot_rho_dip_climb(i) = 0.0_pReal else significantSlipStress rho_dip_distance = 3.0_pReal*prm%mu*prm%b_sl(i)/(16.0_pReal*PI*abs(tau)) - rho_dip_distance = math_clip(rho_dip_distance, right = dst%Lambda_sl(i,me)) + rho_dip_distance = math_clip(rho_dip_distance, right = dst%Lambda_sl(i,en)) rho_dip_distance = math_clip(rho_dip_distance, left = rho_dip_distance_min(i)) dot_rho_dip_formation(i) = 2.0_pReal*(rho_dip_distance-rho_dip_distance_min(i))/prm%b_sl(i) & - * stt%rho_mob(i,me)*abs(dot_gamma_sl(i)) + * stt%rho_mob(i,en)*abs(dot_gamma_sl(i)) if (dEq(rho_dip_distance,rho_dip_distance_min(i))) then dot_rho_dip_climb(i) = 0.0_pReal @@ -698,25 +698,25 @@ module subroutine dislotwin_dotState(Mp,T,ph,me) v_cl = 2.0_pReal*prm%omega*b_d**2.0_pReal*exp(-prm%Q_cl/(kB*T)) & * (exp(abs(sigma_cl)*prm%b_sl(i)**3.0_pReal/(kB*T)) - 1.0_pReal) - dot_rho_dip_climb(i) = 4.0_pReal*v_cl*stt%rho_dip(i,me) & + dot_rho_dip_climb(i) = 4.0_pReal*v_cl*stt%rho_dip(i,en) & / (rho_dip_distance-rho_dip_distance_min(i)) endif endif significantSlipStress enddo slipState - dot%rho_mob(:,me) = abs(dot_gamma_sl)/(prm%b_sl*dst%Lambda_sl(:,me)) & + dot%rho_mob(:,en) = abs(dot_gamma_sl)/(prm%b_sl*dst%Lambda_sl(:,en)) & - dot_rho_dip_formation & - - 2.0_pReal*rho_dip_distance_min/prm%b_sl * stt%rho_mob(:,me)*abs(dot_gamma_sl) + - 2.0_pReal*rho_dip_distance_min/prm%b_sl * stt%rho_mob(:,en)*abs(dot_gamma_sl) - dot%rho_dip(:,me) = dot_rho_dip_formation & - - 2.0_pReal*rho_dip_distance_min/prm%b_sl * stt%rho_dip(:,me)*abs(dot_gamma_sl) & + dot%rho_dip(:,en) = dot_rho_dip_formation & + - 2.0_pReal*rho_dip_distance_min/prm%b_sl * stt%rho_dip(:,en)*abs(dot_gamma_sl) & - dot_rho_dip_climb - call kinetics_twin(Mp,T,dot_gamma_sl,ph,me,dot_gamma_tw) - dot%f_tw(:,me) = f_unrotated*dot_gamma_tw/prm%gamma_char + call kinetics_twin(Mp,T,dot_gamma_sl,ph,en,dot_gamma_tw) + dot%f_tw(:,en) = f_unrotated*dot_gamma_tw/prm%gamma_char - call kinetics_trans(Mp,T,dot_gamma_sl,ph,me,dot_gamma_tr) - dot%f_tr(:,me) = f_unrotated*dot_gamma_tr + call kinetics_trans(Mp,T,dot_gamma_sl,ph,en,dot_gamma_tr) + dot%f_tr(:,en) = f_unrotated*dot_gamma_tr end associate @@ -726,11 +726,11 @@ end subroutine dislotwin_dotState !-------------------------------------------------------------------------------------------------- !> @brief Calculate derived quantities from state. !-------------------------------------------------------------------------------------------------- -module subroutine dislotwin_dependentState(T,ph,me) +module subroutine dislotwin_dependentState(T,ph,en) integer, intent(in) :: & ph, & - me + en real(pReal), intent(in) :: & T @@ -752,50 +752,50 @@ module subroutine dislotwin_dependentState(T,ph,me) stt => state(ph),& dst => dependentState(ph)) - sumf_tw = sum(stt%f_tw(1:prm%sum_N_tw,me)) - sumf_tr = sum(stt%f_tr(1:prm%sum_N_tr,me)) + sumf_tw = sum(stt%f_tw(1:prm%sum_N_tw,en)) + sumf_tr = sum(stt%f_tr(1:prm%sum_N_tr,en)) Gamma = prm%Gamma_sf_0K + prm%dGamma_sf_dT * T !* rescaled volume fraction for topology - f_over_t_tw = stt%f_tw(1:prm%sum_N_tw,me)/prm%t_tw ! this is per system ... + f_over_t_tw = stt%f_tw(1:prm%sum_N_tw,en)/prm%t_tw ! this is per system ... f_over_t_tr = sumf_tr/prm%t_tr ! but this not ! ToDo ...Physically correct, but naming could be adjusted - inv_lambda_sl = sqrt(matmul(prm%forestProjection,stt%rho_mob(:,me)+stt%rho_dip(:,me)))/prm%i_sl + inv_lambda_sl = sqrt(matmul(prm%forestProjection,stt%rho_mob(:,en)+stt%rho_dip(:,en)))/prm%i_sl if (prm%sum_N_tw > 0 .and. prm%sum_N_sl > 0) & inv_lambda_sl = inv_lambda_sl + matmul(prm%h_sl_tw,f_over_t_tw)/(1.0_pReal-sumf_tw) if (prm%sum_N_tr > 0 .and. prm%sum_N_sl > 0) & inv_lambda_sl = inv_lambda_sl + matmul(prm%h_sl_tr,f_over_t_tr)/(1.0_pReal-sumf_tr) - dst%Lambda_sl(:,me) = prm%D / (1.0_pReal+prm%D*inv_lambda_sl) + dst%Lambda_sl(:,en) = prm%D / (1.0_pReal+prm%D*inv_lambda_sl) inv_lambda_tw_tw = matmul(prm%h_tw_tw,f_over_t_tw)/(1.0_pReal-sumf_tw) - dst%Lambda_tw(:,me) = prm%i_tw*prm%D/(1.0_pReal+prm%D*inv_lambda_tw_tw) + dst%Lambda_tw(:,en) = prm%i_tw*prm%D/(1.0_pReal+prm%D*inv_lambda_tw_tw) inv_lambda_tr_tr = matmul(prm%h_tr_tr,f_over_t_tr)/(1.0_pReal-sumf_tr) - dst%Lambda_tr(:,me) = prm%i_tr*prm%D/(1.0_pReal+prm%D*inv_lambda_tr_tr) + dst%Lambda_tr(:,en) = prm%i_tr*prm%D/(1.0_pReal+prm%D*inv_lambda_tr_tr) !* threshold stress for dislocation motion - dst%tau_pass(:,me) = prm%mu*prm%b_sl* sqrt(matmul(prm%h_sl_sl,stt%rho_mob(:,me)+stt%rho_dip(:,me))) + dst%tau_pass(:,en) = prm%mu*prm%b_sl* sqrt(matmul(prm%h_sl_sl,stt%rho_mob(:,en)+stt%rho_dip(:,en))) !* threshold stress for growing twin/martensite if(prm%sum_N_tw == prm%sum_N_sl) & - dst%tau_hat_tw(:,me) = Gamma/(3.0_pReal*prm%b_tw) & + dst%tau_hat_tw(:,en) = Gamma/(3.0_pReal*prm%b_tw) & + 3.0_pReal*prm%b_tw*prm%mu/(prm%L_tw*prm%b_sl) ! slip Burgers here correct? if(prm%sum_N_tr == prm%sum_N_sl) & - dst%tau_hat_tr(:,me) = Gamma/(3.0_pReal*prm%b_tr) & + dst%tau_hat_tr(:,en) = Gamma/(3.0_pReal*prm%b_tr) & + 3.0_pReal*prm%b_tr*prm%mu/(prm%L_tr*prm%b_sl) & ! slip Burgers here correct? + prm%h*prm%delta_G/ (3.0_pReal*prm%b_tr) - dst%V_tw(:,me) = (PI/4.0_pReal)*prm%t_tw*dst%Lambda_tw(:,me)**2.0_pReal - dst%V_tr(:,me) = (PI/4.0_pReal)*prm%t_tr*dst%Lambda_tr(:,me)**2.0_pReal + dst%V_tw(:,en) = (PI/4.0_pReal)*prm%t_tw*dst%Lambda_tw(:,en)**2.0_pReal + dst%V_tr(:,en) = (PI/4.0_pReal)*prm%t_tr*dst%Lambda_tr(:,en)**2.0_pReal x0 = prm%mu*prm%b_tw**2.0_pReal/(Gamma*8.0_pReal*PI)*(2.0_pReal+prm%nu)/(1.0_pReal-prm%nu) ! ToDo: In the paper, this is the Burgers vector for slip and is the same for twin and trans - dst%tau_r_tw(:,me) = prm%mu*prm%b_tw/(2.0_pReal*PI)*(1.0_pReal/(x0+prm%x_c_tw)+cos(pi/3.0_pReal)/x0) + dst%tau_r_tw(:,en) = prm%mu*prm%b_tw/(2.0_pReal*PI)*(1.0_pReal/(x0+prm%x_c_tw)+cos(pi/3.0_pReal)/x0) x0 = prm%mu*prm%b_tr**2.0_pReal/(Gamma*8.0_pReal*PI)*(2.0_pReal+prm%nu)/(1.0_pReal-prm%nu) ! ToDo: In the paper, this is the Burgers vector for slip - dst%tau_r_tr(:,me) = prm%mu*prm%b_tr/(2.0_pReal*PI)*(1.0_pReal/(x0+prm%x_c_tr)+cos(pi/3.0_pReal)/x0) + dst%tau_r_tr(:,en) = prm%mu*prm%b_tr/(2.0_pReal*PI)*(1.0_pReal/(x0+prm%x_c_tr)+cos(pi/3.0_pReal)/x0) end associate @@ -860,7 +860,7 @@ end subroutine plastic_dislotwin_results ! NOTE: Against the common convention, the result (i.e. intent(out)) variables are the last to ! have the optional arguments at the end !-------------------------------------------------------------------------------------------------- -pure subroutine kinetics_slip(Mp,T,ph,me, & +pure subroutine kinetics_slip(Mp,T,ph,en, & dot_gamma_sl,ddot_gamma_dtau_slip,tau_slip) real(pReal), dimension(3,3), intent(in) :: & @@ -869,7 +869,7 @@ pure subroutine kinetics_slip(Mp,T,ph,me, & T !< temperature integer, intent(in) :: & ph, & - me + en real(pReal), dimension(param(ph)%sum_N_sl), intent(out) :: & dot_gamma_sl @@ -898,7 +898,7 @@ pure subroutine kinetics_slip(Mp,T,ph,me, & tau(i) = math_tensordot(Mp,prm%P_sl(1:3,1:3,i)) enddo - tau_eff = abs(tau)-dst%tau_pass(:,me) + tau_eff = abs(tau)-dst%tau_pass(:,en) significantStress: where(tau_eff > tol_math_check) stressRatio = tau_eff/prm%tau_0 @@ -907,7 +907,7 @@ pure subroutine kinetics_slip(Mp,T,ph,me, & v_wait_inverse = prm%v_0**(-1.0_pReal) * exp(BoltzmannRatio*(1.0_pReal-StressRatio_p)** prm%q) v_run_inverse = prm%B/(tau_eff*prm%b_sl) - dot_gamma_sl = sign(stt%rho_mob(:,me)*prm%b_sl/(v_wait_inverse+v_run_inverse),tau) + dot_gamma_sl = sign(stt%rho_mob(:,en)*prm%b_sl/(v_wait_inverse+v_run_inverse),tau) dV_wait_inverse_dTau = -1.0_pReal * v_wait_inverse * prm%p * prm%q * BoltzmannRatio & * (stressRatio**(prm%p-1.0_pReal)) & @@ -916,7 +916,7 @@ pure subroutine kinetics_slip(Mp,T,ph,me, & dV_run_inverse_dTau = -1.0_pReal * v_run_inverse/tau_eff dV_dTau = -1.0_pReal * (dV_wait_inverse_dTau+dV_run_inverse_dTau) & / (v_wait_inverse+v_run_inverse)**2.0_pReal - ddot_gamma_dtau = dV_dTau*stt%rho_mob(:,me)*prm%b_sl + ddot_gamma_dtau = dV_dTau*stt%rho_mob(:,en)*prm%b_sl else where significantStress dot_gamma_sl = 0.0_pReal ddot_gamma_dtau = 0.0_pReal @@ -937,7 +937,7 @@ end subroutine kinetics_slip ! NOTE: Against the common convention, the result (i.e. intent(out)) variables are the last to ! have the optional arguments at the end. !-------------------------------------------------------------------------------------------------- -pure subroutine kinetics_twin(Mp,T,dot_gamma_sl,ph,me,& +pure subroutine kinetics_twin(Mp,T,dot_gamma_sl,ph,en,& dot_gamma_tw,ddot_gamma_dtau_tw) real(pReal), dimension(3,3), intent(in) :: & @@ -946,7 +946,7 @@ pure subroutine kinetics_twin(Mp,T,dot_gamma_sl,ph,me,& T !< temperature integer, intent(in) :: & ph, & - me + en real(pReal), dimension(param(ph)%sum_N_sl), intent(in) :: & dot_gamma_sl @@ -970,11 +970,11 @@ pure subroutine kinetics_twin(Mp,T,dot_gamma_sl,ph,me,& isFCC: if (prm%fccTwinTransNucleation) then s1=prm%fcc_twinNucleationSlipPair(1,i) s2=prm%fcc_twinNucleationSlipPair(2,i) - if (tau(i) < dst%tau_r_tw(i,me)) then ! ToDo: correct? - Ndot0=(abs(dot_gamma_sl(s1))*(stt%rho_mob(s2,me)+stt%rho_dip(s2,me))+& - abs(dot_gamma_sl(s2))*(stt%rho_mob(s1,me)+stt%rho_dip(s1,me)))/& ! ToDo: MD: it would be more consistent to use shearrates from state + if (tau(i) < dst%tau_r_tw(i,en)) then ! ToDo: correct? + Ndot0=(abs(dot_gamma_sl(s1))*(stt%rho_mob(s2,en)+stt%rho_dip(s2,en))+& + abs(dot_gamma_sl(s2))*(stt%rho_mob(s1,en)+stt%rho_dip(s1,en)))/& ! ToDo: MD: it would be more consistent to use shearrates from state (prm%L_tw*prm%b_sl(i))*& - (1.0_pReal-exp(-prm%V_cs/(kB*T)*(dst%tau_r_tw(i,me)-tau(i)))) ! P_ncs + (1.0_pReal-exp(-prm%V_cs/(kB*T)*(dst%tau_r_tw(i,en)-tau(i)))) ! P_ncs else Ndot0=0.0_pReal end if @@ -984,8 +984,8 @@ pure subroutine kinetics_twin(Mp,T,dot_gamma_sl,ph,me,& enddo significantStress: where(tau > tol_math_check) - StressRatio_r = (dst%tau_hat_tw(:,me)/tau)**prm%r - dot_gamma_tw = prm%gamma_char * dst%V_tw(:,me) * Ndot0*exp(-StressRatio_r) + StressRatio_r = (dst%tau_hat_tw(:,en)/tau)**prm%r + dot_gamma_tw = prm%gamma_char * dst%V_tw(:,en) * Ndot0*exp(-StressRatio_r) ddot_gamma_dtau = (dot_gamma_tw*prm%r/tau)*StressRatio_r else where significantStress dot_gamma_tw = 0.0_pReal @@ -1006,7 +1006,7 @@ end subroutine kinetics_twin ! NOTE: Against the common convention, the result (i.e. intent(out)) variables are the last to ! have the optional arguments at the end. !-------------------------------------------------------------------------------------------------- -pure subroutine kinetics_trans(Mp,T,dot_gamma_sl,ph,me,& +pure subroutine kinetics_trans(Mp,T,dot_gamma_sl,ph,en,& dot_gamma_tr,ddot_gamma_dtau_tr) real(pReal), dimension(3,3), intent(in) :: & @@ -1015,7 +1015,7 @@ pure subroutine kinetics_trans(Mp,T,dot_gamma_sl,ph,me,& T !< temperature integer, intent(in) :: & ph, & - me + en real(pReal), dimension(param(ph)%sum_N_sl), intent(in) :: & dot_gamma_sl @@ -1038,11 +1038,11 @@ pure subroutine kinetics_trans(Mp,T,dot_gamma_sl,ph,me,& isFCC: if (prm%fccTwinTransNucleation) then s1=prm%fcc_twinNucleationSlipPair(1,i) s2=prm%fcc_twinNucleationSlipPair(2,i) - if (tau(i) < dst%tau_r_tr(i,me)) then ! ToDo: correct? - Ndot0=(abs(dot_gamma_sl(s1))*(stt%rho_mob(s2,me)+stt%rho_dip(s2,me))+& - abs(dot_gamma_sl(s2))*(stt%rho_mob(s1,me)+stt%rho_dip(s1,me)))/& ! ToDo: MD: it would be more consistent to use shearrates from state + if (tau(i) < dst%tau_r_tr(i,en)) then ! ToDo: correct? + Ndot0=(abs(dot_gamma_sl(s1))*(stt%rho_mob(s2,en)+stt%rho_dip(s2,en))+& + abs(dot_gamma_sl(s2))*(stt%rho_mob(s1,en)+stt%rho_dip(s1,en)))/& ! ToDo: MD: it would be more consistent to use shearrates from state (prm%L_tr*prm%b_sl(i))*& - (1.0_pReal-exp(-prm%V_cs/(kB*T)*(dst%tau_r_tr(i,me)-tau(i)))) ! P_ncs + (1.0_pReal-exp(-prm%V_cs/(kB*T)*(dst%tau_r_tr(i,en)-tau(i)))) ! P_ncs else Ndot0=0.0_pReal end if @@ -1052,8 +1052,8 @@ pure subroutine kinetics_trans(Mp,T,dot_gamma_sl,ph,me,& enddo significantStress: where(tau > tol_math_check) - StressRatio_s = (dst%tau_hat_tr(:,me)/tau)**prm%s - dot_gamma_tr = dst%V_tr(:,me) * Ndot0*exp(-StressRatio_s) + StressRatio_s = (dst%tau_hat_tr(:,en)/tau)**prm%s + dot_gamma_tr = dst%V_tr(:,en) * Ndot0*exp(-StressRatio_s) ddot_gamma_dtau = (dot_gamma_tr*prm%s/tau)*StressRatio_s else where significantStress dot_gamma_tr = 0.0_pReal diff --git a/src/phase_mechanical_plastic_isotropic.f90 b/src/phase_mechanical_plastic_isotropic.f90 index d1f191f28..5ab73895f 100644 --- a/src/phase_mechanical_plastic_isotropic.f90 +++ b/src/phase_mechanical_plastic_isotropic.f90 @@ -156,7 +156,7 @@ end function plastic_isotropic_init !-------------------------------------------------------------------------------------------------- !> @brief Calculate plastic velocity gradient and its tangent. !-------------------------------------------------------------------------------------------------- -module subroutine isotropic_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,me) +module subroutine isotropic_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,en) real(pReal), dimension(3,3), intent(out) :: & Lp !< plastic velocity gradient @@ -167,7 +167,7 @@ module subroutine isotropic_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,me) Mp !< Mandel stress integer, intent(in) :: & ph, & - me + en real(pReal), dimension(3,3) :: & Mp_dev !< deviatoric part of the Mandel stress @@ -185,7 +185,7 @@ module subroutine isotropic_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,me) norm_Mp_dev = sqrt(squarenorm_Mp_dev) if (norm_Mp_dev > 0.0_pReal) then - dot_gamma = prm%dot_gamma_0 * (sqrt(1.5_pReal) * norm_Mp_dev/(prm%M*stt%xi(me))) **prm%n + dot_gamma = prm%dot_gamma_0 * (sqrt(1.5_pReal) * norm_Mp_dev/(prm%M*stt%xi(en))) **prm%n Lp = dot_gamma/prm%M * Mp_dev/norm_Mp_dev forall (k=1:3,l=1:3,m=1:3,n=1:3) & @@ -248,13 +248,13 @@ module subroutine plastic_isotropic_LiAndItsTangent(Li,dLi_dMi,Mi,ph,en) !-------------------------------------------------------------------------------------------------- !> @brief Calculate the rate of change of microstructure. !-------------------------------------------------------------------------------------------------- -module subroutine isotropic_dotState(Mp,ph,me) +module subroutine isotropic_dotState(Mp,ph,en) real(pReal), dimension(3,3), intent(in) :: & Mp !< Mandel stress integer, intent(in) :: & ph, & - me + en real(pReal) :: & dot_gamma, & !< strainrate @@ -270,7 +270,7 @@ module subroutine isotropic_dotState(Mp,ph,me) norm_Mp = sqrt(math_tensordot(math_deviatoric33(Mp),math_deviatoric33(Mp))) endif - dot_gamma = prm%dot_gamma_0 * (sqrt(1.5_pReal) * norm_Mp /(prm%M*stt%xi(me))) **prm%n + dot_gamma = prm%dot_gamma_0 * (sqrt(1.5_pReal) * norm_Mp /(prm%M*stt%xi(en))) **prm%n if (dot_gamma > 1e-12_pReal) then if (dEq0(prm%c_1)) then @@ -280,15 +280,15 @@ module subroutine isotropic_dotState(Mp,ph,me) + asinh( (dot_gamma / prm%c_1)**(1.0_pReal / prm%c_2))**(1.0_pReal / prm%c_3) & / prm%c_4 * (dot_gamma / prm%dot_gamma_0)**(1.0_pReal / prm%n) endif - dot%xi(me) = dot_gamma & + dot%xi(en) = dot_gamma & * ( prm%h_0 + prm%h_ln * log(dot_gamma) ) & - * abs( 1.0_pReal - stt%xi(me)/xi_inf_star )**prm%a & - * sign(1.0_pReal, 1.0_pReal - stt%xi(me)/xi_inf_star) + * abs( 1.0_pReal - stt%xi(en)/xi_inf_star )**prm%a & + * sign(1.0_pReal, 1.0_pReal - stt%xi(en)/xi_inf_star) else - dot%xi(me) = 0.0_pReal + dot%xi(en) = 0.0_pReal endif - dot%gamma(me) = dot_gamma ! ToDo: not really used + dot%gamma(en) = dot_gamma ! ToDo: not really used end associate diff --git a/src/phase_mechanical_plastic_kinehardening.f90 b/src/phase_mechanical_plastic_kinehardening.f90 index 936a7f884..18f467617 100644 --- a/src/phase_mechanical_plastic_kinehardening.f90 +++ b/src/phase_mechanical_plastic_kinehardening.f90 @@ -230,7 +230,7 @@ end function plastic_kinehardening_init !-------------------------------------------------------------------------------------------------- !> @brief Calculate plastic velocity gradient and its tangent. !-------------------------------------------------------------------------------------------------- -pure module subroutine kinehardening_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,me) +pure module subroutine kinehardening_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,en) real(pReal), dimension(3,3), intent(out) :: & Lp !< plastic velocity gradient @@ -241,7 +241,7 @@ pure module subroutine kinehardening_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,me) Mp !< Mandel stress integer, intent(in) :: & ph, & - me + en integer :: & i,k,l,m,n @@ -254,7 +254,7 @@ pure module subroutine kinehardening_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,me) associate(prm => param(ph)) - call kinetics(Mp,ph,me,gdot_pos,gdot_neg,dgdot_dtau_pos,dgdot_dtau_neg) + call kinetics(Mp,ph,en,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) @@ -272,13 +272,13 @@ end subroutine kinehardening_LpAndItsTangent !-------------------------------------------------------------------------------------------------- !> @brief Calculate the rate of change of microstructure. !-------------------------------------------------------------------------------------------------- -module subroutine plastic_kinehardening_dotState(Mp,ph,me) +module subroutine plastic_kinehardening_dotState(Mp,ph,en) real(pReal), dimension(3,3), intent(in) :: & Mp !< Mandel stress integer, intent(in) :: & ph, & - me + en real(pReal) :: & sumGamma @@ -289,22 +289,22 @@ module subroutine plastic_kinehardening_dotState(Mp,ph,me) associate(prm => param(ph), stt => state(ph),& dot => dotState(ph)) - call kinetics(Mp,ph,me,gdot_pos,gdot_neg) - dot%accshear(:,me) = abs(gdot_pos+gdot_neg) - sumGamma = sum(stt%accshear(:,me)) + call kinetics(Mp,ph,en,gdot_pos,gdot_neg) + dot%accshear(:,en) = abs(gdot_pos+gdot_neg) + sumGamma = sum(stt%accshear(:,en)) - dot%crss(:,me) = matmul(prm%interaction_SlipSlip,dot%accshear(:,me)) & + dot%crss(:,en) = matmul(prm%interaction_SlipSlip,dot%accshear(:,en)) & * ( prm%h_inf_f & + (prm%h_0_f - prm%h_inf_f + prm%h_0_f*prm%h_inf_f*sumGamma/prm%xi_inf_f) & * exp(-sumGamma*prm%h_0_f/prm%xi_inf_f) & ) - dot%crss_back(:,me) = stt%sense(:,me)*dot%accshear(:,me) * & + dot%crss_back(:,en) = stt%sense(:,en)*dot%accshear(:,en) * & ( prm%h_inf_b + & (prm%h_0_b - prm%h_inf_b & - + prm%h_0_b*prm%h_inf_b/(prm%xi_inf_b+stt%chi0(:,me))*(stt%accshear(:,me)-stt%gamma0(:,me))& - ) *exp(-(stt%accshear(:,me)-stt%gamma0(:,me)) *prm%h_0_b/(prm%xi_inf_b+stt%chi0(:,me))) & + + prm%h_0_b*prm%h_inf_b/(prm%xi_inf_b+stt%chi0(:,en))*(stt%accshear(:,en)-stt%gamma0(:,en))& + ) *exp(-(stt%accshear(:,en)-stt%gamma0(:,en)) *prm%h_0_b/(prm%xi_inf_b+stt%chi0(:,en))) & ) end associate @@ -315,13 +315,13 @@ end subroutine plastic_kinehardening_dotState !-------------------------------------------------------------------------------------------------- !> @brief Calculate (instantaneous) incremental change of microstructure. !-------------------------------------------------------------------------------------------------- -module subroutine plastic_kinehardening_deltaState(Mp,ph,me) +module subroutine plastic_kinehardening_deltaState(Mp,ph,en) real(pReal), dimension(3,3), intent(in) :: & Mp !< Mandel stress integer, intent(in) :: & ph, & - me + en real(pReal), dimension(param(ph)%sum_N_sl) :: & gdot_pos,gdot_neg, & @@ -329,22 +329,22 @@ module subroutine plastic_kinehardening_deltaState(Mp,ph,me) associate(prm => param(ph), stt => state(ph), dlt => deltaState(ph)) - call kinetics(Mp,ph,me,gdot_pos,gdot_neg) - sense = merge(state(ph)%sense(:,me), & ! keep existing... + call kinetics(Mp,ph,en,gdot_pos,gdot_neg) + sense = merge(state(ph)%sense(:,en), & ! keep existing... sign(1.0_pReal,gdot_pos+gdot_neg), & ! ...or have a defined dEq0(gdot_pos+gdot_neg,1e-10_pReal)) ! current sense of shear direction !-------------------------------------------------------------------------------------------------- -! switch in sense me shear? - where(dNeq(sense,stt%sense(:,me),0.1_pReal)) - dlt%sense (:,me) = sense - stt%sense(:,me) ! switch sense - dlt%chi0 (:,me) = abs(stt%crss_back(:,me)) - stt%chi0(:,me) ! remember current backstress magnitude - dlt%gamma0(:,me) = stt%accshear(:,me) - stt%gamma0(:,me) ! remember current accumulated shear +! switch in sense of shear? + where(dNeq(sense,stt%sense(:,en),0.1_pReal)) + dlt%sense (:,en) = sense - stt%sense(:,en) ! switch sense + dlt%chi0 (:,en) = abs(stt%crss_back(:,en)) - stt%chi0(:,en) ! remember current backstress magnitude + dlt%gamma0(:,en) = stt%accshear(:,en) - stt%gamma0(:,en) ! remember current accumulated shear else where - dlt%sense (:,me) = 0.0_pReal - dlt%chi0 (:,me) = 0.0_pReal - dlt%gamma0(:,me) = 0.0_pReal + dlt%sense (:,en) = 0.0_pReal + dlt%chi0 (:,en) = 0.0_pReal + dlt%gamma0(:,en) = 0.0_pReal end where end associate @@ -397,14 +397,14 @@ end subroutine plastic_kinehardening_results ! NOTE: Against the common convention, the result (i.e. intent(out)) variables are the last to ! have the optional arguments at the end. !-------------------------------------------------------------------------------------------------- -pure subroutine kinetics(Mp,ph,me, & +pure subroutine kinetics(Mp,ph,en, & gdot_pos,gdot_neg,dgdot_dtau_pos,dgdot_dtau_neg) real(pReal), dimension(3,3), intent(in) :: & Mp !< Mandel stress integer, intent(in) :: & ph, & - me + en real(pReal), intent(out), dimension(param(ph)%sum_N_sl) :: & gdot_pos, & @@ -421,21 +421,21 @@ pure subroutine kinetics(Mp,ph,me, & associate(prm => param(ph), stt => state(ph)) do i = 1, prm%sum_N_sl - tau_pos(i) = math_tensordot(Mp,prm%nonSchmid_pos(1:3,1:3,i)) - stt%crss_back(i,me) - tau_neg(i) = merge(math_tensordot(Mp,prm%nonSchmid_neg(1:3,1:3,i)) - stt%crss_back(i,me), & + tau_pos(i) = math_tensordot(Mp,prm%nonSchmid_pos(1:3,1:3,i)) - stt%crss_back(i,en) + tau_neg(i) = merge(math_tensordot(Mp,prm%nonSchmid_neg(1:3,1:3,i)) - stt%crss_back(i,en), & 0.0_pReal, prm%nonSchmidActive) enddo where(dNeq0(tau_pos)) gdot_pos = prm%dot_gamma_0 * merge(0.5_pReal,1.0_pReal, prm%nonSchmidActive) & ! 1/2 if non-Schmid active - * sign(abs(tau_pos/stt%crss(:,me))**prm%n, tau_pos) + * sign(abs(tau_pos/stt%crss(:,en))**prm%n, tau_pos) else where gdot_pos = 0.0_pReal end where where(dNeq0(tau_neg)) gdot_neg = prm%dot_gamma_0 * 0.5_pReal & ! only used if non-Schmid active, always 1/2 - * sign(abs(tau_neg/stt%crss(:,me))**prm%n, tau_neg) + * sign(abs(tau_neg/stt%crss(:,en))**prm%n, tau_neg) else where gdot_neg = 0.0_pReal end where diff --git a/src/phase_mechanical_plastic_nonlocal.f90 b/src/phase_mechanical_plastic_nonlocal.f90 index ba7f76881..05ba371ce 100644 --- a/src/phase_mechanical_plastic_nonlocal.f90 +++ b/src/phase_mechanical_plastic_nonlocal.f90 @@ -49,7 +49,7 @@ submodule(phase:plastic) nonlocal !END DEPRECATED real(pReal), dimension(:,:,:,:,:,:), allocatable :: & - compatibility !< slip system compatibility between me and my neighbors + compatibility !< slip system compatibility between en and my neighbors type :: tInitialParameters !< container type for internal constitutive parameters real(pReal) :: & @@ -558,11 +558,11 @@ end function plastic_nonlocal_init !-------------------------------------------------------------------------------------------------- !> @brief calculates quantities characterizing the microstructure !-------------------------------------------------------------------------------------------------- -module subroutine nonlocal_dependentState(ph, me, ip, el) +module subroutine nonlocal_dependentState(ph, en, ip, el) integer, intent(in) :: & ph, & - me, & + en, & ip, & el @@ -615,9 +615,9 @@ module subroutine nonlocal_dependentState(ph, me, ip, el) associate(prm => param(ph),dst => microstructure(ph), stt => state(ph)) - rho = getRho(ph,me) + rho = getRho(ph,en) - stt%rho_forest(:,me) = matmul(prm%forestProjection_Edge, sum(abs(rho(:,edg)),2)) & + stt%rho_forest(:,en) = matmul(prm%forestProjection_Edge, sum(abs(rho(:,edg)),2)) & + matmul(prm%forestProjection_Screw,sum(abs(rho(:,scr)),2)) @@ -627,13 +627,13 @@ module subroutine nonlocal_dependentState(ph, me, ip, el) myInteractionMatrix = prm%h_sl_sl & * spread(( 1.0_pReal - prm%f_F & + prm%f_F & - * log(0.35_pReal * prm%b_sl * sqrt(max(stt%rho_forest(:,me),prm%rho_significant))) & + * log(0.35_pReal * prm%b_sl * sqrt(max(stt%rho_forest(:,en),prm%rho_significant))) & / log(0.35_pReal * prm%b_sl * 1e6_pReal))** 2.0_pReal,2,prm%sum_N_sl) else myInteractionMatrix = prm%h_sl_sl endif - dst%tau_pass(:,me) = prm%mu * prm%b_sl & + dst%tau_pass(:,en) = prm%mu * prm%b_sl & * sqrt(matmul(myInteractionMatrix,sum(abs(rho),2))) !*** calculate the dislocation stress of the neighboring excess dislocation densities @@ -643,10 +643,10 @@ module subroutine nonlocal_dependentState(ph, me, ip, el) ! ToDo: MD: this is most likely only correct for F_i = I !################################################################################################# - rho0 = getRho0(ph,me) + rho0 = getRho0(ph,en) if (.not. phase_localPlasticity(material_phaseAt(1,el)) .and. prm%shortRangeStressCorrection) then - invFp = math_inv33(phase_mechanical_Fp(ph)%data(1:3,1:3,me)) - invFe = math_inv33(phase_mechanical_Fe(ph)%data(1:3,1:3,me)) + invFp = math_inv33(phase_mechanical_Fp(ph)%data(1:3,1:3,en)) + invFe = math_inv33(phase_mechanical_Fe(ph)%data(1:3,1:3,en)) rho_edg_delta = rho0(:,mob_edg_pos) - rho0(:,mob_edg_neg) rho_scr_delta = rho0(:,mob_scr_pos) - rho0(:,mob_scr_neg) @@ -654,7 +654,7 @@ module subroutine nonlocal_dependentState(ph, me, ip, el) rhoExcess(1,:) = rho_edg_delta rhoExcess(2,:) = rho_scr_delta - FVsize = geom(ph)%V_0(me) ** (1.0_pReal/3.0_pReal) + FVsize = geom(ph)%V_0(en) ** (1.0_pReal/3.0_pReal) !* loop through my neighborhood and get the connection vectors (in lattice frame) and the excess densities @@ -734,7 +734,7 @@ module subroutine nonlocal_dependentState(ph, me, ip, el) where(rhoTotal > 0.0_pReal) rhoExcessGradient_over_rho = rhoExcessGradient / rhoTotal ! ... gives the local stress correction when multiplied with a factor - dst%tau_back(s,me) = - prm%mu * prm%b_sl(s) / (2.0_pReal * PI) & + dst%tau_back(s,en) = - prm%mu * prm%b_sl(s) / (2.0_pReal * PI) & * ( rhoExcessGradient_over_rho(1) / (1.0_pReal - prm%nu) & + rhoExcessGradient_over_rho(2)) enddo @@ -745,9 +745,9 @@ module subroutine nonlocal_dependentState(ph, me, ip, el) .and. ((debugConstitutive%element == el .and. debugConstitutive%ip == ip)& .or. .not. debugConstitutive%selective)) then print'(/,a,i8,1x,i2,1x,i1,/)', '<< CONST >> nonlocal_microstructure at el ip ',el,ip - print'(a,/,12x,12(e10.3,1x))', '<< CONST >> rhoForest', stt%rho_forest(:,me) - print'(a,/,12x,12(f10.5,1x))', '<< CONST >> tauThreshold / MPa', dst%tau_pass(:,me)*1e-6 - print'(a,/,12x,12(f10.5,1x),/)', '<< CONST >> tauBack / MPa', dst%tau_back(:,me)*1e-6 + print'(a,/,12x,12(e10.3,1x))', '<< CONST >> rhoForest', stt%rho_forest(:,en) + print'(a,/,12x,12(f10.5,1x))', '<< CONST >> tauThreshold / MPa', dst%tau_pass(:,en)*1e-6 + print'(a,/,12x,12(f10.5,1x),/)', '<< CONST >> tauBack / MPa', dst%tau_back(:,en)*1e-6 endif #endif @@ -760,14 +760,14 @@ end subroutine nonlocal_dependentState !> @brief calculates plastic velocity gradient and its tangent !-------------------------------------------------------------------------------------------------- module subroutine nonlocal_LpAndItsTangent(Lp,dLp_dMp, & - Mp,Temperature,ph,me) + Mp,Temperature,ph,en) 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) :: & ph, & - me + en real(pReal), intent(in) :: & Temperature !< temperature @@ -800,7 +800,7 @@ module subroutine nonlocal_LpAndItsTangent(Lp,dLp_dMp, & ns = prm%sum_N_sl !*** shortcut to state variables - rho = getRho(ph,me) + rho = getRho(ph,en) rhoSgl = rho(:,sgl) do s = 1,ns @@ -815,12 +815,12 @@ module subroutine nonlocal_LpAndItsTangent(Lp,dLp_dMp, & tauNS(s,4) = math_tensordot(Mp, -prm%nonSchmid_pos(1:3,1:3,s)) endif enddo - tauNS = tauNS + spread(dst%tau_back(:,me),2,4) - tau = tau + dst%tau_back(:,me) + tauNS = tauNS + spread(dst%tau_back(:,en),2,4) + tau = tau + dst%tau_back(:,en) ! edges call kinetics(v(:,1), dv_dtau(:,1), dv_dtauNS(:,1), & - tau, tauNS(:,1), dst%tau_pass(:,me),1,Temperature, ph) + tau, tauNS(:,1), dst%tau_pass(:,en),1,Temperature, ph) v(:,2) = v(:,1) dv_dtau(:,2) = dv_dtau(:,1) dv_dtauNS(:,2) = dv_dtauNS(:,1) @@ -829,7 +829,7 @@ module subroutine nonlocal_LpAndItsTangent(Lp,dLp_dMp, & if (prm%nonSchmidActive) then do t = 3,4 call kinetics(v(:,t), dv_dtau(:,t), dv_dtauNS(:,t), & - tau, tauNS(:,t), dst%tau_pass(:,me),2,Temperature, ph) + tau, tauNS(:,t), dst%tau_pass(:,en),2,Temperature, ph) enddo else v(:,3:4) = spread(v(:,1),2,2) @@ -837,7 +837,7 @@ module subroutine nonlocal_LpAndItsTangent(Lp,dLp_dMp, & dv_dtauNS(:,3:4) = spread(dv_dtauNS(:,1),2,2) endif - stt%v(:,me) = pack(v,.true.) + stt%v(:,en) = pack(v,.true.) !*** Bauschinger effect forall (s = 1:ns, t = 5:8, rhoSgl(s,t) * v(s,t-4) < 0.0_pReal) & @@ -866,13 +866,13 @@ end subroutine nonlocal_LpAndItsTangent !-------------------------------------------------------------------------------------------------- !> @brief (instantaneous) incremental change of microstructure !-------------------------------------------------------------------------------------------------- -module subroutine plastic_nonlocal_deltaState(Mp,ph,me) +module subroutine plastic_nonlocal_deltaState(Mp,ph,en) real(pReal), dimension(3,3), intent(in) :: & Mp !< MandelStress integer, intent(in) :: & ph, & - me + en integer :: & ns, & ! short notation for the total number of active slip systems @@ -898,10 +898,10 @@ module subroutine plastic_nonlocal_deltaState(Mp,ph,me) ns = prm%sum_N_sl !*** shortcut to state variables - forall (s = 1:ns, t = 1:4) v(s,t) = plasticState(ph)%state(iV(s,t,ph),me) - forall (s = 1:ns, c = 1:2) dUpperOld(s,c) = plasticState(ph)%state(iD(s,c,ph),me) + forall (s = 1:ns, t = 1:4) v(s,t) = plasticState(ph)%state(iV(s,t,ph),en) + forall (s = 1:ns, c = 1:2) dUpperOld(s,c) = plasticState(ph)%state(iD(s,c,ph),en) - rho = getRho(ph,me) + rho = getRho(ph,en) rhoDip = rho(:,dip) !**************************************************************************** @@ -922,7 +922,7 @@ module subroutine plastic_nonlocal_deltaState(Mp,ph,me) !*** calculate limits for stable dipole height do s = 1,prm%sum_N_sl - tau(s) = math_tensordot(Mp, prm%Schmid(1:3,1:3,s)) +dst%tau_back(s,me) + tau(s) = math_tensordot(Mp, prm%Schmid(1:3,1:3,s)) +dst%tau_back(s,en) if (abs(tau(s)) < 1.0e-15_pReal) tau(s) = 1.0e-15_pReal enddo @@ -946,10 +946,10 @@ module subroutine plastic_nonlocal_deltaState(Mp,ph,me) / (dUpperOld(s,c) - prm%minDipoleHeight(s,c)) forall (t=1:4) deltaRhoDipole2SingleStress(:,t) = -0.5_pReal * deltaRhoDipole2SingleStress(:,(t-1)/2+9) - forall (s = 1:ns, c = 1:2) plasticState(ph)%state(iD(s,c,ph),me) = dUpper(s,c) + forall (s = 1:ns, c = 1:2) plasticState(ph)%state(iD(s,c,ph),en) = dUpper(s,c) - plasticState(ph)%deltaState(:,me) = 0.0_pReal - del%rho(:,me) = reshape(deltaRhoRemobilization + deltaRhoDipole2SingleStress, [10*ns]) + plasticState(ph)%deltaState(:,en) = 0.0_pReal + del%rho(:,en) = reshape(deltaRhoRemobilization + deltaRhoDipole2SingleStress, [10*ns]) end associate @@ -960,7 +960,7 @@ end subroutine plastic_nonlocal_deltaState !> @brief calculates the rate of change of microstructure !--------------------------------------------------------------------------------------------------- module subroutine nonlocal_dotState(Mp, Temperature,timestep, & - ph,me,ip,el) + ph,en,ip,el) real(pReal), dimension(3,3), intent(in) :: & Mp !< MandelStress @@ -969,7 +969,7 @@ module subroutine nonlocal_dotState(Mp, Temperature,timestep, & timestep !< substepped crystallite time increment integer, intent(in) :: & ph, & - me, & + en, & ip, & !< current integration point el !< current element number @@ -1017,13 +1017,13 @@ module subroutine nonlocal_dotState(Mp, Temperature,timestep, & tau = 0.0_pReal gdot = 0.0_pReal - rho = getRho(ph,me) + rho = getRho(ph,en) rhoSgl = rho(:,sgl) rhoDip = rho(:,dip) - rho0 = getRho0(ph,me) + rho0 = getRho0(ph,en) my_rhoSgl0 = rho0(:,sgl) - forall (s = 1:ns, t = 1:4) v(s,t) = plasticState(ph)%state(iV(s,t,ph),me) + forall (s = 1:ns, t = 1:4) v(s,t) = plasticState(ph)%state(iV(s,t,ph),en) gdot = rhoSgl(:,1:4) * v * spread(prm%b_sl,2,4) #ifdef DEBUG @@ -1038,7 +1038,7 @@ module subroutine nonlocal_dotState(Mp, Temperature,timestep, & !**************************************************************************** !*** limits for stable dipole height do s = 1,ns - tau(s) = math_tensordot(Mp, prm%Schmid(1:3,1:3,s)) + dst%tau_back(s,me) + tau(s) = math_tensordot(Mp, prm%Schmid(1:3,1:3,s)) + dst%tau_back(s,en) if (abs(tau(s)) < 1.0e-15_pReal) tau(s) = 1.0e-15_pReal enddo @@ -1059,20 +1059,20 @@ module subroutine nonlocal_dotState(Mp, Temperature,timestep, & isBCC: if (lattice_structure(ph) == LATTICE_bcc_ID) then forall (s = 1:ns, sum(abs(v(s,1:4))) > 0.0_pReal) rhoDotMultiplication(s,1:2) = sum(abs(gdot(s,3:4))) / prm%b_sl(s) & ! assuming double-cross-slip of screws to be decisive for multiplication - * sqrt(stt%rho_forest(s,me)) / prm%i_sl(s) ! & ! mean free path + * sqrt(stt%rho_forest(s,en)) / prm%i_sl(s) ! & ! mean free path ! * 2.0_pReal * sum(abs(v(s,3:4))) / sum(abs(v(s,1:4))) ! ratio of screw to overall velocity determines edge generation rhoDotMultiplication(s,3:4) = sum(abs(gdot(s,3:4))) /prm%b_sl(s) & ! assuming double-cross-slip of screws to be decisive for multiplication - * sqrt(stt%rho_forest(s,me)) / prm%i_sl(s) ! & ! mean free path + * sqrt(stt%rho_forest(s,en)) / prm%i_sl(s) ! & ! mean free path ! * 2.0_pReal * sum(abs(v(s,1:2))) / sum(abs(v(s,1:4))) ! ratio of edge to overall velocity determines screw generation endforall else isBCC rhoDotMultiplication(:,1:4) = spread( & (sum(abs(gdot(:,1:2)),2) * prm%f_ed_mult + sum(abs(gdot(:,3:4)),2)) & - * sqrt(stt%rho_forest(:,me)) / prm%i_sl / prm%b_sl, 2, 4) + * sqrt(stt%rho_forest(:,en)) / 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,ph),me) + forall (s = 1:ns, t = 1:4) v0(s,t) = plasticState(ph)%state0(iV(s,t,ph),en) !**************************************************************************** @@ -1115,10 +1115,10 @@ module subroutine nonlocal_dotState(Mp, Temperature,timestep, & if (lattice_structure(ph) == LATTICE_fcc_ID) & forall (s = 1:ns, prm%colinearSystem(s) > 0) & rhoDotAthermalAnnihilation(prm%colinearSystem(s),1:2) = - rhoDotAthermalAnnihilation(s,10) & - * 0.25_pReal * sqrt(stt%rho_forest(s,me)) * (dUpper(s,2) + dLower(s,2)) * prm%f_ed + * 0.25_pReal * sqrt(stt%rho_forest(s,en)) * (dUpper(s,2) + dLower(s,2)) * prm%f_ed - !*** thermally activated annihilation me edge dipoles by climb + !*** thermally activated annihilation of edge dipoles by climb rhoDotThermalAnnihilation = 0.0_pReal selfDiffusion = prm%D_0 * exp(-prm%Q_cl / (kB * Temperature)) vClimb = prm%V_at * selfDiffusion * prm%mu & @@ -1128,7 +1128,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, ph,me,ip,el) & + rhoDot = rhoDotFlux(timestep, ph,en,ip,el) & + rhoDotMultiplication & + rhoDotSingle2DipoleGlide & + rhoDotAthermalAnnihilation & @@ -1145,8 +1145,8 @@ module subroutine nonlocal_dotState(Mp, Temperature,timestep, & #endif plasticState(ph)%dotState = IEEE_value(1.0_pReal,IEEE_quiet_NaN) else - dot%rho(:,me) = pack(rhoDot,.true.) - dot%gamma(:,me) = sum(gdot,2) + dot%rho(:,en) = pack(rhoDot,.true.) + dot%gamma(:,en) = sum(gdot,2) endif end associate @@ -1157,13 +1157,13 @@ end subroutine nonlocal_dotState !--------------------------------------------------------------------------------------------------- !> @brief calculates the rate of change of microstructure !--------------------------------------------------------------------------------------------------- -function rhoDotFlux(timestep,ph,me,ip,el) +function rhoDotFlux(timestep,ph,en,ip,el) real(pReal), intent(in) :: & timestep !< substepped crystallite time increment integer, intent(in) :: & ph, & - me, & + en, & ip, & !< current integration point el !< current element number @@ -1174,11 +1174,11 @@ function rhoDotFlux(timestep,ph,me,ip,el) n, & !< index of my current neighbor neighbor_el, & !< element number of my neighbor neighbor_ip, & !< integration point of my neighbor - neighbor_n, & !< neighbor index pointing to me when looking from my neighbor + neighbor_n, & !< neighbor index pointing to en when looking from my neighbor opposite_neighbor, & !< index of my opposite neighbor opposite_ip, & !< ip of my opposite neighbor opposite_el, & !< element index of my opposite neighbor - opposite_n, & !< neighbor index pointing to me when looking from my opposite neighbor + opposite_n, & !< neighbor index pointing to en when looking from my opposite neighbor t, & !< type of dislocation no,& !< neighbor offset shortcut np,& !< neighbor phase shortcut @@ -1204,12 +1204,12 @@ function rhoDotFlux(timestep,ph,me,ip,el) neighbor_F, & !< total deformation gradient of my neighbor my_Fe, & !< my elastic deformation gradient neighbor_Fe, & !< elastic deformation gradient of my neighbor - Favg !< average total deformation gradient of me and my neighbor + Favg !< average total deformation gradient of en and my neighbor real(pReal), dimension(3) :: & - normal_neighbor2me, & !< interface normal pointing from my neighbor to me in neighbor's lattice configuration - normal_neighbor2me_defConf, & !< interface normal pointing from my neighbor to me in shared deformed configuration - normal_me2neighbor, & !< interface normal pointing from me to my neighbor in my lattice configuration - normal_me2neighbor_defConf !< interface normal pointing from me to my neighbor in shared deformed configuration + normal_neighbor2me, & !< interface normal pointing from my neighbor to en in neighbor's lattice configuration + normal_neighbor2me_defConf, & !< interface normal pointing from my neighbor to en in shared deformed configuration + normal_me2neighbor, & !< interface normal pointing from en to my neighbor in my lattice configuration + normal_me2neighbor_defConf !< interface normal pointing from en to my neighbor in shared deformed configuration real(pReal) :: & area, & !< area of the current interface transmissivity, & !< overall transmissivity of dislocation flux to neighboring material point @@ -1224,16 +1224,16 @@ function rhoDotFlux(timestep,ph,me,ip,el) gdot = 0.0_pReal - rho = getRho(ph,me) + rho = getRho(ph,en) rhoSgl = rho(:,sgl) - rho0 = getRho0(ph,me) + rho0 = getRho0(ph,en) my_rhoSgl0 = rho0(:,sgl) - forall (s = 1:ns, t = 1:4) v(s,t) = plasticState(ph)%state(iV(s,t,ph),me) !ToDo: MD: I think we should use state0 here + forall (s = 1:ns, t = 1:4) v(s,t) = plasticState(ph)%state(iV(s,t,ph),en) !ToDo: MD: I think we should use state0 here gdot = rhoSgl(:,1:4) * v * spread(prm%b_sl,2,4) - forall (s = 1:ns, t = 1:4) v0(s,t) = plasticState(ph)%state0(iV(s,t,ph),me) + forall (s = 1:ns, t = 1:4) v0(s,t) = plasticState(ph)%state0(iV(s,t,ph),en) !**************************************************************************** !*** calculate dislocation fluxes (only for nonlocal plasticity) @@ -1268,8 +1268,8 @@ function rhoDotFlux(timestep,ph,me,ip,el) m(1:3,:,3) = -prm%slip_transverse m(1:3,:,4) = prm%slip_transverse - my_F = phase_mechanical_F(ph)%data(1:3,1:3,me) - my_Fe = matmul(my_F, math_inv33(phase_mechanical_Fp(ph)%data(1:3,1:3,me))) + my_F = phase_mechanical_F(ph)%data(1:3,1:3,en) + my_Fe = matmul(my_F, math_inv33(phase_mechanical_Fp(ph)%data(1:3,1:3,en))) neighbors: do n = 1,nIPneighbors @@ -1316,7 +1316,7 @@ function rhoDotFlux(timestep,ph,me,ip,el) .or. neighbor_rhoSgl0 < prm%rho_significant) & neighbor_rhoSgl0 = 0.0_pReal normal_neighbor2me_defConf = math_det33(Favg) * matmul(math_inv33(transpose(Favg)), & - IPareaNormal(1:3,neighbor_n,neighbor_ip,neighbor_el)) ! normal of the interface in (average) deformed configuration (pointing neighbor => me) + IPareaNormal(1:3,neighbor_n,neighbor_ip,neighbor_el)) ! normal of the interface in (average) deformed configuration (pointing neighbor => en) normal_neighbor2me = matmul(transpose(neighbor_Fe), normal_neighbor2me_defConf) & / math_det33(neighbor_Fe) ! interface normal in the lattice configuration of my neighbor area = IParea(neighbor_n,neighbor_ip,neighbor_el) * norm2(normal_neighbor2me) @@ -1325,7 +1325,7 @@ function rhoDotFlux(timestep,ph,me,ip,el) do t = 1,4 c = (t + 1) / 2 topp = t + mod(t,2) - mod(t+1,2) - if (neighbor_v0(s,t) * math_inner(m(1:3,s,t), normal_neighbor2me) > 0.0_pReal & ! flux from my neighbor to me == entering flux for me + if (neighbor_v0(s,t) * math_inner(m(1:3,s,t), normal_neighbor2me) > 0.0_pReal & ! flux from my neighbor to en == entering flux for en .and. v0(s,t) * neighbor_v0(s,t) >= 0.0_pReal ) then ! ... only if no sign change in flux density lineLength = neighbor_rhoSgl0(s,t) * neighbor_v0(s,t) & * math_inner(m(1:3,s,t), normal_neighbor2me) * area ! positive line length that wants to enter through this interface @@ -1344,7 +1344,7 @@ function rhoDotFlux(timestep,ph,me,ip,el) !* FLUX FROM ME TO MY NEIGHBOR !* This is not considered, if my opposite neighbor has a different constitutive law than nonlocal (still considered for nonlocal law with local properties). - !* Then, we assume, that the opposite(!) neighbor sends an equal amount of dislocations to me. + !* Then, we assume, that the opposite(!) neighbor sends an equal amount of dislocations to en. !* So the net flux in the direction of my neighbor is equal to zero: !* leaving flux to neighbor == entering flux from opposite neighbor !* In case of reduced transmissivity, part of the leaving flux is stored as dead dislocation density. @@ -1353,7 +1353,7 @@ function rhoDotFlux(timestep,ph,me,ip,el) if (phase_plasticity(material_phaseAt(1,opposite_el)) == PLASTICITY_NONLOCAL_ID) then normal_me2neighbor_defConf = math_det33(Favg) & - * matmul(math_inv33(transpose(Favg)),IPareaNormal(1:3,n,ip,el)) ! normal of the interface in (average) deformed configuration (pointing me => neighbor) + * matmul(math_inv33(transpose(Favg)),IPareaNormal(1:3,n,ip,el)) ! normal of the interface in (average) deformed configuration (pointing en => neighbor) normal_me2neighbor = matmul(transpose(my_Fe), normal_me2neighbor_defConf) & / math_det33(my_Fe) ! interface normal in my lattice configuration area = IParea(n,ip,el) * norm2(normal_me2neighbor) @@ -1361,7 +1361,7 @@ function rhoDotFlux(timestep,ph,me,ip,el) do s = 1,ns do t = 1,4 c = (t + 1) / 2 - if (v0(s,t) * math_inner(m(1:3,s,t), normal_me2neighbor) > 0.0_pReal ) then ! flux from me to my neighbor == leaving flux for me (might also be a pure flux from my mobile density to dead density if interface not at all transmissive) + if (v0(s,t) * math_inner(m(1:3,s,t), normal_me2neighbor) > 0.0_pReal ) then ! flux from en to my neighbor == leaving flux for en (might also be a pure flux from my mobile density to dead density if interface not at all transmissive) if (v0(s,t) * neighbor_v0(s,t) >= 0.0_pReal) then ! no sign change in flux density transmissivity = sum(compatibility(c,:,s,n,ip,el)**2.0_pReal) ! overall transmissivity from this slip system to my neighbor else ! sign change in flux density means sign change in stress which does not allow for dislocations to arive at the neighbor @@ -1403,13 +1403,13 @@ module subroutine plastic_nonlocal_updateCompatibility(orientation,ph,i,e) integer :: & n, & ! neighbor index - me, & + en, & neighbor_e, & ! element index of my neighbor neighbor_i, & ! integration point index of my neighbor neighbor_me, & neighbor_phase, & ns, & ! number of active slip systems - s1, & ! slip system index (me) + s1, & ! slip system index (en) s2 ! slip system index (my neighbor) real(pReal), dimension(2,param(ph)%sum_N_sl,param(ph)%sum_N_sl,nIPneighbors) :: & my_compatibility ! my_compatibility for current element and ip @@ -1424,7 +1424,7 @@ module subroutine plastic_nonlocal_updateCompatibility(orientation,ph,i,e) associate(prm => param(ph)) ns = prm%sum_N_sl - me = material_phaseMemberAt(1,i,e) + en = material_phaseMemberAt(1,i,e) !*** start out fully compatible my_compatibility = 0.0_pReal forall(s1 = 1:ns) my_compatibility(:,s1,s1,:) = 1.0_pReal @@ -1450,7 +1450,7 @@ module subroutine plastic_nonlocal_updateCompatibility(orientation,ph,i,e) elseif (prm%chi_GB >= 0.0_pReal) then !* GRAIN BOUNDARY ! !* fixed transmissivity for adjacent ips with different texture (only if explicitly given in material.config) - if (any(dNeq(material_orientation0(1,ph,me)%asQuaternion(), & + if (any(dNeq(material_orientation0(1,ph,en)%asQuaternion(), & material_orientation0(1,neighbor_phase,neighbor_me)%asQuaternion())) .and. & (.not. phase_localPlasticity(neighbor_phase))) & forall(s1 = 1:ns) my_compatibility(:,s1,s1,n) = sqrt(prm%chi_GB) @@ -1767,21 +1767,21 @@ end subroutine kinetics !> @brief returns copy of current dislocation densities from state !> @details raw values is rectified !-------------------------------------------------------------------------------------------------- -pure function getRho(ph,me) +pure function getRho(ph,en) - integer, intent(in) :: ph, me + integer, intent(in) :: ph, en real(pReal), dimension(param(ph)%sum_N_sl,10) :: getRho associate(prm => param(ph)) - getRho = reshape(state(ph)%rho(:,me),[prm%sum_N_sl,10]) + getRho = reshape(state(ph)%rho(:,en),[prm%sum_N_sl,10]) ! ensure positive densities (not for imm, they have a sign) getRho(:,mob) = max(getRho(:,mob),0.0_pReal) getRho(:,dip) = max(getRho(:,dip),0.0_pReal) - where(abs(getRho) < max(prm%rho_min/geom(ph)%V_0(me)**(2.0_pReal/3.0_pReal),prm%rho_significant)) & + where(abs(getRho) < max(prm%rho_min/geom(ph)%V_0(en)**(2.0_pReal/3.0_pReal),prm%rho_significant)) & getRho = 0.0_pReal end associate @@ -1793,21 +1793,21 @@ end function getRho !> @brief returns copy of current dislocation densities from state !> @details raw values is rectified !-------------------------------------------------------------------------------------------------- -pure function getRho0(ph,me) +pure function getRho0(ph,en) - integer, intent(in) :: ph, me + integer, intent(in) :: ph, en real(pReal), dimension(param(ph)%sum_N_sl,10) :: getRho0 associate(prm => param(ph)) - getRho0 = reshape(state0(ph)%rho(:,me),[prm%sum_N_sl,10]) + getRho0 = reshape(state0(ph)%rho(:,en),[prm%sum_N_sl,10]) ! ensure positive densities (not for imm, they have a sign) getRho0(:,mob) = max(getRho0(:,mob),0.0_pReal) getRho0(:,dip) = max(getRho0(:,dip),0.0_pReal) - where (abs(getRho0) < max(prm%rho_min/geom(ph)%V_0(me)**(2.0_pReal/3.0_pReal),prm%rho_significant)) & + where (abs(getRho0) < max(prm%rho_min/geom(ph)%V_0(en)**(2.0_pReal/3.0_pReal),prm%rho_significant)) & getRho0 = 0.0_pReal end associate diff --git a/src/phase_mechanical_plastic_phenopowerlaw.f90 b/src/phase_mechanical_plastic_phenopowerlaw.f90 index 77c47907f..b2db9e346 100644 --- a/src/phase_mechanical_plastic_phenopowerlaw.f90 +++ b/src/phase_mechanical_plastic_phenopowerlaw.f90 @@ -283,7 +283,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,ph,me) +pure module subroutine phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,en) real(pReal), dimension(3,3), intent(out) :: & Lp !< plastic velocity gradient @@ -294,7 +294,7 @@ pure module subroutine phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,me) Mp !< Mandel stress integer, intent(in) :: & ph, & - me + en integer :: & i,k,l,m,n @@ -309,7 +309,7 @@ pure module subroutine phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,me) associate(prm => param(ph)) - call kinetics_slip(Mp,ph,me,gdot_slip_pos,gdot_slip_neg,dgdot_dtauslip_pos,dgdot_dtauslip_neg) + call kinetics_slip(Mp,ph,en,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) & @@ -318,7 +318,7 @@ pure module subroutine phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,me) + dgdot_dtauslip_neg(i) * prm%P_sl(k,l,i) * prm%nonSchmid_neg(m,n,i) enddo slipSystems - call kinetics_twin(Mp,ph,me,gdot_twin,dgdot_dtautwin) + call kinetics_twin(Mp,ph,en,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) & @@ -334,13 +334,13 @@ end subroutine phenopowerlaw_LpAndItsTangent !-------------------------------------------------------------------------------------------------- !> @brief Calculate the rate of change of microstructure. !-------------------------------------------------------------------------------------------------- -module subroutine phenopowerlaw_dotState(Mp,ph,me) +module subroutine phenopowerlaw_dotState(Mp,ph,en) real(pReal), dimension(3,3), intent(in) :: & Mp !< Mandel stress integer, intent(in) :: & ph, & - me + en real(pReal) :: & c_SlipSlip,c_TwinSlip,c_TwinTwin, & @@ -353,8 +353,8 @@ module subroutine phenopowerlaw_dotState(Mp,ph,me) associate(prm => param(ph), stt => state(ph), & dot => dotState(ph)) - sumGamma = sum(stt%gamma_slip(:,me)) - sumF = sum(stt%gamma_twin(:,me)/prm%gamma_char) + sumGamma = sum(stt%gamma_slip(:,en)) + sumF = sum(stt%gamma_twin(:,en)/prm%gamma_char) !-------------------------------------------------------------------------------------------------- ! system-independent (nonlinear) prefactors to M_Xx (X influenced by x) matrices @@ -366,23 +366,23 @@ module subroutine phenopowerlaw_dotState(Mp,ph,me) ! calculate left and right vectors left_SlipSlip = 1.0_pReal + prm%h_int xi_slip_sat_offset = prm%f_sat_sl_tw*sqrt(sumF) - right_SlipSlip = abs(1.0_pReal-stt%xi_slip(:,me) / (prm%xi_inf_sl+xi_slip_sat_offset)) **prm%a_sl & - * sign(1.0_pReal,1.0_pReal-stt%xi_slip(:,me) / (prm%xi_inf_sl+xi_slip_sat_offset)) + right_SlipSlip = abs(1.0_pReal-stt%xi_slip(:,en) / (prm%xi_inf_sl+xi_slip_sat_offset)) **prm%a_sl & + * sign(1.0_pReal,1.0_pReal-stt%xi_slip(:,en) / (prm%xi_inf_sl+xi_slip_sat_offset)) !-------------------------------------------------------------------------------------------------- ! shear rates - call kinetics_slip(Mp,ph,me,gdot_slip_pos,gdot_slip_neg) - dot%gamma_slip(:,me) = abs(gdot_slip_pos+gdot_slip_neg) - call kinetics_twin(Mp,ph,me,dot%gamma_twin(:,me)) + call kinetics_slip(Mp,ph,en,gdot_slip_pos,gdot_slip_neg) + dot%gamma_slip(:,en) = abs(gdot_slip_pos+gdot_slip_neg) + call kinetics_twin(Mp,ph,en,dot%gamma_twin(:,en)) !-------------------------------------------------------------------------------------------------- ! hardening - dot%xi_slip(:,me) = c_SlipSlip * left_SlipSlip * & - matmul(prm%h_sl_sl,dot%gamma_slip(:,me)*right_SlipSlip) & - + matmul(prm%h_sl_tw,dot%gamma_twin(:,me)) + dot%xi_slip(:,en) = c_SlipSlip * left_SlipSlip * & + matmul(prm%h_sl_sl,dot%gamma_slip(:,en)*right_SlipSlip) & + + matmul(prm%h_sl_tw,dot%gamma_twin(:,en)) - dot%xi_twin(:,me) = c_TwinSlip * matmul(prm%h_tw_sl,dot%gamma_slip(:,me)) & - + c_TwinTwin * matmul(prm%h_tw_tw,dot%gamma_twin(:,me)) + dot%xi_twin(:,en) = c_TwinSlip * matmul(prm%h_tw_sl,dot%gamma_slip(:,en)) & + + c_TwinTwin * matmul(prm%h_tw_tw,dot%gamma_twin(:,en)) end associate end subroutine phenopowerlaw_dotState @@ -430,14 +430,14 @@ end subroutine plastic_phenopowerlaw_results ! NOTE: Against the common convention, the result (i.e. intent(out)) variables are the last to ! have the optional arguments at the end. !-------------------------------------------------------------------------------------------------- -pure subroutine kinetics_slip(Mp,ph,me, & +pure subroutine kinetics_slip(Mp,ph,en, & gdot_slip_pos,gdot_slip_neg,dgdot_dtau_slip_pos,dgdot_dtau_slip_neg) real(pReal), dimension(3,3), intent(in) :: & Mp !< Mandel stress integer, intent(in) :: & ph, & - me + en real(pReal), intent(out), dimension(param(ph)%sum_N_sl) :: & gdot_slip_pos, & @@ -461,14 +461,14 @@ pure subroutine kinetics_slip(Mp,ph,me, & where(dNeq0(tau_slip_pos)) gdot_slip_pos = prm%dot_gamma_0_sl * merge(0.5_pReal,1.0_pReal, prm%nonSchmidActive) & ! 1/2 if non-Schmid active - * sign(abs(tau_slip_pos/stt%xi_slip(:,me))**prm%n_sl, tau_slip_pos) + * sign(abs(tau_slip_pos/stt%xi_slip(:,en))**prm%n_sl, tau_slip_pos) else where gdot_slip_pos = 0.0_pReal end where where(dNeq0(tau_slip_neg)) gdot_slip_neg = prm%dot_gamma_0_sl * 0.5_pReal & ! only used if non-Schmid active, always 1/2 - * sign(abs(tau_slip_neg/stt%xi_slip(:,me))**prm%n_sl, tau_slip_neg) + * sign(abs(tau_slip_neg/stt%xi_slip(:,en))**prm%n_sl, tau_slip_neg) else where gdot_slip_neg = 0.0_pReal end where @@ -499,14 +499,14 @@ end subroutine kinetics_slip ! NOTE: Against the common convention, the result (i.e. intent(out)) variables are the last to ! have the optional arguments at the end. !-------------------------------------------------------------------------------------------------- -pure subroutine kinetics_twin(Mp,ph,me,& +pure subroutine kinetics_twin(Mp,ph,en,& gdot_twin,dgdot_dtau_twin) real(pReal), dimension(3,3), intent(in) :: & Mp !< Mandel stress integer, intent(in) :: & ph, & - me + en real(pReal), dimension(param(ph)%sum_N_tw), intent(out) :: & gdot_twin @@ -524,8 +524,8 @@ pure subroutine kinetics_twin(Mp,ph,me,& enddo where(tau_twin > 0.0_pReal) - gdot_twin = (1.0_pReal-sum(stt%gamma_twin(:,me)/prm%gamma_char)) & ! only twin in untwinned volume fraction - * prm%dot_gamma_0_tw*(abs(tau_twin)/stt%xi_twin(:,me))**prm%n_tw + gdot_twin = (1.0_pReal-sum(stt%gamma_twin(:,en)/prm%gamma_char)) & ! only twin in untwinned volume fraction + * prm%dot_gamma_0_tw*(abs(tau_twin)/stt%xi_twin(:,en))**prm%n_tw else where gdot_twin = 0.0_pReal end where diff --git a/src/phase_thermal.f90 b/src/phase_thermal.f90 index b51ac74f5..8d6ccaf6f 100644 --- a/src/phase_thermal.f90 +++ b/src/phase_thermal.f90 @@ -26,7 +26,7 @@ submodule(phase) thermal integer(kind(THERMAL_UNDEFINED_ID)), dimension(:,:), allocatable :: & thermal_source - type(tDataContainer), dimension(:), allocatable :: current ! ?? not very telling name. Better: "field" ?? MD: current(ho)%T(me) reads quite good + type(tDataContainer), dimension(:), allocatable :: current ! ?? not very telling name. Better: "field" ?? MD: current(ho)%T(en) reads quite good type(tThermalParameters), dimension(:), allocatable :: param @@ -46,23 +46,23 @@ submodule(phase) thermal end function externalheat_init - module subroutine externalheat_dotState(ph, me) + module subroutine externalheat_dotState(ph, en) integer, intent(in) :: & ph, & - me + en end subroutine externalheat_dotState - module function dissipation_f_T(ph,me) result(f_T) + module function dissipation_f_T(ph,en) result(f_T) integer, intent(in) :: & ph, & - me + en real(pReal) :: f_T end function dissipation_f_T - module function externalheat_f_T(ph,me) result(f_T) + module function externalheat_f_T(ph,en) result(f_T) integer, intent(in) :: & ph, & - me + en real(pReal) :: f_T end function externalheat_f_T @@ -137,9 +137,9 @@ end subroutine thermal_init !---------------------------------------------------------------------------------------------- !< @brief calculates thermal dissipation rate !---------------------------------------------------------------------------------------------- -module function phase_f_T(ph,me) result(f) +module function phase_f_T(ph,en) result(f) - integer, intent(in) :: ph, me + integer, intent(in) :: ph, en real(pReal) :: f @@ -152,10 +152,10 @@ module function phase_f_T(ph,me) result(f) select case(thermal_source(so,ph)) case (THERMAL_DISSIPATION_ID) - f = f + dissipation_f_T(ph,me) + f = f + dissipation_f_T(ph,en) case (THERMAL_EXTERNALHEAT_ID) - f = f + externalheat_f_T(ph,me) + f = f + externalheat_f_T(ph,en) end select @@ -167,9 +167,9 @@ end function phase_f_T !-------------------------------------------------------------------------------------------------- !> @brief contains the constitutive equation for calculating the rate of change of microstructure !-------------------------------------------------------------------------------------------------- -function phase_thermal_collectDotState(ph,me) result(broken) +function phase_thermal_collectDotState(ph,en) result(broken) - integer, intent(in) :: ph, me + integer, intent(in) :: ph, en logical :: broken integer :: i @@ -180,9 +180,9 @@ function phase_thermal_collectDotState(ph,me) result(broken) SourceLoop: do i = 1, thermal_Nsources(ph) if (thermal_source(i,ph) == THERMAL_EXTERNALHEAT_ID) & - call externalheat_dotState(ph,me) + call externalheat_dotState(ph,en) - broken = broken .or. any(IEEE_is_NaN(thermalState(ph)%p(i)%dotState(:,me))) + broken = broken .or. any(IEEE_is_NaN(thermalState(ph)%p(i)%dotState(:,en))) enddo SourceLoop @@ -218,14 +218,14 @@ module function phase_K_T(co,ce) result(K) end function phase_K_T -module function thermal_stress(Delta_t,ph,me) result(converged_) ! ?? why is this called "stress" when it seems closer to "updateState" ?? +module function thermal_stress(Delta_t,ph,en) result(converged_) ! ?? why is this called "stress" when it seems closer to "updateState" ?? real(pReal), intent(in) :: Delta_t - integer, intent(in) :: ph, me + integer, intent(in) :: ph, en logical :: converged_ - converged_ = .not. integrateThermalState(Delta_t,ph,me) + converged_ = .not. integrateThermalState(Delta_t,ph,en) end function thermal_stress @@ -233,10 +233,10 @@ end function thermal_stress !-------------------------------------------------------------------------------------------------- !> @brief integrate state with 1st order explicit Euler method !-------------------------------------------------------------------------------------------------- -function integrateThermalState(Delta_t, ph,me) result(broken) +function integrateThermalState(Delta_t, ph,en) result(broken) real(pReal), intent(in) :: Delta_t - integer, intent(in) :: ph, me + integer, intent(in) :: ph, en logical :: & broken @@ -244,13 +244,13 @@ function integrateThermalState(Delta_t, ph,me) result(broken) so, & sizeDotState - broken = phase_thermal_collectDotState(ph,me) + broken = phase_thermal_collectDotState(ph,en) if (broken) return do so = 1, thermal_Nsources(ph) sizeDotState = thermalState(ph)%p(so)%sizeDotState - thermalState(ph)%p(so)%state(1:sizeDotState,me) = thermalState(ph)%p(so)%state0(1:sizeDotState,me) & - + thermalState(ph)%p(so)%dotState(1:sizeDotState,me) * Delta_t + thermalState(ph)%p(so)%state(1:sizeDotState,en) = thermalState(ph)%p(so)%state0(1:sizeDotState,en) & + + thermalState(ph)%p(so)%dotState(1:sizeDotState,en) * Delta_t enddo end function integrateThermalState @@ -273,13 +273,13 @@ end subroutine thermal_forward !---------------------------------------------------------------------------------------------- !< @brief Get temperature (for use by non-thermal physics) !---------------------------------------------------------------------------------------------- -module function thermal_T(ph,me) result(T) +module function thermal_T(ph,en) result(T) - integer, intent(in) :: ph, me + integer, intent(in) :: ph, en real(pReal) :: T - T = current(ph)%T(me) + T = current(ph)%T(en) end function thermal_T @@ -287,13 +287,13 @@ end function thermal_T !---------------------------------------------------------------------------------------------- !< @brief Get rate of temperature (for use by non-thermal physics) !---------------------------------------------------------------------------------------------- -module function thermal_dot_T(ph,me) result(dot_T) +module function thermal_dot_T(ph,en) result(dot_T) - integer, intent(in) :: ph, me + integer, intent(in) :: ph, en real(pReal) :: dot_T - dot_T = current(ph)%dot_T(me) + dot_T = current(ph)%dot_T(en) end function thermal_dot_T diff --git a/src/phase_thermal_dissipation.f90 b/src/phase_thermal_dissipation.f90 index 3a4ee651a..cbdf5bb99 100644 --- a/src/phase_thermal_dissipation.f90 +++ b/src/phase_thermal_dissipation.f90 @@ -69,15 +69,15 @@ end function dissipation_init !-------------------------------------------------------------------------------------------------- !> @brief Ninstancess dissipation rate !-------------------------------------------------------------------------------------------------- -module function dissipation_f_T(ph,me) result(f_T) +module function dissipation_f_T(ph,en) result(f_T) - integer, intent(in) :: ph, me + integer, intent(in) :: ph, en real(pReal) :: & f_T associate(prm => param(ph)) - f_T = prm%kappa*sum(abs(mechanical_S(ph,me)*mechanical_L_p(ph,me))) + f_T = prm%kappa*sum(abs(mechanical_S(ph,en)*mechanical_L_p(ph,en))) end associate end function dissipation_f_T diff --git a/src/phase_thermal_externalheat.f90 b/src/phase_thermal_externalheat.f90 index 6d4403ab8..46da127f3 100644 --- a/src/phase_thermal_externalheat.f90 +++ b/src/phase_thermal_externalheat.f90 @@ -81,18 +81,18 @@ end function externalheat_init !> @brief rate of change of state !> @details state only contains current time to linearly interpolate given heat powers !-------------------------------------------------------------------------------------------------- -module subroutine externalheat_dotState(ph, me) +module subroutine externalheat_dotState(ph, en) integer, intent(in) :: & ph, & - me + en integer :: & so so = source_thermal_externalheat_offset(ph) - thermalState(ph)%p(so)%dotState(1,me) = 1.0_pReal ! state is current time + thermalState(ph)%p(so)%dotState(1,en) = 1.0_pReal ! state is current time end subroutine externalheat_dotState @@ -100,11 +100,11 @@ end subroutine externalheat_dotState !-------------------------------------------------------------------------------------------------- !> @brief returns local heat generation rate !-------------------------------------------------------------------------------------------------- -module function externalheat_f_T(ph,me) result(f_T) +module function externalheat_f_T(ph,en) result(f_T) integer, intent(in) :: & ph, & - me + en real(pReal) :: & f_T @@ -117,14 +117,14 @@ module function externalheat_f_T(ph,me) result(f_T) associate(prm => param(ph)) do interval = 1, prm%nIntervals ! scan through all rate segments - frac_time = (thermalState(ph)%p(so)%state(1,me) - prm%t_n(interval)) & + frac_time = (thermalState(ph)%p(so)%state(1,en) - prm%t_n(interval)) & / (prm%t_n(interval+1) - prm%t_n(interval)) ! fractional time within segment if ( (frac_time < 0.0_pReal .and. interval == 1) & .or. (frac_time >= 1.0_pReal .and. interval == prm%nIntervals) & .or. (frac_time >= 0.0_pReal .and. frac_time < 1.0_pReal) ) & f_T = prm%f_T(interval ) * (1.0_pReal - frac_time) + & prm%f_T(interval+1) * frac_time ! interpolate heat rate between segment boundaries... - ! ...or extrapolate if outside me bounds + ! ...or extrapolate if outside of bounds enddo end associate