From 5af6cc288b6e892f594c611689d1fc42d5f52181 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 19 Dec 2021 21:46:10 +0100 Subject: [PATCH 1/3] whitespace adjustments --- src/phase_mechanical_plastic_nonlocal.f90 | 313 +++++++++++----------- 1 file changed, 158 insertions(+), 155 deletions(-) diff --git a/src/phase_mechanical_plastic_nonlocal.f90 b/src/phase_mechanical_plastic_nonlocal.f90 index bf0dc3b19..7ff016514 100644 --- a/src/phase_mechanical_plastic_nonlocal.f90 +++ b/src/phase_mechanical_plastic_nonlocal.f90 @@ -772,66 +772,67 @@ module subroutine nonlocal_LpAndItsTangent(Lp,dLp_dMp, & tau, & !< resolved shear stress including backstress terms dot_gamma !< shear rate + associate(prm => param(ph),dst=>dependentState(ph),stt=>state(ph)) - !*** shortcut to state variables - rho = getRho(ph,en) - rhoSgl = rho(:,sgl) + !*** shortcut to state variables + rho = getRho(ph,en) + rhoSgl = rho(:,sgl) - do s = 1,prm%sum_N_sl - tau(s) = math_tensordot(Mp, prm%P_sl(1:3,1:3,s)) - tauNS(s,1) = tau(s) - tauNS(s,2) = tau(s) - if (tau(s) > 0.0_pReal) then - tauNS(s,3) = math_tensordot(Mp, +prm%P_nS_pos(1:3,1:3,s)) - tauNS(s,4) = math_tensordot(Mp, -prm%P_nS_neg(1:3,1:3,s)) - else - tauNS(s,3) = math_tensordot(Mp, +prm%P_nS_neg(1:3,1:3,s)) - tauNS(s,4) = math_tensordot(Mp, -prm%P_nS_pos(1:3,1:3,s)) - end if - end do - 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(:,en),1,Temperature, ph) - v(:,2) = v(:,1) - dv_dtau(:,2) = dv_dtau(:,1) - dv_dtauNS(:,2) = dv_dtauNS(:,1) - - !screws - if (prm%nonSchmidActive) then - do t = 3,4 - call kinetics(v(:,t), dv_dtau(:,t), dv_dtauNS(:,t), & - tau, tauNS(:,t), dst%tau_pass(:,en),2,Temperature, ph) + do s = 1,prm%sum_N_sl + tau(s) = math_tensordot(Mp, prm%P_sl(1:3,1:3,s)) + tauNS(s,1) = tau(s) + tauNS(s,2) = tau(s) + if (tau(s) > 0.0_pReal) then + tauNS(s,3) = math_tensordot(Mp, +prm%P_nS_pos(1:3,1:3,s)) + tauNS(s,4) = math_tensordot(Mp, -prm%P_nS_neg(1:3,1:3,s)) + else + tauNS(s,3) = math_tensordot(Mp, +prm%P_nS_neg(1:3,1:3,s)) + tauNS(s,4) = math_tensordot(Mp, -prm%P_nS_pos(1:3,1:3,s)) + end if end do - else - v(:,3:4) = spread(v(:,1),2,2) - dv_dtau(:,3:4) = spread(dv_dtau(:,1),2,2) - dv_dtauNS(:,3:4) = spread(dv_dtauNS(:,1),2,2) - end if + tauNS = tauNS + spread(dst%tau_back(:,en),2,4) + tau = tau + dst%tau_back(:,en) - stt%v(:,en) = pack(v,.true.) + ! edges + call kinetics(v(:,1), dv_dtau(:,1), dv_dtauNS(:,1), & + 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) - !*** Bauschinger effect - forall (s = 1:prm%sum_N_sl, t = 5:8, rhoSgl(s,t) * v(s,t-4) < 0.0_pReal) & - rhoSgl(s,t-4) = rhoSgl(s,t-4) + abs(rhoSgl(s,t)) + !screws + if (prm%nonSchmidActive) then + do t = 3,4 + call kinetics(v(:,t), dv_dtau(:,t), dv_dtauNS(:,t), & + tau, tauNS(:,t), dst%tau_pass(:,en),2,Temperature, ph) + end do + else + v(:,3:4) = spread(v(:,1),2,2) + dv_dtau(:,3:4) = spread(dv_dtau(:,1),2,2) + dv_dtauNS(:,3:4) = spread(dv_dtauNS(:,1),2,2) + end if - dot_gamma = sum(rhoSgl(:,1:4) * v, 2) * prm%b_sl + stt%v(:,en) = pack(v,.true.) - Lp = 0.0_pReal - dLp_dMp = 0.0_pReal - do s = 1,prm%sum_N_sl - Lp = Lp + dot_gamma(s) * prm%P_sl(1:3,1:3,s) - forall (i=1:3,j=1:3,k=1:3,l=1:3) & - dLp_dMp(i,j,k,l) = dLp_dMp(i,j,k,l) & - + prm%P_sl(i,j,s) * prm%P_sl(k,l,s) & - * sum(rhoSgl(s,1:4) * dv_dtau(s,1:4)) * prm%b_sl(s) & - + prm%P_sl(i,j,s) & - * (+ prm%P_nS_pos(k,l,s) * rhoSgl(s,3) * dv_dtauNS(s,3) & - - prm%P_nS_neg(k,l,s) * rhoSgl(s,4) * dv_dtauNS(s,4)) * prm%b_sl(s) - end do + !*** Bauschinger effect + forall (s = 1:prm%sum_N_sl, t = 5:8, rhoSgl(s,t) * v(s,t-4) < 0.0_pReal) & + rhoSgl(s,t-4) = rhoSgl(s,t-4) + abs(rhoSgl(s,t)) + + dot_gamma = sum(rhoSgl(:,1:4) * v, 2) * prm%b_sl + + Lp = 0.0_pReal + dLp_dMp = 0.0_pReal + do s = 1,prm%sum_N_sl + Lp = Lp + dot_gamma(s) * prm%P_sl(1:3,1:3,s) + forall (i=1:3,j=1:3,k=1:3,l=1:3) & + dLp_dMp(i,j,k,l) = dLp_dMp(i,j,k,l) & + + prm%P_sl(i,j,s) * prm%P_sl(k,l,s) & + * sum(rhoSgl(s,1:4) * dv_dtau(s,1:4)) * prm%b_sl(s) & + + prm%P_sl(i,j,s) & + * (+ prm%P_nS_pos(k,l,s) * rhoSgl(s,3) * dv_dtauNS(s,3) & + - prm%P_nS_neg(k,l,s) * rhoSgl(s,4) * dv_dtauNS(s,4)) * prm%b_sl(s) + end do end associate @@ -870,7 +871,8 @@ module subroutine plastic_nonlocal_deltaState(Mp,ph,en) dUpper, & ! current maximum stable dipole distance for edges and screws dUpperOld, & ! old maximum stable dipole distance for edges and screws deltaDUpper ! change in maximum stable dipole distance for edges and screws - + + associate(prm => param(ph),dst => dependentState(ph),del => deltaState(ph)) mu = elastic_mu(ph,en) @@ -1394,78 +1396,79 @@ module subroutine plastic_nonlocal_updateCompatibility(orientation,ph,i,e) belowThreshold type(rotation) :: mis + associate(prm => param(ph)) - ns = prm%sum_N_sl + ns = prm%sum_N_sl - 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 + 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 - neighbors: do n = 1,nIPneighbors - neighbor_e = IPneighborhood(1,n,i,e) - neighbor_i = IPneighborhood(2,n,i,e) - neighbor_me = material_phaseMemberAt(1,neighbor_i,neighbor_e) - neighbor_phase = material_phaseAt(1,neighbor_e) + neighbors: do n = 1,nIPneighbors + neighbor_e = IPneighborhood(1,n,i,e) + neighbor_i = IPneighborhood(2,n,i,e) + neighbor_me = material_phaseMemberAt(1,neighbor_i,neighbor_e) + neighbor_phase = material_phaseAt(1,neighbor_e) - if (neighbor_e <= 0 .or. neighbor_i <= 0) then - !* FREE SURFACE - forall(s1 = 1:ns) my_compatibility(:,s1,s1,n) = sqrt(prm%chi_surface) - elseif (neighbor_phase /= ph) then - !* PHASE BOUNDARY - if (plasticState(neighbor_phase)%nonlocal .and. plasticState(ph)%nonlocal) & - forall(s1 = 1:ns) my_compatibility(:,s1,s1,n) = 0.0_pReal - elseif (prm%chi_GB >= 0.0_pReal) then - !* GRAIN BOUNDARY - if (any(dNeq(phase_O_0(ph)%data(en)%asQuaternion(), & - phase_O_0(neighbor_phase)%data(neighbor_me)%asQuaternion())) .and. & - plasticState(neighbor_phase)%nonlocal) & - forall(s1 = 1:ns) my_compatibility(:,s1,s1,n) = sqrt(prm%chi_GB) - else - !* GRAIN BOUNDARY ? - !* Compatibility defined by relative orientation of slip systems: - !* The my_compatibility value is defined as the product of the slip normal projection and the slip direction projection. - !* Its sign is always positive for screws, for edges it has the same sign as the slip normal projection. - !* Since the sum for each slip system can easily exceed one (which would result in a transmissivity larger than one), - !* only values above or equal to a certain threshold value are considered. This threshold value is chosen, such that - !* the number of compatible slip systems is minimized with the sum of the original compatibility values exceeding one. - !* Finally the smallest compatibility value is decreased until the sum is exactly equal to one. - !* All values below the threshold are set to zero. - mis = orientation(ph)%data(en)%misorientation(orientation(neighbor_phase)%data(neighbor_me)) - mySlipSystems: do s1 = 1,ns - neighborSlipSystems: do s2 = 1,ns - my_compatibility(1,s2,s1,n) = math_inner(prm%slip_normal(1:3,s1), & - mis%rotate(prm%slip_normal(1:3,s2))) & - * abs(math_inner(prm%slip_direction(1:3,s1), & - mis%rotate(prm%slip_direction(1:3,s2)))) - my_compatibility(2,s2,s1,n) = abs(math_inner(prm%slip_normal(1:3,s1), & - mis%rotate(prm%slip_normal(1:3,s2)))) & - * abs(math_inner(prm%slip_direction(1:3,s1), & - mis%rotate(prm%slip_direction(1:3,s2)))) - end do neighborSlipSystems + if (neighbor_e <= 0 .or. neighbor_i <= 0) then + !* FREE SURFACE + forall(s1 = 1:ns) my_compatibility(:,s1,s1,n) = sqrt(prm%chi_surface) + elseif (neighbor_phase /= ph) then + !* PHASE BOUNDARY + if (plasticState(neighbor_phase)%nonlocal .and. plasticState(ph)%nonlocal) & + forall(s1 = 1:ns) my_compatibility(:,s1,s1,n) = 0.0_pReal + elseif (prm%chi_GB >= 0.0_pReal) then + !* GRAIN BOUNDARY + if (any(dNeq(phase_O_0(ph)%data(en)%asQuaternion(), & + phase_O_0(neighbor_phase)%data(neighbor_me)%asQuaternion())) .and. & + plasticState(neighbor_phase)%nonlocal) & + forall(s1 = 1:ns) my_compatibility(:,s1,s1,n) = sqrt(prm%chi_GB) + else + !* GRAIN BOUNDARY ? + !* Compatibility defined by relative orientation of slip systems: + !* The my_compatibility value is defined as the product of the slip normal projection and the slip direction projection. + !* Its sign is always positive for screws, for edges it has the same sign as the slip normal projection. + !* Since the sum for each slip system can easily exceed one (which would result in a transmissivity larger than one), + !* only values above or equal to a certain threshold value are considered. This threshold value is chosen, such that + !* the number of compatible slip systems is minimized with the sum of the original compatibility values exceeding one. + !* Finally the smallest compatibility value is decreased until the sum is exactly equal to one. + !* All values below the threshold are set to zero. + mis = orientation(ph)%data(en)%misorientation(orientation(neighbor_phase)%data(neighbor_me)) + mySlipSystems: do s1 = 1,ns + neighborSlipSystems: do s2 = 1,ns + my_compatibility(1,s2,s1,n) = math_inner(prm%slip_normal(1:3,s1), & + mis%rotate(prm%slip_normal(1:3,s2))) & + * abs(math_inner(prm%slip_direction(1:3,s1), & + mis%rotate(prm%slip_direction(1:3,s2)))) + my_compatibility(2,s2,s1,n) = abs(math_inner(prm%slip_normal(1:3,s1), & + mis%rotate(prm%slip_normal(1:3,s2)))) & + * abs(math_inner(prm%slip_direction(1:3,s1), & + mis%rotate(prm%slip_direction(1:3,s2)))) + end do neighborSlipSystems - my_compatibilitySum = 0.0_pReal - belowThreshold = .true. - do while (my_compatibilitySum < 1.0_pReal .and. any(belowThreshold)) - thresholdValue = maxval(my_compatibility(2,:,s1,n), belowThreshold) ! screws always positive - nThresholdValues = real(count(my_compatibility(2,:,s1,n) >= thresholdValue),pReal) - where (my_compatibility(2,:,s1,n) >= thresholdValue) belowThreshold = .false. - if (my_compatibilitySum + thresholdValue * nThresholdValues > 1.0_pReal) & - where (abs(my_compatibility(:,:,s1,n)) >= thresholdValue) & - my_compatibility(:,:,s1,n) = sign((1.0_pReal - my_compatibilitySum)/nThresholdValues,& - my_compatibility(:,:,s1,n)) - my_compatibilitySum = my_compatibilitySum + nThresholdValues * thresholdValue - end do + my_compatibilitySum = 0.0_pReal + belowThreshold = .true. + do while (my_compatibilitySum < 1.0_pReal .and. any(belowThreshold)) + thresholdValue = maxval(my_compatibility(2,:,s1,n), belowThreshold) ! screws always positive + nThresholdValues = real(count(my_compatibility(2,:,s1,n) >= thresholdValue),pReal) + where (my_compatibility(2,:,s1,n) >= thresholdValue) belowThreshold = .false. + if (my_compatibilitySum + thresholdValue * nThresholdValues > 1.0_pReal) & + where (abs(my_compatibility(:,:,s1,n)) >= thresholdValue) & + my_compatibility(:,:,s1,n) = sign((1.0_pReal - my_compatibilitySum)/nThresholdValues,& + my_compatibility(:,:,s1,n)) + my_compatibilitySum = my_compatibilitySum + nThresholdValues * thresholdValue + end do - where(belowThreshold) my_compatibility(1,:,s1,n) = 0.0_pReal - where(belowThreshold) my_compatibility(2,:,s1,n) = 0.0_pReal + where(belowThreshold) my_compatibility(1,:,s1,n) = 0.0_pReal + where(belowThreshold) my_compatibility(2,:,s1,n) = 0.0_pReal - end do mySlipSystems - end if + end do mySlipSystems + end if - end do neighbors + end do neighbors - compatibility(:,:,:,:,i,e) = my_compatibility + compatibility(:,:,:,:,i,e) = my_compatibility end associate @@ -1647,52 +1650,52 @@ pure subroutine kinetics(v, dv_dtau, dv_dtauNS, tau, tauNS, tauThreshold, c, T, criticalStress_S !< maximum obstacle strength associate(prm => param(ph)) - v = 0.0_pReal - dv_dtau = 0.0_pReal - dv_dtauNS = 0.0_pReal + v = 0.0_pReal + dv_dtau = 0.0_pReal + dv_dtauNS = 0.0_pReal - do s = 1,prm%sum_N_sl - if (abs(tau(s)) > tauThreshold(s)) then + do s = 1,prm%sum_N_sl + if (abs(tau(s)) > tauThreshold(s)) then - !* Peierls contribution - tauEff = max(0.0_pReal, abs(tauNS(s)) - tauThreshold(s)) - lambda_P = prm%b_sl(s) - activationVolume_P = prm%w *prm%b_sl(s)**3 - criticalStress_P = prm%peierlsStress(s,c) - activationEnergy_P = criticalStress_P * activationVolume_P - tauRel_P = min(1.0_pReal, tauEff / criticalStress_P) - tPeierls = 1.0_pReal / prm%nu_a & - * exp(activationEnergy_P / (K_B * T) & - * (1.0_pReal - tauRel_P**prm%p)**prm%q) - dtPeierls_dtau = merge(tPeierls * prm%p * prm%q * activationVolume_P / (K_B * T) & - * (1.0_pReal - tauRel_P**prm%p)**(prm%q-1.0_pReal) * tauRel_P**(prm%p-1.0_pReal), & - 0.0_pReal, & - tauEff < criticalStress_P) + !* Peierls contribution + tauEff = max(0.0_pReal, abs(tauNS(s)) - tauThreshold(s)) + lambda_P = prm%b_sl(s) + activationVolume_P = prm%w *prm%b_sl(s)**3 + criticalStress_P = prm%peierlsStress(s,c) + activationEnergy_P = criticalStress_P * activationVolume_P + tauRel_P = min(1.0_pReal, tauEff / criticalStress_P) + tPeierls = 1.0_pReal / prm%nu_a & + * exp(activationEnergy_P / (K_B * T) & + * (1.0_pReal - tauRel_P**prm%p)**prm%q) + dtPeierls_dtau = merge(tPeierls * prm%p * prm%q * activationVolume_P / (K_B * T) & + * (1.0_pReal - tauRel_P**prm%p)**(prm%q-1.0_pReal) * tauRel_P**(prm%p-1.0_pReal), & + 0.0_pReal, & + tauEff < criticalStress_P) - ! Contribution from solid solution strengthening - tauEff = abs(tau(s)) - tauThreshold(s) - lambda_S = prm%b_sl(s) / sqrt(prm%c_sol) - activationVolume_S = prm%f_sol * prm%b_sl(s)**3 / sqrt(prm%c_sol) - criticalStress_S = prm%Q_sol / activationVolume_S - tauRel_S = min(1.0_pReal, tauEff / criticalStress_S) - tSolidSolution = 1.0_pReal / prm%nu_a & - * exp(prm%Q_sol / (K_B * T)* (1.0_pReal - tauRel_S**prm%p)**prm%q) - dtSolidSolution_dtau = merge(tSolidSolution * prm%p * prm%q * activationVolume_S / (K_B * T) & - * (1.0_pReal - tauRel_S**prm%p)**(prm%q-1.0_pReal)* tauRel_S**(prm%p-1.0_pReal), & - 0.0_pReal, & - tauEff < criticalStress_S) + ! Contribution from solid solution strengthening + tauEff = abs(tau(s)) - tauThreshold(s) + lambda_S = prm%b_sl(s) / sqrt(prm%c_sol) + activationVolume_S = prm%f_sol * prm%b_sl(s)**3 / sqrt(prm%c_sol) + criticalStress_S = prm%Q_sol / activationVolume_S + tauRel_S = min(1.0_pReal, tauEff / criticalStress_S) + tSolidSolution = 1.0_pReal / prm%nu_a & + * exp(prm%Q_sol / (K_B * T)* (1.0_pReal - tauRel_S**prm%p)**prm%q) + dtSolidSolution_dtau = merge(tSolidSolution * prm%p * prm%q * activationVolume_S / (K_B * T) & + * (1.0_pReal - tauRel_S**prm%p)**(prm%q-1.0_pReal)* tauRel_S**(prm%p-1.0_pReal), & + 0.0_pReal, & + tauEff < criticalStress_S) - !* viscous glide velocity - tauEff = abs(tau(s)) - tauThreshold(s) + !* viscous glide velocity + tauEff = abs(tau(s)) - tauThreshold(s) - v(s) = sign(1.0_pReal,tau(s)) & - / (tPeierls / lambda_P + tSolidSolution / lambda_S + prm%B /(prm%b_sl(s) * tauEff)) - dv_dtau(s) = v(s)**2 * (dtSolidSolution_dtau / lambda_S + prm%B / (prm%b_sl(s) * tauEff**2)) - dv_dtauNS(s) = v(s)**2 * dtPeierls_dtau / lambda_P + v(s) = sign(1.0_pReal,tau(s)) & + / (tPeierls / lambda_P + tSolidSolution / lambda_S + prm%B /(prm%b_sl(s) * tauEff)) + dv_dtau(s) = v(s)**2 * (dtSolidSolution_dtau / lambda_S + prm%B / (prm%b_sl(s) * tauEff**2)) + dv_dtauNS(s) = v(s)**2 * dtPeierls_dtau / lambda_P - end if - end do + end if + end do end associate From 00230d482f30ce8155852d1e201613bbf6b37a06 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 19 Dec 2021 22:07:23 +0100 Subject: [PATCH 2/3] use data from other physics directly more clear code, simplified interfaces --- src/phase.f90 | 2 +- src/phase_mechanical_plastic.f90 | 22 +--- ...phase_mechanical_plastic_dislotungsten.f90 | 8 +- src/phase_mechanical_plastic_dislotwin.f90 | 112 +++++++++--------- src/phase_mechanical_plastic_nonlocal.f90 | 21 ++-- src/phase_thermal.f90 | 2 +- 6 files changed, 81 insertions(+), 86 deletions(-) diff --git a/src/phase.f90 b/src/phase.f90 index 214b0f5fa..6035b4491 100644 --- a/src/phase.f90 +++ b/src/phase.f90 @@ -155,7 +155,7 @@ module phase real(pReal), dimension(3,3) :: P end function phase_P - module function thermal_T(ph,en) result(T) + pure module function thermal_T(ph,en) result(T) integer, intent(in) :: ph,en real(pReal) :: T end function thermal_T diff --git a/src/phase_mechanical_plastic.f90 b/src/phase_mechanical_plastic.f90 index 4951712fd..0a24b0cbf 100644 --- a/src/phase_mechanical_plastic.f90 +++ b/src/phase_mechanical_plastic.f90 @@ -73,47 +73,37 @@ submodule(phase:mechanical) plastic en end subroutine kinehardening_LpAndItsTangent - module subroutine dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,ph,en) + module subroutine dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,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 - real(pReal), intent(in) :: & - T integer, intent(in) :: & ph, & en end subroutine dislotwin_LpAndItsTangent - pure module subroutine dislotungsten_LpAndItsTangent(Lp,dLp_dMp,Mp,T,ph,en) + pure module subroutine dislotungsten_LpAndItsTangent(Lp,dLp_dMp,Mp,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 - real(pReal), intent(in) :: & - T integer, intent(in) :: & ph, & en end subroutine dislotungsten_LpAndItsTangent - module subroutine nonlocal_LpAndItsTangent(Lp,dLp_dMp, & - Mp,Temperature,ph,en) + module subroutine nonlocal_LpAndItsTangent(Lp,dLp_dMp,Mp,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 !< Mandel stress - real(pReal), intent(in) :: & - Temperature integer, intent(in) :: & ph, & en @@ -282,13 +272,13 @@ module subroutine plastic_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, & call kinehardening_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,en) case (PLASTIC_NONLOCAL_ID) plasticType - call nonlocal_LpAndItsTangent(Lp,dLp_dMp,Mp, thermal_T(ph,en),ph,en) + call nonlocal_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,en) case (PLASTIC_DISLOTWIN_ID) plasticType - call dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp, thermal_T(ph,en),ph,en) + call dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,en) case (PLASTIC_DISLOTUNGSTEN_ID) plasticType - call dislotungsten_LpAndItsTangent(Lp,dLp_dMp,Mp, thermal_T(ph,en),ph,en) + call dislotungsten_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,en) end select plasticType diff --git a/src/phase_mechanical_plastic_dislotungsten.f90 b/src/phase_mechanical_plastic_dislotungsten.f90 index c2a584ea1..cd71a7fd7 100644 --- a/src/phase_mechanical_plastic_dislotungsten.f90 +++ b/src/phase_mechanical_plastic_dislotungsten.f90 @@ -257,27 +257,27 @@ end function plastic_dislotungsten_init !> @brief Calculate plastic velocity gradient and its tangent. !-------------------------------------------------------------------------------------------------- pure module subroutine dislotungsten_LpAndItsTangent(Lp,dLp_dMp, & - Mp,T,ph,en) + Mp,ph,en) real(pReal), dimension(3,3), intent(out) :: & Lp !< plastic velocity gradient real(pReal), dimension(3,3,3,3), intent(out) :: & dLp_dMp !< derivative of Lp with respect to the Mandel stress - real(pReal), dimension(3,3), intent(in) :: & Mp !< Mandel stress - real(pReal), intent(in) :: & - T !< temperature integer, intent(in) :: & ph, & en integer :: & i,k,l,m,n + real(pReal) :: & + T !< temperature real(pReal), dimension(param(ph)%sum_N_sl) :: & dot_gamma_pos,dot_gamma_neg, & ddot_gamma_dtau_pos,ddot_gamma_dtau_neg + T = thermal_T(ph,en) Lp = 0.0_pReal dLp_dMp = 0.0_pReal diff --git a/src/phase_mechanical_plastic_dislotwin.f90 b/src/phase_mechanical_plastic_dislotwin.f90 index c5e043e3f..cfc9a002f 100644 --- a/src/phase_mechanical_plastic_dislotwin.f90 +++ b/src/phase_mechanical_plastic_dislotwin.f90 @@ -513,20 +513,20 @@ end function plastic_dislotwin_homogenizedC !-------------------------------------------------------------------------------------------------- !> @brief Calculate plastic velocity gradient and its tangent. !-------------------------------------------------------------------------------------------------- -module subroutine dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,ph,en) +module subroutine dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,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,en - real(pReal), intent(in) :: T integer :: i,k,l,m,n real(pReal) :: & - f_unrotated,StressRatio_p,& - E_kB_T, & - ddot_gamma_dtau, & - tau + f_unrotated,StressRatio_p,& + E_kB_T, & + ddot_gamma_dtau, & + tau, & + T real(pReal), dimension(param(ph)%sum_N_sl) :: & dot_gamma_sl,ddot_gamma_dtau_sl real(pReal), dimension(param(ph)%sum_N_tw) :: & @@ -556,69 +556,71 @@ module subroutine dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,ph,en) 0, 1, 1 & ],pReal),[ 3,6]) - associate(prm => param(ph), stt => state(ph)) - - f_unrotated = 1.0_pReal & - - sum(stt%f_tw(1:prm%sum_N_tw,en)) & - - sum(stt%f_tr(1:prm%sum_N_tr,en)) + T = thermal_T(ph,en) Lp = 0.0_pReal dLp_dMp = 0.0_pReal - call kinetics_sl(Mp,T,ph,en,dot_gamma_sl,ddot_gamma_dtau_sl) - 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) & - dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) & - + ddot_gamma_dtau_sl(i) * prm%P_sl(k,l,i) * prm%P_sl(m,n,i) - end do slipContribution + associate(prm => param(ph), stt => state(ph)) - call kinetics_tw(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) & - dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) & - + ddot_gamma_dtau_tw(i)* prm%P_tw(k,l,i)*prm%P_tw(m,n,i) - end do twinContibution + f_unrotated = 1.0_pReal & + - sum(stt%f_tw(1:prm%sum_N_tw,en)) & + - sum(stt%f_tr(1:prm%sum_N_tr,en)) - call kinetics_tr(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) & - dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) & - + ddot_gamma_dtau_tr(i)* prm%P_tr(k,l,i)*prm%P_tr(m,n,i) - end do transContibution + call kinetics_sl(Mp,T,ph,en,dot_gamma_sl,ddot_gamma_dtau_sl) + 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) & + dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) & + + ddot_gamma_dtau_sl(i) * prm%P_sl(k,l,i) * prm%P_sl(m,n,i) + end do slipContribution - Lp = Lp * f_unrotated - dLp_dMp = dLp_dMp * f_unrotated + call kinetics_tw(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) & + dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) & + + ddot_gamma_dtau_tw(i)* prm%P_tw(k,l,i)*prm%P_tw(m,n,i) + end do twinContibution - shearBandingContribution: if (dNeq0(prm%v_sb)) then + call kinetics_tr(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) & + dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) & + + ddot_gamma_dtau_tr(i)* prm%P_tr(k,l,i)*prm%P_tr(m,n,i) + end do transContibution - E_kB_T = prm%E_sb/(K_B*T) - call math_eigh33(eigValues,eigVectors,Mp) ! is Mp symmetric by design? + Lp = Lp * f_unrotated + dLp_dMp = dLp_dMp * f_unrotated - do i = 1,6 - P_sb = 0.5_pReal * math_outer(matmul(eigVectors,sb_sComposition(1:3,i)),& - matmul(eigVectors,sb_mComposition(1:3,i))) - tau = math_tensordot(Mp,P_sb) + shearBandingContribution: if (dNeq0(prm%v_sb)) then - significantShearBandStress: if (abs(tau) > tol_math_check) then - StressRatio_p = (abs(tau)/prm%xi_sb)**prm%p_sb - dot_gamma_sb = sign(prm%v_sb*exp(-E_kB_T*(1-StressRatio_p)**prm%q_sb), tau) - ddot_gamma_dtau = abs(dot_gamma_sb)*E_kB_T*prm%p_sb*prm%q_sb/prm%xi_sb & - * (abs(tau)/prm%xi_sb)**(prm%p_sb-1.0_pReal) & - * (1.0_pReal-StressRatio_p)**(prm%q_sb-1.0_pReal) + E_kB_T = prm%E_sb/(K_B*T) + call math_eigh33(eigValues,eigVectors,Mp) ! is Mp symmetric by design? - Lp = Lp + dot_gamma_sb * P_sb - forall (k=1:3,l=1:3,m=1:3,n=1:3) & - dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) & - + ddot_gamma_dtau * P_sb(k,l) * P_sb(m,n) - end if significantShearBandStress - end do + do i = 1,6 + P_sb = 0.5_pReal * math_outer(matmul(eigVectors,sb_sComposition(1:3,i)),& + matmul(eigVectors,sb_mComposition(1:3,i))) + tau = math_tensordot(Mp,P_sb) - end if shearBandingContribution + significantShearBandStress: if (abs(tau) > tol_math_check) then + StressRatio_p = (abs(tau)/prm%xi_sb)**prm%p_sb + dot_gamma_sb = sign(prm%v_sb*exp(-E_kB_T*(1-StressRatio_p)**prm%q_sb), tau) + ddot_gamma_dtau = abs(dot_gamma_sb)*E_kB_T*prm%p_sb*prm%q_sb/prm%xi_sb & + * (abs(tau)/prm%xi_sb)**(prm%p_sb-1.0_pReal) & + * (1.0_pReal-StressRatio_p)**(prm%q_sb-1.0_pReal) - end associate + Lp = Lp + dot_gamma_sb * P_sb + forall (k=1:3,l=1:3,m=1:3,n=1:3) & + dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) & + + ddot_gamma_dtau * P_sb(k,l) * P_sb(m,n) + end if significantShearBandStress + end do + + end if shearBandingContribution + + end associate end subroutine dislotwin_LpAndItsTangent diff --git a/src/phase_mechanical_plastic_nonlocal.f90 b/src/phase_mechanical_plastic_nonlocal.f90 index 7ff016514..7118055fe 100644 --- a/src/phase_mechanical_plastic_nonlocal.f90 +++ b/src/phase_mechanical_plastic_nonlocal.f90 @@ -741,7 +741,7 @@ end subroutine nonlocal_dependentState !> @brief calculates plastic velocity gradient and its tangent !-------------------------------------------------------------------------------------------------- module subroutine nonlocal_LpAndItsTangent(Lp,dLp_dMp, & - Mp,Temperature,ph,en) + Mp,ph,en) real(pReal), dimension(3,3), intent(out) :: & Lp !< plastic velocity gradient real(pReal), dimension(3,3,3,3), intent(out) :: & @@ -749,9 +749,6 @@ module subroutine nonlocal_LpAndItsTangent(Lp,dLp_dMp, & integer, intent(in) :: & ph, & en - real(pReal), intent(in) :: & - Temperature !< temperature - real(pReal), dimension(3,3), intent(in) :: & Mp !< derivative of Lp with respect to Mp @@ -771,8 +768,14 @@ module subroutine nonlocal_LpAndItsTangent(Lp,dLp_dMp, & real(pReal), dimension(param(ph)%sum_N_sl) :: & tau, & !< resolved shear stress including backstress terms dot_gamma !< shear rate + real(pReal) :: & + Temperature !< temperature + Temperature = thermal_T(ph,en) + Lp = 0.0_pReal + dLp_dMp = 0.0_pReal + associate(prm => param(ph),dst=>dependentState(ph),stt=>state(ph)) !*** shortcut to state variables @@ -821,8 +824,6 @@ module subroutine nonlocal_LpAndItsTangent(Lp,dLp_dMp, & dot_gamma = sum(rhoSgl(:,1:4) * v, 2) * prm%b_sl - Lp = 0.0_pReal - dLp_dMp = 0.0_pReal do s = 1,prm%sum_N_sl Lp = Lp + dot_gamma(s) * prm%P_sl(1:3,1:3,s) forall (i=1:3,j=1:3,k=1:3,l=1:3) & @@ -1649,10 +1650,12 @@ pure subroutine kinetics(v, dv_dtau, dv_dtauNS, tau, tauNS, tauThreshold, c, T, criticalStress_P, & !< maximum obstacle strength criticalStress_S !< maximum obstacle strength + + v = 0.0_pReal + dv_dtau = 0.0_pReal + dv_dtauNS = 0.0_pReal + associate(prm => param(ph)) - v = 0.0_pReal - dv_dtau = 0.0_pReal - dv_dtauNS = 0.0_pReal do s = 1,prm%sum_N_sl if (abs(tau(s)) > tauThreshold(s)) then diff --git a/src/phase_thermal.f90 b/src/phase_thermal.f90 index 4ea21d039..a3e4dd628 100644 --- a/src/phase_thermal.f90 +++ b/src/phase_thermal.f90 @@ -271,7 +271,7 @@ end subroutine thermal_forward !---------------------------------------------------------------------------------------------- !< @brief Get temperature (for use by non-thermal physics) !---------------------------------------------------------------------------------------------- -module function thermal_T(ph,en) result(T) +pure module function thermal_T(ph,en) result(T) integer, intent(in) :: ph, en real(pReal) :: T From e10dea5b6cefeaa668a3e6f608cd1b64306c5391 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 19 Dec 2021 22:53:48 +0100 Subject: [PATCH 3/3] easier to understand --- src/phase_mechanical_plastic_dislotwin.f90 | 34 +++++++++++----------- 1 file changed, 17 insertions(+), 17 deletions(-) diff --git a/src/phase_mechanical_plastic_dislotwin.f90 b/src/phase_mechanical_plastic_dislotwin.f90 index cfc9a002f..78681fa24 100644 --- a/src/phase_mechanical_plastic_dislotwin.f90 +++ b/src/phase_mechanical_plastic_dislotwin.f90 @@ -476,18 +476,18 @@ module function plastic_dislotwin_homogenizedC(ph,en) result(homogenizedC) C66_tw, & C66_tr integer :: i - real(pReal) :: f_unrotated + real(pReal) :: f_matrix C = elastic_C66(ph,en) associate(prm => param(ph), stt => state(ph)) - f_unrotated = 1.0_pReal & - - sum(stt%f_tw(1:prm%sum_N_tw,en)) & - - sum(stt%f_tr(1:prm%sum_N_tr,en)) + f_matrix = 1.0_pReal & + - sum(stt%f_tw(1:prm%sum_N_tw,en)) & + - sum(stt%f_tr(1:prm%sum_N_tr,en)) - homogenizedC = f_unrotated * C + homogenizedC = f_matrix * C twinActive: if (prm%sum_N_tw > 0) then C66_tw = lattice_C66_twin(prm%N_tw,C,phase_lattice(ph),phase_cOverA(ph)) @@ -522,7 +522,7 @@ module subroutine dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,en) integer :: i,k,l,m,n real(pReal) :: & - f_unrotated,StressRatio_p,& + f_matrix,StressRatio_p,& E_kB_T, & ddot_gamma_dtau, & tau, & @@ -563,9 +563,9 @@ module subroutine dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,en) associate(prm => param(ph), stt => state(ph)) - f_unrotated = 1.0_pReal & - - sum(stt%f_tw(1:prm%sum_N_tw,en)) & - - sum(stt%f_tr(1:prm%sum_N_tr,en)) + f_matrix = 1.0_pReal & + - sum(stt%f_tw(1:prm%sum_N_tw,en)) & + - sum(stt%f_tr(1:prm%sum_N_tr,en)) call kinetics_sl(Mp,T,ph,en,dot_gamma_sl,ddot_gamma_dtau_sl) slipContribution: do i = 1, prm%sum_N_sl @@ -591,8 +591,8 @@ module subroutine dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,en) + ddot_gamma_dtau_tr(i)* prm%P_tr(k,l,i)*prm%P_tr(m,n,i) end do transContibution - Lp = Lp * f_unrotated - dLp_dMp = dLp_dMp * f_unrotated + Lp = Lp * f_matrix + dLp_dMp = dLp_dMp * f_matrix shearBandingContribution: if (dNeq0(prm%v_sb)) then @@ -640,7 +640,7 @@ module subroutine dislotwin_dotState(Mp,T,ph,en) integer :: i real(pReal) :: & - f_unrotated, & + f_matrix, & d_hat, & v_cl, & !< climb velocity tau, & @@ -663,9 +663,9 @@ module subroutine dislotwin_dotState(Mp,T,ph,en) mu = elastic_mu(ph,en) nu = elastic_nu(ph,en) - f_unrotated = 1.0_pReal & - - sum(stt%f_tw(1:prm%sum_N_tw,en)) & - - sum(stt%f_tr(1:prm%sum_N_tr,en)) + f_matrix = 1.0_pReal & + - sum(stt%f_tw(1:prm%sum_N_tw,en)) & + - sum(stt%f_tr(1:prm%sum_N_tr,en)) call kinetics_sl(Mp,T,ph,en,dot_gamma_sl) dot%gamma_sl(:,en) = abs(dot_gamma_sl) @@ -711,10 +711,10 @@ module subroutine dislotwin_dotState(Mp,T,ph,en) - dot_rho_dip_climb call kinetics_tw(Mp,T,dot_gamma_sl,ph,en,dot_gamma_tw) - dot%f_tw(:,en) = f_unrotated*dot_gamma_tw/prm%gamma_char + dot%f_tw(:,en) = f_matrix*dot_gamma_tw/prm%gamma_char call kinetics_tr(Mp,T,dot_gamma_sl,ph,en,dot_gamma_tr) - dot%f_tr(:,en) = f_unrotated*dot_gamma_tr + dot%f_tr(:,en) = f_matrix*dot_gamma_tr end associate