223031012
/
DAMASK_EICMD
mirror of https://github.com/achalhp/DAMASK_EICMD.git
1
0
Fork 0
Code Issues Projects Releases Wiki Activity

new names (as in HDF5 out)

This commit is contained in:
Martin Diehl 2021-04-12 21:21:15 +02:00
parent 887524bcc1
commit 74548d5f51
6 changed files with 213 additions and 215 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

16
src/phase.f90
View File

@ -135,18 +135,18 @@ module phase
end subroutine mechanical_restartRead end subroutine mechanical_restartRead
module function mechanical_S(ph,me) result(S) module function mechanical_S(ph,en) result(S)
integer, intent(in) :: ph,me integer, intent(in) :: ph,en
real(pReal), dimension(3,3) :: S real(pReal), dimension(3,3) :: S
end function mechanical_S end function mechanical_S
module function mechanical_L_p(ph,me) result(L_p) module function mechanical_L_p(ph,en) result(L_p)
integer, intent(in) :: ph,me integer, intent(in) :: ph,en
real(pReal), dimension(3,3) :: L_p real(pReal), dimension(3,3) :: L_p
end function mechanical_L_p end function mechanical_L_p
module function mechanical_F_e(ph,me) result(F_e) module function mechanical_F_e(ph,en) result(F_e)
integer, intent(in) :: ph,me integer, intent(in) :: ph,en
real(pReal), dimension(3,3) :: F_e real(pReal), dimension(3,3) :: F_e
end function mechanical_F_e end function mechanical_F_e
@ -239,8 +239,8 @@ module phase
logical :: converged_ logical :: converged_
end function crystallite_stress end function crystallite_stress
module function phase_homogenizedC(ph,me) result(C) module function phase_homogenizedC(ph,en) result(C)
integer, intent(in) :: ph, me integer, intent(in) :: ph, en
real(pReal), dimension(6,6) :: C real(pReal), dimension(6,6) :: C
end function phase_homogenizedC end function phase_homogenizedC

322
src/phase_mechanical.f90
View File

