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

switched order of Lp-Li nesting in stress integration loop for better convergence. temperature is integrated during stress integration so it need not be a dot state

This commit is contained in:
Pratheek Shanthraj 2014-11-12 16:40:50 +00:00
parent e89199d119
commit fc9f4835c3
3 changed files with 236 additions and 223 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

25
code/constitutive.f90
View File

@ -796,15 +796,19 @@ end function constitutive_getFi
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief contains the constitutive equation for calculating the intermediate deformation gradient !> @brief contains the constitutive equation for calculating the intermediate deformation gradient
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine constitutive_putFi(Tstar_v, dt, ipc, ip, el) subroutine constitutive_putFi(Tstar_v, Lp, dt, ipc, ip, el)
use prec, only: & use prec, only: &
pReal pReal
use material, only: & use material, only: &
phase_damage, & phase_damage, &
phase_thermal, &
material_phase, & material_phase, &
LOCAL_DAMAGE_anisoBrittle_ID LOCAL_DAMAGE_anisoBrittle_ID, &
LOCAL_THERMAL_adiabatic_ID
use damage_anisoBrittle, only: & use damage_anisoBrittle, only: &
damage_anisoBrittle_putFd damage_anisoBrittle_putFd
use thermal_adiabatic, only: &
thermal_adiabatic_putFT
implicit none implicit none
integer(pInt), intent(in) :: & integer(pInt), intent(in) :: &
@ -813,6 +817,8 @@ subroutine constitutive_putFi(Tstar_v, dt, ipc, ip, el)
el !< element number el !< element number
real(pReal), intent(in), dimension(6) :: & real(pReal), intent(in), dimension(6) :: &
Tstar_v !< 2nd Piola-Kirchhoff stress Tstar_v !< 2nd Piola-Kirchhoff stress
real(pReal), intent(in), dimension(3,3) :: &
Lp !< plastic velocity gradient
real(pReal), intent(in) :: & real(pReal), intent(in) :: &
dt dt
@ -821,6 +827,12 @@ subroutine constitutive_putFi(Tstar_v, dt, ipc, ip, el)
call damage_anisoBrittle_putFd (Tstar_v, dt, ipc, ip, el) call damage_anisoBrittle_putFd (Tstar_v, dt, ipc, ip, el)
end select end select
select case (phase_thermal(material_phase(ipc,ip,el)))
case (LOCAL_THERMAL_adiabatic_ID)
call thermal_adiabatic_putFT (Tstar_v, Lp, dt, ipc, ip, el)
end select
end subroutine constitutive_putFi end subroutine constitutive_putFi
@ -1010,7 +1022,6 @@ subroutine constitutive_collectDotState(Tstar_v, Lp, FeArray, FpArray, subdt, su
use material, only: & use material, only: &
phase_plasticity, & phase_plasticity, &
phase_damage, & phase_damage, &
phase_thermal, &
phase_vacancy, & phase_vacancy, &
material_phase, & material_phase, &
homogenization_maxNgrains, & homogenization_maxNgrains, &
@ -1026,7 +1037,6 @@ subroutine constitutive_collectDotState(Tstar_v, Lp, FeArray, FpArray, subdt, su
LOCAL_DAMAGE_anisoDuctile_ID, & LOCAL_DAMAGE_anisoDuctile_ID, &
LOCAL_DAMAGE_anisoBrittle_ID, & LOCAL_DAMAGE_anisoBrittle_ID, &
LOCAL_DAMAGE_gurson_ID, & LOCAL_DAMAGE_gurson_ID, &
LOCAL_THERMAL_adiabatic_ID, &
LOCAL_VACANCY_generation_ID LOCAL_VACANCY_generation_ID
use constitutive_j2, only: & use constitutive_j2, only: &
constitutive_j2_dotState constitutive_j2_dotState
@ -1050,8 +1060,6 @@ subroutine constitutive_collectDotState(Tstar_v, Lp, FeArray, FpArray, subdt, su
damage_anisoDuctile_dotState damage_anisoDuctile_dotState
use damage_gurson, only: & use damage_gurson, only: &
damage_gurson_dotState damage_gurson_dotState
use thermal_adiabatic, only: &
thermal_adiabatic_dotState
use vacancy_generation, only: & use vacancy_generation, only: &
vacancy_generation_dotState vacancy_generation_dotState
@ -1115,11 +1123,6 @@ subroutine constitutive_collectDotState(Tstar_v, Lp, FeArray, FpArray, subdt, su
call damage_gurson_dotState(Tstar_v, Lp, ipc, ip, el) call damage_gurson_dotState(Tstar_v, Lp, ipc, ip, el)
end select end select
select case (phase_thermal(material_phase(ipc,ip,el)))
case (LOCAL_THERMAL_adiabatic_ID)
call thermal_adiabatic_dotState(Tstar_v, Lp, ipc, ip, el)
end select
select case (phase_vacancy(material_phase(ipc,ip,el))) select case (phase_vacancy(material_phase(ipc,ip,el)))
case (LOCAL_VACANCY_generation_ID) case (LOCAL_VACANCY_generation_ID)
call constitutive_getAccumulatedSlip(nSlip,accumulatedSlip,ipc,ip,el) call constitutive_getAccumulatedSlip(nSlip,accumulatedSlip,ipc,ip,el)

350
code/crystallite.f90
View File

@ -3418,6 +3418,7 @@ logical function crystallite_integrateStress(&
dRI_dLp, & ! partial derivative of residuumI (Jacobian for NEwton-Raphson scheme) dRI_dLp, & ! partial derivative of residuumI (Jacobian for NEwton-Raphson scheme)
dRI_dLp2 ! working copy of dRIdLp dRI_dLp2 ! working copy of dRIdLp
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, &
dFe_dLp3333, & ! partial derivative of elastic deformation gradient dFe_dLp3333, & ! partial derivative of elastic deformation gradient
dInvFi_dLp3333, & dInvFi_dLp3333, &
dFe_dInvFi3333, & dFe_dInvFi3333, &
@ -3477,7 +3478,6 @@ 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_old = crystallite_Lp(1:3,1:3,g,i,e) ! consider present Lp good (i.e. worth remembering) ...
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
@ -3497,6 +3497,7 @@ 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
Liguess = 0.0_pReal 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
invFi_current = math_inv33(Fi_current) invFi_current = math_inv33(Fi_current)
@ -3514,229 +3515,236 @@ logical function crystallite_integrateStress(&
!* start LpLoop with normal step length !* start LpLoop with normal step length
NiterationStress = 0_pInt NiterationStressI = 0_pInt
jacoCounter = 0_pInt jacoICounter = 0_pInt
steplength0 = 1.0_pReal steplengthI0 = 1.0_pReal
steplength = steplength0 steplengthI = steplengthI0
residuum_old = 0.0_pReal residuumI_old = 0.0_pReal
LpLoop: do LiLoop: do
NiterationStress = NiterationStress + 1_pInt NiterationStressI = NiterationStressI + 1_pInt
loopsExeced: if (NiterationStress > nStress) then IloopsExeced: if (NiterationStressI > 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 inelastic loop limit',nStress, &
' at el (elFE) ip g ', e,mesh_element(1,e),i,g ' at el (elFE) ip g ', e,mesh_element(1,e),i,g
#endif #endif
return return
endif loopsExeced endif IloopsExeced
B = math_I3 - dt*Lpguess
NiterationStressI = 0_pInt invFi = math_mul33x33(invFi_current,math_I3 - dt*Liguess)
Liguess_old = Liguess detInvFi = math_det33(invFi)
jacoICounter = 0_pInt Fi = math_inv33(invFi)
steplengthI0 = 1.0_pReal
steplengthI = steplengthI0 NiterationStress = 0_pInt
residuumI_old = 0.0_pReal jacoCounter = 0_pInt
LiLoop: do ! inner stress integration loop for consistency with Fi steplength0 = 1.0_pReal
NiterationStressI = NiterationStressI + 1_pInt steplength = steplength0
IloopsExeced: if (NiterationStressI > nStress) then residuum_old = 0.0_pReal
Lpguess_old = Lpguess
LpLoop: do ! inner stress integration loop for consistency with Fi
NiterationStress = NiterationStress + 1_pInt
loopsExeced: if (NiterationStress > 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 loop limit',nStress, &
' at el (elFE) ip g ', e,mesh_element(1,e),i,g ' at el (elFE) ip g ', e,mesh_element(1,e),i,g
#endif #endif
return return
endif IloopsExeced endif loopsExeced
!* calculate (elastic) 2nd Piola--Kirchhoff stress tensor and its tangent from constitutive law !* calculate (elastic) 2nd Piola--Kirchhoff stress tensor and its tangent from constitutive law
invFi = math_mul33x33(invFi_current,math_I3 - dt*Liguess) B = math_I3 - dt*Lpguess
detInvFi = math_det33(invFi)
Fi = math_inv33(invFi)
Fe = math_mul33x33(math_mul33x33(A,B),invFi) ! current elastic deformation tensor Fe = math_mul33x33(math_mul33x33(A,B),invFi) ! current elastic deformation tensor
call constitutive_TandItsTangent(Tstar_unloaded, dT_dFe3333, Fe, g,i,e) ! call constitutive law to calculate 2nd Piola-Kirchhoff stress and its derivative in unloaded configuration call constitutive_TandItsTangent(Tstar_unloaded, dT_dFe3333_unloaded, Fe, g,i,e) ! call constitutive law to calculate 2nd Piola-Kirchhoff stress and its derivative in unloaded configuration
Tstar = math_mul33x33(invFi, & Tstar = math_mul33x33(invFi, &
math_mul33x33(Tstar_unloaded,math_transpose33(invFi)))/detInvFi ! push Tstar forward from unloaded to plastic (lattice) configuration math_mul33x33(Tstar_unloaded,math_transpose33(invFi)))/detInvFi ! push Tstar forward from unloaded to plastic (lattice) configuration
do o=1_pInt,3_pInt; do p=1_pInt,3_pInt
dT_dFe3333(1:3,1:3,o,p) = math_mul33x33(invFi, &
math_mul33x33(dT_dFe3333_unloaded(1:3,1:3,o,p),math_transpose33(invFi)))/detInvFi
enddo; enddo
Tstar_v = math_Mandel33to6(Tstar) Tstar_v = math_Mandel33to6(Tstar)
!* calculate intermediate velocity gradient and its tangent from constitutive law
call constitutive_LiAndItsTangent(Li_constitutive, dLi_dT_constitutive, Tstar_v, Lpguess, & !* calculate plastic velocity gradient and its tangent from constitutive law
if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) then
call system_clock(count=tick,count_rate=tickrate,count_max=maxticks)
endif
call constitutive_LpAndItsTangent(Lp_constitutive, dLp_dT_constitutive, Tstar_v, &
g, i, e) g, i, e)
if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) then
call system_clock(count=tock,count_rate=tickrate,count_max=maxticks)
!$OMP CRITICAL (debugTimingLpTangent)
debug_cumLpCalls = debug_cumLpCalls + 1_pInt
debug_cumLpTicks = debug_cumLpTicks + tock-tick
!$OMP FLUSH (debug_cumLpTicks)
if (tock < tick) debug_cumLpTicks = debug_cumLpTicks + maxticks
!$OMP END CRITICAL (debugTimingLpTangent)
endif
#ifndef _OPENMP
if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt &
.and. ((e == debug_e .and. i == debug_i .and. g == debug_g) &
.or. .not. iand(debug_level(debug_crystallite), debug_levelSelective) /= 0_pInt)) then
write(6,'(a,i3,/)') '<< CRYST >> iteration ', NiterationStress
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)
endif
#endif
!* 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 aTol = 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
residuumI = Liguess - Li_constitutive residuum = Lpguess - Lp_constitutive
if (any(residuumI /= residuumI)) then ! NaN in residuum...
return ! ...me = .false. to inform integrator about problem if (any(residuum /= residuum)) then ! NaN in residuum...
elseif (math_norm33(residuumI) < aTol) then ! converged if below absolute tolerance #ifndef _OPENMP
exit LiLoop ! ...leave iteration loop if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) &
elseif (math_norm33(residuumI) < math_norm33(residuumI_old) .or. NiterationStress == 1_pInt ) then ! not converged, but improved norm of residuum (always proceed in first iteration)... write(6,'(a,i8,1x,a,i8,a,1x,i2,1x,i3,a,i3,a)') '<< CRYST >> integrateStress encountered NaN at el (elFE) ip g ', &
residuumI_old = residuumI ! ...remember old values and... e,mesh_element(1,e),i,g, &
Liguess_old = Liguess ' ; iteration ', NiterationStress,&
steplengthI = steplengthI0 ! ...proceed with normal step length (calculate new search direction) ' >> returning..!'
else ! not converged and residuum not improved... #endif
steplengthI = 0.5_pReal * steplengthI ! ...try with smaller step length in same direction return ! ...me = .false. to inform integrator about problem
Liguess = Liguess_old + steplengthI * deltaLi elseif (math_norm33(residuum) < aTol) then ! converged if below absolute tolerance
cycle LiLoop 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)...
residuum_old = residuum ! ...remember old values and...
Lpguess_old = Lpguess
steplength = steplength0 ! ...proceed with normal step length (calculate new search direction)
else ! not converged and residuum not improved...
steplength = 0.5_pReal * steplength ! ...try with smaller step length in same direction
Lpguess = Lpguess_old + steplength * deltaLp
cycle LpLoop
endif endif
!* calculate Jacobian for correction term !* calculate Jacobian for correction term
if (mod(jacoICounter, iJacoLpresiduum) == 0_pInt) then if (mod(jacoCounter, iJacoLpresiduum) == 0_pInt) then
dInvFi_dLp3333 = 0.0_pReal dFe_dLp3333 = 0.0_pReal
do o=1_pInt,3_pInt
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_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
temp_33 = math_mul33x33(invFi,Tstar_unloaded)/detInvFi dFe_dLp3333 = -dt * dFe_dLp3333
do o=1_pInt,3_pInt dFe_dLp = math_Plain3333to99(dFe_dLp3333)
dT_dInvFi3333(1:3,o,o,1:3) = dT_dInvFi3333(1:3,o,o,1:3) + temp_33 dT_dFe_constitutive = math_Plain3333to99(dT_dFe3333)
dT_dInvFi3333(o,1:3,o,1:3) = dT_dInvFi3333(o,1:3,o,1:3) + temp_33 dR_dLp = math_identity2nd(9_pInt) - &
enddo math_mul99x99(dLp_dT_constitutive, math_mul99x99(dT_dFe_constitutive , dFe_dLp))
temp_3333 = math_mul3333xx3333(dT_dFe3333,dFe_dInvFi3333) dR_dLp2 = dR_dLp ! will be overwritten in first call to LAPACK routine
do o=1_pInt,3_pInt; do p=1_pInt,3_pInt work = math_plain33to9(residuum)
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,dR_dLp2,9,ipiv,work,9,ierr) ! solve dR/dLp * delta Lp = -res for dR/dLp
#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,dR_dLp2,9,ipiv,work,9,ierr) ! solve dR/dLp * delta Lp = -res for dR/dLp
#endif #endif
if (ierr /= 0_pInt) then if (ierr /= 0_pInt) then
#ifndef _OPENMP
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 dR/dLp inversion at el ip g ', &
e,mesh_element(1,e),i,g
if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt &
.and. ((e == debug_e .and. i == debug_i .and. g == debug_g)&
.or. .not. iand(debug_level(debug_crystallite), debug_levelSelective) /= 0_pInt)) then
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 >> dFe_dLp',transpose(dFe_dLp)
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 >> dLp_dT_constitutive',transpose(dLp_dT_constitutive)
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 >> Lp_constitutive',math_transpose33(Lp_constitutive)
write(6,'(a,/,3(12x,3(e20.7,1x)/))') '<< CRYST >> Lpguess',math_transpose33(Lpguess)
endif
endif
#endif
return return
endif endif
deltaLi = - math_plain9to33(work) deltaLp = - math_plain9to33(work)
endif endif
jacoICounter = jacoICounter + 1_pInt ! increase counter for jaco update jacoCounter = jacoCounter + 1_pInt ! increase counter for jaco update
Liguess = Liguess + steplengthI * deltaLi Lpguess = Lpguess + steplength * deltaLp
enddo LiLoop
!* calculate plastic velocity gradient and its tangent from constitutive law enddo LpLoop
if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) then !* calculate intermediate velocity gradient and its tangent from constitutive law
call system_clock(count=tick,count_rate=tickrate,count_max=maxticks)
endif
call constitutive_LpAndItsTangent(Lp_constitutive, dLp_dT_constitutive, Tstar_v, & call constitutive_LiAndItsTangent(Li_constitutive, dLi_dT_constitutive, Tstar_v, Lpguess, &
g, i, e) g, i, e)
if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) then
call system_clock(count=tock,count_rate=tickrate,count_max=maxticks)
!$OMP CRITICAL (debugTimingLpTangent)
debug_cumLpCalls = debug_cumLpCalls + 1_pInt
debug_cumLpTicks = debug_cumLpTicks + tock-tick
!$OMP FLUSH (debug_cumLpTicks)
if (tock < tick) debug_cumLpTicks = debug_cumLpTicks + maxticks
!$OMP END CRITICAL (debugTimingLpTangent)
endif
#ifndef _OPENMP
if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt &
.and. ((e == debug_e .and. i == debug_i .and. g == debug_g) &
.or. .not. iand(debug_level(debug_crystallite), debug_levelSelective) /= 0_pInt)) then
write(6,'(a,i3,/)') '<< CRYST >> iteration ', NiterationStress
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)
endif
#endif
!* 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 aTol = 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
residuum = Lpguess - Lp_constitutive residuumI = Liguess - Li_constitutive
if (any(residuumI /= residuumI)) then ! NaN in residuum...
if (any(residuum /= residuum)) then ! NaN in residuum... return ! ...me = .false. to inform integrator about problem
#ifndef _OPENMP elseif (math_norm33(residuumI) < aTol) then ! converged if below absolute tolerance
if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) & exit LiLoop ! ...leave iteration loop
write(6,'(a,i8,1x,a,i8,a,1x,i2,1x,i3,a,i3,a)') '<< CRYST >> integrateStress encountered NaN at el (elFE) ip g ', & 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)...
e,mesh_element(1,e),i,g, & residuumI_old = residuumI ! ...remember old values and...
' ; iteration ', NiterationStress,& Liguess_old = Liguess
' >> returning..!' steplengthI = steplengthI0 ! ...proceed with normal step length (calculate new search direction)
#endif else ! not converged and residuum not improved...
return ! ...me = .false. to inform integrator about problem steplengthI = 0.5_pReal * steplengthI ! ...try with smaller step length in same direction
elseif (math_norm33(residuum) < aTol) then ! converged if below absolute tolerance Liguess = Liguess_old + steplengthI * deltaLi
exit LpLoop ! ...leave iteration loop cycle LiLoop
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)...
residuum_old = residuum ! ...remember old values and...
Lpguess_old = Lpguess
steplength = steplength0 ! ...proceed with normal step length (calculate new search direction)
else ! not converged and residuum not improved...
steplength = 0.5_pReal * steplength ! ...try with smaller step length in same direction
Lpguess = Lpguess_old + steplength * deltaLp
cycle LpLoop
endif endif
!* calculate Jacobian for correction term !* calculate Jacobian for correction term
if (mod(jacoCounter, iJacoLpresiduum) == 0_pInt) then if (mod(jacoICounter, iJacoLpresiduum) == 0_pInt) then
dFe_dLp3333 = 0.0_pReal dInvFi_dLp3333 = 0.0_pReal
do o=1_pInt,3_pInt
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
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) dT_dInvFi3333(1:3,1:3,p,o) = -Tstar * Fi(o,p)
enddo; enddo enddo; enddo
dFe_dLp3333 = -dt * dFe_dLp3333 temp_33 = math_mul33x33(invFi,Tstar_unloaded)/detInvFi
dFe_dLp = math_Plain3333to99(dFe_dLp3333) do o=1_pInt,3_pInt
dT_dFe_constitutive = math_Plain3333to99(dT_dFe3333) dT_dInvFi3333(1:3,o,o,1:3) = dT_dInvFi3333(1:3,o,o,1:3) + temp_33
dR_dLp = math_identity2nd(9_pInt) - & dT_dInvFi3333(o,1:3,o,1:3) = dT_dInvFi3333(o,1:3,o,1:3) + temp_33
math_mul99x99(dLp_dT_constitutive, math_mul99x99(dT_dFe_constitutive , dFe_dLp)) enddo
dR_dLp2 = dR_dLp ! will be overwritten in first call to LAPACK routine temp_3333 = math_mul3333xx3333(dT_dFe3333,dFe_dInvFi3333)
work = math_plain33to9(residuum) 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,dR_dLp2,9,ipiv,work,9,ierr) ! solve dR/dLp * delta Lp = -res for dR/dLp call dgesv(9,1,dRI_dLp2,9,ipiv,work,9,ierr) ! solve dR/dLp * delta Lp = -res for dR/dLp
#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,dRI_dLp2,9,ipiv,work,9,ierr) ! solve dR/dLp * delta Lp = -res for dR/dLp
#endif #endif
if (ierr /= 0_pInt) then if (ierr /= 0_pInt) then
#ifndef _OPENMP
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 dR/dLp inversion at el ip g ', &
e,mesh_element(1,e),i,g
if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt &
.and. ((e == debug_e .and. i == debug_i .and. g == debug_g)&
.or. .not. iand(debug_level(debug_crystallite), debug_levelSelective) /= 0_pInt)) then
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 >> dFe_dLp',transpose(dFe_dLp)
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 >> dLp_dT_constitutive',transpose(dLp_dT_constitutive)
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 >> Lp_constitutive',math_transpose33(Lp_constitutive)
write(6,'(a,/,3(12x,3(e20.7,1x)/))') '<< CRYST >> Lpguess',math_transpose33(Lpguess)
endif
endif
#endif
return return
endif endif
deltaLp = - math_plain9to33(work) deltaLi = - math_plain9to33(work)
endif endif
jacoCounter = jacoCounter + 1_pInt ! increase counter for jaco update jacoICounter = jacoICounter + 1_pInt ! increase counter for jaco update
Lpguess = Lpguess + steplength * deltaLp
enddo LpLoop
Liguess = Liguess + steplengthI * deltaLi
enddo LiLoop
!* calculate new plastic and elastic deformation gradient !* calculate new plastic and elastic deformation gradient
@ -3758,22 +3766,20 @@ logical function crystallite_integrateStress(&
endif endif
Fe_new = math_mul33x33(math_mul33x33(Fg_new,invFp_new),invFi) ! calc resulting Fe Fe_new = math_mul33x33(math_mul33x33(Fg_new,invFp_new),invFi) ! calc resulting Fe
!* calculate 1st Piola-Kirchhoff stress !* calculate 1st Piola-Kirchhoff stress
crystallite_P(1:3,1:3,g,i,e) = math_mul33x33(math_mul33x33(Fe_new,Fi), & crystallite_P(1:3,1:3,g,i,e) = math_mul33x33(math_mul33x33(Fe_new,Fi), &
math_mul33x33(math_Mandel6to33(Tstar_v), & math_mul33x33(math_Mandel6to33(Tstar_v), &
math_transpose33(invFp_new)))/detInvFi math_transpose33(invFp_new)))/detInvFi
!* store local values in global variables !* store local values in global variables
crystallite_Lp(1:3,1:3,g,i,e) = Lpguess crystallite_Lp(1:3,1:3,g,i,e) = Lpguess
crystallite_Tstar_v(1:6,g,i,e) = Tstar_v crystallite_Tstar_v(1:6,g,i,e) = Tstar_v
crystallite_Fp(1:3,1:3,g,i,e) = Fp_new crystallite_Fp(1:3,1:3,g,i,e) = Fp_new
crystallite_Fe(1:3,1:3,g,i,e) = Fe_new crystallite_Fe(1:3,1:3,g,i,e) = Fe_new
crystallite_invFp(1:3,1:3,g,i,e) = invFp_new crystallite_invFp(1:3,1:3,g,i,e) = invFp_new
call constitutive_putFi(Tstar_v, dt, g, i, e) call constitutive_putFi(Tstar_v, Lpguess, dt, g, i, e)
!* set return flag to true !* set return flag to true

84
code/thermal_adiabatic.f90
View File

@ -39,9 +39,9 @@ module thermal_adiabatic
thermal_adiabatic_init, & thermal_adiabatic_init, &
thermal_adiabatic_stateInit, & thermal_adiabatic_stateInit, &
thermal_adiabatic_aTolState, & thermal_adiabatic_aTolState, &
thermal_adiabatic_dotState, &
thermal_adiabatic_LTAndItsTangent, & thermal_adiabatic_LTAndItsTangent, &
thermal_adiabatic_getFT, & thermal_adiabatic_getFT, &
thermal_adiabatic_putFT, &
thermal_adiabatic_getFT0, & thermal_adiabatic_getFT0, &
thermal_adiabatic_getPartionedFT0, & thermal_adiabatic_getPartionedFT0, &
thermal_adiabatic_getTemperature, & thermal_adiabatic_getTemperature, &
@ -183,7 +183,7 @@ subroutine thermal_adiabatic_init(fileUnit,temperature_init)
endif endif
enddo outputsLoop enddo outputsLoop
! Determine size of state array ! Determine size of state array
sizeDotState = 1_pInt sizeDotState = 0_pInt
sizeState = 1_pInt sizeState = 1_pInt
thermalState(phase)%sizeState = sizeState thermalState(phase)%sizeState = sizeState
thermalState(phase)%sizeDotState = sizeDotState thermalState(phase)%sizeDotState = sizeDotState
@ -250,43 +250,6 @@ subroutine thermal_adiabatic_aTolState(phase,instance)
thermalState(phase)%aTolState = tempTol thermalState(phase)%aTolState = tempTol
end subroutine thermal_adiabatic_aTolState end subroutine thermal_adiabatic_aTolState
!--------------------------------------------------------------------------------------------------
!> @brief calculates derived quantities from state
!--------------------------------------------------------------------------------------------------
subroutine thermal_adiabatic_dotState(Tstar_v, Lp, ipc, ip, el)
use lattice, only: &
lattice_massDensity, &
lattice_specificHeat
use material, only: &
mappingConstitutive, &
phase_thermalInstance, &
thermalState
use math, only: &
math_Mandel6to33
implicit none
integer(pInt), intent(in) :: &
ipc, & !< component-ID of integration point
ip, & !< integration point
el !< element
real(pReal), intent(in), dimension(6) :: &
Tstar_v !< 2nd Piola Kirchhoff stress tensor (Mandel)
real(pReal), intent(in), dimension(3,3) :: &
Lp
integer(pInt) :: &
instance, phase, constituent
phase = mappingConstitutive(2,ipc,ip,el)
constituent = mappingConstitutive(1,ipc,ip,el)
instance = phase_thermalInstance(phase)
thermalState(phase)%dotState(1,constituent) = &
0.95_pReal &
* sum(abs(math_Mandel6to33(Tstar_v)*Lp)) &
/ (lattice_massDensity(phase)*lattice_specificHeat(phase))
end subroutine thermal_adiabatic_dotState
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief contains the constitutive equation for calculating the velocity gradient !> @brief contains the constitutive equation for calculating the velocity gradient
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -335,7 +298,7 @@ subroutine thermal_adiabatic_LTAndItsTangent(LT, dLT_dTstar, Tstar_v, Lp, ipc, i
LT = Tdot*lattice_thermalExpansion33(1:3,1:3,phase) LT = Tdot*lattice_thermalExpansion33(1:3,1:3,phase)
dLT_dTstar3333 = 0.0_pReal dLT_dTstar3333 = 0.0_pReal
forall (i=1_pInt:3_pInt,j=1_pInt:3_pInt,k=1_pInt:3_pInt,l=1_pInt:3_pInt) & forall (i=1_pInt:3_pInt,j=1_pInt:3_pInt,k=1_pInt:3_pInt,l=1_pInt:3_pInt) &
dLT_dTstar3333(i,j,k,l) = dLT_dTstar3333(i,j,k,l) + Lp(k,l)*lattice_thermalExpansion33(i,j,phase) dLT_dTstar3333(i,j,k,l) = Lp(k,l)*lattice_thermalExpansion33(i,j,phase)
dLT_dTstar3333 = 0.95_pReal*dLT_dTstar3333/(lattice_massDensity(phase)*lattice_specificHeat(phase)) dLT_dTstar3333 = 0.95_pReal*dLT_dTstar3333/(lattice_massDensity(phase)*lattice_specificHeat(phase))
dLT_dTstar = math_Plain3333to99(dLT_dTstar3333) dLT_dTstar = math_Plain3333to99(dLT_dTstar3333)
@ -375,6 +338,47 @@ pure function thermal_adiabatic_getFT(ipc, ip, el)
end function thermal_adiabatic_getFT end function thermal_adiabatic_getFT
!--------------------------------------------------------------------------------------------------
!> @brief returns local thermal deformation gradient
!--------------------------------------------------------------------------------------------------
subroutine thermal_adiabatic_putFT(Tstar_v, Lp, dt, ipc, ip, el)
use material, only: &
mappingConstitutive, &
thermalState
use math, only: &
math_Mandel6to33
use lattice, only: &
lattice_massDensity, &
lattice_specificHeat, &
lattice_thermalExpansion33
implicit none
integer(pInt), intent(in) :: &
ipc, & !< grain number
ip, & !< integration point number
el !< element number
real(pReal), intent(in), dimension(6) :: &
Tstar_v !< 2nd Piola-Kirchhoff stress
real(pReal), intent(in), dimension(3,3) :: &
Lp !< plastic velocity gradient
real(pReal), intent(in) :: &
dt
integer(pInt) :: &
phase, &
constituent
real(pReal) :: &
Tdot
phase = mappingConstitutive(2,ipc,ip,el)
constituent = mappingConstitutive(1,ipc,ip,el)
Tdot = 0.95_pReal &
* sum(abs(math_Mandel6to33(Tstar_v))*Lp) &
/ (lattice_massDensity(phase)*lattice_specificHeat(phase))
thermalState(phase)%state(1,constituent) = thermalState(phase)%subState0(1,constituent) + Tdot*dt
end subroutine thermal_adiabatic_putFT
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief returns local thermal deformation gradient !> @brief returns local thermal deformation gradient
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------