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

In crystallite_stressAndItsTangent: 

state is now correctly collected during perturbation method
This commit is contained in:
Luc Hantcherli 2009-08-13 10:04:14 +00:00
parent 86211bf0ce
commit 93543c21e9
1 changed files with 42 additions and 59 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

101
code/crystallite.f90
View File

@ -270,6 +270,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
Tstar, & ! 2nd Piola-Kirchhoff stress tensor
myF, & ! local copy of the deformation gradient
myFp, & ! local copy of the plastic deformation gradient
myInvFp, & ! local copy of the invert of plastic deformation gradient
myFe, & ! local copy of the elastic deformation gradient
myLp, & ! local copy of the plastic velocity gradient
myP ! local copy of the 1st Piola-Kirchhoff stress tensor
@ -410,7 +411,6 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
.and. crystallite_onTrack &
.and. .not. crystallite_converged)
! --+>> preguess for state <<+--
!
! incrementing by crystallite_subdt
@ -450,7 +450,6 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
enddo; enddo; enddo
!$OMPEND PARALLEL DO
! --+>> state loop <<+--
NiterationState = 0_pInt
@ -460,7 +459,6 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
NiterationState = NiterationState + 1_pInt
! --+>> stress integration <<+--
!
! incrementing by crystallite_subdt
@ -472,15 +470,16 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
do e = FEsolving_execElem(1),FEsolving_execElem(2) ! iterate over elements to be processed
myNgrains = homogenization_Ngrains(mesh_element(3,e))
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed
do g = 1,myNgrains
! debugger = (e == 1 .and. i == 1 .and. g == 1)
if (crystallite_nonfinished(g,i,e)) & ! all undone crystallites
crystallite_onTrack(g,i,e) = crystallite_integrateStress(g,i,e)
enddo
do g = 1,myNgrains
debugger = (e == 1 .and. i == 1 .and. g == 1)
if (crystallite_nonfinished(g,i,e)) & ! all undone crystallites
crystallite_onTrack(g,i,e) = crystallite_integrateStress(g,i,e)
enddo
enddo
enddo
!$OMPEND PARALLEL DO
crystallite_nonfinished = crystallite_nonfinished .and. crystallite_onTrack
! --+>> state integration <<+--
!
@ -490,8 +489,6 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
! first loop for collection of state evolution based on old state
! second loop for updating to new state
crystallite_nonfinished = crystallite_nonfinished .and. crystallite_onTrack
!$OMP PARALLEL DO
do e = FEsolving_execElem(1),FEsolving_execElem(2) ! iterate over elements to be processed
myNgrains = homogenization_Ngrains(mesh_element(3,e))
@ -505,7 +502,6 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
myNgrains = homogenization_Ngrains(mesh_element(3,e))
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed
do g = 1,myNgrains
! debugger = (e == 1 .and. i == 1 .and. g == 1)
if (crystallite_nonfinished(g,i,e)) & ! all undone crystallites
call constitutive_collectDotState(crystallite_Tstar_v(:,g,i,e), crystallite_Fp(:,:,g,i,e), &
crystallite_invFp(:,:,g,i,e), crystallite_Temperature(g,i,e), g, i, e)
@ -516,11 +512,12 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
myNgrains = homogenization_Ngrains(mesh_element(3,e))
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed
do g = 1,myNgrains
! debugger = (e == 1 .and. i == 1 .and. g == 1)
debugger = (e == 1 .and. i == 1 .and. g == 1)
if (crystallite_nonfinished(g,i,e)) then ! all undone crystallites
crystallite_stateConverged(g,i,e) = crystallite_updateState(g,i,e) ! update state
crystallite_temperatureConverged(g,i,e) = crystallite_updateTemperature(g,i,e) ! update temperature
crystallite_converged(g,i,e) = crystallite_stateConverged(g,i,e) .and. crystallite_temperatureConverged(g,i,e)
if (debugger) write (6,*) g,i,e,'converged after updState',crystallite_converged(g,i,e)
if (crystallite_converged(g,i,e)) then
!$OMP CRITICAL (distributionState)
debug_CrystalliteStateLoopDistribution(NiterationState) = &
@ -573,7 +570,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
myNgrains = homogenization_Ngrains(mesh_element(3,e))
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed
do g = 1,myNgrains
! debugger = (g == 1 .and. i == 1 .and. e == 1)
! debugger = (g == 1 .and. i == 1 .and. e == 1)
if (crystallite_converged(g,i,e)) then ! grain converged in above iteration
mySizeState = constitutive_sizeState(g,i,e) ! number of state variables for this grain
mySizeDotState = constitutive_sizeDotState(g,i,e) ! number of dotStates for this grain
@ -582,6 +579,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
myTemperature = crystallite_Temperature(g,i,e) ! ... Temperature, ...
myF = crystallite_subF(:,:,g,i,e) ! ... and kinematics
myFp = crystallite_Fp(:,:,g,i,e)
myInvFp = crystallite_invFp(:,:,g,i,e)
myFe = crystallite_Fe(:,:,g,i,e)
myLp = crystallite_Lp(:,:,g,i,e)
myTstar_v = crystallite_Tstar_v(:,g,i,e)
@ -616,6 +614,9 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
NiterationState = NiterationState + 1_pInt
onTrack = crystallite_integrateStress(g,i,e) ! stress of perturbed situation (overwrites _P,_Tstar_v,_Fp,_Lp,_Fe)
if (onTrack) then
call constitutive_collectDotState(crystallite_Tstar_v(:,g,i,e), crystallite_Fp(:,:,g,i,e), &
crystallite_invFp(:,:,g,i,e), crystallite_Temperature(g,i,e), g, i, e)
stateConverged = crystallite_updateState(g,i,e) ! update state
temperatureConverged = crystallite_updateTemperature(g,i,e) ! update temperature
converged = stateConverged .and. temperatureConverged
@ -637,13 +638,15 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
constitutive_dotState(g,i,e)%p = myDotState ! ... dotState, ...
crystallite_Temperature(g,i,e) = myTemperature ! ... temperature, ...
crystallite_Fp(:,:,g,i,e) = myFp ! ... and kinematics
crystallite_invFp(:,:,g,i,e) = myInvFp
crystallite_Fe(:,:,g,i,e) = myFe
crystallite_Lp(:,:,g,i,e) = myLp
crystallite_Tstar_v(:,g,i,e) = myTstar_v
crystallite_P(:,:,g,i,e) = myP
!$OMP CRITICAL (out)
debug_StiffnessStateLoopDistribution(NiterationState) = &
debug_StiffnessstateLoopDistribution(NiterationState) + 1
debug_StiffnessstateLoopDistribution(NiterationState) + 1
if (nState < NiterationState) write(6,*) 'ohh shit!! stiffenss state loop debugging exceeded',NiterationState
!$OMPEND CRITICAL (out)
enddo
enddo
@ -702,7 +705,7 @@ endsubroutine
maxticks
mySize = constitutive_sizeDotState(g,i,e)
! calculate the residuum
call system_clock(count=tick,count_rate=tickrate,count_max=maxticks)
residuum = constitutive_state(g,i,e)%p(1:mySize) - constitutive_subState0(g,i,e)%p(1:mySize) - &
@ -715,33 +718,20 @@ endsubroutine
! if NaN occured then return without changing the state
if (any(residuum/=residuum)) then
crystallite_updateState = .false. ! indicate state update failed
if (debugger) then
!$OMP CRITICAL (write2out)
write(6,*) '::: updateState encountered NaN',g,i,e
write(6,*)
!$OMPEND CRITICAL (write2out)
endif
!$OMP CRITICAL (write2out)
write(6,*) '::: updateState encountered NaN',e,i,g
!$OMPEND CRITICAL (write2out)
return
endif
! update the microstructure
constitutive_state(g,i,e)%p(1:mySize) = constitutive_state(g,i,e)%p(1:mySize) - residuum
! setting flag to true if state is below relative tolerance, otherwise set it to false
! setting flag to true if state is below relative tolerance, otherwise set it to false <<<updated 31.07.2009>>>
crystallite_updateState = all(constitutive_state(g,i,e)%p(1:mySize) == 0.0_pReal .or. &
abs(residuum) < rTol_crystalliteState*abs(constitutive_state(g,i,e)%p(1:mySize)))
if (debugger) then
!$OMP CRITICAL (write2out)
if (crystallite_updateState) then
write(6,*) '::: updateState did not converge',g,i,e
write(6,*)
else
write(6,*) '::: updateState converged',g,i,e
write(6,*)
endif
write(6,'(a,/,12(f10.5,x))') 'resid tolerance',abs(residuum/rTol_crystalliteState/constitutive_state(g,i,e)%p(1:mySize))
write(6,*)
!$OMPEND CRITICAL (write2out)
write(6,'(a,/,12(f10.5,x))') 'resid tolerance',abs(residuum/rTol_crystalliteState/constitutive_state(g,i,e)%p(1:mySize))
endif
return
@ -797,7 +787,7 @@ endsubroutine
if (residuum/=residuum) then
crystallite_updateTemperature = .false. ! indicate update failed
!$OMP CRITICAL (write2out)
write(6,*) '::: updateTemperature encountered NaN',g,i,e
write(6,*) '::: updateTemperature encountered NaN',e,i,g
!$OMPEND CRITICAL (write2out)
return
endif
@ -805,7 +795,7 @@ endsubroutine
! update the microstructure
crystallite_Temperature(g,i,e) = crystallite_Temperature(g,i,e) - residuum
! setting flag to true if residuum is below relative tolerance (or zero Kelvin), otherwise set it to false
! setting flag to true if residuum is below relative tolerance (or zero Kelvin), otherwise set it to false <<<updated 31.07.2009>>>
crystallite_updateTemperature = crystallite_Temperature(g,i,e) == 0.0_pReal .or. &
abs(residuum) < rTol_crystalliteTemperature*crystallite_Temperature(g,i,e)
@ -920,13 +910,9 @@ endsubroutine
! inversion of Fp_current...
invFp_current = math_inv3x3(Fp_current)
if (all(invFp_current == 0.0_pReal)) then ! ... failed?
if (debugger) then
!$OMP CRITICAL (write2out)
write(6,*) '::: integrateStress failed on invFp_current inversion',g,i,e
write(6,*)
write(6,'(a11,3(i3,x),/,3(3(f12.7,x)/))') 'invFp_new at ',g,i,e,invFp_new
!$OMPEND CRITICAL (write2out)
endif
!$OMP CRITICAL (write2out)
write(6,*) '::: integrateStress failed on invFp_current inversion',e,i,g
!$OMPEND CRITICAL (write2out)
return
endif
@ -987,11 +973,9 @@ LpLoop: do
! NaN occured at regular speed?
if (any(residuum/=residuum) .and. leapfrog == 1.0) then
if (debugger) then
!$OMP CRITICAL (write2out)
write(6,*) '::: integrateStress encountered NaN at iteration', NiterationStress,'at',g,i,e
!$OMPEND CRITICAL (write2out)
endif
!$OMP CRITICAL (write2out)
write(6,*) '::: integrateStress encountered NaN at iteration', NiterationStress,'at',e,i,g
!$OMPEND CRITICAL (write2out)
return
! something went wrong at accelerated speed?
@ -1024,11 +1008,11 @@ LpLoop: do
if (error) then
if (debugger) then
!$OMP CRITICAL (write2out)
write(6,*) '::: integrateStress failed on dR/dLp inversion at iteration', NiterationStress
write(6,*)
write(6,'(a9,3(i3,x),/,9(9(f12.7,x)/))') 'dRdLp at ',g,i,e,dRdLp
write(6,'(a20,3(i3,x),/,9(9(e12.2,x)/))') 'dLp_constitutive at ',g,i,e,dLp_constitutive
write(6,'(a11,3(i3,x),/,3(3(f12.7,x)/))') 'Lpguess at ',g,i,e,Lpguess
write(6,*) '::: integrateStress failed on dR/dLp inversion at iteration', NiterationStress
write(6,'(a9,3(i3,x),/,9(9(e12.2,x)/))') 'dTdLp at ',g,i,e,dTdLp
write(6,'(a20,3(i3,x),/,9(9(e12.2,x)/))') 'dLp_constitutive at ',g,i,e,dLp_constitutive
write(6,'(a9,3(i3,x),/,9(9(f12.7,x)/))') 'dRdLp at ',g,i,e,dRdLp
write(6,'(a11,3(i3,x),/,3(3(f12.7,x)/))') 'Lpguess at ',g,i,e,Lpguess
!$OMPEND CRITICAL (write2out)
endif
return
@ -1056,9 +1040,7 @@ LpLoop: do
if (error) then
if (debugger) then
!$OMP CRITICAL (write2out)
write(6,*) '::: integrateStress failed on invFp_new inversion at iteration', NiterationStress
write(6,*)
write(6,'(a11,3(i3,x),/,3(3(f12.7,x)/))') 'invFp_new at ',g,i,e,invFp_new
write(6,*) '::: integrateStress failed on invFp_new inversion at iteration', NiterationStress
!$OMPEND CRITICAL (write2out)
endif
return
@ -1082,15 +1064,16 @@ LpLoop: do
crystallite_integrateStress = .true.
if (debugger) then
!$OMP CRITICAL (write2out)
write(6,*) '::: integrateStress converged at iteration', NiterationStress
write(6,*)
write(6,'(a,/,3(3(f12.7,x)/))') 'P / MPa',crystallite_P(:,:,g,i,e)/1e6
write(6,'(a,/,3(3(f12.7,x)/))') 'Lp',crystallite_Lp(:,:,g,i,e)
write(6,*) '::: integrateStress converged at iteration', NiterationStress
write(6,*)
write(6,'(a,/,3(3(f12.7,x)/))') 'P / MPa',crystallite_P(:,:,g,i,e)/1e6
write(6,'(a,/,3(3(f12.7,x)/))') 'Lp',crystallite_Lp(:,:,g,i,e)
!$OMP CRITICAL (write2out)
endif
!$OMP CRITICAL (distributionStress)
debug_StressLoopDistribution(NiterationStress) = debug_StressLoopDistribution(NiterationStress) + 1
if (nStress < NiterationStress) write(6,*) 'ohh shit!! debug loop of stress exceeded',NiterationStress
!$OMPEND CRITICAL (distributionStress)
return