@ -60,7 +60,7 @@ submodule(phase) mechanical
module subroutine plastic_init module subroutine plastic_init
end subroutine plastic_init end subroutine plastic_init
module subroutine plastic_isotropic_LiAndItsTangent(Li,dLi_dMi,Mi,ph,me) module subroutine plastic_isotropic_LiAndItsTangent(Li,dLi_dMi,Mi,ph,en)
real(pReal), dimension(3,3), intent(out) :: & real(pReal), dimension(3,3), intent(out) :: &
Li !< inleastic velocity gradient Li !< inleastic velocity gradient
real(pReal), dimension(3,3,3,3), intent(out) :: & real(pReal), dimension(3,3,3,3), intent(out) :: &
@ -69,34 +69,34 @@ submodule(phase) mechanical
Mi !< Mandel stress Mi !< Mandel stress
integer, intent(in) :: & integer, intent(in) :: &
ph, & ph, &
me en
end subroutine plastic_isotropic_LiAndItsTangent end subroutine plastic_isotropic_LiAndItsTangent
module function plastic_dotState(subdt,co,ip,el,ph,me) result(broken) module function plastic_dotState(subdt,co,ip,el,ph,en) result(broken)
integer, intent(in) :: & integer, intent(in) :: &
co, & !< component-ID of integration point co, & !< component-ID of integration point
ip, & !< integration point ip, & !< integration point
el, & !< element el, & !< element
ph, & ph, &
me en
real(pReal), intent(in) :: & real(pReal), intent(in) :: &
subdt !< timestep subdt !< timestep
logical :: broken logical :: broken
end function plastic_dotState end function plastic_dotState
module function plastic_deltaState(ph, me) result(broken) module function plastic_deltaState(ph, en) result(broken)
integer, intent(in) :: & integer, intent(in) :: &
ph, & ph, &
me en
logical :: & logical :: &
broken broken
end function plastic_deltaState end function plastic_deltaState
module subroutine phase_LiAndItsTangents(Li, dLi_dS, dLi_dFi, & module subroutine phase_LiAndItsTangents(Li, dLi_dS, dLi_dFi, &
S, Fi, ph,me) S, Fi, ph,en)
integer, intent(in) :: & integer, intent(in) :: &
ph,me ph,en
real(pReal), intent(in), dimension(3,3) :: & real(pReal), intent(in), dimension(3,3) :: &
S !< 2nd Piola-Kirchhoff stress S !< 2nd Piola-Kirchhoff stress
real(pReal), intent(in), dimension(3,3) :: & real(pReal), intent(in), dimension(3,3) :: &
@ -111,9 +111,9 @@ submodule(phase) mechanical
module subroutine plastic_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, & module subroutine plastic_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, &
S, Fi, ph,me) S, Fi, ph,en)
integer, intent(in) :: & integer, intent(in) :: &
ph,me ph,en
real(pReal), intent(in), dimension(3,3) :: & real(pReal), intent(in), dimension(3,3) :: &
S, & !< 2nd Piola-Kirchhoff stress S, & !< 2nd Piola-Kirchhoff stress
Fi !< intermediate deformation gradient Fi !< intermediate deformation gradient
@ -155,9 +155,9 @@ submodule(phase) mechanical
character(len=*), intent(in) :: group character(len=*), intent(in) :: group
end subroutine plastic_nonlocal_results end subroutine plastic_nonlocal_results
module function plastic_dislotwin_homogenizedC(ph,me) result(homogenizedC) module function plastic_dislotwin_homogenizedC(ph,en) result(homogenizedC)
real(pReal), dimension(6,6) :: homogenizedC real(pReal), dimension(6,6) :: homogenizedC
integer, intent(in) :: ph,me integer, intent(in) :: ph,en
end function plastic_dislotwin_homogenizedC end function plastic_dislotwin_homogenizedC
@ -189,7 +189,7 @@ module subroutine mechanical_init(materials,phases)
co, & co, &
ce, & ce, &
ph, & ph, &
me, & en, &
stiffDegradationCtr, & stiffDegradationCtr, &
Nmembers Nmembers
class(tNode), pointer :: & class(tNode), pointer :: &
@ -286,21 +286,21 @@ module subroutine mechanical_init(materials,phases)
constituent => constituents%get(co) constituent => constituents%get(co)
ph = material_phaseAt(co,el) ph = material_phaseAt(co,el)
me = material_phaseMemberAt(co,ip,el) en = material_phaseMemberAt(co,ip,el)
call material_orientation0(co,ph,me)%fromQuaternion(constituent%get_as1dFloat('O',requiredSize=4)) call material_orientation0(co,ph,en)%fromQuaternion(constituent%get_as1dFloat('O',requiredSize=4))
phase_mechanical_Fp0(ph)%data(1:3,1:3,me) = material_orientation0(co,ph,me)%asMatrix() ! Fp reflects initial orientation (see 10.1016/j.actamat.2006.01.005) phase_mechanical_Fp0(ph)%data(1:3,1:3,en) = material_orientation0(co,ph,en)%asMatrix() ! Fp reflects initial orientation (see 10.1016/j.actamat.2006.01.005)
phase_mechanical_Fp0(ph)%data(1:3,1:3,me) = phase_mechanical_Fp0(ph)%data(1:3,1:3,me) & phase_mechanical_Fp0(ph)%data(1:3,1:3,en) = phase_mechanical_Fp0(ph)%data(1:3,1:3,en) &
/ math_det33(phase_mechanical_Fp0(ph)%data(1:3,1:3,me))**(1.0_pReal/3.0_pReal) / math_det33(phase_mechanical_Fp0(ph)%data(1:3,1:3,en))**(1.0_pReal/3.0_pReal)
phase_mechanical_Fi0(ph)%data(1:3,1:3,me) = math_I3 phase_mechanical_Fi0(ph)%data(1:3,1:3,en) = math_I3
phase_mechanical_F0(ph)%data(1:3,1:3,me) = math_I3 phase_mechanical_F0(ph)%data(1:3,1:3,en) = math_I3
phase_mechanical_Fe(ph)%data(1:3,1:3,me) = math_inv33(matmul(phase_mechanical_Fi0(ph)%data(1:3,1:3,me), & phase_mechanical_Fe(ph)%data(1:3,1:3,en) = math_inv33(matmul(phase_mechanical_Fi0(ph)%data(1:3,1:3,en), &
phase_mechanical_Fp0(ph)%data(1:3,1:3,me))) ! assuming that euler angles are given in internal strain free configuration phase_mechanical_Fp0(ph)%data(1:3,1:3,en))) ! assuming that euler angles are given in internal strain free configuration
phase_mechanical_Fp(ph)%data(1:3,1:3,me) = phase_mechanical_Fp0(ph)%data(1:3,1:3,me) phase_mechanical_Fp(ph)%data(1:3,1:3,en) = phase_mechanical_Fp0(ph)%data(1:3,1:3,en)
phase_mechanical_Fi(ph)%data(1:3,1:3,me) = phase_mechanical_Fi0(ph)%data(1:3,1:3,me) phase_mechanical_Fi(ph)%data(1:3,1:3,en) = phase_mechanical_Fi0(ph)%data(1:3,1:3,en)
phase_mechanical_F(ph)%data(1:3,1:3,me) = phase_mechanical_F0(ph)%data(1:3,1:3,me) phase_mechanical_F(ph)%data(1:3,1:3,en) = phase_mechanical_F0(ph)%data(1:3,1:3,en)
enddo enddo
enddo; enddo enddo; enddo
@ -353,11 +353,11 @@ end subroutine mechanical_init
!> the elastic and intermediate deformation gradients using Hooke's law !> the elastic and intermediate deformation gradients using Hooke's law
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine phase_hooke_SandItsTangents(S, dS_dFe, dS_dFi, & subroutine phase_hooke_SandItsTangents(S, dS_dFe, dS_dFi, &
Fe, Fi, ph, me) Fe, Fi, ph, en)
integer, intent(in) :: & integer, intent(in) :: &
ph, & ph, &
me en
real(pReal), intent(in), dimension(3,3) :: & real(pReal), intent(in), dimension(3,3) :: &
Fe, & !< elastic deformation gradient Fe, & !< elastic deformation gradient
Fi !< intermediate deformation gradient Fi !< intermediate deformation gradient
@ -373,12 +373,12 @@ subroutine phase_hooke_SandItsTangents(S, dS_dFe, dS_dFi, &
d, & !< counter in degradation loop d, & !< counter in degradation loop
i, j i, j
C = math_66toSym3333(phase_homogenizedC(ph,me)) C = math_66toSym3333(phase_homogenizedC(ph,en))
DegradationLoop: do d = 1, phase_NstiffnessDegradations(ph) DegradationLoop: do d = 1, phase_NstiffnessDegradations(ph)
degradationType: select case(phase_stiffnessDegradation(d,ph)) degradationType: select case(phase_stiffnessDegradation(d,ph))
case (STIFFNESS_DEGRADATION_damage_ID) degradationType case (STIFFNESS_DEGRADATION_damage_ID) degradationType
C = C * damage_phi(ph,me)**2 C = C * damage_phi(ph,en)**2
end select degradationType end select degradationType
enddo DegradationLoop enddo DegradationLoop
@ -486,7 +486,7 @@ function integrateStress(F,subFp0,subFi0,Delta_t,co,ip,el) result(broken)
o, & o, &
p, & p, &
ph, & ph, &
me, & en, &
jacoCounterLp, & jacoCounterLp, &
jacoCounterLi ! counters to check for Jacobian update jacoCounterLi ! counters to check for Jacobian update
logical :: error,broken logical :: error,broken
@ -495,12 +495,12 @@ function integrateStress(F,subFp0,subFi0,Delta_t,co,ip,el) result(broken)
broken = .true. broken = .true.
ph = material_phaseAt(co,el) ph = material_phaseAt(co,el)
me = material_phaseMemberAt(co,ip,el) en = material_phaseMemberAt(co,ip,el)
call plastic_dependentState(co,ip,el) call plastic_dependentState(co,ip,el)
Lpguess = phase_mechanical_Lp(ph)%data(1:3,1:3,me) ! take as first guess Lpguess = phase_mechanical_Lp(ph)%data(1:3,1:3,en) ! take as first guess
Liguess = phase_mechanical_Li(ph)%data(1:3,1:3,me) ! take as first guess Liguess = phase_mechanical_Li(ph)%data(1:3,1:3,en) ! take as first guess
call math_invert33(invFp_current,devNull,error,subFp0) call math_invert33(invFp_current,devNull,error,subFp0)
if (error) return ! error if (error) return ! error
@ -535,10 +535,10 @@ function integrateStress(F,subFp0,subFi0,Delta_t,co,ip,el) result(broken)
B = math_I3 - Delta_t*Lpguess B = math_I3 - Delta_t*Lpguess
Fe = matmul(matmul(A,B), invFi_new) Fe = matmul(matmul(A,B), invFi_new)
call phase_hooke_SandItsTangents(S, dS_dFe, dS_dFi, & call phase_hooke_SandItsTangents(S, dS_dFe, dS_dFi, &
Fe, Fi_new, ph, me) Fe, Fi_new, ph, en)
call plastic_LpAndItsTangents(Lp_constitutive, dLp_dS, dLp_dFi, & call plastic_LpAndItsTangents(Lp_constitutive, dLp_dS, dLp_dFi, &
S, Fi_new, ph,me) S, Fi_new, ph,en)
!* update current residuum and check for convergence of loop !* update current residuum and check for convergence of loop
atol_Lp = max(num%rtol_crystalliteStress * max(norm2(Lpguess),norm2(Lp_constitutive)), & ! absolute tolerance from largest acceptable relative error atol_Lp = max(num%rtol_crystalliteStress * max(norm2(Lpguess),norm2(Lp_constitutive)), & ! absolute tolerance from largest acceptable relative error
@ -579,7 +579,7 @@ function integrateStress(F,subFp0,subFi0,Delta_t,co,ip,el) result(broken)
enddo LpLoop enddo LpLoop
call phase_LiAndItsTangents(Li_constitutive, dLi_dS, dLi_dFi, & call phase_LiAndItsTangents(Li_constitutive, dLi_dS, dLi_dFi, &
S, Fi_new, ph,me) S, Fi_new, ph,en)
!* update current residuum and check for convergence of loop !* update current residuum and check for convergence of loop
atol_Li = max(num%rtol_crystalliteStress * max(norm2(Liguess),norm2(Li_constitutive)), & ! absolute tolerance from largest acceptable relative error atol_Li = max(num%rtol_crystalliteStress * max(norm2(Liguess),norm2(Li_constitutive)), & ! absolute tolerance from largest acceptable relative error
@ -629,13 +629,13 @@ function integrateStress(F,subFp0,subFi0,Delta_t,co,ip,el) result(broken)
call math_invert33(Fp_new,devNull,error,invFp_new) call math_invert33(Fp_new,devNull,error,invFp_new)
if (error) return ! error if (error) return ! error
phase_mechanical_P(ph)%data(1:3,1:3,me) = matmul(matmul(F,invFp_new),matmul(S,transpose(invFp_new))) phase_mechanical_P(ph)%data(1:3,1:3,en) = matmul(matmul(F,invFp_new),matmul(S,transpose(invFp_new)))
phase_mechanical_S(ph)%data(1:3,1:3,me) = S phase_mechanical_S(ph)%data(1:3,1:3,en) = S
phase_mechanical_Lp(ph)%data(1:3,1:3,me) = Lpguess phase_mechanical_Lp(ph)%data(1:3,1:3,en) = Lpguess
phase_mechanical_Li(ph)%data(1:3,1:3,me) = Liguess phase_mechanical_Li(ph)%data(1:3,1:3,en) = Liguess
phase_mechanical_Fp(ph)%data(1:3,1:3,me) = Fp_new / math_det33(Fp_new)**(1.0_pReal/3.0_pReal) ! regularize phase_mechanical_Fp(ph)%data(1:3,1:3,en) = Fp_new / math_det33(Fp_new)**(1.0_pReal/3.0_pReal) ! regularize
phase_mechanical_Fi(ph)%data(1:3,1:3,me) = Fi_new phase_mechanical_Fi(ph)%data(1:3,1:3,en) = Fi_new
phase_mechanical_Fe(ph)%data(1:3,1:3,me) = matmul(matmul(F,invFp_new),invFi_new) phase_mechanical_Fe(ph)%data(1:3,1:3,en) = matmul(matmul(F,invFp_new),invFi_new)
broken = .false. broken = .false.
end function integrateStress end function integrateStress
@ -660,7 +660,7 @@ function integrateStateFPI(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el) resul
integer :: & integer :: &
NiterationState, & !< number of iterations in state loop NiterationState, & !< number of iterations in state loop
ph, & ph, &
me, & en, &
sizeDotState sizeDotState
real(pReal) :: & real(pReal) :: &
zeta zeta
@ -671,38 +671,37 @@ function integrateStateFPI(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el) resul
ph = material_phaseAt(co,el) ph = material_phaseAt(co,el)
me = material_phaseMemberAt(co,ip,el) en = material_phaseMemberAt(co,ip,el)
broken = plastic_dotState(Delta_t, co,ip,el,ph,me) broken = plastic_dotState(Delta_t, co,ip,el,ph,en)
if(broken) return if(broken) return
sizeDotState = plasticState(ph)%sizeDotState sizeDotState = plasticState(ph)%sizeDotState
plasticState(ph)%state(1:sizeDotState,me) = subState0 & plasticState(ph)%state(1:sizeDotState,en) = subState0 &
+ plasticState(ph)%dotState (1:sizeDotState,me) * Delta_t + plasticState(ph)%dotState (1:sizeDotState,en) * Delta_t
dotState(1:sizeDotState,2) = 0.0_pReal
iteration: do NiterationState = 1, num%nState iteration: do NiterationState = 1, num%nState
if(nIterationState > 1) dotState(1:sizeDotState,2) = dotState(1:sizeDotState,1) dotState(1:sizeDotState,2) = merge(dotState(1:sizeDotState,1),0.0, nIterationState > 1)
dotState(1:sizeDotState,1) = plasticState(ph)%dotState(:,me) dotState(1:sizeDotState,1) = plasticState(ph)%dotState(:,en)
broken = integrateStress(F,subFp0,subFi0,Delta_t,co,ip,el) broken = integrateStress(F,subFp0,subFi0,Delta_t,co,ip,el)
if(broken) exit iteration if(broken) exit iteration
broken = plastic_dotState(Delta_t, co,ip,el,ph,me) broken = plastic_dotState(Delta_t, co,ip,el,ph,en)
if(broken) exit iteration if(broken) exit iteration
zeta = damper(plasticState(ph)%dotState(:,me),dotState(1:sizeDotState,1),& zeta = damper(plasticState(ph)%dotState(:,en),dotState(1:sizeDotState,1),&
dotState(1:sizeDotState,2)) dotState(1:sizeDotState,2))
plasticState(ph)%dotState(:,me) = plasticState(ph)%dotState(:,me) * zeta & plasticState(ph)%dotState(:,en) = plasticState(ph)%dotState(:,en) * zeta &
+ dotState(1:sizeDotState,1) * (1.0_pReal - zeta) + dotState(1:sizeDotState,1) * (1.0_pReal - zeta)
r(1:sizeDotState) = plasticState(ph)%state (1:sizeDotState,me) & r(1:sizeDotState) = plasticState(ph)%state(1:sizeDotState,en) &
- subState0 & - subState0 &
- plasticState(ph)%dotState (1:sizeDotState,me) * Delta_t - plasticState(ph)%dotState(1:sizeDotState,en) * Delta_t
plasticState(ph)%state(1:sizeDotState,me) = plasticState(ph)%state(1:sizeDotState,me) & plasticState(ph)%state(1:sizeDotState,en) = plasticState(ph)%state(1:sizeDotState,en) &
- r(1:sizeDotState) - r(1:sizeDotState)
if (converged(r(1:sizeDotState),plasticState(ph)%state(1:sizeDotState,me),plasticState(ph)%atol(1:sizeDotState))) then if (converged(r(1:sizeDotState),plasticState(ph)%state(1:sizeDotState,en),plasticState(ph)%atol(1:sizeDotState))) then
broken = plastic_deltaState(ph,me) broken = plastic_deltaState(ph,en)
exit iteration exit iteration
endif endif
@ -714,20 +713,19 @@ function integrateStateFPI(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el) resul
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief calculate the damping for correction of state and dot state !> @brief calculate the damping for correction of state and dot state
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
real(pReal) pure function damper(current,previous,previous2) real(pReal) pure function damper(omega_0,omega_1,omega_2)
real(pReal), dimension(:), intent(in) ::& real(pReal), dimension(:), intent(in) :: &
current, previous, previous2 omega_0, omega_1, omega_2
real(pReal) :: dot_prod12, dot_prod22 real(pReal) :: dot_prod12, dot_prod22
dot_prod12 = dot_product(current - previous, previous - previous2) dot_prod12 = dot_product(omega_0 - omega_1, omega_1 - omega_2)
dot_prod22 = dot_product(previous - previous2, previous - previous2) dot_prod22 = dot_product(omega_1 - omega_2, omega_1 - omega_2)
if ((dot_product(current,previous) < 0.0_pReal .or. dot_prod12 < 0.0_pReal) .and. dot_prod22 > 0.0_pReal) then
damper = 0.75_pReal + 0.25_pReal * tanh(2.0_pReal + 4.0_pReal * dot_prod12 / dot_prod22) damper = merge(0.75_pReal + 0.25_pReal * tanh(2.0_pReal + 4.0_pReal * dot_prod12 / dot_prod22), &
else 1.0_pReal, &
damper = 1.0_pReal (min(dot_product(omega_0,omega_1), dot_prod12) < 0.0_pReal) .and. dot_prod22 > 0.0_pReal)
endif
end function damper end function damper
@ -751,21 +749,21 @@ function integrateStateEuler(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el) res
integer :: & integer :: &
ph, & ph, &
me, & en, &
sizeDotState sizeDotState
ph = material_phaseAt(co,el) ph = material_phaseAt(co,el)
me = material_phaseMemberAt(co,ip,el) en = material_phaseMemberAt(co,ip,el)
broken = plastic_dotState(Delta_t, co,ip,el,ph,me) broken = plastic_dotState(Delta_t, co,ip,el,ph,en)
if(broken) return if(broken) return
sizeDotState = plasticState(ph)%sizeDotState sizeDotState = plasticState(ph)%sizeDotState
plasticState(ph)%state(1:sizeDotState,me) = subState0 & plasticState(ph)%state(1:sizeDotState,en) = subState0 &
+ plasticState(ph)%dotState(1:sizeDotState,me) * Delta_t + plasticState(ph)%dotState(1:sizeDotState,en) * Delta_t
broken = plastic_deltaState(ph,me) broken = plastic_deltaState(ph,en)
if(broken) return if(broken) return
broken = integrateStress(F,subFp0,subFi0,Delta_t,co,ip,el) broken = integrateStress(F,subFp0,subFi0,Delta_t,co,ip,el)
@ -790,34 +788,34 @@ function integrateStateAdaptiveEuler(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip
integer :: & integer :: &
ph, & ph, &
me, & en, &
sizeDotState sizeDotState
real(pReal), dimension(phase_plasticity_maxSizeDotState) :: residuum_plastic real(pReal), dimension(phase_plasticity_maxSizeDotState) :: residuum_plastic
ph = material_phaseAt(co,el) ph = material_phaseAt(co,el)
me = material_phaseMemberAt(co,ip,el) en = material_phaseMemberAt(co,ip,el)
broken = plastic_dotState(Delta_t, co,ip,el,ph,me) broken = plastic_dotState(Delta_t, co,ip,el,ph,en)
if(broken) return if(broken) return
sizeDotState = plasticState(ph)%sizeDotState sizeDotState = plasticState(ph)%sizeDotState
residuum_plastic(1:sizeDotState) = - plasticState(ph)%dotstate(1:sizeDotState,me) * 0.5_pReal * Delta_t residuum_plastic(1:sizeDotState) = - plasticState(ph)%dotstate(1:sizeDotState,en) * 0.5_pReal * Delta_t
plasticState(ph)%state(1:sizeDotState,me) = subState0 & plasticState(ph)%state(1:sizeDotState,en) = subState0 &
+ plasticState(ph)%dotstate(1:sizeDotState,me) * Delta_t + plasticState(ph)%dotstate(1:sizeDotState,en) * Delta_t
broken = plastic_deltaState(ph,me) broken = plastic_deltaState(ph,en)
if(broken) return if(broken) return
broken = integrateStress(F,subFp0,subFi0,Delta_t,co,ip,el) broken = integrateStress(F,subFp0,subFi0,Delta_t,co,ip,el)
if(broken) return if(broken) return
broken = plastic_dotState(Delta_t, co,ip,el,ph,me) broken = plastic_dotState(Delta_t, co,ip,el,ph,en)
if(broken) return if(broken) return
broken = .not. converged(residuum_plastic(1:sizeDotState) + 0.5_pReal * plasticState(ph)%dotState(:,me) * Delta_t, & broken = .not. converged(residuum_plastic(1:sizeDotState) + 0.5_pReal * plasticState(ph)%dotState(:,en) * Delta_t, &
plasticState(ph)%state(1:sizeDotState,me), & plasticState(ph)%state(1:sizeDotState,en), &
plasticState(ph)%atol(1:sizeDotState)) plasticState(ph)%atol(1:sizeDotState))
end function integrateStateAdaptiveEuler end function integrateStateAdaptiveEuler
@ -908,55 +906,55 @@ function integrateStateRK(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el,A,B,C,D
stage, & ! stage index in integration stage loop stage, & ! stage index in integration stage loop
n, & n, &
ph, & ph, &
me, & en, &
sizeDotState sizeDotState
real(pReal), dimension(phase_plasticity_maxSizeDotState,size(B)) :: plastic_RKdotState real(pReal), dimension(phase_plasticity_maxSizeDotState,size(B)) :: plastic_RKdotState
ph = material_phaseAt(co,el) ph = material_phaseAt(co,el)
me = material_phaseMemberAt(co,ip,el) en = material_phaseMemberAt(co,ip,el)
broken = plastic_dotState(Delta_t,co,ip,el,ph,me) broken = plastic_dotState(Delta_t,co,ip,el,ph,en)
if(broken) return if(broken) return
sizeDotState = plasticState(ph)%sizeDotState sizeDotState = plasticState(ph)%sizeDotState
do stage = 1, size(A,1) do stage = 1, size(A,1)
plastic_RKdotState(1:sizeDotState,stage) = plasticState(ph)%dotState(:,me) plastic_RKdotState(1:sizeDotState,stage) = plasticState(ph)%dotState(:,en)
plasticState(ph)%dotState(:,me) = A(1,stage) * plastic_RKdotState(1:sizeDotState,1) plasticState(ph)%dotState(:,en) = A(1,stage) * plastic_RKdotState(1:sizeDotState,1)
do n = 2, stage do n = 2, stage
plasticState(ph)%dotState(:,me) = plasticState(ph)%dotState(:,me) & plasticState(ph)%dotState(:,en) = plasticState(ph)%dotState(:,en) &
+ A(n,stage) * plastic_RKdotState(1:sizeDotState,n) + A(n,stage) * plastic_RKdotState(1:sizeDotState,n)
enddo enddo
plasticState(ph)%state(1:sizeDotState,me) = subState0 & plasticState(ph)%state(1:sizeDotState,en) = subState0 &
+ plasticState(ph)%dotState (1:sizeDotState,me) * Delta_t + plasticState(ph)%dotState (1:sizeDotState,en) * Delta_t
broken = integrateStress(F_0 + (F - F_0) * Delta_t * C(stage),subFp0,subFi0,Delta_t * C(stage),co,ip,el) broken = integrateStress(F_0 + (F - F_0) * Delta_t * C(stage),subFp0,subFi0,Delta_t * C(stage),co,ip,el)
if(broken) exit if(broken) exit
broken = plastic_dotState(Delta_t*C(stage),co,ip,el,ph,me) broken = plastic_dotState(Delta_t*C(stage),co,ip,el,ph,en)
if(broken) exit if(broken) exit
enddo enddo
if(broken) return if(broken) return
plastic_RKdotState(1:sizeDotState,size(B)) = plasticState (ph)%dotState(:,me) plastic_RKdotState(1:sizeDotState,size(B)) = plasticState (ph)%dotState(:,en)
plasticState(ph)%dotState(:,me) = matmul(plastic_RKdotState(1:sizeDotState,1:size(B)),B) plasticState(ph)%dotState(:,en) = matmul(plastic_RKdotState(1:sizeDotState,1:size(B)),B)
plasticState(ph)%state(1:sizeDotState,me) = subState0 & plasticState(ph)%state(1:sizeDotState,en) = subState0 &
+ plasticState(ph)%dotState (1:sizeDotState,me) * Delta_t + plasticState(ph)%dotState (1:sizeDotState,en) * Delta_t
if(present(DB)) & if(present(DB)) &
broken = .not. converged(matmul(plastic_RKdotState(1:sizeDotState,1:size(DB)),DB) * Delta_t, & broken = .not. converged(matmul(plastic_RKdotState(1:sizeDotState,1:size(DB)),DB) * Delta_t, &
plasticState(ph)%state(1:sizeDotState,me), & plasticState(ph)%state(1:sizeDotState,en), &
plasticState(ph)%atol(1:sizeDotState)) plasticState(ph)%atol(1:sizeDotState))
if(broken) return if(broken) return
broken = plastic_deltaState(ph,me) broken = plastic_deltaState(ph,en)
if(broken) return if(broken) return
broken = integrateStress(F,subFp0,subFi0,Delta_t,co,ip,el) broken = integrateStress(F,subFp0,subFi0,Delta_t,co,ip,el)
@ -1087,14 +1085,14 @@ end subroutine mechanical_forward
!> @brief returns the homogenize elasticity matrix !> @brief returns the homogenize elasticity matrix
!> ToDo: homogenizedC66 would be more consistent !> ToDo: homogenizedC66 would be more consistent
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module function phase_homogenizedC(ph,me) result(C) module function phase_homogenizedC(ph,en) result(C)
real(pReal), dimension(6,6) :: C real(pReal), dimension(6,6) :: C
integer, intent(in) :: ph, me integer, intent(in) :: ph, en
plasticType: select case (phase_plasticity(ph)) plasticType: select case (phase_plasticity(ph))
case (PLASTICITY_DISLOTWIN_ID) plasticType case (PLASTICITY_DISLOTWIN_ID) plasticType
C = plastic_dislotwin_homogenizedC(ph,me) C = plastic_dislotwin_homogenizedC(ph,en)
case default plasticType case default plasticType
C = lattice_C66(1:6,1:6,ph) C = lattice_C66(1:6,1:6,ph)
end select plasticType end select plasticType
@ -1117,7 +1115,7 @@ module function crystallite_stress(dt,co,ip,el) result(converged_)
real(pReal) :: & real(pReal) :: &
formerSubStep formerSubStep
integer :: & integer :: &
ph, me, sizeDotState ph, en, sizeDotState
logical :: todo logical :: todo
real(pReal) :: subFrac,subStep real(pReal) :: subFrac,subStep
real(pReal), dimension(3,3) :: & real(pReal), dimension(3,3) :: &
@ -1131,19 +1129,19 @@ module function crystallite_stress(dt,co,ip,el) result(converged_)
ph = material_phaseAt(co,el) ph = material_phaseAt(co,el)
me = material_phaseMemberAt(co,ip,el) en = material_phaseMemberAt(co,ip,el)
sizeDotState = plasticState(ph)%sizeDotState sizeDotState = plasticState(ph)%sizeDotState
subLi0 = phase_mechanical_Li0(ph)%data(1:3,1:3,me) subLi0 = phase_mechanical_Li0(ph)%data(1:3,1:3,en)
subLp0 = phase_mechanical_Lp0(ph)%data(1:3,1:3,me) subLp0 = phase_mechanical_Lp0(ph)%data(1:3,1:3,en)
subState0 = plasticState(ph)%State0(:,me) subState0 = plasticState(ph)%State0(:,en)
if (damageState(ph)%sizeState > 0) & if (damageState(ph)%sizeState > 0) &
damageState(ph)%subState0(:,me) = damageState(ph)%state0(:,me) damageState(ph)%subState0(:,en) = damageState(ph)%state0(:,en)
subFp0 = phase_mechanical_Fp0(ph)%data(1:3,1:3,me) subFp0 = phase_mechanical_Fp0(ph)%data(1:3,1:3,en)
subFi0 = phase_mechanical_Fi0(ph)%data(1:3,1:3,me) subFi0 = phase_mechanical_Fi0(ph)%data(1:3,1:3,en)
subF0 = phase_mechanical_F0(ph)%data(1:3,1:3,me) subF0 = phase_mechanical_F0(ph)%data(1:3,1:3,en)
subFrac = 0.0_pReal subFrac = 0.0_pReal
subStep = 1.0_pReal/num%subStepSizeCryst subStep = 1.0_pReal/num%subStepSizeCryst
todo = .true. todo = .true.
@ -1161,29 +1159,29 @@ module function crystallite_stress(dt,co,ip,el) result(converged_)
if (todo) then if (todo) then
subF0 = subF subF0 = subF
subLp0 = phase_mechanical_Lp(ph)%data(1:3,1:3,me) subLp0 = phase_mechanical_Lp(ph)%data(1:3,1:3,en)
subLi0 = phase_mechanical_Li(ph)%data(1:3,1:3,me) subLi0 = phase_mechanical_Li(ph)%data(1:3,1:3,en)
subFp0 = phase_mechanical_Fp(ph)%data(1:3,1:3,me) subFp0 = phase_mechanical_Fp(ph)%data(1:3,1:3,en)
subFi0 = phase_mechanical_Fi(ph)%data(1:3,1:3,me) subFi0 = phase_mechanical_Fi(ph)%data(1:3,1:3,en)
subState0 = plasticState(ph)%state(:,me) subState0 = plasticState(ph)%state(:,en)
if (damageState(ph)%sizeState > 0) & if (damageState(ph)%sizeState > 0) &
damageState(ph)%subState0(:,me) = damageState(ph)%state(:,me) damageState(ph)%subState0(:,en) = damageState(ph)%state(:,en)
endif endif
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! cut back (reduced time and restore) ! cut back (reduced time and restore)
else else
subStep = num%subStepSizeCryst * subStep subStep = num%subStepSizeCryst * subStep
phase_mechanical_Fp(ph)%data(1:3,1:3,me) = subFp0 phase_mechanical_Fp(ph)%data(1:3,1:3,en) = subFp0
phase_mechanical_Fi(ph)%data(1:3,1:3,me) = subFi0 phase_mechanical_Fi(ph)%data(1:3,1:3,en) = subFi0
phase_mechanical_S(ph)%data(1:3,1:3,me) = phase_mechanical_S0(ph)%data(1:3,1:3,me) ! why no subS0 ? is S0 of any use? phase_mechanical_S(ph)%data(1:3,1:3,en) = phase_mechanical_S0(ph)%data(1:3,1:3,en) ! why no subS0 ? is S0 of any use?
if (subStep < 1.0_pReal) then ! actual (not initial) cutback if (subStep < 1.0_pReal) then ! actual (not initial) cutback
phase_mechanical_Lp(ph)%data(1:3,1:3,me) = subLp0 phase_mechanical_Lp(ph)%data(1:3,1:3,en) = subLp0
phase_mechanical_Li(ph)%data(1:3,1:3,me) = subLi0 phase_mechanical_Li(ph)%data(1:3,1:3,en) = subLi0
endif endif
plasticState(ph)%state(:,me) = subState0 plasticState(ph)%state(:,en) = subState0
if (damageState(ph)%sizeState > 0) & if (damageState(ph)%sizeState > 0) &
damageState(ph)%state(:,me) = damageState(ph)%subState0(:,me) damageState(ph)%state(:,en) = damageState(ph)%subState0(:,en)
todo = subStep > num%subStepMinCryst ! still on track or already done (beyond repair) todo = subStep > num%subStepMinCryst ! still on track or already done (beyond repair)
endif endif
@ -1192,7 +1190,7 @@ module function crystallite_stress(dt,co,ip,el) result(converged_)
! prepare for integration ! prepare for integration
if (todo) then if (todo) then
subF = subF0 & subF = subF0 &
+ subStep * (phase_mechanical_F(ph)%data(1:3,1:3,me) - phase_mechanical_F0(ph)%data(1:3,1:3,me)) + subStep * (phase_mechanical_F(ph)%data(1:3,1:3,en) - phase_mechanical_F0(ph)%data(1:3,1:3,en))
converged_ = .not. integrateState(subF0,subF,subFp0,subFi0,subState0(1:sizeDotState),subStep * dt,co,ip,el) converged_ = .not. integrateState(subF0,subF,subFp0,subFi0,subState0(1:sizeDotState),subStep * dt,co,ip,el)
converged_ = converged_ .and. .not. integrateDamageState(subStep * dt,co,ip,el) converged_ = converged_ .and. .not. integrateDamageState(subStep * dt,co,ip,el)
endif endif
@ -1212,22 +1210,22 @@ module subroutine mechanical_restore(ce,includeL)
includeL !< protect agains fake cutback includeL !< protect agains fake cutback
integer :: & integer :: &
co, ph, me co, ph, en
do co = 1,homogenization_Nconstituents(material_homogenizationID(ce)) do co = 1,homogenization_Nconstituents(material_homogenizationID(ce))
ph = material_phaseID(co,ce) ph = material_phaseID(co,ce)
me = material_phaseEntry(co,ce) en = material_phaseEntry(co,ce)
if (includeL) then if (includeL) then
phase_mechanical_Lp(ph)%data(1:3,1:3,me) = phase_mechanical_Lp0(ph)%data(1:3,1:3,me) phase_mechanical_Lp(ph)%data(1:3,1:3,en) = phase_mechanical_Lp0(ph)%data(1:3,1:3,en)
phase_mechanical_Li(ph)%data(1:3,1:3,me) = phase_mechanical_Li0(ph)%data(1:3,1:3,me) phase_mechanical_Li(ph)%data(1:3,1:3,en) = phase_mechanical_Li0(ph)%data(1:3,1:3,en)
endif ! maybe protecting everything from overwriting makes more sense endif ! maybe protecting everything from overwriting makes more sense
phase_mechanical_Fp(ph)%data(1:3,1:3,me) = phase_mechanical_Fp0(ph)%data(1:3,1:3,me) phase_mechanical_Fp(ph)%data(1:3,1:3,en) = phase_mechanical_Fp0(ph)%data(1:3,1:3,en)
phase_mechanical_Fi(ph)%data(1:3,1:3,me) = phase_mechanical_Fi0(ph)%data(1:3,1:3,me) phase_mechanical_Fi(ph)%data(1:3,1:3,en) = phase_mechanical_Fi0(ph)%data(1:3,1:3,en)
phase_mechanical_S(ph)%data(1:3,1:3,me) = phase_mechanical_S0(ph)%data(1:3,1:3,me) phase_mechanical_S(ph)%data(1:3,1:3,en) = phase_mechanical_S0(ph)%data(1:3,1:3,en)
plasticState(ph)%state(:,me) = plasticState(ph)%State0(:,me) plasticState(ph)%state(:,en) = plasticState(ph)%State0(:,en)
enddo enddo
end subroutine mechanical_restore end subroutine mechanical_restore
@ -1245,7 +1243,7 @@ module function phase_mechanical_dPdF(dt,co,ce) result(dPdF)
integer :: & integer :: &
o, & o, &
p, ph, me p, ph, en
real(pReal), dimension(3,3) :: devNull, & real(pReal), dimension(3,3) :: devNull, &
invSubFp0,invSubFi0,invFp,invFi, & invSubFp0,invSubFi0,invFp,invFi, &
temp_33_1, temp_33_2, temp_33_3 temp_33_1, temp_33_2, temp_33_3
@ -1266,20 +1264,20 @@ module function phase_mechanical_dPdF(dt,co,ce) result(dPdF)
ph = material_phaseID(co,ce) ph = material_phaseID(co,ce)
me = material_phaseEntry(co,ce) en = material_phaseEntry(co,ce)
call phase_hooke_SandItsTangents(devNull,dSdFe,dSdFi, & call phase_hooke_SandItsTangents(devNull,dSdFe,dSdFi, &
phase_mechanical_Fe(ph)%data(1:3,1:3,me), & phase_mechanical_Fe(ph)%data(1:3,1:3,en), &
phase_mechanical_Fi(ph)%data(1:3,1:3,me),ph,me) phase_mechanical_Fi(ph)%data(1:3,1:3,en),ph,en)
call phase_LiAndItsTangents(devNull,dLidS,dLidFi, & call phase_LiAndItsTangents(devNull,dLidS,dLidFi, &
phase_mechanical_S(ph)%data(1:3,1:3,me), & phase_mechanical_S(ph)%data(1:3,1:3,en), &
phase_mechanical_Fi(ph)%data(1:3,1:3,me), & phase_mechanical_Fi(ph)%data(1:3,1:3,en), &
ph,me) ph,en)
invFp = math_inv33(phase_mechanical_Fp(ph)%data(1:3,1:3,me)) invFp = math_inv33(phase_mechanical_Fp(ph)%data(1:3,1:3,en))
invFi = math_inv33(phase_mechanical_Fi(ph)%data(1:3,1:3,me)) invFi = math_inv33(phase_mechanical_Fi(ph)%data(1:3,1:3,en))
invSubFp0 = math_inv33(phase_mechanical_Fp0(ph)%data(1:3,1:3,me)) invSubFp0 = math_inv33(phase_mechanical_Fp0(ph)%data(1:3,1:3,en))
invSubFi0 = math_inv33(phase_mechanical_Fi0(ph)%data(1:3,1:3,me)) invSubFi0 = math_inv33(phase_mechanical_Fi0(ph)%data(1:3,1:3,en))
if (sum(abs(dLidS)) < tol_math_check) then if (sum(abs(dLidS)) < tol_math_check) then
dFidS = 0.0_pReal dFidS = 0.0_pReal
@ -1305,15 +1303,15 @@ module function phase_mechanical_dPdF(dt,co,ce) result(dPdF)
endif endif
call plastic_LpAndItsTangents(devNull,dLpdS,dLpdFi, & call plastic_LpAndItsTangents(devNull,dLpdS,dLpdFi, &
phase_mechanical_S(ph)%data(1:3,1:3,me), & phase_mechanical_S(ph)%data(1:3,1:3,en), &
phase_mechanical_Fi(ph)%data(1:3,1:3,me),ph,me) phase_mechanical_Fi(ph)%data(1:3,1:3,en),ph,en)
dLpdS = math_mul3333xx3333(dLpdFi,dFidS) + dLpdS dLpdS = math_mul3333xx3333(dLpdFi,dFidS) + dLpdS
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! calculate dSdF ! calculate dSdF
temp_33_1 = transpose(matmul(invFp,invFi)) temp_33_1 = transpose(matmul(invFp,invFi))
temp_33_2 = matmul(phase_mechanical_F(ph)%data(1:3,1:3,me),invSubFp0) temp_33_2 = matmul(phase_mechanical_F(ph)%data(1:3,1:3,en),invSubFp0)
temp_33_3 = matmul(matmul(phase_mechanical_F(ph)%data(1:3,1:3,me),invFp), invSubFi0) temp_33_3 = matmul(matmul(phase_mechanical_F(ph)%data(1:3,1:3,en),invFp), invSubFi0)
do o=1,3; do p=1,3 do o=1,3; do p=1,3
rhs_3333(p,o,1:3,1:3) = matmul(dSdFe(p,o,1:3,1:3),temp_33_1) rhs_3333(p,o,1:3,1:3) = matmul(dSdFe(p,o,1:3,1:3),temp_33_1)
@ -1341,9 +1339,9 @@ module function phase_mechanical_dPdF(dt,co,ce) result(dPdF)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! assemble dPdF ! assemble dPdF
temp_33_1 = matmul(phase_mechanical_S(ph)%data(1:3,1:3,me),transpose(invFp)) temp_33_1 = matmul(phase_mechanical_S(ph)%data(1:3,1:3,en),transpose(invFp))
temp_33_2 = matmul(phase_mechanical_F(ph)%data(1:3,1:3,me),invFp) temp_33_2 = matmul(phase_mechanical_F(ph)%data(1:3,1:3,en),invFp)
temp_33_3 = matmul(temp_33_2,phase_mechanical_S(ph)%data(1:3,1:3,me)) temp_33_3 = matmul(temp_33_2,phase_mechanical_S(ph)%data(1:3,1:3,en))
dPdF = 0.0_pReal dPdF = 0.0_pReal
do p=1,3 do p=1,3
@ -1351,7 +1349,7 @@ module function phase_mechanical_dPdF(dt,co,ce) result(dPdF)
enddo enddo
do o=1,3; do p=1,3 do o=1,3; do p=1,3
dPdF(1:3,1:3,p,o) = dPdF(1:3,1:3,p,o) & dPdF(1:3,1:3,p,o) = dPdF(1:3,1:3,p,o) &
+ matmul(matmul(phase_mechanical_F(ph)%data(1:3,1:3,me),dFpinvdF(1:3,1:3,p,o)),temp_33_1) & + matmul(matmul(phase_mechanical_F(ph)%data(1:3,1:3,en),dFpinvdF(1:3,1:3,p,o)),temp_33_1) &
+ matmul(matmul(temp_33_2,dSdF(1:3,1:3,p,o)),transpose(invFp)) & + matmul(matmul(temp_33_2,dSdF(1:3,1:3,p,o)),transpose(invFp)) &
+ matmul(temp_33_3,transpose(dFpinvdF(1:3,1:3,p,o))) + matmul(temp_33_3,transpose(dFpinvdF(1:3,1:3,p,o)))
enddo; enddo enddo; enddo
@ -1396,13 +1394,13 @@ end subroutine mechanical_restartRead
!---------------------------------------------------------------------------------------------- !----------------------------------------------------------------------------------------------
!< @brief Get first Piola-Kichhoff stress (for use by non-mech physics) !< @brief Get first Piola-Kichhoff stress (for use by non-mech physics)
!---------------------------------------------------------------------------------------------- !----------------------------------------------------------------------------------------------
module function mechanical_S(ph,me) result(S) module function mechanical_S(ph,en) result(S)
integer, intent(in) :: ph,me integer, intent(in) :: ph,en
real(pReal), dimension(3,3) :: S real(pReal), dimension(3,3) :: S
S = phase_mechanical_S(ph)%data(1:3,1:3,me) S = phase_mechanical_S(ph)%data(1:3,1:3,en)
end function mechanical_S end function mechanical_S
@ -1410,13 +1408,13 @@ end function mechanical_S
!---------------------------------------------------------------------------------------------- !----------------------------------------------------------------------------------------------
!< @brief Get plastic velocity gradient (for use by non-mech physics) !< @brief Get plastic velocity gradient (for use by non-mech physics)
!---------------------------------------------------------------------------------------------- !----------------------------------------------------------------------------------------------
module function mechanical_L_p(ph,me) result(L_p) module function mechanical_L_p(ph,en) result(L_p)
integer, intent(in) :: ph,me integer, intent(in) :: ph,en
real(pReal), dimension(3,3) :: L_p real(pReal), dimension(3,3) :: L_p
L_p = phase_mechanical_Lp(ph)%data(1:3,1:3,me) L_p = phase_mechanical_Lp(ph)%data(1:3,1:3,en)
end function mechanical_L_p end function mechanical_L_p
@ -1424,13 +1422,13 @@ end function mechanical_L_p
!---------------------------------------------------------------------------------------------- !----------------------------------------------------------------------------------------------
!< @brief Get elastic deformation gradient (for use by non-mech physics) !< @brief Get elastic deformation gradient (for use by non-mech physics)
!---------------------------------------------------------------------------------------------- !----------------------------------------------------------------------------------------------
module function mechanical_F_e(ph,me) result(F_e) module function mechanical_F_e(ph,en) result(F_e)
integer, intent(in) :: ph,me integer, intent(in) :: ph,en
real(pReal), dimension(3,3) :: F_e real(pReal), dimension(3,3) :: F_e
F_e = phase_mechanical_Fe(ph)%data(1:3,1:3,me) F_e = phase_mechanical_Fe(ph)%data(1:3,1:3,en)
end function mechanical_F_e end function mechanical_F_e

12
src/phase_mechanical_eigen.f90
View File

@ -145,10 +145,10 @@ end function kinematics_active2
! ToDo: MD: S is Mi? ! ToDo: MD: S is Mi?
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module subroutine phase_LiAndItsTangents(Li, dLi_dS, dLi_dFi, & module subroutine phase_LiAndItsTangents(Li, dLi_dS, dLi_dFi, &
S, Fi, ph,me) S, Fi, ph,en)
integer, intent(in) :: & integer, intent(in) :: &
ph,me ph,en
real(pReal), intent(in), dimension(3,3) :: & real(pReal), intent(in), dimension(3,3) :: &
S !< 2nd Piola-Kirchhoff stress S !< 2nd Piola-Kirchhoff stress
real(pReal), intent(in), dimension(3,3) :: & real(pReal), intent(in), dimension(3,3) :: &
@ -179,7 +179,7 @@ module subroutine phase_LiAndItsTangents(Li, dLi_dS, dLi_dFi, &
plasticType: select case (phase_plasticity(ph)) plasticType: select case (phase_plasticity(ph))
case (PLASTICITY_isotropic_ID) plasticType case (PLASTICITY_isotropic_ID) plasticType
call plastic_isotropic_LiAndItsTangent(my_Li, my_dLi_dS, S ,ph,me) call plastic_isotropic_LiAndItsTangent(my_Li, my_dLi_dS, S ,ph,en)
Li = Li + my_Li Li = Li + my_Li
dLi_dS = dLi_dS + my_dLi_dS dLi_dS = dLi_dS + my_dLi_dS
active = .true. active = .true.
@ -188,7 +188,7 @@ module subroutine phase_LiAndItsTangents(Li, dLi_dS, dLi_dFi, &
KinematicsLoop: do k = 1, Nmodels(ph) KinematicsLoop: do k = 1, Nmodels(ph)
kinematicsType: select case (model(k,ph)) kinematicsType: select case (model(k,ph))
case (KINEMATICS_thermal_expansion_ID) kinematicsType case (KINEMATICS_thermal_expansion_ID) kinematicsType
call thermalexpansion_LiAndItsTangent(my_Li, my_dLi_dS, ph,me) call thermalexpansion_LiAndItsTangent(my_Li, my_dLi_dS, ph,en)
Li = Li + my_Li Li = Li + my_Li
dLi_dS = dLi_dS + my_dLi_dS dLi_dS = dLi_dS + my_dLi_dS
active = .true. active = .true.
@ -197,12 +197,12 @@ module subroutine phase_LiAndItsTangents(Li, dLi_dS, dLi_dFi, &
select case (model_damage(ph)) select case (model_damage(ph))
case (KINEMATICS_cleavage_opening_ID) case (KINEMATICS_cleavage_opening_ID)
call damage_anisobrittle_LiAndItsTangent(my_Li, my_dLi_dS, S, ph, me) call damage_anisobrittle_LiAndItsTangent(my_Li, my_dLi_dS, S, ph, en)
Li = Li + my_Li Li = Li + my_Li
dLi_dS = dLi_dS + my_dLi_dS dLi_dS = dLi_dS + my_dLi_dS
active = .true. active = .true.
case (KINEMATICS_slipplane_opening_ID) case (KINEMATICS_slipplane_opening_ID)
call damage_isoductile_LiAndItsTangent(my_Li, my_dLi_dS, S, ph, me) call damage_isoductile_LiAndItsTangent(my_Li, my_dLi_dS, S, ph, en)
Li = Li + my_Li Li = Li + my_Li
dLi_dS = dLi_dS + my_dLi_dS dLi_dS = dLi_dS + my_dLi_dS
active = .true. active = .true.

60
src/phase_mechanical_plastic.f90
View File

@ -240,9 +240,9 @@ end subroutine plastic_init
! Mp in, dLp_dMp out ! Mp in, dLp_dMp out
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module subroutine plastic_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, & module subroutine plastic_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, &
S, Fi, ph,me) S, Fi, ph,en)
integer, intent(in) :: & integer, intent(in) :: &
ph,me ph,en
real(pReal), intent(in), dimension(3,3) :: & real(pReal), intent(in), dimension(3,3) :: &
S, & !< 2nd Piola-Kirchhoff stress S, & !< 2nd Piola-Kirchhoff stress
Fi !< intermediate deformation gradient Fi !< intermediate deformation gradient
@ -250,7 +250,7 @@ module subroutine plastic_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, &
Lp !< plastic velocity gradient Lp !< plastic velocity gradient
real(pReal), intent(out), dimension(3,3,3,3) :: & real(pReal), intent(out), dimension(3,3,3,3) :: &
dLp_dS, & dLp_dS, &
dLp_dFi !< derivative me Lp with respect to Fi dLp_dFi !< derivative en Lp with respect to Fi
real(pReal), dimension(3,3,3,3) :: & real(pReal), dimension(3,3,3,3) :: &
dLp_dMp !< derivative of Lp with respect to Mandel stress dLp_dMp !< derivative of Lp with respect to Mandel stress
@ -270,22 +270,22 @@ module subroutine plastic_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, &
dLp_dMp = 0.0_pReal dLp_dMp = 0.0_pReal
case (PLASTICITY_ISOTROPIC_ID) plasticType case (PLASTICITY_ISOTROPIC_ID) plasticType
call isotropic_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,me) call isotropic_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,en)
case (PLASTICITY_PHENOPOWERLAW_ID) plasticType case (PLASTICITY_PHENOPOWERLAW_ID) plasticType
call phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,me) call phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,en)
case (PLASTICITY_KINEHARDENING_ID) plasticType case (PLASTICITY_KINEHARDENING_ID) plasticType
call kinehardening_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,me) call kinehardening_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,en)
case (PLASTICITY_NONLOCAL_ID) plasticType case (PLASTICITY_NONLOCAL_ID) plasticType
call nonlocal_LpAndItsTangent(Lp,dLp_dMp,Mp, thermal_T(ph,me),ph,me) call nonlocal_LpAndItsTangent(Lp,dLp_dMp,Mp, thermal_T(ph,en),ph,en)
case (PLASTICITY_DISLOTWIN_ID) plasticType case (PLASTICITY_DISLOTWIN_ID) plasticType
call dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp, thermal_T(ph,me),ph,me) call dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp, thermal_T(ph,en),ph,en)
case (PLASTICITY_DISLOTUNGSTEN_ID) plasticType case (PLASTICITY_DISLOTUNGSTEN_ID) plasticType
call dislotungsten_LpAndItsTangent(Lp,dLp_dMp,Mp, thermal_T(ph,me),ph,me) call dislotungsten_LpAndItsTangent(Lp,dLp_dMp,Mp, thermal_T(ph,en),ph,en)
end select plasticType end select plasticType
@ -301,14 +301,14 @@ end subroutine plastic_LpAndItsTangents
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief contains the constitutive equation for calculating the rate of change of microstructure !> @brief contains the constitutive equation for calculating the rate of change of microstructure
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module function plastic_dotState(subdt,co,ip,el,ph,me) result(broken) module function plastic_dotState(subdt,co,ip,el,ph,en) result(broken)
integer, intent(in) :: & integer, intent(in) :: &
co, & !< component-ID of integration point co, & !< component-ID of integration point
ip, & !< integration point ip, & !< integration point
el, & !< element el, & !< element
ph, & ph, &
me en
real(pReal), intent(in) :: & real(pReal), intent(in) :: &
subdt !< timestep subdt !< timestep
real(pReal), dimension(3,3) :: & real(pReal), dimension(3,3) :: &
@ -316,30 +316,30 @@ module function plastic_dotState(subdt,co,ip,el,ph,me) result(broken)
logical :: broken logical :: broken
Mp = matmul(matmul(transpose(phase_mechanical_Fi(ph)%data(1:3,1:3,me)),& Mp = matmul(matmul(transpose(phase_mechanical_Fi(ph)%data(1:3,1:3,en)),&
phase_mechanical_Fi(ph)%data(1:3,1:3,me)),phase_mechanical_S(ph)%data(1:3,1:3,me)) phase_mechanical_Fi(ph)%data(1:3,1:3,en)),phase_mechanical_S(ph)%data(1:3,1:3,en))
plasticType: select case (phase_plasticity(ph)) plasticType: select case (phase_plasticity(ph))
case (PLASTICITY_ISOTROPIC_ID) plasticType case (PLASTICITY_ISOTROPIC_ID) plasticType
call isotropic_dotState(Mp,ph,me) call isotropic_dotState(Mp,ph,en)
case (PLASTICITY_PHENOPOWERLAW_ID) plasticType case (PLASTICITY_PHENOPOWERLAW_ID) plasticType
call phenopowerlaw_dotState(Mp,ph,me) call phenopowerlaw_dotState(Mp,ph,en)
case (PLASTICITY_KINEHARDENING_ID) plasticType case (PLASTICITY_KINEHARDENING_ID) plasticType
call plastic_kinehardening_dotState(Mp,ph,me) call plastic_kinehardening_dotState(Mp,ph,en)
case (PLASTICITY_DISLOTWIN_ID) plasticType case (PLASTICITY_DISLOTWIN_ID) plasticType
call dislotwin_dotState(Mp,thermal_T(ph,me),ph,me) call dislotwin_dotState(Mp,thermal_T(ph,en),ph,en)
case (PLASTICITY_DISLOTUNGSTEN_ID) plasticType case (PLASTICITY_DISLOTUNGSTEN_ID) plasticType
call dislotungsten_dotState(Mp,thermal_T(ph,me),ph,me) call dislotungsten_dotState(Mp,thermal_T(ph,en),ph,en)
case (PLASTICITY_NONLOCAL_ID) plasticType case (PLASTICITY_NONLOCAL_ID) plasticType
call nonlocal_dotState(Mp,thermal_T(ph,me),subdt,ph,me,ip,el) call nonlocal_dotState(Mp,thermal_T(ph,en),subdt,ph,en,ip,el)
end select plasticType end select plasticType
broken = any(IEEE_is_NaN(plasticState(ph)%dotState(:,me))) broken = any(IEEE_is_NaN(plasticState(ph)%dotState(:,en)))
end function plastic_dotState end function plastic_dotState
@ -383,11 +383,11 @@ end subroutine plastic_dependentState
!> @brief for constitutive models having an instantaneous change of state !> @brief for constitutive models having an instantaneous change of state
!> will return false if delta state is not needed/supported by the constitutive model !> will return false if delta state is not needed/supported by the constitutive model
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module function plastic_deltaState(ph, me) result(broken) module function plastic_deltaState(ph, en) result(broken)
integer, intent(in) :: & integer, intent(in) :: &
ph, & ph, &
me en
logical :: & logical :: &
broken broken
@ -398,18 +398,18 @@ module function plastic_deltaState(ph, me) result(broken)
mySize mySize
Mp = matmul(matmul(transpose(phase_mechanical_Fi(ph)%data(1:3,1:3,me)),& Mp = matmul(matmul(transpose(phase_mechanical_Fi(ph)%data(1:3,1:3,en)),&
phase_mechanical_Fi(ph)%data(1:3,1:3,me)),phase_mechanical_S(ph)%data(1:3,1:3,me)) phase_mechanical_Fi(ph)%data(1:3,1:3,en)),phase_mechanical_S(ph)%data(1:3,1:3,en))
plasticType: select case (phase_plasticity(ph)) plasticType: select case (phase_plasticity(ph))
case (PLASTICITY_KINEHARDENING_ID) plasticType case (PLASTICITY_KINEHARDENING_ID) plasticType
call plastic_kinehardening_deltaState(Mp,ph,me) call plastic_kinehardening_deltaState(Mp,ph,en)
broken = any(IEEE_is_NaN(plasticState(ph)%deltaState(:,me))) broken = any(IEEE_is_NaN(plasticState(ph)%deltaState(:,en)))
case (PLASTICITY_NONLOCAL_ID) plasticType case (PLASTICITY_NONLOCAL_ID) plasticType
call plastic_nonlocal_deltaState(Mp,ph,me) call plastic_nonlocal_deltaState(Mp,ph,en)
broken = any(IEEE_is_NaN(plasticState(ph)%deltaState(:,me))) broken = any(IEEE_is_NaN(plasticState(ph)%deltaState(:,en)))
case default case default
broken = .false. broken = .false.
@ -422,8 +422,8 @@ module function plastic_deltaState(ph, me) result(broken)
myOffset = plasticState(ph)%offsetDeltaState myOffset = plasticState(ph)%offsetDeltaState
mySize = plasticState(ph)%sizeDeltaState mySize = plasticState(ph)%sizeDeltaState
plasticState(ph)%state(myOffset + 1:myOffset + mySize,me) = & plasticState(ph)%state(myOffset + 1:myOffset + mySize,en) = &
plasticState(ph)%state(myOffset + 1:myOffset + mySize,me) + plasticState(ph)%deltaState(1:mySize,me) plasticState(ph)%state(myOffset + 1:myOffset + mySize,en) + plasticState(ph)%deltaState(1:mySize,en)
end select end select
endif endif

