simplified

This commit is contained in:
Martin Diehl 2020-02-11 17:47:48 +01:00
parent e212f91fac
commit 5b72110d0a
1 changed files with 25 additions and 47 deletions

View File

@ -818,7 +818,7 @@ logical function integrateStress(ipc,ip,el,timeFraction)
Fe, & ! elastic deformation gradient Fe, & ! elastic deformation gradient
temp_33 temp_33
real(pReal), dimension(9) :: work ! needed for matrix inversion by LAPACK real(pReal), dimension(9) :: work ! needed for matrix inversion by LAPACK
integer, dimension(9) :: devNull ! needed for matrix inversion by LAPACK integer, dimension(9) :: devNull_9 ! needed for matrix inversion by LAPACK
real(pReal), dimension(9,9) :: dRLp_dLp, & ! partial derivative of residuum (Jacobian for Newton-Raphson scheme) real(pReal), dimension(9,9) :: dRLp_dLp, & ! partial derivative of residuum (Jacobian for Newton-Raphson scheme)
dRLp_dLp2, & ! working copy of dRdLp dRLp_dLp2, & ! working copy of dRdLp
dRLi_dLi ! partial derivative of residuumI (Jacobian for Newton-Raphson scheme) dRLi_dLi ! partial derivative of residuumI (Jacobian for Newton-Raphson scheme)
@ -836,7 +836,8 @@ logical function integrateStress(ipc,ip,el,timeFraction)
steplengthLi, & steplengthLi, &
dt, & ! time increment dt, & ! time increment
aTolLp, & aTolLp, &
aTolLi aTolLi, &
devNull
integer NiterationStressLp, & ! number of stress integrations integer NiterationStressLp, & ! number of stress integrations
NiterationStressLi, & ! number of inner stress integrations NiterationStressLi, & ! number of inner stress integrations
ierr, & ! error indicator for LAPACK ierr, & ! error indicator for LAPACK
@ -844,6 +845,7 @@ logical function integrateStress(ipc,ip,el,timeFraction)
p, & p, &
jacoCounterLp, & jacoCounterLp, &
jacoCounterLi ! counters to check for Jacobian update jacoCounterLi ! counters to check for Jacobian update
logical :: error
external :: & external :: &
dgesv dgesv
@ -861,53 +863,38 @@ logical function integrateStress(ipc,ip,el,timeFraction)
Lpguess = crystallite_Lp(1:3,1:3,ipc,ip,el) ! take as first guess Lpguess = crystallite_Lp(1:3,1:3,ipc,ip,el) ! take as first guess
Liguess = crystallite_Li(1:3,1:3,ipc,ip,el) ! take as first guess Liguess = crystallite_Li(1:3,1:3,ipc,ip,el) ! take as first guess
Liguess_old = Liguess
invFp_current = math_inv33(crystallite_subFp0(1:3,1:3,ipc,ip,el)) call math_invert33(invFp_current,devNull,error,crystallite_subFp0(1:3,1:3,ipc,ip,el))
failedInversionFp: if (all(dEq0(invFp_current))) then if (error) return
return call math_invert33(invFi_current,devNull,error,crystallite_subFi0(1:3,1:3,ipc,ip,el))
endif failedInversionFp if (error) return
A = matmul(Fg_new,invFp_current) ! intermediate tensor needed later to calculate dFe_dLp
A = matmul(Fg_new,invFp_current) ! intermediate tensor needed later to calculate dFe_dLp
invFi_current = math_inv33(crystallite_subFi0(1:3,1:3,ipc,ip,el))
failedInversionFi: if (all(dEq0(invFi_current))) then
return
endif failedInversionFi
!* start Li loop with normal step length !* start Li loop with normal step length
NiterationStressLi = 0
jacoCounterLi = 0 jacoCounterLi = 0
steplengthLi = 1.0_pReal steplengthLi = 1.0_pReal
residuumLi_old = 0.0_pReal residuumLi_old = 0.0_pReal
Liguess_old = Liguess
LiLoop: do LiLoop: do NiterationStressLi=1,num%nStress
NiterationStressLi = NiterationStressLi + 1
LiLoopLimit: if (NiterationStressLi > num%nStress) then
return
endif LiLoopLimit
invFi_new = matmul(invFi_current,math_I3 - dt*Liguess) invFi_new = matmul(invFi_current,math_I3 - dt*Liguess)
Fi_new = math_inv33(invFi_new) Fi_new = math_inv33(invFi_new)
detInvFi = math_det33(invFi_new) detInvFi = math_det33(invFi_new)
!* start Lp loop with normal step length !* start Lp loop with normal step length
NiterationStressLp = 0
jacoCounterLp = 0 jacoCounterLp = 0
steplengthLp = 1.0_pReal steplengthLp = 1.0_pReal
residuumLp_old = 0.0_pReal residuumLp_old = 0.0_pReal
Lpguess_old = Lpguess Lpguess_old = Lpguess
LpLoop: do LpLoop: do NiterationStressLp=1,num%nStress
NiterationStressLp = NiterationStressLp + 1
LpLoopLimit: if (NiterationStressLp > num%nStress) then
return
endif LpLoopLimit
B = math_I3 - dt*Lpguess B = math_I3 - dt*Lpguess
Fe = matmul(matmul(A,B), invFi_new) Fe = matmul(matmul(A,B), invFi_new)
call constitutive_SandItsTangents(S, dS_dFe, dS_dFi, & call constitutive_SandItsTangents(S, dS_dFe, dS_dFi, &
Fe, Fi_new, ipc, ip, el) ! call constitutive law to calculate 2nd Piola-Kirchhoff stress and its derivative in unloaded configuration Fe, Fi_new, ipc, ip, el)
call constitutive_LpAndItsTangents(Lp_constitutive, dLp_dS, dLp_dFi, & call constitutive_LpAndItsTangents(Lp_constitutive, dLp_dS, dLp_dFi, &
S, Fi_new, ipc, ip, el) S, Fi_new, ipc, ip, el)
@ -921,8 +908,7 @@ logical function integrateStress(ipc,ip,el,timeFraction)
return ! ...me = .false. to inform integrator about problem return ! ...me = .false. to inform integrator about problem
elseif (norm2(residuumLp) < aTolLp) then ! converged if below absolute tolerance elseif (norm2(residuumLp) < aTolLp) then ! converged if below absolute tolerance
exit LpLoop ! ...leave iteration loop exit LpLoop ! ...leave iteration loop
elseif ( NiterationStressLp == 1 & elseif (NiterationStressLp == 1 .or. norm2(residuumLp) < norm2(residuumLp_old)) then ! not converged, but improved norm of residuum (always proceed in first iteration)...
.or. norm2(residuumLp) < norm2(residuumLp_old)) then ! not converged, but improved norm of residuum (always proceed in first iteration)...
residuumLp_old = residuumLp ! ...remember old values and... residuumLp_old = residuumLp ! ...remember old values and...
Lpguess_old = Lpguess Lpguess_old = Lpguess
steplengthLp = 1.0_pReal ! ...proceed with normal step length (calculate new search direction) steplengthLp = 1.0_pReal ! ...proceed with normal step length (calculate new search direction)
@ -942,10 +928,8 @@ logical function integrateStress(ipc,ip,el,timeFraction)
- math_3333to99(math_mul3333xx3333(math_mul3333xx3333(dLp_dS,dS_dFe),dFe_dLp)) - math_3333to99(math_mul3333xx3333(math_mul3333xx3333(dLp_dS,dS_dFe),dFe_dLp))
dRLp_dLp2 = dRLp_dLp ! will be overwritten in first call to LAPACK routine dRLp_dLp2 = dRLp_dLp ! will be overwritten in first call to LAPACK routine
work = math_33to9(residuumLp) work = math_33to9(residuumLp)
call dgesv(9,1,dRLp_dLp2,9,devNull,work,9,ierr) ! solve dRLp/dLp * delta Lp = -res for delta Lp call dgesv(9,1,dRLp_dLp2,9,devNull_9,work,9,ierr) ! solve dRLp/dLp * delta Lp = -res for delta Lp
if (ierr /= 0) then if (ierr /= 0) return
return
endif
deltaLp = - math_9to33(work) deltaLp = - math_9to33(work)
endif endif
jacoCounterLp = jacoCounterLp + 1 jacoCounterLp = jacoCounterLp + 1
@ -965,8 +949,7 @@ logical function integrateStress(ipc,ip,el,timeFraction)
return ! ...me = .false. to inform integrator about problem return ! ...me = .false. to inform integrator about problem
elseif (norm2(residuumLi) < aTolLi) then ! converged if below absolute tolerance elseif (norm2(residuumLi) < aTolLi) then ! converged if below absolute tolerance
exit LiLoop ! ...leave iteration loop exit LiLoop ! ...leave iteration loop
elseif ( NiterationStressLi == 1 & elseif (NiterationStressLi == 1 .or. norm2(residuumLi) < norm2(residuumLi_old)) then ! not converged, but improved norm of residuum (always proceed in first iteration)...
.or. norm2(residuumLi) < norm2(residuumLi_old)) then ! not converged, but improved norm of residuum (always proceed in first iteration)...
residuumLi_old = residuumLi ! ...remember old values and... residuumLi_old = residuumLi ! ...remember old values and...
Liguess_old = Liguess Liguess_old = Liguess
steplengthLi = 1.0_pReal ! ...proceed with normal step length (calculate new search direction) steplengthLi = 1.0_pReal ! ...proceed with normal step length (calculate new search direction)
@ -987,14 +970,12 @@ logical function integrateStress(ipc,ip,el,timeFraction)
dFi_dLi(1:3,1:3,o,p) = matmul(matmul(Fi_new,dFi_dLi(1:3,1:3,o,p)),Fi_new) dFi_dLi(1:3,1:3,o,p) = matmul(matmul(Fi_new,dFi_dLi(1:3,1:3,o,p)),Fi_new)
enddo; enddo enddo; enddo
dRLi_dLi = math_identity2nd(9) & dRLi_dLi = math_identity2nd(9) &
- math_3333to99(math_mul3333xx3333(dLi_dS, math_mul3333xx3333(dS_dFe, dFe_dLi) + & - math_3333to99(math_mul3333xx3333(dLi_dS, math_mul3333xx3333(dS_dFe, dFe_dLi) &
math_mul3333xx3333(dS_dFi, dFi_dLi))) & + math_mul3333xx3333(dS_dFi, dFi_dLi))) &
- math_3333to99(math_mul3333xx3333(dLi_dFi, dFi_dLi)) - math_3333to99(math_mul3333xx3333(dLi_dFi, dFi_dLi))
work = math_33to9(residuumLi) work = math_33to9(residuumLi)
call dgesv(9,1,dRLi_dLi,9,devNull,work,9,ierr) ! solve dRLi/dLp * delta Li = -res for delta Li call dgesv(9,1,dRLi_dLi,9,devNull_9,work,9,ierr) ! solve dRLi/dLp * delta Li = -res for delta Li
if (ierr /= 0) then if (ierr /= 0) return
return
endif
deltaLi = - math_9to33(work) deltaLi = - math_9to33(work)
endif endif
@ -1006,10 +987,8 @@ logical function integrateStress(ipc,ip,el,timeFraction)
!* calculate new plastic and elastic deformation gradient !* calculate new plastic and elastic deformation gradient
invFp_new = matmul(invFp_current,B) invFp_new = matmul(invFp_current,B)
invFp_new = invFp_new / math_det33(invFp_new)**(1.0_pReal/3.0_pReal) ! regularize invFp_new = invFp_new / math_det33(invFp_new)**(1.0_pReal/3.0_pReal) ! regularize
Fp_new = math_inv33(invFp_new) call math_invert33(Fp_new,devNull,error,invFp_new)
failedInversionInvFp: if (all(dEq0(Fp_new))) then if (error) return
return
endif failedInversionInvFp
Fe_new = matmul(matmul(Fg_new,invFp_new),invFi_new) Fe_new = matmul(matmul(Fg_new,invFp_new),invFi_new)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -1025,7 +1004,6 @@ logical function integrateStress(ipc,ip,el,timeFraction)
crystallite_invFp(1:3,1:3,ipc,ip,el) = invFp_new crystallite_invFp(1:3,1:3,ipc,ip,el) = invFp_new
crystallite_invFi(1:3,1:3,ipc,ip,el) = invFi_new crystallite_invFi(1:3,1:3,ipc,ip,el) = invFi_new
end function integrateStress end function integrateStress