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

naming reflect change from Piola-Kirchhoff to Mandel stress 

not using the symmetrized stress anymore to avoid handling of
symmetrized Schmid matrizes. The time saved when calculating the double
contraction is probably anyway lost during the conversion from (3,3) to
(6) of Mstar
This commit is contained in:
Martin Diehl 2018-08-03 14:53:40 +02:00
parent 0a6bde70c5
commit 5cf2973715
2 changed files with 39 additions and 27 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

8
src/constitutive.f90
View File

@ -482,7 +482,8 @@ subroutine constitutive_LpAndItsTangent(Lp, dLp_dTstar3333, dLp_dFi3333, Tstar_v
real(pReal), dimension(9,9) :: & real(pReal), dimension(9,9) :: &
dLp_dMstar !< derivative of Lp with respect to Mstar (4th-order tensor) dLp_dMstar !< derivative of Lp with respect to Mstar (4th-order tensor)
real(pReal), dimension(3,3) :: & real(pReal), dimension(3,3) :: &
temp_33 temp_33, &
Mstar
integer(pInt) :: & integer(pInt) :: &
ho, & !< homogenization ho, & !< homogenization
tme !< thermal member position tme !< thermal member position
@ -492,7 +493,8 @@ subroutine constitutive_LpAndItsTangent(Lp, dLp_dTstar3333, dLp_dFi3333, Tstar_v
ho = material_homog(ip,el) ho = material_homog(ip,el)
tme = thermalMapping(ho)%p(ip,el) tme = thermalMapping(ho)%p(ip,el)
Mstar_v = math_Mandel33to6(math_mul33x33(math_mul33x33(transpose(Fi),Fi),math_Mandel6to33(Tstar_v))) Mstar = math_mul33x33(math_mul33x33(transpose(Fi),Fi),math_Mandel6to33(Tstar_v))
Mstar_v = math_Mandel33to6(Mstar)
plasticityType: select case (phase_plasticity(material_phase(ipc,ip,el))) plasticityType: select case (phase_plasticity(material_phase(ipc,ip,el)))
case (PLASTICITY_NONE_ID) plasticityType case (PLASTICITY_NONE_ID) plasticityType
@ -501,7 +503,7 @@ subroutine constitutive_LpAndItsTangent(Lp, dLp_dTstar3333, dLp_dFi3333, Tstar_v
case (PLASTICITY_ISOTROPIC_ID) plasticityType case (PLASTICITY_ISOTROPIC_ID) plasticityType
call plastic_isotropic_LpAndItsTangent (Lp,dLp_dMstar,Mstar_v,ipc,ip,el) call plastic_isotropic_LpAndItsTangent (Lp,dLp_dMstar,Mstar_v,ipc,ip,el)
case (PLASTICITY_PHENOPOWERLAW_ID) plasticityType case (PLASTICITY_PHENOPOWERLAW_ID) plasticityType
call plastic_phenopowerlaw_LpAndItsTangent (Lp,dLp_dMstar,Mstar_v,ipc,ip,el) call plastic_phenopowerlaw_LpAndItsTangent (Lp,dLp_dMstar,Mstar,ipc,ip,el)
case (PLASTICITY_KINEHARDENING_ID) plasticityType case (PLASTICITY_KINEHARDENING_ID) plasticityType
call plastic_kinehardening_LpAndItsTangent (Lp,dLp_dMstar,Mstar_v,ipc,ip,el) call plastic_kinehardening_LpAndItsTangent (Lp,dLp_dMstar,Mstar_v,ipc,ip,el)
case (PLASTICITY_NONLOCAL_ID) plasticityType case (PLASTICITY_NONLOCAL_ID) plasticityType

56
src/plastic_phenopowerlaw.f90
View File

@ -320,7 +320,7 @@ subroutine plastic_phenopowerlaw_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! allocate state arrays ! allocate state arrays
NipcMyPhase = count(material_phase == p) ! number of IPCs containing my phase NipcMyPhase = count(material_phase == p) ! number of IPCs containing my phase
sizeState = size(['tau_slip ','accshear_slip']) * prm%TotalNslip & sizeState = size(['tau_slip ','accshear_slip']) * prm%TotalNslip &
+ size(['tau_twin ','accshear_twin']) * prm%TotalNtwin & + size(['tau_twin ','accshear_twin']) * prm%TotalNtwin &
+ size(['sum(gamma)', 'sum(f) ']) + size(['sum(gamma)', 'sum(f) '])
@ -357,6 +357,16 @@ subroutine plastic_phenopowerlaw_init
index_myFamily = sum(prm%Nslip(1:f-1_pInt)) index_myFamily = sum(prm%Nslip(1:f-1_pInt))
mySlipSystems: do j = 1_pInt,prm%Nslip(f) mySlipSystems: do j = 1_pInt,prm%Nslip(f)
!prm%Schmid_pos(1:3,1:3,index_myFamily+j) = lattice_Sslip(1:3,1:3,1,index_myFamily+j,p)
!prm%Schmid_neg(1:3,1:3,index_myFamily+j) = lattice_Sslip(1:3,1:3,2,index_myFamily+j,p)
!do k = 1,size(prm%nonSchmidCoeff)
! prm%nonSchmid_pos(1:3,1:3,k,index_myFamily+j) = lattice_Sslip(1:3,1:3,2*k, index_myFamily+j,p)
! prm%nonSchmid_neg(1:3,1:3,k,index_myFamily+j) = lattice_Sslip(1:3,1:3,2*k+1,index_myFamily+j,p)
! !nonSchmid_tensor(1:3,1:3,1) = nonSchmid_tensor(1:3,1:3,1) + prm%nonSchmidCoeff(k)*&
! ! lattice_Sslip(1:3,1:3,2*k,index_myFamily+i,ph)
! !nonSchmid_tensor(1:3,1:3,2) = nonSchmid_tensor(1:3,1:3,2) + prm%nonSchmidCoeff(k)*&
! ! lattice_Sslip(1:3,1:3,2*k+1,index_myFamily+i,ph)
!enddo
otherSlipFamilies: do o = 1_pInt,size(prm%Nslip,1) otherSlipFamilies: do o = 1_pInt,size(prm%Nslip,1)
index_otherFamily = sum(prm%Nslip(1:o-1_pInt)) index_otherFamily = sum(prm%Nslip(1:o-1_pInt))
otherSlipSystems: do k = 1_pInt,prm%Nslip(o) otherSlipSystems: do k = 1_pInt,prm%Nslip(o)
@ -467,17 +477,16 @@ end subroutine plastic_phenopowerlaw_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief calculates plastic velocity gradient and its tangent !> @brief calculates plastic velocity gradient and its tangent
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,ipc,ip,el) subroutine plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dMstar99,Mstar,ipc,ip,el)
use prec, only: & use prec, only: &
dNeq0 dNeq0
use math, only: & use math, only: &
math_mul33xx33,&
math_Plain3333to99, & math_Plain3333to99, &
math_Mandel6to33 math_Mandel6to33
use lattice, only: & use lattice, only: &
lattice_Sslip, & lattice_Sslip, &
lattice_Sslip_v, &
lattice_Stwin, & lattice_Stwin, &
lattice_Stwin_v, &
lattice_NslipSystem, & lattice_NslipSystem, &
lattice_NtwinSystem lattice_NtwinSystem
use material, only: & use material, only: &
@ -489,14 +498,14 @@ subroutine plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,ipc,ip,
real(pReal), dimension(3,3), intent(out) :: & real(pReal), dimension(3,3), intent(out) :: &
Lp !< plastic velocity gradient Lp !< plastic velocity gradient
real(pReal), dimension(9,9), intent(out) :: & real(pReal), dimension(9,9), intent(out) :: &
dLp_dTstar99 !< derivative of Lp with respect to 2nd Piola Kirchhoff stress dLp_dMstar99 !< derivative of Lp with respect to the Mandel stress
integer(pInt), intent(in) :: & integer(pInt), intent(in) :: &
ipc, & !< component-ID of integration point ipc, & !< component-ID of integration point
ip, & !< integration point ip, & !< integration point
el !< element el !< element
real(pReal), dimension(6), intent(in) :: & real(pReal), dimension(3,3), intent(in) :: &
Tstar_v !< 2nd Piola Kirchhoff stress tensor in Mandel notation Mstar !< Mandel stress
integer(pInt) :: & integer(pInt) :: &
index_myFamily, & index_myFamily, &
@ -509,7 +518,7 @@ subroutine plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,ipc,ip,
dgdot_dtauslip_pos,dgdot_dtauslip_neg, & dgdot_dtauslip_pos,dgdot_dtauslip_neg, &
gdot_twin,dgdot_dtautwin,tau_twin gdot_twin,dgdot_dtautwin,tau_twin
real(pReal), dimension(3,3,3,3) :: & real(pReal), dimension(3,3,3,3) :: &
dLp_dTstar3333 !< derivative of Lp with respect to Tstar as 4th order tensor dLp_dMstar3333 !< derivative of Lp with respect to Mstar as 4th order tensor
real(pReal), dimension(3,3,2) :: & real(pReal), dimension(3,3,2) :: &
nonSchmid_tensor nonSchmid_tensor
type(tParameters) :: prm type(tParameters) :: prm
@ -518,30 +527,31 @@ subroutine plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,ipc,ip,
of = phasememberAt(ipc,ip,el) of = phasememberAt(ipc,ip,el)
ph = material_phase(ipc,ip,el) ph = material_phase(ipc,ip,el)
associate(prm => param(phase_plasticityInstance(ph)), stt => state(phase_plasticityInstance(ph))) associate(prm => param(phase_plasticityInstance(ph)),&
stt => state(phase_plasticityInstance(ph)))
Lp = 0.0_pReal Lp = 0.0_pReal
dLp_dTstar3333 = 0.0_pReal dLp_dMstar3333 = 0.0_pReal
dLp_dTstar99 = 0.0_pReal dLp_dMstar99 = 0.0_pReal
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! Slip part ! Slip part
j = 0_pInt j = 0_pInt
slipFamilies: do f = 1_pInt,size(prm%Nslip,1) slipFamilies: do f = 1_pInt,size(prm%Nslip,1)
index_myFamily = sum(lattice_NslipSystem(1:f-1_pInt,ph)) ! at which index starts my family index_myFamily = sum(lattice_NslipSystem(1:f-1_pInt,ph)) ! at which index starts my family
slipSystems: do i = 1_pInt,prm%Nslip(f) slipSystems: do i = 1_pInt,prm%Nslip(f)
j = j+1_pInt j = j+1_pInt
! Calculation of Lp ! Calculation of Lp
tau_slip_pos = dot_product(Tstar_v,lattice_Sslip_v(1:6,1,index_myFamily+i,ph)) tau_slip_pos = math_mul33xx33(Mstar,lattice_Sslip(1:3,1:3,1,index_myFamily+i,ph))
tau_slip_neg = tau_slip_pos tau_slip_neg = tau_slip_pos
nonSchmid_tensor(1:3,1:3,1) = lattice_Sslip(1:3,1:3,1,index_myFamily+i,ph) nonSchmid_tensor(1:3,1:3,1) = lattice_Sslip(1:3,1:3,1,index_myFamily+i,ph)
nonSchmid_tensor(1:3,1:3,2) = nonSchmid_tensor(1:3,1:3,1) nonSchmid_tensor(1:3,1:3,2) = nonSchmid_tensor(1:3,1:3,1)
do k = 1,size(prm%nonSchmidCoeff) do k = 1,size(prm%nonSchmidCoeff)
tau_slip_pos = tau_slip_pos + prm%nonSchmidCoeff(k)* & tau_slip_pos = tau_slip_pos + prm%nonSchmidCoeff(k)* &
dot_product(Tstar_v,lattice_Sslip_v(1:6,2*k,index_myFamily+i,ph)) math_mul33xx33(Mstar,lattice_Sslip(1:3,1:3,2*k,index_myFamily+i,ph))
tau_slip_neg = tau_slip_neg + prm%nonSchmidCoeff(k)* & tau_slip_neg = tau_slip_neg + prm%nonSchmidCoeff(k)* &
dot_product(Tstar_v,lattice_Sslip_v(1:6,2*k+1,index_myFamily+i,ph)) math_mul33xx33(Mstar,lattice_Sslip(1:3,1:3,2*k+1,index_myFamily+i,ph))
nonSchmid_tensor(1:3,1:3,1) = nonSchmid_tensor(1:3,1:3,1) + prm%nonSchmidCoeff(k)*& nonSchmid_tensor(1:3,1:3,1) = nonSchmid_tensor(1:3,1:3,1) + prm%nonSchmidCoeff(k)*&
lattice_Sslip(1:3,1:3,2*k,index_myFamily+i,ph) lattice_Sslip(1:3,1:3,2*k,index_myFamily+i,ph)
nonSchmid_tensor(1:3,1:3,2) = nonSchmid_tensor(1:3,1:3,2) + prm%nonSchmidCoeff(k)*& nonSchmid_tensor(1:3,1:3,2) = nonSchmid_tensor(1:3,1:3,2) + prm%nonSchmidCoeff(k)*&
@ -560,7 +570,7 @@ subroutine plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,ipc,ip,
if (dNeq0(tau_slip_pos)) then if (dNeq0(tau_slip_pos)) then
dgdot_dtauslip_pos = gdot_slip_pos*prm%n_slip/tau_slip_pos dgdot_dtauslip_pos = gdot_slip_pos*prm%n_slip/tau_slip_pos
forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) & forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) &
dLp_dTstar3333(k,l,m,n) = dLp_dTstar3333(k,l,m,n) + & dLp_dMstar3333(k,l,m,n) = dLp_dMstar3333(k,l,m,n) + &
dgdot_dtauslip_pos*lattice_Sslip(k,l,1,index_myFamily+i,ph)* & dgdot_dtauslip_pos*lattice_Sslip(k,l,1,index_myFamily+i,ph)* &
nonSchmid_tensor(m,n,1) nonSchmid_tensor(m,n,1)
endif endif
@ -568,7 +578,7 @@ subroutine plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,ipc,ip,
if (dNeq0(tau_slip_neg)) then if (dNeq0(tau_slip_neg)) then
dgdot_dtauslip_neg = gdot_slip_neg*prm%n_slip/tau_slip_neg dgdot_dtauslip_neg = gdot_slip_neg*prm%n_slip/tau_slip_neg
forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) & forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) &
dLp_dTstar3333(k,l,m,n) = dLp_dTstar3333(k,l,m,n) + & dLp_dMstar3333(k,l,m,n) = dLp_dMstar3333(k,l,m,n) + &
dgdot_dtauslip_neg*lattice_Sslip(k,l,1,index_myFamily+i,ph)* & dgdot_dtauslip_neg*lattice_Sslip(k,l,1,index_myFamily+i,ph)* &
nonSchmid_tensor(m,n,2) nonSchmid_tensor(m,n,2)
endif endif
@ -584,7 +594,7 @@ subroutine plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,ipc,ip,
j = j+1_pInt j = j+1_pInt
! Calculation of Lp ! Calculation of Lp
tau_twin = dot_product(Tstar_v,lattice_Stwin_v(1:6,index_myFamily+i,ph)) tau_twin = math_mul33xx33(Mstar,lattice_Stwin(1:3,1:3,index_myFamily+i,ph))
gdot_twin = (1.0_pReal-stt%sumF(of))*prm%gdot0_twin*& gdot_twin = (1.0_pReal-stt%sumF(of))*prm%gdot0_twin*&
(abs(tau_twin)/stt%s_twin(j,of))**& (abs(tau_twin)/stt%s_twin(j,of))**&
prm%n_twin*max(0.0_pReal,sign(1.0_pReal,tau_twin)) prm%n_twin*max(0.0_pReal,sign(1.0_pReal,tau_twin))
@ -594,14 +604,14 @@ subroutine plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,ipc,ip,
if (dNeq0(gdot_twin)) then if (dNeq0(gdot_twin)) then
dgdot_dtautwin = gdot_twin*prm%n_twin/tau_twin dgdot_dtautwin = gdot_twin*prm%n_twin/tau_twin
forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) & forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) &
dLp_dTstar3333(k,l,m,n) = dLp_dTstar3333(k,l,m,n) + & dLp_dMstar3333(k,l,m,n) = dLp_dMstar3333(k,l,m,n) + &
dgdot_dtautwin*lattice_Stwin(k,l,index_myFamily+i,ph)* & dgdot_dtautwin*lattice_Stwin(k,l,index_myFamily+i,ph)* &
lattice_Stwin(m,n,index_myFamily+i,ph) lattice_Stwin(m,n,index_myFamily+i,ph)
endif endif
enddo twinSystems enddo twinSystems
enddo twinFamilies enddo twinFamilies
dLp_dTstar99 = math_Plain3333to99(dLp_dTstar3333) dLp_dMstar99 = math_Plain3333to99(dLp_dMstar3333)
end associate end associate
end subroutine plastic_phenopowerlaw_LpAndItsTangent end subroutine plastic_phenopowerlaw_LpAndItsTangent
@ -650,9 +660,9 @@ subroutine plastic_phenopowerlaw_dotState(Tstar_v,ipc,ip,el)
of = phasememberAt(ipc,ip,el) of = phasememberAt(ipc,ip,el)
ph = material_phase(ipc,ip,el) ph = material_phase(ipc,ip,el)
associate( prm => param(phase_plasticityInstance(ph)), & associate(prm => param(phase_plasticityInstance(ph)), &
stt => state(phase_plasticityInstance(ph)), & stt => state(phase_plasticityInstance(ph)), &
dst => dotState(phase_plasticityInstance(ph))) dst => dotState(phase_plasticityInstance(ph)))
dst%whole(:,of) = 0.0_pReal dst%whole(:,of) = 0.0_pReal