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

restructured stress integration. Double loop for stress integration (Lp and stress) with cutback step in IntegrateStress is replaced by single loop for Lp without cutback. loop counter and limits are accordingly renamed and implemented

This commit is contained in:
Christoph Kords 2009-05-28 16:38:40 +00:00
parent aead2b8d20
commit ab57a40bab
5 changed files with 356 additions and 332 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

40
trunk/IO.f90
View File

@ -837,19 +837,19 @@ END FUNCTION
!$OMP CRITICAL (write2out) !$OMP CRITICAL (write2out)
write(6,*) write(6,*)
write(6,*) '+------------------------------+' write(6,'(a38)') '+------------------------------------+'
write(6,*) '+ ERROR +' write(6,'(a38)') '+ error +'
write(6,*) '+ +' write(6,'(a38)') '+ +'
write(6,*) msg write(6,*) '+ ',msg
if (present(ext_msg)) write(6,*) ext_msg if (present(ext_msg)) write(6,*) '+ ',ext_msg
if (present(e)) then if (present(e)) then
if (present(i) .and. present(g)) then if (present(i) .and. present(g)) then
write(6,'(a10,x,i6,x,a2,x,i2,x,a5,x,i4)') 'at element',e,'IP',i,'grain',g write(6,'(a12,x,i6,x,a2,x,i2,x,a5,x,i4,a2)') '+ at element',e,'IP',i,'grain',g,' +'
else else
write(6,'(a2,x,i6)') 'at',e write(6,'(a17,i6,a14)') '+ at ',e,' +'
endif endif
endif endif
write(6,*) '+------------------------------+' write(6,'(a38)') '+------------------------------------+'
call debug_info() call debug_info()
call flush(6) call flush(6)
@ -877,29 +877,29 @@ END FUNCTION
character(len=80) msg character(len=80) msg
select case (ID) select case (ID)
case (650)
msg = 'Polar decomposition failed'
case (600) case (600)
msg = 'Crystallite responds elastically' msg = '+ Crystallite responds elastically +'
case (650)
msg = '+ Polar decomposition failed +'
case default case default
msg = 'Unknown warning number...' msg = '+ Unknown warning number... +'
end select end select
!$OMP CRITICAL (write2out) !$OMP CRITICAL (write2out)
write(6,*) write(6,*)
write(6,*) '+------------------------------+' write(6,'(a38)') '+------------------------------------+'
write(6,*) '+ warning +' write(6,'(a38)') '+ warning +'
write(6,*) '+ +' write(6,'(a38)') '+ +'
write(6,*) msg write(6,'(a38)') msg
if (present(ext_msg)) write(6,*) ext_msg if (present(ext_msg)) write(6,*) '+ ',ext_msg
if (present(e)) then if (present(e)) then
if (present(i) .and. present(g)) then if (present(i) .and. present(g)) then
write(6,'(a10,x,i6,x,a2,x,i2,x,a5,x,i4)') 'at element',e,'IP',i,'grain',g write(6,'(a12,x,i6,x,a2,x,i2,x,a5,x,i4,a2)') '+ at element',e,'IP',i,'grain',g,' +'
else else
write(6,'(a2,x,i6)') 'at',e write(6,'(a17,i6,a14)') '+ at ',e,' +'
endif endif
endif endif
write(6,*) '+------------------------------+' write(6,'(a38)') '+------------------------------------+'
END SUBROUTINE END SUBROUTINE

561
trunk/crystallite.f90
View File