12
src/phase_mechanical_plastic_dislotwin.f90
View File

@ -483,10 +483,10 @@ end function plastic_dislotwin_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Return the homogenized elasticity matrix. !> @brief Return the homogenized elasticity matrix.
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module function plastic_dislotwin_homogenizedC(ph,me) result(homogenizedC) module function plastic_dislotwin_homogenizedC(ph,en) result(homogenizedC)
integer, intent(in) :: & integer, intent(in) :: &
ph, me ph, en
real(pReal), dimension(6,6) :: & real(pReal), dimension(6,6) :: &
homogenizedC homogenizedC
@ -498,17 +498,17 @@ module function plastic_dislotwin_homogenizedC(ph,me) result(homogenizedC)
stt => state(ph)) stt => state(ph))
f_unrotated = 1.0_pReal & f_unrotated = 1.0_pReal &
- sum(stt%f_tw(1:prm%sum_N_tw,me)) & - sum(stt%f_tw(1:prm%sum_N_tw,en)) &
- sum(stt%f_tr(1:prm%sum_N_tr,me)) - sum(stt%f_tr(1:prm%sum_N_tr,en))
homogenizedC = f_unrotated * prm%C66 homogenizedC = f_unrotated * prm%C66
do i=1,prm%sum_N_tw do i=1,prm%sum_N_tw
homogenizedC = homogenizedC & homogenizedC = homogenizedC &
+ stt%f_tw(i,me)*prm%C66_tw(1:6,1:6,i) + stt%f_tw(i,en)*prm%C66_tw(1:6,1:6,i)
enddo enddo
do i=1,prm%sum_N_tr do i=1,prm%sum_N_tr
homogenizedC = homogenizedC & homogenizedC = homogenizedC &
+ stt%f_tr(i,me)*prm%C66_tr(1:6,1:6,i) + stt%f_tr(i,en)*prm%C66_tr(1:6,1:6,i)
enddo enddo
end associate end associate

