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

cleaned up stress integrations and some documentation

This commit is contained in:
Pratheek Shanthraj 2014-11-13 12:53:20 +00:00
parent c84797c584
commit 23d68114e3
1 changed files with 112 additions and 125 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

237
code/crystallite.f90
View File

@ -3389,14 +3389,14 @@ logical function crystallite_integrateStress(&
Lpguess, & ! current guess for plastic velocity gradient Lpguess, & ! current guess for plastic velocity gradient
Lpguess_old, & ! known last good guess for plastic velocity gradient Lpguess_old, & ! known last good guess for plastic velocity gradient
Lp_constitutive, & ! plastic velocity gradient resulting from constitutive law Lp_constitutive, & ! plastic velocity gradient resulting from constitutive law
residuum, & ! current residuum of plastic velocity gradient residuumLp, & ! current residuum of plastic velocity gradient
residuum_old, & ! last residuum of plastic velocity gradient residuumLp_old, & ! last residuum of plastic velocity gradient
deltaLp, & ! direction of next guess deltaLp, & ! direction of next guess
Liguess, & ! current guess for plastic velocity gradient Liguess, & ! current guess for intermediate velocity gradient
Liguess_old, & ! known last good guess for plastic velocity gradient Liguess_old, & ! known last good guess for intermediate velocity gradient
Li_constitutive, & ! plastic velocity gradient resulting from constitutive law Li_constitutive, & ! intermediate velocity gradient resulting from constitutive law
residuumI, & ! current residuum of plastic velocity gradient residuumLi, & ! current residuum of intermediate velocity gradient
residuumI_old, & ! last residuum of plastic velocity gradient residuumLi_old, & ! last residuum of intermediate velocity gradient
deltaLi, & ! direction of next guess deltaLi, & ! direction of next guess
Tstar, & ! 2nd Piola-Kirchhoff Stress in plastic (lattice) configuration Tstar, & ! 2nd Piola-Kirchhoff Stress in plastic (lattice) configuration
Tstar_unloaded, & ! 2nd Piola-Kirchhoff Stress in unloaded configuration Tstar_unloaded, & ! 2nd Piola-Kirchhoff Stress in unloaded configuration
@ -3409,37 +3409,35 @@ logical function crystallite_integrateStress(&
real(pReal), dimension(6):: Tstar_v ! 2nd Piola-Kirchhoff Stress in Mandel-Notation real(pReal), dimension(6):: Tstar_v ! 2nd Piola-Kirchhoff Stress in Mandel-Notation
real(pReal), dimension(9):: work ! needed for matrix inversion by LAPACK real(pReal), dimension(9):: work ! needed for matrix inversion by LAPACK
integer(pInt), dimension(9) :: ipiv ! needed for matrix inversion by LAPACK integer(pInt), dimension(9) :: ipiv ! needed for matrix inversion by LAPACK
real(pReal), dimension(9,9) :: dLp_dT_constitutive, & ! partial derivative of plastic velocity gradient calculated by constitutive law real(pReal), dimension(9,9) :: dLp_dT_constitutive99, & ! partial derivative of plastic velocity gradient calculated by constitutive law
dLi_dT_constitutive, & ! partial derivative of intermediate velocity gradient calculated by constitutive law dLi_dT_constitutive99, & ! partial derivative of intermediate velocity gradient calculated by constitutive law
dT_dFe_constitutive, & ! partial derivative of 2nd Piola-Kirchhoff stress calculated by constitutive law dT_dFe99, & ! partial derivative of 2nd Piola-Kirchhoff stress calculated by constitutive law
dFe_dLp, & ! partial derivative of elastic deformation gradient dFe_dLp99, & ! partial derivative of elastic deformation gradient
dR_dLp, & ! partial derivative of residuum (Jacobian for NEwton-Raphson scheme) dFe_dLi99, &
dR_dLp2, & ! working copy of dRdLp dRLp_dLp, & ! partial derivative of residuum (Jacobian for NEwton-Raphson scheme)
dRI_dLp, & ! partial derivative of residuumI (Jacobian for NEwton-Raphson scheme) dRLp_dLp2, & ! working copy of dRdLp
dRI_dLp2 ! working copy of dRIdLp dRLi_dLi ! partial derivative of residuumI (Jacobian for NEwton-Raphson scheme)
real(pReal), dimension(3,3,3,3):: dT_dFe3333, & ! partial derivative of 2nd Piola-Kirchhoff stress real(pReal), dimension(3,3,3,3):: dT_dFe3333, & ! partial derivative of 2nd Piola-Kirchhoff stress
dT_dFe3333_unloaded, & dT_dFe3333_unloaded, &
dFe_dLp3333, & ! partial derivative of elastic deformation gradient dFe_dLp3333, & ! partial derivative of elastic deformation gradient
dInvFi_dLp3333, & dFe_dLi3333
dFe_dInvFi3333, &
dT_dInvFi3333, &
temp_3333
real(pReal) det, & ! determinant real(pReal) det, & ! determinant
detInvFi, & detInvFi, &
steplength0, & steplengthLp0, &
steplength, & steplengthLp, &
steplengthI0, & steplengthLi0, &
steplengthI, & steplengthLi, &
dt, & ! time increment dt, & ! time increment
aTol aTolLp, &
aTolLi
logical error ! flag indicating an error logical error ! flag indicating an error
integer(pInt) NiterationStress, & ! number of stress integrations integer(pInt) NiterationStressLp, & ! number of stress integrations
NiterationStressI, & ! number of inner stress integrations NiterationStressLi, & ! number of inner stress integrations
ierr, & ! error indicator for LAPACK ierr, & ! error indicator for LAPACK
o, & o, &
p, & p, &
jacoCounter, & jacoCounterLp, &
jacoICounter ! counters to check for Jacobian update jacoCounterLi ! counters to check for Jacobian update
integer(pLongInt) tick, & integer(pLongInt) tick, &
tock, & tock, &
tickrate, & tickrate, &
@ -3478,7 +3476,7 @@ logical function crystallite_integrateStress(&
!* feed local variables !* feed local variables
Fp_current = crystallite_subFp0(1:3,1:3,g,i,e) ! "Fp_current" is only used as temp var here... Fp_current = crystallite_subFp0(1:3,1:3,g,i,e) ! "Fp_current" is only used as temp var here...
Lpguess = crystallite_Lp(1:3,1:3,g,i,e) ! ... and take it as first guess Lpguess = crystallite_Lp (1:3,1:3,g,i,e) ! ... and take it as first guess
!* inversion of Fp_current... !* inversion of Fp_current...
@ -3497,9 +3495,15 @@ logical function crystallite_integrateStress(&
endif endif
A = math_mul33x33(Fg_new,invFp_current) ! intermediate tensor needed later to calculate dFe_dLp A = math_mul33x33(Fg_new,invFp_current) ! intermediate tensor needed later to calculate dFe_dLp
Liguess_old = 0.0_pReal !* feed local variables
Liguess = 0.0_pReal
Fi_current = constitutive_getFi0(g,i,e) ! intermediate configuration, assume decomposition as F = Fe Fi Fp Fi_current = constitutive_getFi0(g,i,e) ! intermediate configuration, assume decomposition as F = Fe Fi Fp
call constitutive_LiAndItsTangent(Liguess, dLi_dT_constitutive99, &
crystallite_Tstar_v(1:6,g,i,e), Lpguess, g, i, e)
Liguess_old = Liguess
!* inversion of Fi_current...
invFi_current = math_inv33(Fi_current) invFi_current = math_inv33(Fi_current)
if (all(invFi_current == 0.0_pReal)) then ! ... failed? if (all(invFi_current == 0.0_pReal)) then ! ... failed?
#ifndef _OPENMP #ifndef _OPENMP
@ -3515,15 +3519,15 @@ logical function crystallite_integrateStress(&
!* start LpLoop with normal step length !* start LpLoop with normal step length
NiterationStressI = 0_pInt NiterationStressLi = 0_pInt
jacoICounter = 0_pInt jacoCounterLi = 0_pInt
steplengthI0 = 1.0_pReal steplengthLi0 = 1.0_pReal
steplengthI = steplengthI0 steplengthLi = steplengthLi0
residuumI_old = 0.0_pReal residuumLi_old = 0.0_pReal
LiLoop: do LiLoop: do
NiterationStressI = NiterationStressI + 1_pInt NiterationStressLi = NiterationStressLi + 1_pInt
IloopsExeced: if (NiterationStressI > nStress) then IloopsExeced: if (NiterationStressLi > nStress) then
#ifndef _OPENMP #ifndef _OPENMP
if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) & if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) &
write(6,'(a,i3,a,i8,1x,a,i8,a,1x,i2,1x,i3,/)') '<< CRYST >> integrateStress reached inelastic loop limit',nStress, & write(6,'(a,i3,a,i8,1x,a,i8,a,1x,i2,1x,i3,/)') '<< CRYST >> integrateStress reached inelastic loop limit',nStress, &
@ -3532,20 +3536,20 @@ logical function crystallite_integrateStress(&
return return
endif IloopsExeced endif IloopsExeced
invFi = math_mul33x33(invFi_current,math_I3 - dt*Liguess) invFi = math_mul33x33(invFi_current,math_I3 - dt*Liguess)
detInvFi = math_det33(invFi) detInvFi = math_det33(invFi)
Fi = math_inv33(invFi) Fi = math_inv33(invFi)
NiterationStress = 0_pInt NiterationStressLp = 0_pInt
jacoCounter = 0_pInt jacoCounterLp = 0_pInt
steplength0 = 1.0_pReal steplengthLp0 = 1.0_pReal
steplength = steplength0 steplengthLp = steplengthLp0
residuum_old = 0.0_pReal residuumLp_old = 0.0_pReal
Lpguess_old = Lpguess Lpguess_old = Lpguess
LpLoop: do ! inner stress integration loop for consistency with Fi LpLoop: do ! inner stress integration loop for consistency with Fi
NiterationStress = NiterationStress + 1_pInt NiterationStressLp = NiterationStressLp + 1_pInt
loopsExeced: if (NiterationStress > nStress) then loopsExeced: if (NiterationStressLp > nStress) then
#ifndef _OPENMP #ifndef _OPENMP
if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) & if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) &
write(6,'(a,i3,a,i8,1x,a,i8,a,1x,i2,1x,i3,/)') '<< CRYST >> integrateStress reached loop limit',nStress, & write(6,'(a,i3,a,i8,1x,a,i8,a,1x,i2,1x,i3,/)') '<< CRYST >> integrateStress reached loop limit',nStress, &
@ -3574,7 +3578,7 @@ logical function crystallite_integrateStress(&
call system_clock(count=tick,count_rate=tickrate,count_max=maxticks) call system_clock(count=tick,count_rate=tickrate,count_max=maxticks)
endif endif
call constitutive_LpAndItsTangent(Lp_constitutive, dLp_dT_constitutive, Tstar_v, & call constitutive_LpAndItsTangent(Lp_constitutive, dLp_dT_constitutive99, Tstar_v, &
g, i, e) g, i, e)
if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) then if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) then
@ -3591,7 +3595,7 @@ logical function crystallite_integrateStress(&
if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt & if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt &
.and. ((e == debug_e .and. i == debug_i .and. g == debug_g) & .and. ((e == debug_e .and. i == debug_i .and. g == debug_g) &
.or. .not. iand(debug_level(debug_crystallite), debug_levelSelective) /= 0_pInt)) then .or. .not. iand(debug_level(debug_crystallite), debug_levelSelective) /= 0_pInt)) then
write(6,'(a,i3,/)') '<< CRYST >> iteration ', NiterationStress write(6,'(a,i3,/)') '<< CRYST >> iteration ', NiterationStressLp
write(6,'(a,/,3(12x,3(e20.7,1x)/))') '<< CRYST >> Lp_constitutive', math_transpose33(Lp_constitutive) write(6,'(a,/,3(12x,3(e20.7,1x)/))') '<< CRYST >> Lp_constitutive', math_transpose33(Lp_constitutive)
write(6,'(a,/,3(12x,3(e20.7,1x)/))') '<< CRYST >> Lpguess', math_transpose33(Lpguess) write(6,'(a,/,3(12x,3(e20.7,1x)/))') '<< CRYST >> Lpguess', math_transpose33(Lpguess)
endif endif
@ -3600,50 +3604,51 @@ logical function crystallite_integrateStress(&
!* update current residuum and check for convergence of loop !* update current residuum and check for convergence of loop
aTol = max(rTol_crystalliteStress * max(math_norm33(Lpguess),math_norm33(Lp_constitutive)), & ! absolute tolerance from largest acceptable relative error aTolLp = max(rTol_crystalliteStress * max(math_norm33(Lpguess),math_norm33(Lp_constitutive)), & ! absolute tolerance from largest acceptable relative error
aTol_crystalliteStress) ! minimum lower cutoff aTol_crystalliteStress) ! minimum lower cutoff
residuum = Lpguess - Lp_constitutive residuumLp = Lpguess - Lp_constitutive
if (any(residuum /= residuum)) then ! NaN in residuum... if (any(residuumLp /= residuumLp)) then ! NaN in residuum...
#ifndef _OPENMP #ifndef _OPENMP
if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) & if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) &
write(6,'(a,i8,1x,a,i8,a,1x,i2,1x,i3,a,i3,a)') '<< CRYST >> integrateStress encountered NaN at el (elFE) ip g ', & write(6,'(a,i8,1x,a,i8,a,1x,i2,1x,i3,a,i3,a)') '<< CRYST >> integrateStress encountered NaN at el (elFE) ip g ', &
e,mesh_element(1,e),i,g, & e,mesh_element(1,e),i,g, &
' ; iteration ', NiterationStress,& ' ; iteration ', NiterationStressLp,&
' >> returning..!' ' >> returning..!'
#endif #endif
return ! ...me = .false. to inform integrator about problem return ! ...me = .false. to inform integrator about problem
elseif (math_norm33(residuum) < aTol) then ! converged if below absolute tolerance elseif (math_norm33(residuumLp) < aTolLp) then ! converged if below absolute tolerance
exit LpLoop ! ...leave iteration loop exit LpLoop ! ...leave iteration loop
elseif (math_norm33(residuum) < math_norm33(residuum_old) .or. NiterationStress == 1_pInt ) then ! not converged, but improved norm of residuum (always proceed in first iteration)... elseif ( NiterationStressLp == 1_pInt &
residuum_old = residuum ! ...remember old values and... .or. math_norm33(residuumLp) < math_norm33(residuumLp_old)) then ! not converged, but improved norm of residuum (always proceed in first iteration)...
Lpguess_old = Lpguess residuumLp_old = residuumLp ! ...remember old values and...
steplength = steplength0 ! ...proceed with normal step length (calculate new search direction) Lpguess_old = Lpguess
steplengthLp = steplengthLp0 ! ...proceed with normal step length (calculate new search direction)
else ! not converged and residuum not improved... else ! not converged and residuum not improved...
steplength = 0.5_pReal * steplength ! ...try with smaller step length in same direction steplengthLp = 0.5_pReal * steplengthLp ! ...try with smaller step length in same direction
Lpguess = Lpguess_old + steplength * deltaLp Lpguess = Lpguess_old + steplengthLp * deltaLp
cycle LpLoop cycle LpLoop
endif endif
!* calculate Jacobian for correction term !* calculate Jacobian for correction term
if (mod(jacoCounter, iJacoLpresiduum) == 0_pInt) then if (mod(jacoCounterLp, iJacoLpresiduum) == 0_pInt) then
dFe_dLp3333 = 0.0_pReal dFe_dLp3333 = 0.0_pReal
do o=1_pInt,3_pInt; do p=1_pInt,3_pInt do o=1_pInt,3_pInt; do p=1_pInt,3_pInt
dFe_dLp3333(o,1:3,p,1:3) = A(o,p)*math_transpose33(invFi) ! dFe_dLp(i,j,k,l) = -dt * A(i,k) invFi(l,j) dFe_dLp3333(o,1:3,p,1:3) = A(o,p)*math_transpose33(invFi) ! dFe_dLp(i,j,k,l) = -dt * A(i,k) invFi(l,j)
enddo; enddo enddo; enddo
dFe_dLp3333 = -dt * dFe_dLp3333 dFe_dLp3333 = - dt * dFe_dLp3333
dFe_dLp = math_Plain3333to99(dFe_dLp3333) dFe_dLp99 = math_Plain3333to99(dFe_dLp3333)
dT_dFe_constitutive = math_Plain3333to99(dT_dFe3333) dT_dFe99 = math_Plain3333to99(dT_dFe3333)
dR_dLp = math_identity2nd(9_pInt) - & dRLp_dLp = math_identity2nd(9_pInt) &
math_mul99x99(dLp_dT_constitutive, math_mul99x99(dT_dFe_constitutive , dFe_dLp)) - math_mul99x99(dLp_dT_constitutive99, math_mul99x99(dT_dFe99 , dFe_dLp99))
dR_dLp2 = dR_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_plain33to9(residuum) work = math_plain33to9(residuumLp)
#if(FLOAT==8) #if(FLOAT==8)
call dgesv(9,1,dR_dLp2,9,ipiv,work,9,ierr) ! solve dR/dLp * delta Lp = -res for dR/dLp call dgesv(9,1,dRLp_dLp2,9,ipiv,work,9,ierr) ! solve dRLp/dLp * delta Lp = -res for delta Lp
#elif(FLOAT==4) #elif(FLOAT==4)
call sgesv(9,1,dR_dLp2,9,ipiv,work,9,ierr) ! solve dR/dLp * delta Lp = -res for dR/dLp call sgesv(9,1,dRLp_dLp2,9,ipiv,work,9,ierr) ! solve dRLp/dLp * delta Lp = -res for delta Lp
#endif #endif
if (ierr /= 0_pInt) then if (ierr /= 0_pInt) then
#ifndef _OPENMP #ifndef _OPENMP
@ -3654,10 +3659,10 @@ logical function crystallite_integrateStress(&
.and. ((e == debug_e .and. i == debug_i .and. g == debug_g)& .and. ((e == debug_e .and. i == debug_i .and. g == debug_g)&
.or. .not. iand(debug_level(debug_crystallite), debug_levelSelective) /= 0_pInt)) then .or. .not. iand(debug_level(debug_crystallite), debug_levelSelective) /= 0_pInt)) then
write(6,*) write(6,*)
write(6,'(a,/,9(12x,9(e15.3,1x)/))') '<< CRYST >> dR_dLp',transpose(dR_dLp) write(6,'(a,/,9(12x,9(e15.3,1x)/))') '<< CRYST >> dR_dLp',transpose(dRLp_dLp)
write(6,'(a,/,9(12x,9(e15.3,1x)/))') '<< CRYST >> dFe_dLp',transpose(dFe_dLp) write(6,'(a,/,9(12x,9(e15.3,1x)/))') '<< CRYST >> dFe_dLp',transpose(dFe_dLp99)
write(6,'(a,/,9(12x,9(e15.3,1x)/))') '<< CRYST >> dT_dFe_constitutive',transpose(dT_dFe_constitutive) write(6,'(a,/,9(12x,9(e15.3,1x)/))') '<< CRYST >> dT_dFe_constitutive',transpose(dT_dFe99)
write(6,'(a,/,9(12x,9(e15.3,1x)/))') '<< CRYST >> dLp_dT_constitutive',transpose(dLp_dT_constitutive) write(6,'(a,/,9(12x,9(e15.3,1x)/))') '<< CRYST >> dLp_dT_constitutive',transpose(dLp_dT_constitutive99)
write(6,'(a,/,3(12x,3(e20.7,1x)/))') '<< CRYST >> A',math_transpose33(A) write(6,'(a,/,3(12x,3(e20.7,1x)/))') '<< CRYST >> A',math_transpose33(A)
write(6,'(a,/,3(12x,3(e20.7,1x)/))') '<< CRYST >> B',math_transpose33(B) write(6,'(a,/,3(12x,3(e20.7,1x)/))') '<< CRYST >> B',math_transpose33(B)
write(6,'(a,/,3(12x,3(e20.7,1x)/))') '<< CRYST >> Lp_constitutive',math_transpose33(Lp_constitutive) write(6,'(a,/,3(12x,3(e20.7,1x)/))') '<< CRYST >> Lp_constitutive',math_transpose33(Lp_constitutive)
@ -3669,81 +3674,63 @@ logical function crystallite_integrateStress(&
endif endif
deltaLp = - math_plain9to33(work) deltaLp = - math_plain9to33(work)
endif endif
jacoCounter = jacoCounter + 1_pInt ! increase counter for jaco update jacoCounterLp = jacoCounterLp + 1_pInt ! increase counter for jaco update
Lpguess = Lpguess + steplength * deltaLp Lpguess = Lpguess + steplengthLp * deltaLp
enddo LpLoop enddo LpLoop
!* calculate intermediate velocity gradient and its tangent from constitutive law !* calculate intermediate velocity gradient and its tangent from constitutive law
call constitutive_LiAndItsTangent(Li_constitutive, dLi_dT_constitutive, Tstar_v, Lpguess, & call constitutive_LiAndItsTangent(Li_constitutive, dLi_dT_constitutive99, Tstar_v, Lpguess, &
g, i, e) g, i, e)
!* update current residuum and check for convergence of loop !* update current residuum and check for convergence of loop
aTol = max(rTol_crystalliteStress * max(math_norm33(Liguess),math_norm33(Li_constitutive)), & ! absolute tolerance from largest acceptable relative error aTolLi = max(rTol_crystalliteStress * max(math_norm33(Liguess),math_norm33(Li_constitutive)), & ! absolute tolerance from largest acceptable relative error
aTol_crystalliteStress) ! minimum lower cutoff aTol_crystalliteStress) ! minimum lower cutoff
residuumI = Liguess - Li_constitutive residuumLi = Liguess - Li_constitutive
if (any(residuumI /= residuumI)) then ! NaN in residuum... if (any(residuumLi /= residuumLi)) then ! NaN in residuum...
return ! ...me = .false. to inform integrator about problem return ! ...me = .false. to inform integrator about problem
elseif (math_norm33(residuumI) < aTol) then ! converged if below absolute tolerance elseif (math_norm33(residuumLi) < aTolLi) then ! converged if below absolute tolerance
exit LiLoop ! ...leave iteration loop exit LiLoop ! ...leave iteration loop
elseif (math_norm33(residuumI) < math_norm33(residuumI_old) .or. NiterationStressI == 1_pInt ) then! not converged, but improved norm of residuum (always proceed in first iteration)... elseif ( NiterationStressLi == 1_pInt &
residuumI_old = residuumI ! ...remember old values and... .or. math_norm33(residuumLi) < math_norm33(residuumLi_old)) then ! not converged, but improved norm of residuum (always proceed in first iteration)...
Liguess_old = Liguess residuumLi_old = residuumLi ! ...remember old values and...
steplengthI = steplengthI0 ! ...proceed with normal step length (calculate new search direction) Liguess_old = Liguess
steplengthLi = steplengthLi0 ! ...proceed with normal step length (calculate new search direction)
else ! not converged and residuum not improved... else ! not converged and residuum not improved...
steplengthI = 0.5_pReal * steplengthI ! ...try with smaller step length in same direction steplengthLi = 0.5_pReal * steplengthLi ! ...try with smaller step length in same direction
Liguess = Liguess_old + steplengthI * deltaLi Liguess = Liguess_old + steplengthLi * deltaLi
cycle LiLoop cycle LiLoop
endif endif
!* calculate Jacobian for correction term !* calculate Jacobian for correction term
if (mod(jacoICounter, iJacoLpresiduum) == 0_pInt) then if (mod(jacoCounterLi, iJacoLpresiduum) == 0_pInt) then
dInvFi_dLp3333 = 0.0_pReal temp_33 = math_mul33x33(math_mul33x33(A,B),invFi_current)
do o=1_pInt,3_pInt dFe_dLi3333 = 0.0_pReal
dInvFi_dLp3333(1:3,o,1:3,o) = invFi_current
enddo
dInvFi_dLp3333 = -dt * dInvFi_dLp3333
dFe_dInvFi3333 = 0.0_pReal
temp_33 = math_mul33x33(A,B)
do o=1_pInt,3_pInt
dFe_dInvFi3333(1:3,o,1:3,o) = temp_33
enddo
dT_dInvFi3333 = 0.0_pReal
do o=1_pInt,3_pInt; do p=1_pInt,3_pInt do o=1_pInt,3_pInt; do p=1_pInt,3_pInt
dT_dInvFi3333(1:3,1:3,p,o) = -Tstar * Fi(o,p) dFe_dLi3333(o,1:3,p,1:3) = temp_33 ! dFe_dLp(i,j,k,l) = -dt * A(i,k) invFi(l,j)
enddo; enddo enddo; enddo
temp_33 = math_mul33x33(invFi,Tstar_unloaded)/detInvFi dFe_dLi3333 = - dt * dFe_dLi3333
do o=1_pInt,3_pInt dFe_dLi99 = math_Plain3333to99(dFe_dLi3333)
dT_dInvFi3333(1:3,o,o,1:3) = dT_dInvFi3333(1:3,o,o,1:3) + temp_33 dRLi_dLi = math_identity2nd(9_pInt) &
dT_dInvFi3333(o,1:3,o,1:3) = dT_dInvFi3333(o,1:3,o,1:3) + temp_33 - math_mul99x99(dLi_dT_constitutive99, math_mul99x99(dT_dFe99, dFe_dLi99))
enddo work = math_plain33to9(residuumLi)
temp_3333 = math_mul3333xx3333(dT_dFe3333,dFe_dInvFi3333)
do o=1_pInt,3_pInt; do p=1_pInt,3_pInt
dT_dInvFi3333(1:3,1:3,o,p) = dT_dInvFi3333(1:3,1:3,o,p) + &
math_mul33x33(math_mul33x33(invFi,temp_3333(1:3,1:3,o,p)), &
math_transpose33(invFi))/detInvFi
enddo; enddo
dRI_dLp = math_identity2nd(9_pInt) - &
math_mul99x99(dLi_dT_constitutive,math_Plain3333to99(math_mul3333xx3333(dT_dInvFi3333,dInvFi_dLp3333)))
dRI_dLp2 = dRI_dLp ! will be overwritten in first call to LAPACK routine
work = math_plain33to9(residuumI)
#if(FLOAT==8) #if(FLOAT==8)
call dgesv(9,1,dRI_dLp2,9,ipiv,work,9,ierr) ! solve dR/dLp * delta Lp = -res for dR/dLp call dgesv(9,1,dRLi_dLi,9,ipiv,work,9,ierr) ! solve dRLi/dLp * delta Li = -res for delta Li
#elif(FLOAT==4) #elif(FLOAT==4)
call sgesv(9,1,dRI_dLp2,9,ipiv,work,9,ierr) ! solve dR/dLp * delta Lp = -res for dR/dLp call sgesv(9,1,dRLi_dLi,9,ipiv,work,9,ierr) ! solve dRLi/dLp * delta Li = -res for delta Li
#endif #endif
if (ierr /= 0_pInt) then if (ierr /= 0_pInt) then
return return
endif endif
deltaLi = - math_plain9to33(work) deltaLi = - math_plain9to33(work)
endif endif
jacoICounter = jacoICounter + 1_pInt ! increase counter for jaco update jacoCounterLi = jacoCounterLi + 1_pInt ! increase counter for jaco update
Liguess = Liguess + steplengthI * deltaLi Liguess = Liguess + steplengthLi * deltaLi
enddo LiLoop enddo LiLoop
!* calculate new plastic and elastic deformation gradient !* calculate new plastic and elastic deformation gradient
@ -3755,7 +3742,7 @@ logical function crystallite_integrateStress(&
#ifndef _OPENMP #ifndef _OPENMP
if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) then if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) then
write(6,'(a,i8,1x,a,i8,a,1x,i2,1x,i3,a,i3)') '<< CRYST >> integrateStress failed on invFp_new inversion at el ip g ',& write(6,'(a,i8,1x,a,i8,a,1x,i2,1x,i3,a,i3)') '<< CRYST >> integrateStress failed on invFp_new inversion at el ip g ',&
e,mesh_element(1,e),i,g, ' ; iteration ', NiterationStress e,mesh_element(1,e),i,g, ' ; iteration ', NiterationStressLp
if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt & if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt &
.and. ((e == debug_e .and. i == debug_i .and. g == debug_g) & .and. ((e == debug_e .and. i == debug_i .and. g == debug_g) &
.or. .not. iand(debug_level(debug_crystallite), debug_levelSelective) /= 0_pInt)) & .or. .not. iand(debug_level(debug_crystallite), debug_levelSelective) /= 0_pInt)) &
@ -3799,8 +3786,8 @@ logical function crystallite_integrateStress(&
if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) then if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) then
!$OMP CRITICAL (distributionStress) !$OMP CRITICAL (distributionStress)
debug_StressLoopDistribution(NiterationStress,numerics_integrationMode) = & debug_StressLoopDistribution(NiterationStressLp,numerics_integrationMode) = &
debug_StressLoopDistribution(NiterationStress,numerics_integrationMode) + 1_pInt debug_StressLoopDistribution(NiterationStressLp,numerics_integrationMode) + 1_pInt
!$OMP END CRITICAL (distributionStress) !$OMP END CRITICAL (distributionStress)
endif endif