@ -164,7 +164,7 @@ MODULE crystallite
return return
end subroutine endsubroutine
!******************************************************************** !********************************************************************
@ -176,7 +176,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
use debug use debug
use IO, only: IO_warning use IO, only: IO_warning
use math use math
use FEsolving, only: FEsolving_execElem, FEsolving_execIP use FEsolving, only: FEsolving_execElem, FEsolving_execIP, theInc
use mesh, only: mesh_element use mesh, only: mesh_element
use material, only: homogenization_Ngrains use material, only: homogenization_Ngrains
use constitutive use constitutive
@ -185,18 +185,18 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
logical, intent(in) :: updateJaco logical, intent(in) :: updateJaco
real(pReal), dimension(3,3) :: invFp,Fe_guess,PK2,myF,myFp,myFe,myLp,myP real(pReal), dimension(3,3) :: invFp,Fe_guess,PK2,myF,myFp,myFe,myLp,myP
real(pReal), dimension(constitutive_maxSizeState) :: myState real(pReal), dimension(constitutive_maxSizeState) :: myState
integer(pInt) crystallite_Niteration integer(pInt) NiterationCrystallite, NiterationState
integer(pInt) g,i,e,k,l, myNgrains, mySizeState integer(pInt) g,i,e,k,l, myNgrains, mySizeState
logical converged logical onTrack,converged
! ------ initialize to starting condition ------ ! ------ initialize to starting condition ------
write (6,*) write (6,*)
write (6,*) 'Crystallite request from Materialpoint' write (6,*) 'Crystallite request from Materialpoint'
write (6,'(a,/,3(3(f12.7,x)/))') 'crystallite_partionedF0 of 1 8 1',crystallite_partionedF0(1:3,:,1,8,1) write (6,'(a,/,3(3(f12.7,x)/))') 'crystallite_partionedF0 of 1 1 1',crystallite_partionedF0(1:3,:,1,1,1)
write (6,'(a,/,3(3(f12.7,x)/))') 'crystallite_partionedFp0 of 1 8 1',crystallite_partionedFp0(1:3,:,1,8,1) write (6,'(a,/,3(3(f12.7,x)/))') 'crystallite_partionedFp0 of 1 1 1',crystallite_partionedFp0(1:3,:,1,1,1)
write (6,'(a,/,3(3(f12.7,x)/))') 'crystallite_partionedF of 1 8 1',crystallite_partionedF(1:3,:,1,8,1) write (6,'(a,/,3(3(f12.7,x)/))') 'crystallite_partionedF of 1 1 1',crystallite_partionedF(1:3,:,1,1,1)
write (6,'(a,/,3(3(f12.7,x)/))') 'crystallite_partionedLp0 of 1 8 1',crystallite_partionedLp0(1:3,:,1,8,1) write (6,'(a,/,3(3(f12.7,x)/))') 'crystallite_partionedLp0 of 1 1 1',crystallite_partionedLp0(1:3,:,1,1,1)
!$OMP PARALLEL DO !$OMP PARALLEL DO
@ -222,19 +222,24 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
! ------ cutback loop ------ ! ------ cutback loop ------
NiterationCrystallite = 0_pInt
do while (any(crystallite_subStep(:,:,FEsolving_execELem(1):FEsolving_execElem(2)) > subStepMin)) ! cutback loop for crystallites do while (any(crystallite_subStep(:,:,FEsolving_execELem(1):FEsolving_execElem(2)) > subStepMin)) ! cutback loop for crystallites
write(6,*)
write(6,*) 'entering cutback at crystallite_stress' NiterationCrystallite = NiterationCrystallite + 1
if (any(.not. crystallite_converged(:,:,FEsolving_execELem(1):FEsolving_execElem(2)) .and. & ! any non-converged grain if (any(.not. crystallite_converged(:,:,FEsolving_execELem(1):FEsolving_execElem(2)) .and. & ! any non-converged grain
.not. crystallite_localConstitution(:,:,FEsolving_execELem(1):FEsolving_execElem(2))) ) & ! has non-local constitution? .not. crystallite_localConstitution(:,:,FEsolving_execELem(1):FEsolving_execElem(2))) ) & ! has non-local constitution?
crystallite_converged(:,:,FEsolving_execELem(1):FEsolving_execElem(2)) = & crystallite_converged(:,:,FEsolving_execELem(1):FEsolving_execElem(2)) = &
crystallite_converged(:,:,FEsolving_execELem(1):FEsolving_execElem(2)) .and. & crystallite_converged(:,:,FEsolving_execELem(1):FEsolving_execElem(2)) .and. &
crystallite_localConstitution(:,:,FEsolving_execELem(1):FEsolving_execElem(2)) ! reset non-local grains' convergence status crystallite_localConstitution(:,:,FEsolving_execELem(1):FEsolving_execElem(2)) ! reset non-local grains' convergence status
!$OMP PARALLEL DO !$OMP PARALLEL DO
do e = FEsolving_execElem(1),FEsolving_execElem(2) ! iterate over elements to be processed do e = FEsolving_execElem(1),FEsolving_execElem(2) ! iterate over elements to be processed
myNgrains = homogenization_Ngrains(mesh_element(3,e)) 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 i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed
do g = 1,myNgrains do g = 1,myNgrains
debugger = (g == 1 .and. i == 1 .and. e == 1)
if (crystallite_converged(g,i,e)) then if (crystallite_converged(g,i,e)) then
crystallite_subFrac(g,i,e) = crystallite_subFrac(g,i,e) + crystallite_subStep(g,i,e) crystallite_subFrac(g,i,e) = crystallite_subFrac(g,i,e) + crystallite_subStep(g,i,e)
crystallite_subStep(g,i,e) = min(1.0_pReal-crystallite_subFrac(g,i,e), 2.0_pReal * crystallite_subStep(g,i,e)) crystallite_subStep(g,i,e) = min(1.0_pReal-crystallite_subFrac(g,i,e), 2.0_pReal * crystallite_subStep(g,i,e))
@ -244,11 +249,26 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
crystallite_subLp0(:,:,g,i,e) = crystallite_Lp(:,:,g,i,e) ! ...plastic velocity gradient crystallite_subLp0(:,:,g,i,e) = crystallite_Lp(:,:,g,i,e) ! ...plastic velocity gradient
constitutive_subState0(g,i,e)%p = constitutive_state(g,i,e)%p ! ...microstructure constitutive_subState0(g,i,e)%p = constitutive_state(g,i,e)%p ! ...microstructure
endif endif
if (debugger) then
!$OMP CRITICAL (write2out)
write(6,'(a21,f6.4,a28,f6.4,a35)') 'winding forward from ', &
crystallite_subFrac(g,i,e)-crystallite_subStep(g,i,e),' to new crystallite_subfrac ', &
crystallite_subFrac(g,i,e),' in crystallite_stressAndItsTangent'
write(6,*)
!$OMP END CRITICAL (write2out)
endif
else else
crystallite_subStep(g,i,e) = 0.5_pReal * crystallite_subStep(g,i,e) ! cut step in half and restore... crystallite_subStep(g,i,e) = 0.5_pReal * crystallite_subStep(g,i,e) ! cut step in half and restore...
crystallite_Fp(:,:,g,i,e) = crystallite_subFp0(:,:,g,i,e) ! ...plastic def grad crystallite_Fp(:,:,g,i,e) = crystallite_subFp0(:,:,g,i,e) ! ...plastic def grad
crystallite_Lp(:,:,g,i,e) = crystallite_subLp0(:,:,g,i,e) ! ...plastic velocity grad crystallite_Lp(:,:,g,i,e) = crystallite_subLp0(:,:,g,i,e) ! ...plastic velocity grad
constitutive_state(g,i,e)%p = constitutive_subState0(g,i,e)%p ! ...microstructure constitutive_state(g,i,e)%p = constitutive_subState0(g,i,e)%p ! ...microstructure
if (debugger) then
!$OMP CRITICAL (write2out)
write(6,'(a78,f6.4)') 'cutback step in crystallite_stressAndItsTangent with new crystallite_subStep: ',&
crystallite_subStep(g,i,e)
write(6,*)
!$OMP END CRITICAL (write2out)
endif
endif endif
crystallite_onTrack(g,i,e) = crystallite_subStep(g,i,e) > subStepMin ! still on track or already done (beyond repair) crystallite_onTrack(g,i,e) = crystallite_subStep(g,i,e) > subStepMin ! still on track or already done (beyond repair)
@ -257,6 +277,12 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
crystallite_subStep(g,i,e) * & crystallite_subStep(g,i,e) * &
(crystallite_partionedF(:,:,g,i,e) - crystallite_partionedF0(:,:,g,i,e)) (crystallite_partionedF(:,:,g,i,e) - crystallite_partionedF0(:,:,g,i,e))
crystallite_subdt(g,i,e) = crystallite_subStep(g,i,e) * crystallite_dt(g,i,e) crystallite_subdt(g,i,e) = crystallite_subStep(g,i,e) * crystallite_dt(g,i,e)
if (debugger) then
!$OMP CRITICAL (write2out)
write(6,'(a36,e8.3)') 'current timestep crystallite_subdt: ',crystallite_subdt(g,i,e)
write(6,*)
!$OMP END CRITICAL (write2out)
endif
crystallite_converged(g,i,e) = .false. ! start out non-converged crystallite_converged(g,i,e) = .false. ! start out non-converged
endif endif
enddo enddo
@ -266,14 +292,15 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
! ------ convergence loop for stress and state ------ ! ------ convergence loop for stress and state ------
crystallite_Niteration = 0_pInt NiterationState = 0_pInt
if (debugger) write(6,*) 'state integration started'
do while (any( crystallite_requested(:,:,FEsolving_execELem(1):FEsolving_execElem(2)) & do while (any( crystallite_requested(:,:,FEsolving_execELem(1):FEsolving_execElem(2)) &
.and. crystallite_onTrack(:,:,FEsolving_execELem(1):FEsolving_execElem(2)) & .and. crystallite_onTrack(:,:,FEsolving_execELem(1):FEsolving_execElem(2)) &
.and. .not. crystallite_converged(:,:,FEsolving_execELem(1):FEsolving_execElem(2)) & .and. .not. crystallite_converged(:,:,FEsolving_execELem(1):FEsolving_execElem(2)) &
) .and. crystallite_Niteration < nCryst) ! convergence loop for crystallite ) .and. NiterationState < nState) ! convergence loop for crystallite
crystallite_Niteration = crystallite_Niteration + 1
NiterationState = NiterationState + 1
! --+>> stress integration <<+-- ! --+>> stress integration <<+--
! !
@ -283,6 +310,8 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
! to account for substepping within _integrateStress ! to account for substepping within _integrateStress
! results in crystallite_Fp,.._Lp ! results in crystallite_Fp,.._Lp
if (debugger) write(6,*) 'stress integration started'
!$OMP PARALLEL DO !$OMP PARALLEL DO
do e = FEsolving_execElem(1),FEsolving_execElem(2) ! iterate over elements to be processed do e = FEsolving_execElem(1),FEsolving_execElem(2) ! iterate over elements to be processed
myNgrains = homogenization_Ngrains(mesh_element(3,e)) myNgrains = homogenization_Ngrains(mesh_element(3,e))
@ -296,8 +325,11 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
enddo enddo
!$OMP END PARALLEL DO !$OMP END PARALLEL DO
write(6,'(i2,x,a10,x,16(l,x))') crystallite_Niteration,'cryst_onT',crystallite_onTrack if (crystallite_requested(1,1,1) .and. crystallite_onTrack(1,1,1)) then
write (6,'(a,/,3(3(f12.7,x)/))') 'Lp of 1 8 1',crystallite_Lp(1:3,:,1,8,1) write(6,*) 'stress integration converged'
write(6,'(a,/,3(3(e15.7,x)/))') 'P of 1 1 1',crystallite_P(:,:,1,1,1)
write(6,'(a,/,3(3(f12.7,x)/))') 'Lp of 1 1 1',crystallite_Lp(:,:,1,1,1)
endif
! --+>> state integration <<+-- ! --+>> state integration <<+--
! !
@ -315,9 +347,12 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
crystallite_converged(g,i,e) = crystallite_updateState(g,i,e) crystallite_converged(g,i,e) = crystallite_updateState(g,i,e)
if (crystallite_converged(g,i,e)) then if (crystallite_converged(g,i,e)) then
!$OMP CRITICAL (distributionState) !$OMP CRITICAL (distributionState)
debug_StateLoopDistribution(crystallite_Niteration) = & debug_StateLoopDistribution(NiterationState) = debug_StateLoopDistribution(NiterationState) + 1
debug_StateLoopDistribution(crystallite_Niteration) + 1
!$OMP END CRITICAL (distributionState) !$OMP END CRITICAL (distributionState)
!$OMP CRITICAL (distributionCrystallite)
debug_CrystalliteLoopDistribution(NiterationCrystallite) = &
debug_CrystalliteLoopDistribution(NiterationCrystallite) + 1
!$OMP END CRITICAL (distributionCrystallite)
endif endif
endif endif
enddo enddo
@ -325,15 +360,14 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
enddo enddo
!$OMP END PARALLEL DO !$OMP END PARALLEL DO
if (crystallite_requested(1,1,1) .and. crystallite_onTrack(1,1,1) .and. crystallite_converged(1,1,1)) then
write(6,*) 'state integration converged'
write(6,'(a20,e8.3)') 'state of 1 1 1: ', constitutive_state(1,1,1)%p(1)
write(6,*)
endif
enddo ! crystallite convergence loop enddo ! crystallite convergence loop
write(6,*)
write(6,'(a10,x,16(f6.4,x))') 'cryst_frac',crystallite_subFrac
write(6,'(a10,x,16(f6.4,x))') 'cryst_step',crystallite_subStep
write(6,'(a10,x,16(l,x))') 'cryst_req',crystallite_requested
write(6,'(a10,x,16(l,x))') 'cryst_onT',crystallite_onTrack
write(6,'(a10,x,16(l,x))') 'cryst_cvg',crystallite_converged
write(6,'(a10,x,16(e8.3,x))') 'cryst_dt',crystallite_subdt
enddo ! cutback loop enddo ! cutback loop
@ -368,9 +402,9 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
if(updateJaco) then ! Jacobian required if(updateJaco) then ! Jacobian required
if (debugger) then if (debugger) then
write (6,*) 'Stiffness calculation started'
!$OMP CRITICAL (write2out) !$OMP CRITICAL (write2out)
write (6,*) 'Jacobian calc' ! write(6,'(a10,x,16(f6.4,x))') 'cryst_dt',crystallite_subdt
write(6,'(a10,x,16(f6.4,x))') 'cryst_dt',crystallite_subdt
!$OMP END CRITICAL (write2out) !$OMP END CRITICAL (write2out)
endif endif
@ -387,30 +421,63 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
myFe = crystallite_Fe(:,:,g,i,e) myFe = crystallite_Fe(:,:,g,i,e)
myLp = crystallite_Lp(:,:,g,i,e) myLp = crystallite_Lp(:,:,g,i,e)
myP = crystallite_P(:,:,g,i,e) myP = crystallite_P(:,:,g,i,e)
if (debugger) then
write (6,*) '#############'
write (6,*) 'central solution'
write (6,*) '#############'
write (6,'(a,/,3(3(f12.4,x)/))') ' P of 1 1 1',myP(1:3,:)/1e6
write (6,'(a,/,3(3(f12.8,x)/))') ' Fp of 1 1 1',myFp(1:3,:)
write (6,'(a,/,3(3(f12.8,x)/))') ' Lp of 1 1 1',myLp(1:3,:)
write (6,'(a,/,f12.4)') 'state of 1 1 1',myState/1e6
endif
do k = 1,3 ! perturbation... do k = 1,3 ! perturbation...
do l = 1,3 ! ...components do l = 1,3 ! ...components
crystallite_subF(:,:,g,i,e) = myF ! initialize perturbed F to match converged crystallite_subF(:,:,g,i,e) = myF ! initialize perturbed F to match converged
crystallite_subF(k,l,g,i,e) = crystallite_subF(k,l,g,i,e) + pert_Fg ! perturb single component crystallite_subF(k,l,g,i,e) = crystallite_subF(k,l,g,i,e) + pert_Fg ! perturb single component
converged = .false. onTrack = .true.
crystallite_Niteration = 0_pInt converged = .false.
do while(.not. converged .and. crystallite_Niteration < nCryst) ! keep cycling until done (potentially non-converged) NiterationState = 0_pInt
crystallite_Niteration = crystallite_Niteration + 1_pInt if (debugger) then
if(crystallite_integrateStress(g,i,e)) & ! stress of perturbed situation (overwrites_P,_Tstar_v,_Fp,_Lp,_Fe) write (6,*) '============='
converged = crystallite_updateState(g,i,e) write (6,'(i1,x,i1)') k,l
end do write (6,*) '============='
if (converged) & ! converged state warrants stiffness update write (6,'(a,/,3(3(f12.6,x)/))') 'pertF of 1 1 1',crystallite_subF(1:3,:,g,i,e)
crystallite_dPdF(:,:,k,l,g,i,e) =(crystallite_p(:,:,g,i,e) - myP)/pert_Fg ! tangent dP_ij/dFg_kl endif
do while(.not. converged .and. onTrack .and. NiterationState < nState) ! keep cycling until done (potentially non-converged)
NiterationState = NiterationState + 1_pInt
if (debugger) write (6,'(a4,x,i6)') 'loop',NiterationState
onTrack = crystallite_integrateStress(g,i,e) ! stress of perturbed situation (overwrites _P,_Tstar_v,_Fp,_Lp,_Fe)
if(onTrack) converged = crystallite_updateState(g,i,e)
if (debugger) then
write (6,*) '-------------'
write (6,'(l,x,l)') onTrack,converged
write (6,'(a,/,3(3(f12.4,x)/))') 'pertP of 1 1 1',crystallite_P(1:3,:,g,i,e)/1e6
write (6,'(a,/,3(3(f12.4,x)/))') 'DP of 1 1 1',(crystallite_P(1:3,:,g,i,e)-myP(1:3,:))/1e6
write (6,'(a,/,f12.4)') 'state of 1 1 1',constitutive_state(g,i,e)%p/1e6
write (6,'(a,/,f12.4)') 'Dstate of 1 1 1',(constitutive_state(g,i,e)%p-myState)/1e6
endif
enddo
if (converged) & ! converged state warrants stiffness update
crystallite_dPdF(:,:,k,l,g,i,e) = (crystallite_P(:,:,g,i,e) - myP)/pert_Fg ! tangent dP_ij/dFg_kl
constitutive_state(g,i,e)%p = myState ! restore unperturbed, converged state...
crystallite_Fp(:,:,g,i,e) = myFp ! ... and kinematics
crystallite_Fe(:,:,g,i,e) = myFe
crystallite_Lp(:,:,g,i,e) = myLp
crystallite_P(:,:,g,i,e) = myP
!$OMP CRITICAL (out) !$OMP CRITICAL (out)
debug_StiffnessStateLoopDistribution(crystallite_Niteration) = & debug_StiffnessStateLoopDistribution(NiterationState) = &
debug_StiffnessstateLoopDistribution(crystallite_Niteration) + 1 debug_StiffnessstateLoopDistribution(NiterationState) + 1
!$OMP END CRITICAL (out) !$OMP END CRITICAL (out)
end do enddo
end do enddo
constitutive_state(g,i,e)%p = myState ! restore unperturbed, converged state... constitutive_state(g,i,e)%p = myState ! restore unperturbed, converged state...
crystallite_Fp(:,:,g,i,e) = myFp ! ... and kinematics crystallite_Fp(:,:,g,i,e) = myFp ! ... and kinematics
crystallite_Fe(:,:,g,i,e) = myFe crystallite_Fe(:,:,g,i,e) = myFe
crystallite_Lp(:,:,g,i,e) = myLp crystallite_Lp(:,:,g,i,e) = myLp
crystallite_P(:,:,g,i,e) = myP crystallite_P(:,:,g,i,e) = myP
if (e == 1 .and. i == 1 .and. g == 1) then
write (6,'(a,/,9(9(f12.4,x)/))') 'dPdF/GPa',crystallite_dPdF(:,:,:,:,1,1,1)/1e9
endif
else ! grain has not converged else ! grain has not converged
crystallite_dPdF(:,:,:,:,g,i,e) = crystallite_fallbackdPdF(:,:,:,:,g,i,e) ! use fallback crystallite_dPdF(:,:,:,:,g,i,e) = crystallite_fallbackdPdF(:,:,:,:,g,i,e) ! use fallback
endif endif
@ -420,7 +487,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
!$OMP END PARALLEL DO !$OMP END PARALLEL DO
endif endif
end subroutine endsubroutine
@ -434,9 +501,13 @@ end subroutine
i,& ! integration point number i,& ! integration point number
e & ! element number e & ! element number
) )
use prec, only: pReal,pInt,rTol_crystalliteState use prec, only: pReal, &
use constitutive, only: constitutive_dotState,constitutive_sizeDotState,& pInt, &
constitutive_subState0,constitutive_state rTol_crystalliteState
use constitutive, only: constitutive_dotState, &
constitutive_sizeDotState, &
constitutive_subState0, &
constitutive_state
use debug use debug
logical crystallite_updateState logical crystallite_updateState
@ -455,18 +526,21 @@ end subroutine
debug_cumDotStateCalls = debug_cumDotStateCalls + 1_pInt debug_cumDotStateCalls = debug_cumDotStateCalls + 1_pInt
debug_cumDotStateTicks = debug_cumDotStateTicks + tock-tick debug_cumDotStateTicks = debug_cumDotStateTicks + tock-tick
if (tock < tick) debug_cumDotStateTicks = debug_cumDotStateTicks + maxticks if (tock < tick) debug_cumDotStateTicks = debug_cumDotStateTicks + maxticks
if (any(constitutive_state(g,i,e)%p(1:mySize)/=constitutive_state(g,i,e)%p(1:mySize))) return ! NaN occured?
constitutive_state(g,i,e)%p(1:mySize) = constitutive_state(g,i,e)%p(1:mySize) - residuum ! update of microstructure constitutive_state(g,i,e)%p(1:mySize) = constitutive_state(g,i,e)%p(1:mySize) - residuum ! update of microstructure
crystallite_updateState = maxval(abs(residuum/constitutive_state(g,i,e)%p(1:mySize)),& crystallite_updateState = maxval(abs(residuum/constitutive_state(g,i,e)%p(1:mySize)),&
constitutive_state(g,i,e)%p(1:mySize) /= 0.0_pReal) < rTol_crystalliteState constitutive_state(g,i,e)%p(1:mySize) /= 0.0_pReal) < rTol_crystalliteState
return return
end function endfunction
!*********************************************************************** !***********************************************************************
!*** calculation of stress (P), stiffness (dPdF), *** !*** calculation of stress (P) with time integration ***
!*** and announcement of any *** !*** based on a residuum in Lp and intermediate ***
!*** acceleration of the Newton-Raphson correction *** !*** acceleration of the Newton-Raphson correction ***
!*********************************************************************** !***********************************************************************
function crystallite_integrateStress(& function crystallite_integrateStress(&
@ -474,264 +548,224 @@ end subroutine
i,& ! integration point number i,& ! integration point number
e) ! element number e) ! element number
use prec, only: pReal,pInt,subStepMin
use debug
use constitutive, only: constitutive_state,constitutive_subState0,constitutive_sizeState
use math
implicit none
integer(pInt), intent(in) :: e,i,g
integer(pInt) mySize, Niteration
real(pReal) dt_aim,subFrac,subStep,det
real(pReal), dimension(3,3) :: inv
real(pReal), dimension(3,3) :: Fg_current,Fg_aim,deltaFg
real(pReal), dimension(3,3) :: Fp_current,Fp_new
real(pReal), dimension(3,3) :: Fe_current,Fe_new
real(pReal), dimension(constitutive_sizeState(g,i,e)) :: interpolatedState
logical crystallite_integrateStress ! still on track?
mySize = constitutive_sizeState(g,i,e)
deltaFg = crystallite_subF(:,:,g,i,e) - crystallite_subF0(:,:,g,i,e)
Fg_current = crystallite_subF0(:,:,g,i,e) ! initialize to start of inc
Fp_current = crystallite_subFp0(:,:,g,i,e)
Fe_current = math_mul33x33(Fg_current,math_inv3x3(Fp_current))
subFrac = 0.0_pReal
subStep = 1.0_pReal
Niteration = 0_pInt
crystallite_integrateStress = .false. ! be pessimisitc
! begin the cutback loop
do while (subStep > subStepMin) ! continue until finished or too much cut backing
Niteration = Niteration + 1
Fg_aim = Fg_current + subStep*deltaFg ! aim for Fg
dt_aim = subStep*crystallite_subdt(g,i,e) ! aim for dt
debugger = (g == 1 .and. i == 8 .and. e == 1)
call TimeIntegration(crystallite_integrateStress,&
crystallite_Lp(:,:,g,i,e),crystallite_Fp(:,:,g,i,e),crystallite_Fe(:,:,g,i,e),&
crystallite_Tstar_v(:,g,i,e),crystallite_P(:,:,g,i,e), &
Fg_aim,Fp_current,crystallite_Temperature(g,i,e),&
( subFrac+subStep)*constitutive_state(g,i,e)%p(1:mySize) + &
(1.0_pReal-subFrac-subStep)*constitutive_subState0(g,i,e)%p(1:mySize),& ! interpolated state
dt_aim,g,i,e)
if (crystallite_integrateStress) then ! happy with time integration
if (e == 1 .and. i == 8 .and. g == 1) then
write(6,*) '*** winding forward in IntegrateStress ***'
write(6,*) subFrac, subFrac+subStep
write (6,'(a,/,3(3(f12.7,x)/))') 'Lp of 1 8 1',crystallite_Lp(1:3,:,1,8,1)
write (6,'(a,/,3(3(f12.7,x)/))') 'Fp of 1 8 1',crystallite_Fp(1:3,:,1,8,1)
endif
Fg_current = Fg_aim ! wind forward
Fe_current = crystallite_Fe(:,:,g,i,e)
Fp_current = crystallite_Fp(:,:,g,i,e)
subFrac = subFrac + subStep
subStep = min(1.0_pReal-subFrac, subStep*2.0_pReal) ! accelerate
else ! time integration encountered trouble
subStep = 0.5_pReal * subStep ! cut time step in half
crystallite_Lp(:,:,g,i,e) = 0.5_pReal*(crystallite_Lp(:,:,g,i,e) + & ! interpolate Lp and L
(math_I3 - math_mul33x33(Fp_current,&
math_mul33x33(math_inv3x3(Fg_aim),Fe_current)))/dt_aim)
endif
enddo ! potential substepping
!$OMP CRITICAL (distributionStress)
debug_StressLoopDistribution(Niteration) = debug_StressLoopDistribution(Niteration) + 1
!$OMP END CRITICAL (distributionStress)
return ! with final happyness
end function
!***********************************************************************
!*** fully-implicit two-level time integration ***
!*** based on a residuum in Lp and intermediate ***
!*** acceleration of the Newton-Raphson correction ***
!***********************************************************************
subroutine TimeIntegration(&
happy,& ! return status
Lpguess,& ! guess of plastic velocity gradient
Fp_new,& ! new plastic deformation gradient
Fe_new,& ! new "elastic" deformation gradient
Tstar_v,& ! Stress vector
P,& ! 1st PK stress (taken as initial guess if /= 0)
!
Fg_new,& ! new total def gradient
Fp_old,& ! former plastic def gradient
Temperature,& ! temperature
state,& ! microstructural state
dt,& ! time increment
grain,& ! grain number
ip,& ! integration point number
cp_en & ! element number
)
use prec use prec
use debug use debug
use mesh, only: mesh_element use constitutive, only: constitutive_microstructure, &
use constitutive, only: constitutive_microstructure,constitutive_homogenizedC,constitutive_LpAndItsTangent,& constitutive_homogenizedC, &
constitutive_sizeState constitutive_LpAndItsTangent
use math use math
use IO
implicit none implicit none
logical, intent(out) :: happy !*** input variables ***!
real(pReal), dimension(3,3), intent(inout) :: Lpguess integer(pInt), intent(in) :: e, & ! element index
real(pReal), dimension(3,3), intent(out) :: Fp_new,Fe_new,P i, & ! integration point index
real(pReal), dimension(6), intent(out) :: Tstar_v g ! grain index
real(pReal), dimension(3,3), intent(in) :: Fg_new,Fp_old
real(pReal), intent(in) :: Temperature,dt
integer(pInt), intent(in) :: cp_en, ip, grain
real(pReal), dimension(constitutive_sizeState(grain,ip,cp_en)), intent(in) :: state
logical error !*** output variables ***!
integer(pInt) Niteration,dummy, i,j,k,l,m,n logical crystallite_integrateStress ! flag indicating if integration suceeded
integer(pLongInt) tick,tock,tickrate,maxticks
real(pReal) p_hydro,det, leapfrog,maxleap
real(pReal), dimension(3,3,3,3) :: C
real(pReal), dimension(9,9) :: dLp,dTdLp,dRdLp,invdRdLp,eye2
real(pReal), dimension(6,6) :: C_66
real(pReal), dimension(3,3) :: invFp_new,invFp_old,Lp,Lpguess_old,Rinner,Rinner_old,A,B,BT,AB,BTA
happy = .false. ! be pessimistic !*** internal variables ***!
eye2 = math_identity2nd(9) real(pReal), dimension(3,3) :: Fg_current, & ! deformation gradient at start of timestep
Fg_new, & ! deformation gradient at end of timestep
Fp_current, & ! plastic deformation gradient at start of timestep
Fp_new, & ! plastic deformation gradient at end of timestep
Fe_current, & ! elastic deformation gradient at start of timestep
Fe_new, & ! elastic deformation gradient at end of timestep
invFp_new, & ! inverse of Fp_new
invFp_current, & ! inverse of Fp_current
Lpguess, & ! current guess for plastic velocity gradient
Lpguess_old, & ! known last good guess for plastic velocity gradient
Lp_constitutive, & ! plastic velocity gradient resulting from constitutive law
residuum, & ! current residuum of plastic velocity gradient
residuum_old, & ! last residuum of plastic velocity gradient
A, &
B, &
BT, &
AB, &
BTA
real(pReal), dimension(6) :: Tstar_v ! 2nd Piola-Kirchhoff Stress in Mandel-Notation
real(pReal), dimension(9,9) :: dLp_constitutive, & ! partial derivative of plastic velocity gradient calculated by constitutive law
dTdLp, & ! partial derivative of 2nd Piola-Kirchhoff stress
dRdLp, & ! partial derivative of residuum (Jacobian for NEwton-Raphson scheme)
invdRdLp ! inverse of dRdLp
real(pReal), dimension(3,3,3,3) :: C ! 4th rank elasticity tensor
real(pReal), dimension(6,6) :: C_66 ! simplified 2nd rank elasticity tensor
real(pReal) p_hydro, & ! volumetric part of 2nd Piola-Kirchhoff Stress
det, & ! determinant
leapfrog, & ! acceleration factor for Newton-Raphson scheme
maxleap ! maximum acceleration factor
logical error ! flag indicating an error
integer(pInt) NiterationStress, & ! number of stress integrations
dummy, &
h, &
j, &
k, &
l, &
m, &
n
integer(pLongInt) tick, &
tock, &
tickrate, &
maxticks
invFp_old = math_inv3x3(Fp_old) ! inversion of Fp_old
if (all(invFp_old == 0.0_pReal)) return ! failed
! write (6,'(a,/,3(3(f12.7,x)/))') 'Fp old',Fp_old
! write (6,'(a,/,3(3(f12.7,x)/))') 'Fp old inv',invFp_old
A = math_mul33x33(transpose(invFp_old), math_mul33x33(transpose(Fg_new),math_mul33x33(Fg_new,invFp_old))) ! be pessimistic
crystallite_integrateStress = .false.
!$OMP CRITICAL (write2out) ! feed local variables
! if (debugger) write (6,'(a,/,3(3(f12.7,x)/))') 'Fg to be calculated',Fg_new Fg_current = crystallite_subF0(:,:,g,i,e)
!$OMP END CRITICAL (write2out) Fg_new = crystallite_subF(:,:,g,i,e)
Fp_current = crystallite_subFp0(:,:,g,i,e)
Fe_current = math_mul33x33(Fg_current,math_inv3x3(Fp_current))
Tstar_v = crystallite_Tstar_v(:,g,i,e)
Lpguess_old = crystallite_Lp(:,:,g,i,e) ! consider present Lp good (i.e. worth remembering) ...
Lpguess = crystallite_Lp(:,:,g,i,e) ! ... and take it as first guess
call constitutive_microstructure(Temperature,grain,ip,cp_en) ! inversion of Fp_current...
C_66 = constitutive_homogenizedC(grain,ip,cp_en) invFp_current = math_inv3x3(Fp_current)
C = math_Mandel66to3333(C_66) ! 4th rank elasticity tensor if (all(invFp_current == 0.0_pReal)) then ! ... failed?
if (debugger) write(6,*) '::: integrateStress failed on invFp_current inversion'
return
endif
Niteration = 0_pInt A = math_mul33x33(transpose(invFp_current), math_mul33x33(transpose(Fg_new),math_mul33x33(Fg_new,invFp_current)))
leapfrog = 1.0_pReal ! correction as suggested by invdRdLp-step
maxleap = 1024.0_pReal ! preassign maximum acceleration level
Lpguess_old = Lpguess ! consider present Lpguess good (i.e. worth remembering) ! update microstructure
call constitutive_microstructure(crystallite_Temperature(g,i,e),g,i,e)
Inner: do ! inner iteration: Lp ! get elasticity tensor
Niteration = Niteration + 1 C_66 = constitutive_homogenizedC(g,i,e)
if (Niteration > nLp) then ! too many loops required C = math_Mandel66to3333(C_66)
Lpguess = Lpguess_old ! do not trust the last update but resort to former one
! start LpLoop with no acceleration
NiterationStress = 0_pInt
leapfrog = 1.0_pReal
maxleap = 1024.0_pReal
LpLoop: do
! increase loop counter
NiterationStress = NiterationStress + 1
! too many loops required ?
if (NiterationStress > nStress) then
if (debugger) write(6,*) '::: integrateStress exceeded nStress loopcount'
return return
endif endif
! write(6,*) 'iteration',Niteration B = math_i3 - crystallite_subdt(g,i,e)*Lpguess
! write (6,'(a,/,3(3(f12.7,x)/))') 'Lpguess',Lpguess
B = math_i3 - dt*Lpguess
BT = transpose(B) BT = transpose(B)
AB = math_mul33x33(A,B) AB = math_mul33x33(A,B)
BTA = math_mul33x33(BT,A) BTA = math_mul33x33(BT,A)
! calculate 2nd Piola-Kirchhoff stress tensor
Tstar_v = 0.5_pReal*math_mul66x6(C_66,math_mandel33to6(math_mul33x33(BT,AB)-math_I3)) Tstar_v = 0.5_pReal*math_mul66x6(C_66,math_mandel33to6(math_mul33x33(BT,AB)-math_I3))
p_hydro = (Tstar_v(1) + Tstar_v(2) + Tstar_v(3))/3.0_pReal p_hydro = sum(Tstar_v(1:3))/3.0_pReal
forall(i=1:3) Tstar_v(i) = Tstar_v(i) - p_hydro ! subtract hydrostatic pressure forall(n=1:3) Tstar_v(n) = Tstar_v(n) - p_hydro ! get deviatoric stress tensor
! write (6,'(a,/,6(f12.7,x))') 'Tstar',Tstar_v
! calculate plastic velocity gradient and its tangent according to constitutive law
call system_clock(count=tick,count_rate=tickrate,count_max=maxticks) call system_clock(count=tick,count_rate=tickrate,count_max=maxticks)
call constitutive_LpAndItsTangent(Lp,dLp, Tstar_v,Temperature,grain,ip,cp_en) call constitutive_LpAndItsTangent(Lp_constitutive,dLp_constitutive,Tstar_v,crystallite_Temperature(g,i,e),g,i,e)
call system_clock(count=tock,count_rate=tickrate,count_max=maxticks) call system_clock(count=tock,count_rate=tickrate,count_max=maxticks)
debug_cumLpCalls = debug_cumLpCalls + 1_pInt debug_cumLpCalls = debug_cumLpCalls + 1_pInt
debug_cumLpTicks = debug_cumLpTicks + tock-tick debug_cumLpTicks = debug_cumLpTicks + tock-tick
if (tock < tick) debug_cumLpTicks = debug_cumLpTicks + maxticks if (tock < tick) debug_cumLpTicks = debug_cumLpTicks + maxticks
Rinner = Lpguess - Lp ! update current residuum ! update current residuum
! write (6,'(a,/,3(3(f12.7,x)/))') 'Lp',Lp residuum = Lpguess - Lp_constitutive
! write (6,'(a,/,3(3(f12.7,x)/))') 'Residuum',Rinner
if (.not.(any(Rinner/=Rinner)) .and. & ! exclude any NaN in residuum ! Check for convergence of loop
( ( maxval(abs(Rinner)) < aTol_crystalliteStress) .or. & ! below abs tol .or. if (.not.(any(residuum/=residuum)) .and. & ! exclude any NaN in residuum
( any(abs(dt*Lpguess) > relevantStrain) .and. & ! worth checking? .and. ( maxval(abs(residuum)) < aTol_crystalliteStress .or. & ! below absolute tolerance .or.
maxval(abs(Rinner/Lpguess),abs(dt*Lpguess) > relevantStrain) < rTol_crystalliteStress & ! below rel tol ( any(abs(crystallite_subdt(g,i,e)*Lpguess) > relevantStrain) .and. & ! worth checking? .and.
maxval(abs(residuum/Lpguess), abs(crystallite_subdt(g,i,e)*Lpguess) > relevantStrain) < rTol_crystalliteStress & ! below relative tolerance
) & ) &
) & ) &
) & ) &
exit Inner ! convergence exit LpLoop
!
! check for acceleration/deceleration in Newton--Raphson correction ! NaN occured at regular speed?
! if (any(residuum/=residuum) .and. leapfrog == 1.0) then
if (any(Rinner/=Rinner) .and. & ! NaN occured at regular speed if (debugger) write(6,*) '::: integrateStress encountered NaN at iteration', NiterationStress
leapfrog == 1.0) then
Lpguess = Lpguess_old ! restore known good guess and croak for cutback
return return
elseif (leapfrog > 1.0_pReal .and. & ! at fast pace ? ! something went wrong at accelerated speed?
(sum(Rinner*Rinner) > sum(Rinner_old*Rinner_old) .or. & ! worse residuum elseif (leapfrog > 1.0_pReal .and. & ! at fast pace .and.
sum(Rinner*Rinner_old) < 0.0_pReal) .or. & ! residuum changed sign (overshoot) ( sum(residuum*residuum) > sum(residuum_old*residuum_old) .or. & ! worse residuum .or.
any(Rinner/=Rinner) ) then ! NaN sum(residuum*residuum_old) < 0.0_pReal .or. & ! residuum changed sign (overshoot) .or.
maxleap = 0.5_pReal * leapfrog ! limit next acceleration any(residuum/=residuum) & ! NaN occured
leapfrog = 1.0_pReal ! grinding halt ) &
) then
maxleap = 0.5_pReal * leapfrog ! limit next acceleration
leapfrog = 1.0_pReal ! grinding halt
else ! better residuum ! restore old residuum and Lp
dTdLp = 0.0_pReal ! calc dT/dLp Lpguess = Lpguess_old
forall (i=1:3,j=1:3,k=1:3,l=1:3,m=1:3,n=1:3) & residuum = residuum_old
dTdLp(3*(i-1)+j,3*(k-1)+l) = dTdLp(3*(i-1)+j,3*(k-1)+l) + &
C(i,j,l,n)*AB(k,n)+C(i,j,m,l)*BTA(m,k) ! residuum got better
dTdLp = -0.5_pReal*dt*dTdLp else
dRdLp = eye2 - math_mul99x99(dLp,dTdLp) ! calc dR/dLp ! calculate Jacobian for correction term
dTdLp = 0.0_pReal
forall (h=1:3,j=1:3,k=1:3,l=1:3,m=1:3,n=1:3) &
dTdLp(3*(h-1)+j,3*(k-1)+l) = dTdLp(3*(h-1)+j,3*(k-1)+l) + &
C(h,j,l,n)*AB(k,n)+C(h,j,m,l)*BTA(m,k)
dTdLp = -0.5_pReal*crystallite_subdt(g,i,e)*dTdLp
dRdLp = math_identity2nd(9) - math_mul99x99(dLp_constitutive,dTdLp)
invdRdLp = 0.0_pReal invdRdLp = 0.0_pReal
call math_invert(9,dRdLp,invdRdLp,dummy,error) ! invert dR/dLp --> dLp/dR call math_invert(9,dRdLp,invdRdLp,dummy,error) ! invert dR/dLp --> dLp/dR
if (error) then if (error) then
if (debugger) then if (debugger) write(6,*) '::: integrateStress failed on dR/dLp inversion at iteration', NiterationStress
!$OMP CRITICAL (write2out)
write (6,*) 'inversion dR/dLp failed',grain,ip,cp_en
! write (6,'(a,/,9(9(e9.3,x)/))') 'dRdLp', dRdLp(1:9,:)
! write (6,'(a,/,3(4(f9.3,x)/))') 'state_new / MPa',state/1e6_pReal
write (6,'(a,/,3(3(f12.7,x)/))') 'Lpguess',Lpguess(1:3,:)
write (6,'(a,/,3(3(e12.7,x)/))') 'Lp',Lp(1:3,:)
write (6,'(a,/,6(f9.3,x))') 'Tstar / MPa',Tstar_v/1e6_pReal
!$OMP END CRITICAL (write2out)
endif
return return
endif endif
Rinner_old = Rinner ! remember current residuum ! remember current residuum and Lpguess
Lpguess_old = Lpguess ! remember current Lp guess residuum_old = residuum
if (Niteration > 1 .and. leapfrog < maxleap) leapfrog = 2.0_pReal * leapfrog ! accelerate if ok Lpguess_old = Lpguess
! accelerate?
if (NiterationStress > 1 .and. leapfrog < maxleap) leapfrog = 2.0_pReal * leapfrog
endif endif
Lpguess = Lpguess_old ! start from current guess ! leapfrog to updated Lp
Rinner = Rinner_old ! use current residuum forall (k=1:3,l=1:3,m=1:3,n=1:3) &
forall (i=1:3,j=1:3,k=1:3,l=1:3) & ! leapfrog to updated Lpguess Lpguess(k,l) = Lpguess(k,l) - leapfrog*invdRdLp(3*(k-1)+l,3*(m-1)+n)*residuum(m,n)
Lpguess(i,j) = Lpguess(i,j) - leapfrog*invdRdLp(3*(i-1)+j,3*(k-1)+l)*Rinner(k,l) enddo LpLoop
enddo Inner
!$OMP CRITICAL (distributionLp) ! calculate new plastic and elastic deformation gradient
debug_LpLoopDistribution(Niteration) = debug_LpLoopDistribution(Niteration) + 1 invFp_new = math_mul33x33(invFp_current,B)
!$OMP END CRITICAL (distributionLp) invFp_new = invFp_new/math_det3x3(invFp_new)**(1.0_pReal/3.0_pReal) ! regularize by det
invFp_new = math_mul33x33(invFp_old,B)
if (debugger) then
write (6,'(a,x,f10.6,/,3(3(f12.7,x)/))') 'Lp(guess)',dt,Lpguess(1:3,:)
write (6,'(a,/,3(3(f12.7,x)/))') 'invFp_old',invFp_old(1:3,:)
write (6,'(a,/,3(3(f12.7,x)/))') 'B',B(1:3,:)
write (6,'(a,/,3(3(f12.7,x)/))') 'invFp_new',invFp_new(1:3,:)
endif
call math_invert3x3(invFp_new,Fp_new,det,error) call math_invert3x3(invFp_new,Fp_new,det,error)
if (debugger) then if (error) then
write (6,'(a,/,3(3(f12.7,x)/))') 'invFp_new',invFp_new(1:3,:) if (debugger) write(6,*) '::: integrateStress failed on invFp_new inversion at iteration', NiterationStress
write (6,'(a,/,3(3(f12.7,x)/))') 'Fp_new',Fp_new(1:3,:) return
write (6,'(a,x,l,x,a,f10.6)') 'with inversion error:',error,'and determinant:',det
endif endif
if (error) return ! inversion failed Fe_new = math_mul33x33(Fg_new,invFp_new) ! calc resulting Fe
Fp_new = Fp_new*det**(1.0_pReal/3.0_pReal) ! regularize Fp by det = det(InvFp_new) !! ! add volumetric component to 2nd Piola-Kirchhoff stress
forall (i=1:3) Tstar_v(i) = Tstar_v(i) + p_hydro ! add hydrostatic component back forall (n=1:3) Tstar_v(n) = Tstar_v(n) + p_hydro
Fe_new = math_mul33x33(Fg_new,invFp_new) ! calc resulting Fe
P = math_mul33x33(Fe_new,math_mul33x33(math_Mandel6to33(Tstar_v),transpose(invFp_new))) ! first PK stress
happy = .true. ! now smile... ! calculate 1st Piola-Kirchhoff stress
crystallite_P(:,:,g,i,e) = math_mul33x33(Fe_new,math_mul33x33(math_Mandel6to33(Tstar_v),transpose(invFp_new)))
! store local values in global variables
crystallite_Lp(:,:,g,i,e) = Lpguess
crystallite_Tstar_v(:,g,i,e) = Tstar_v
crystallite_Fp(:,:,g,i,e) = Fp_new
crystallite_Fe(:,:,g,i,e) = Fe_new
! set return flag to true
crystallite_integrateStress = .true.
if (debugger) write(6,*) '::: integrateStress finished at iteration', NiterationStress
!$OMP CRITICAL (distributionStress)
debug_StressLoopDistribution(NiterationStress) = debug_StressLoopDistribution(NiterationStress) + 1
!$OMP END CRITICAL (distributionStress)
return return
end subroutine endfunction
@ -780,8 +814,9 @@ function crystallite_postResults(&
constitutive_postResults(Tstar_v,Temperature,dt,g,i,e) constitutive_postResults(Tstar_v,Temperature,dt,g,i,e)
return return
end function endfunction
END MODULE END MODULE
!############################################################## !##############################################################

48
trunk/debug.f90
View File

@ -6,10 +6,10 @@
implicit none implicit none
integer(pInt), dimension(nLp) :: debug_LpLoopDistribution = 0_pInt
integer(pInt), dimension(nStress) :: debug_StressLoopDistribution = 0_pInt integer(pInt), dimension(nStress) :: debug_StressLoopDistribution = 0_pInt
integer(pInt), dimension(nCryst) :: debug_StateLoopDistribution = 0_pInt integer(pInt), dimension(nState) :: debug_StateLoopDistribution = 0_pInt
integer(pInt), dimension(nCryst) :: debug_StiffnessStateLoopDistribution = 0_pInt integer(pInt), dimension(nState) :: debug_StiffnessStateLoopDistribution = 0_pInt
integer(pInt), dimension(nCryst) :: debug_CrystalliteLoopDistribution = 0_pInt
integer(pLongInt) :: debug_cumLpTicks = 0_pInt integer(pLongInt) :: debug_cumLpTicks = 0_pInt
integer(pLongInt) :: debug_cumDotStateTicks = 0_pInt integer(pLongInt) :: debug_cumDotStateTicks = 0_pInt
integer(pInt) :: debug_cumLpCalls = 0_pInt integer(pInt) :: debug_cumLpCalls = 0_pInt
@ -27,10 +27,10 @@ SUBROUTINE debug_reset()
use prec use prec
implicit none implicit none
debug_LpLoopDistribution = 0_pInt ! initialize debugging data debug_StressLoopDistribution = 0_pInt ! initialize debugging data
debug_StressLoopDistribution = 0_pInt
debug_StateLoopDistribution = 0_pInt debug_StateLoopDistribution = 0_pInt
debug_StiffnessStateLoopDistribution = 0_pInt debug_StiffnessStateLoopDistribution = 0_pInt
debug_CrystalliteLoopDistribution = 0_pInt
debug_cumLpTicks = 0_pInt debug_cumLpTicks = 0_pInt
debug_cumDotStateTicks = 0_pInt debug_cumDotStateTicks = 0_pInt
debug_cumLpCalls = 0_pInt debug_cumLpCalls = 0_pInt
@ -67,49 +67,49 @@ END SUBROUTINE
write(6,'(a33,x,i12)') 'total CPU ticks :',debug_cumDotStateTicks write(6,'(a33,x,i12)') 'total CPU ticks :',debug_cumDotStateTicks
endif endif
integral = 0_pInt
write(6,*)
write(6,*) 'distribution_LpLoop :'
do i=1,nLp
if (debug_LpLoopDistribution(i) /= 0) then
integral = integral + i*debug_LpLoopDistribution(i)
write(6,*) i,debug_LpLoopDistribution(i)
endif
enddo
write(6,*) 'total',sum(debug_LpLoopDistribution),integral
integral = 0_pInt integral = 0_pInt
write(6,*) write(6,*)
write(6,*) 'distribution_StressLoop :' write(6,*) 'distribution_StressLoop :'
do i=1,nStress do i=1,nStress
if (debug_StressLoopDistribution(i) /= 0) then if (debug_StressLoopDistribution(i) /= 0) then
integral = integral + i*debug_StressLoopDistribution(i) integral = integral + i*debug_StressLoopDistribution(i)
write(6,*) i,debug_StressLoopDistribution(i) write(6,'(i25,i10)') i,debug_StressLoopDistribution(i)
endif endif
enddo enddo
write(6,*) 'total',sum(debug_StressLoopDistribution),integral write(6,'(a15,i10,i10)') ' total',sum(debug_StressLoopDistribution),integral
integral = 0_pInt integral = 0_pInt
write(6,*) write(6,*)
write(6,*) 'distribution_StateLoop :' write(6,*) 'distribution_StateLoop :'
do i=1,nCryst do i=1,nState
if (debug_StateLoopDistribution(i) /= 0) then if (debug_StateLoopDistribution(i) /= 0) then
integral = integral + i*debug_StateLoopDistribution(i) integral = integral + i*debug_StateLoopDistribution(i)
write(6,*) i,debug_StateLoopDistribution(i) write(6,'(i25,i10)') i,debug_StateLoopDistribution(i)
endif endif
enddo enddo
write(6,*) 'total',sum(debug_StateLoopDistribution),integral write(6,'(a15,i10,i10)') ' total',sum(debug_StateLoopDistribution),integral
integral = 0_pInt integral = 0_pInt
write(6,*) write(6,*)
write(6,*) 'distribution_StiffnessStateLoop :' write(6,*) 'distribution_StiffnessStateLoop :'
do i=1,nCryst do i=1,nState
if (debug_StiffnessStateLoopDistribution(i) /= 0) then if (debug_StiffnessStateLoopDistribution(i) /= 0) then
integral = integral + i*debug_StiffnessStateLoopDistribution(i) integral = integral + i*debug_StiffnessStateLoopDistribution(i)
write(6,*) i,debug_StiffnessStateLoopDistribution(i) write(6,'(i25,i10)') i,debug_StiffnessStateLoopDistribution(i)
endif endif
enddo enddo
write(6,*) 'total',sum(debug_StiffnessStateLoopDistribution),integral write(6,'(a15,i10,i10)') ' total',sum(debug_StiffnessStateLoopDistribution),integral
integral = 0_pInt
write(6,*)
write(6,*) 'distribution_CrystalliteLoop :'
do i=1,nCryst
if (debug_CrystalliteLoopDistribution(i) /= 0) then
integral = integral + i*debug_CrystalliteLoopDistribution(i)
write(6,'(i25,i10)') i,debug_CrystalliteLoopDistribution(i)
endif
enddo
write(6,'(a15,i10,i10)') ' total',sum(debug_CrystalliteLoopDistribution),integral
write(6,*) write(6,*)
END SUBROUTINE END SUBROUTINE

11
trunk/homogenization.f90
View File

@ -257,17 +257,6 @@ subroutine materialpoint_stressAndItsTangent(&
enddo enddo
!$OMP END PARALLEL DO !$OMP END PARALLEL DO
! tell what is requested
write(6,'(a14,x,16(l,x))') 'matpnt_req',materialpoint_requested
write(6,'(a14,x,16(l,x))') 'matpnt_don',materialpoint_doneAndHappy(1,:,:)
write(6,'(a14,x,16(l,x))') 'matpnt_hpy',materialpoint_doneAndHappy(1,:,:)
write(6,'(a14,x,16(f6.4,x))') 'matpnt_frac',materialpoint_subFrac
write(6,'(a14,x,16(f6.4,x))') 'matpnt_step',materialpoint_subStep
write(6,'(a10,x,16(e8.3,x))') 'matpnt_dt',materialpoint_subdt
write(6,*)
write (6,'(a,/,3(3(f12.7,x)/))') 'subF0 of 8 1',materialpoint_subF0(1:3,:,8,1)
write (6,'(a,/,3(3(f12.7,x)/))') 'subF of 8 1',materialpoint_subF(1:3,:,8,1)
! ------ convergence loop material point homogenization ------ ! ------ convergence loop material point homogenization ------
homogenization_Niteration = 0_pInt homogenization_Niteration = 0_pInt

6
trunk/prec.f90
View File

@ -23,9 +23,9 @@
real(pReal), parameter :: pert_Fg = 1.0e-6_pReal ! strain perturbation for FEM Jacobi real(pReal), parameter :: pert_Fg = 1.0e-6_pReal ! strain perturbation for FEM Jacobi
integer(pInt), parameter :: nReg = 1_pInt ! regularization attempts for Jacobi inversion integer(pInt), parameter :: nReg = 1_pInt ! regularization attempts for Jacobi inversion
integer(pInt), parameter :: nHomog = 10_pInt ! homogenization loop limit integer(pInt), parameter :: nHomog = 10_pInt ! homogenization loop limit
integer(pInt), parameter :: nCryst = 10_pInt ! crystallite loop limit (state update) integer(pInt), parameter :: nCryst = 20_pInt ! crystallite loop limit (only for debugging info, real loop limit is "subStepMin")
integer(pInt), parameter :: nStress = 20_pInt ! stress loop limit integer(pInt), parameter :: nState = 10_pInt ! state loop limit
integer(pInt), parameter :: nLp = 40_pInt ! stress loop limit integer(pInt), parameter :: nStress = 40_pInt ! stress loop limit
real(pReal), parameter :: rTol_crystalliteState = 1.0e-5_pReal ! relative tolerance in crystallite state loop real(pReal), parameter :: rTol_crystalliteState = 1.0e-5_pReal ! relative tolerance in crystallite state loop
real(pReal), parameter :: rTol_crystalliteStress = 1.0e-6_pReal ! relative tolerance in crystallite stress loop real(pReal), parameter :: rTol_crystalliteStress = 1.0e-6_pReal ! relative tolerance in crystallite stress loop
real(pReal), parameter :: aTol_crystalliteStress = 1.0e-8_pReal ! absolute tolerance in crystallite stress loop real(pReal), parameter :: aTol_crystalliteStress = 1.0e-8_pReal ! absolute tolerance in crystallite stress loop