6
src/phase_mechanical_plastic_isotropic.f90
View File

@ -208,7 +208,7 @@ end subroutine isotropic_LpAndItsTangent
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Calculate inelastic velocity gradient and its tangent. !> @brief Calculate inelastic velocity gradient and its tangent.
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module subroutine plastic_isotropic_LiAndItsTangent(Li,dLi_dMi,Mi,ph,me) module subroutine plastic_isotropic_LiAndItsTangent(Li,dLi_dMi,Mi,ph,en)
real(pReal), dimension(3,3), intent(out) :: & real(pReal), dimension(3,3), intent(out) :: &
Li !< inleastic velocity gradient Li !< inleastic velocity gradient
@ -219,7 +219,7 @@ module subroutine plastic_isotropic_LiAndItsTangent(Li,dLi_dMi,Mi,ph,me)
Mi !< Mandel stress Mi !< Mandel stress
integer, intent(in) :: & integer, intent(in) :: &
ph, & ph, &
me en
real(pReal) :: & real(pReal) :: &
tr !< trace of spherical part of Mandel stress (= 3 x pressure) tr !< trace of spherical part of Mandel stress (= 3 x pressure)
@ -232,7 +232,7 @@ module subroutine plastic_isotropic_LiAndItsTangent(Li,dLi_dMi,Mi,ph,me)
if (prm%dilatation .and. abs(tr) > 0.0_pReal) then ! no stress or J2 plasticity --> Li and its derivative are zero if (prm%dilatation .and. abs(tr) > 0.0_pReal) then ! no stress or J2 plasticity --> Li and its derivative are zero
Li = math_I3 & Li = math_I3 &
* prm%dot_gamma_0/prm%M * (3.0_pReal*prm%M*stt%xi(me))**(-prm%n) & * prm%dot_gamma_0/prm%M * (3.0_pReal*prm%M*stt%xi(en))**(-prm%n) &
* tr * abs(tr)**(prm%n-1.0_pReal) * tr * abs(tr)**(prm%n-1.0_pReal)
forall (k=1:3,l=1:3,m=1:3,n=1:3) dLi_dMi(k,l,m,n) = prm%n / tr * Li(k,l) * math_I3(m,n) forall (k=1:3,l=1:3,m=1:3,n=1:3) dLi_dMi(k,l,m,n) = prm%n / tr * Li(k,l) * math_I3(m,n)
else else