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

crystallite: 

added full check whether residuum in inner loop is NaN
SingleCrystallite now takes/stores Lp guess

CPFEM_Taylor (missing in GIA!!):
former fully plastic guess for Lp has been exchanged by keeping the last converged Lp (global array CPFEM_Lp) to serve as new best guess for the next time step. This speeds up the inner loop of TimeIntegration.
This commit is contained in:
Philip Eisenlohr 2008-04-16 17:00:28 +00:00
parent fa2d6b9b6d
commit cdb2dd8808
3 changed files with 56 additions and 80 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

30
trunk/CPFEM_Taylor.f90
View File

@ -18,6 +18,7 @@
real(pReal), dimension (:,:,:,:), allocatable :: CPFEM_jaco_bar
real(pReal), dimension (:,:,:,:), allocatable :: CPFEM_jaco_knownGood
real(pReal), dimension (:,:,:,:), allocatable :: CPFEM_results
real(pReal), dimension (:,:,:,:,:), allocatable :: CPFEM_Lp
real(pReal), dimension (:,:,:,:,:), allocatable :: CPFEM_Fp_old
real(pReal), dimension (:,:,:,:,:), allocatable :: CPFEM_Fp_new
real(pReal), parameter :: CPFEM_odd_stress = 1e15_pReal, CPFEM_odd_jacobian = 1e50_pReal
@ -58,6 +59,9 @@
allocate(CPFEM_results(CPFEM_Nresults+constitutive_maxNresults,constitutive_maxNgrains,mesh_maxNips,mesh_NcpElems))
CPFEM_results = 0.0_pReal
!
! *** Plastic velocity gradient ***
allocate(CPFEM_Lp(3,3,constitutive_maxNgrains,mesh_maxNips,mesh_NcpElems)) ; CPFEM_Lp = 0.0_pReal
! *** Plastic deformation gradient at (t=t0) and (t=t1) ***
allocate(CPFEM_Fp_new(3,3,constitutive_maxNgrains,mesh_maxNips,mesh_NcpElems)) ; CPFEM_Fp_new = 0.0_pReal
allocate(CPFEM_Fp_old(3,3,constitutive_maxNgrains,mesh_maxNips,mesh_NcpElems))
@ -77,6 +81,7 @@
write(6,*) 'CPFEM_jaco_bar: ', shape(CPFEM_jaco_bar)
write(6,*) 'CPFEM_jaco_knownGood: ', shape(CPFEM_jaco_knownGood)
write(6,*) 'CPFEM_results: ', shape(CPFEM_results)
write(6,*) 'CPFEM_Lp: ', shape(CPFEM_Lp)
write(6,*) 'CPFEM_Fp_old: ', shape(CPFEM_Fp_old)
write(6,*) 'CPFEM_Fp_new: ', shape(CPFEM_Fp_new)
write(6,*)
@ -140,19 +145,16 @@
!
cp_en = mesh_FEasCP('elem',CPFEM_en)
if (cp_en == 1 .and. CPFEM_in == 1) &
write(6,'(a6,x,i4,x,a4,x,i4,x,a10,x,f8.4,x,a10,x,i2,x,a10,x,i2,x,a10,x,i2,x,a10,x,i2)') &
'elem',cp_en,'IP',CPFEM_in,&
write(6,'(a10,x,f8.4,x,a10,x,i2,x,a10,x,i2,x,a10,x,i2,x,a10,x,i2)') &
'theTime',theTime,'theInc',theInc,'theCycle',theCycle,'theLovl',theLovl,&
'mode',CPFEM_mode
!
select case (CPFEM_mode)
case (2,1) ! regular computation (with aging of results)
if (.not. CPFEM_calc_done) then ! puuh, me needs doing all the work...
write (6,*) 'puuh me needs doing all the work', cp_en
if (CPFEM_mode == 1) then ! age results at start of new increment
CPFEM_Fp_old = CPFEM_Fp_new
constitutive_state_old = constitutive_state_new
write (6,*) '#### aged results'
endif
debug_cutbackDistribution = 0_pInt ! initialize debugging data
debug_InnerLoopDistribution = 0_pInt
@ -199,8 +201,6 @@
! return the local stress and the jacobian from storage
CPFEM_stress(1:CPFEM_ngens) = CPFEM_stress_bar(1:CPFEM_ngens,CPFEM_in,cp_en)
CPFEM_jaco(1:CPFEM_ngens,1:CPFEM_ngens) = CPFEM_jaco_bar(1:CPFEM_ngens,1:CPFEM_ngens,CPFEM_in,cp_en)
! if (cp_en == 1 .and. CPFEM_in == 1) write (6,*) 'stress',CPFEM_stress
! if (cp_en == 1 .and. CPFEM_in == 1 .and. CPFEM_updateJaco) write (6,*) 'stiffness',CPFEM_jaco
!
return
!
@ -235,26 +235,26 @@
real(pReal), dimension(3,3,3,3) :: dPdF,dPdF_bar_old
!
CPFEM_PK1_bar(:,:,CPFEM_in,cp_en) = 0.0_pReal ! zero out average first PK stress
if (updateJaco) then
if (updateJaco) then
dPdF_bar_old = CPFEM_dPdF_bar(:,:,:,:,CPFEM_in,cp_en) ! remember former average consistent tangent
CPFEM_dPdF_bar(:,:,:,:,CPFEM_in,cp_en) = 0.0_pReal ! zero out avg consistent tangent for later assembly
endif
CPFEM_dPdF_bar(:,:,:,:,CPFEM_in,cp_en) = 0.0_pReal ! zero out avg consistent tangent for later assembly
endif
do grain = 1,texture_Ngrains(mesh_element(4,cp_en))
do grain = 1,texture_Ngrains(mesh_element(4,cp_en))
dPdF = dPdF_bar_old ! preguess consistent tangent of grain with avg
call SingleCrystallite(msg,PK1,dPdF,&
CPFEM_results(5:4+constitutive_Nresults(grain,CPFEM_in,cp_en),grain,CPFEM_in,cp_en),&
CPFEM_Lp(:,:,grain,CPFEM_in,cp_en),&
CPFEM_Fp_new(:,:,grain,CPFEM_in,cp_en),Fe1,constitutive_state_new(:,grain,CPFEM_in,cp_en),& ! output up to here
CPFEM_dt,cp_en,CPFEM_in,grain,.true.,&
CPFEM_Temperature(CPFEM_in,cp_en),&
CPFEM_ffn1_bar(:,:,CPFEM_in,cp_en),CPFEM_ffn_bar(:,:,CPFEM_in,cp_en),&
CPFEM_Fp_old(:,:,grain,CPFEM_in,cp_en),constitutive_state_old(:,grain,CPFEM_in,cp_en))
volfrac = constitutive_matVolFrac(grain,CPFEM_in,cp_en)*constitutive_texVolFrac(grain,CPFEM_in,cp_en)
volfrac = constitutive_matVolFrac(grain,CPFEM_in,cp_en)*constitutive_texVolFrac(grain,CPFEM_in,cp_en)
CPFEM_PK1_bar(:,:,CPFEM_in,cp_en) = CPFEM_PK1_bar(:,:,CPFEM_in,cp_en) + volfrac*PK1
if (updateJaco) CPFEM_dPdF_bar(:,:,:,:,CPFEM_in,cp_en) = &
CPFEM_dPdF_bar(:,:,:,:,CPFEM_in,cp_en) + volfrac*dPdF ! add up crystallite stiffnesses
! (may have "holes" corresponding
if (updateJaco) CPFEM_dPdF_bar(:,:,:,:,CPFEM_in,cp_en) = &
CPFEM_dPdF_bar(:,:,:,:,CPFEM_in,cp_en) + volfrac*dPdF ! add up crystallite stiffnesses
! (may have "holes" corresponding
! to former avg tangent)
!
! update results plotted in MENTAT

100
trunk/crystallite.f90
View File

@ -19,6 +19,7 @@ CONTAINS
P,& ! first PK stress
dPdF,& ! consistent tangent
post_results,& ! plot results from constitutive model
Lp,& ! plastic velocity gradient
Fp_new,& ! new plastic deformation gradient
Fe_new,& ! new "elastic" deformation gradient
state_new,& ! new state variable array
@ -48,13 +49,13 @@ CONTAINS
real(pReal) Temperature
real(pReal) dt,dt_aim,subFrac,subStep,invJ,det
real(pReal), dimension(3,3) :: Lp,Lp_pert,inv
real(pReal), dimension(3,3) :: Fg_old,Fg_current,Fg_aim,Fg_new,Fg_pert,deltaFg
real(pReal), dimension(3,3) :: Fg_old,Fg_current,Fg_new,Fg_pert,Fg_aim,deltaFg
real(pReal), dimension(3,3) :: Fp_old,Fp_current,Fp_new,Fp_pert
real(pReal), dimension(3,3) :: Fe_old,Fe_current,Fe_new,Fe_pert
real(pReal), dimension(3,3) :: Tstar,tau,P,P_pert
real(pReal), dimension(3,3,3,3) :: dPdF
real(pReal), dimension(constitutive_Nstatevars(grain,ip,cp_en)) :: state_old,state_new
real(pReal), dimension(constitutive_Nstatevars(grain,ip,cp_en)) :: state_current,state_bestguess,state_pert
real(pReal), dimension(constitutive_Nstatevars(grain,ip,cp_en)) :: state_old,state_new
real(pReal), dimension(constitutive_Nstatevars(grain,ip,cp_en)) :: state_current,state_bestguess,state_pert
real(pReal), dimension(constitutive_Nresults(grain,ip,cp_en)) :: post_results
!
deltaFg = Fg_new-Fg_old
@ -65,7 +66,7 @@ CONTAINS
Fg_aim = Fg_old ! make "new", "aim" a synonym for "old"
Fp_new = Fp_old
call math_invert3x3(Fp_old,inv,det,error)
Fe_new = matmul(Fg_old,inv)
Fe_new = matmul(Fg_old,inv)
state_bestguess = state_new ! remember potentially available state guess
state_new = state_old
!
@ -78,23 +79,17 @@ CONTAINS
Fg_current = Fg_aim ! wind forward
Fp_current = Fp_new
Fe_current = Fe_new
state_current = state_new
state_current = state_new
endif
!
Fg_aim = Fg_current + subStep*deltaFg ! aim for Fg
dt_aim = subStep*dt ! aim for dt
msg = '' ! error free so far
if (guessNew) then ! calculate new Lp guess when cutted back
if (dt_aim /= 0.0_pReal) then
call math_invert3x3(Fg_aim,inv,det,error)
Lp = (math_I3-matmul(Fp_current,matmul(inv,Fe_current)))/dt ! fully plastic initial guess
else
Lp = 0.0_pReal ! fully elastic guess
endif
state_new = state_bestguess ! try best available guess for state
endif
if (guessNew) &
state_new = state_bestguess ! try best available guess for state
call TimeIntegration(msg,Lp,Fp_new,Fe_new,P,state_new,post_results,.true., & ! def gradients and PK2 at end of time step
dt_aim,cp_en,ip,grain,Temperature,Fg_aim,Fp_current,state_current)
dt_aim,cp_en,ip,grain,Temperature,Fg_aim,Fp_current,state_current,0)
!
if (msg == 'ok') then
cuttedBack = .false. ! no cut back required
@ -105,8 +100,8 @@ CONTAINS
nCutbacks = nCutbacks + 1 ! record additional cutback
cuttedBack = .true. ! encountered problems -->
guessNew = .true. ! redo plastic Lp guess
subStep = subStep / 2.0_pReal ! cut time step in half
state_bestguess = state_current ! current state is then best guess
subStep = subStep / 2.0_pReal ! cut time step in half
state_bestguess = state_current ! current state is then best guess
endif
enddo ! potential substepping
!
@ -120,10 +115,11 @@ CONTAINS
Fg_pert(k,l) = Fg_pert(k,l) + pert_Fg ! perturb single component
Lp_pert = Lp
state_pert = state_new ! initial guess from end of time step
call TimeIntegration(msg,Lp,Fp_pert,Fe_pert,P_pert,state_pert,post_results,.false., & ! def gradients and PK2 at end of time step
dt_aim,cp_en,ip,grain,Temperature,Fg_pert,Fp_current,state_current)
call TimeIntegration(msg,Lp_pert,Fp_pert,Fe_pert,P_pert,state_pert,post_results,.false., & ! def gradients and PK2 at end of time step
dt_aim,cp_en,ip,grain,Temperature,Fg_pert,Fp_current,state_current,1)
if (msg == 'ok') &
dPdF(:,:,k,l) = (P_pert-P)/pert_Fg ! constructing tangent dP_ij/dFg_kl only if valid forward difference
! otherwise leave component unchanged
enddo
enddo
endif
@ -154,8 +150,8 @@ CONTAINS
Temperature,& ! temperature
Fg_new,& ! new total def gradient
Fp_old,& ! former plastic def gradient
state_old) ! former microstructure
!
state_old,& ! former microstructure
flag)
use prec
use debug
use mesh, only: mesh_element
@ -168,7 +164,7 @@ CONTAINS
character(len=*) msg
logical failed,wantsConstitutiveResults
integer(pInt) cp_en, ip, grain
integer(pInt) iOuter,iInner,dummy, i,j,k,l,m,n
integer(pInt) iOuter,iInner,dummy, i,j,k,l,m,n,flag
real(pReal) dt, Temperature, det, p_hydro, leapfrog,maxleap
real(pReal), dimension(6) :: Tstar_v
real(pReal), dimension(9,9) :: dLp,dTdLp,dRdLp,invdRdLp,eye2
@ -216,14 +212,7 @@ Inner: do ! inner iteration: Lp
msg = 'limit Inner iteration'
debug_InnerLoopDistribution(nInner) = debug_InnerLoopDistribution(nInner)+1
return
endif
! if (debugger) then
! write (6,*) iInner,iOuter
!write (6,'(a,/,9(9(e9.3,x)/))') 'dRdLp', dRdLp(1:9,:)
! write (6,'(a,/,3(3(e9.3,x)/))') 'Lpguess',Lpguess(1:3,:)
! write (6,*) 'state',state
! endif
endif
B = math_i3 - dt*Lpguess
BT = transpose(B)
@ -232,33 +221,28 @@ Inner: do ! inner iteration: Lp
Tstar_v = 0.5_pReal*matmul(C_66,math_mandel33to6(matmul(BT,AB)-math_I3))
Tstar = math_Mandel6to33(Tstar_v)
p_hydro=(Tstar_v(1)+Tstar_v(2)+Tstar_v(3))/3.0_pReal
forall(i=1:3) Tstar_v(i) = Tstar_v(i)-p_hydro ! subtract hydrostatic pressure
forall(i=1:3) Tstar_v(i) = Tstar_v(i)-p_hydro ! subtract hydrostatic pressure
call constitutive_LpAndItsTangent(Lp,dLp, &
Tstar_v,state,Temperature,grain,ip,cp_en)
! if (debugger) then
! write (6,'(a,/,6(f9.3,x))') 'Tstar[MPa]',1.0e-6_pReal*Tstar_v
! write (6,'(a,/,3(3(e9.3,x)/))') 'Lp',Lp(1:3,:)
! write (6,'(a,/,9(9(e9.3,x)/))') 'dLp', dLp(1:9,:)
! endif
Tstar_v,state,Temperature,grain,ip,cp_en)
Rinner = Lpguess - Lp ! update current residuum
! if (debugger) then
! write (6,'(a,/,3(3(e9.3,x)/))') 'Rinner',Rinner(1:3,:)
! endif
Rinner = Lpguess - Lp ! update current residuum
if ((maxval(abs(Rinner)) < abstol_Inner) .or. &
(any(abs(dt*Lpguess) > relevantStrain) .and. &
maxval(abs(Rinner/Lpguess),abs(dt*Lpguess) > relevantStrain) < reltol_Inner)) &
exit Inner
if (not(any(Rinner.NE.Rinner)) .and. & ! exclude all NaN in residuum
( (maxval(abs(Rinner)) < abstol_Inner) .or. & ! below abs tol .or.
( any(abs(dt*Lpguess) > relevantStrain) .and. & ! worth checking? .and.
maxval(abs(Rinner/Lpguess),abs(dt*Lpguess) > relevantStrain) < reltol_Inner & ! below rel tol
) &
) &
) &
exit Inner ! convergence
!
! check for acceleration/deceleration in Newton--Raphson correction
if (leapfrog > 1.0_pReal .and. &
(sum(Rinner*Rinner) > sum(Rinner_old*Rinner_old) .or. & ! worse residuum
sum(Rinner*Rinner_old) < 0.0_pReal) .or. & ! residuum changed sign (overshoot)
any(Rinner /= Rinner) ) then ! NaN
sum(Rinner*Rinner_old) < 0.0_pReal) .or. & ! residuum changed sign (overshoot)
any(Rinner.NE.Rinner) ) then ! NaN
maxleap = 0.5_pReal * leapfrog ! limit next acceleration
leapfrog = 1.0_pReal ! grinding halt
leapfrog = 1.0_pReal ! grinding halt
! if (debugger) write(6,*) 'grinding HALT'
else ! better residuum
dTdLp = 0.0_pReal ! calc dT/dLp
@ -271,20 +255,12 @@ Inner: do ! inner iteration: Lp
call math_invert(9,dRdLp,invdRdLp,dummy,failed) ! invert dR/dLp --> dLp/dR
if (failed) then
msg = 'inversion dR/dLp'
! if (debugger) then
! write (6,*) msg
! write (6,'(a,/,9(9(e9.3,x)/))') 'dRdLp', dRdLp(1:9,:)
! write (6,*) 'state',state
! write (6,'(a,/,3(3(f12.7,x)/))') 'Lpguess',Lpguess(1:3,:)
! write (6,*) 'Tstar',Tstar_v
! endif
return
endif
!
Rinner_old = Rinner ! remember current residuum
Lpguess_old = Lpguess ! remember current Lp guess
if (iInner > 1 .and. leapfrog < maxleap) leapfrog = 2.0_pReal * leapfrog ! accelerate
! if (debugger) write(6,*) 'leapfrogging at',leapfrog
if (iInner > 1 .and. leapfrog < maxleap) leapfrog = 2.0_pReal * leapfrog ! accelerate
endif
!
Lpguess = Lpguess_old ! start from current guess
@ -296,9 +272,9 @@ Inner: do ! inner iteration: Lp
debug_InnerLoopDistribution(iInner) = debug_InnerLoopDistribution(iInner)+1
ROuter = state - state_old - &
dt*constitutive_dotState(Tstar_v,state,Temperature,&
grain,ip,cp_en) ! residuum from evolution of microstructure
state = state - ROuter ! update of microstructure
if (maxval(abs(Router/state),state /= 0.0_pReal) < reltol_Outer) exit Outer
grain,ip,cp_en) ! residuum from evolution of microstructure
state = state - ROuter ! update of microstructure
if (maxval(abs(Router/state),state /= 0.0_pReal) < reltol_Outer) exit Outer ! convergence ?
enddo Outer
!
debug_OuterLoopDistribution(iOuter) = debug_OuterLoopDistribution(iOuter)+1

6
trunk/prec.f90
View File

@ -17,9 +17,9 @@
integer(pInt), parameter :: nCutback = 20_pInt ! max cutbacks accounted for in debug distribution
integer(pInt), parameter :: nReg = 1_pInt ! regularization attempts for Jacobi inversion
real(pReal), parameter :: pert_Fg = 1.0e-6_pReal ! strain perturbation for FEM Jacobi
integer(pInt), parameter :: nOuter = 20_pInt ! outer loop limit
integer(pInt), parameter :: nInner = 200_pInt ! inner loop limit
real(pReal), parameter :: reltol_Outer = 1.0e-4_pReal ! relative tolerance in outer loop (state)
integer(pInt), parameter :: nOuter = 20_pInt ! outer loop limit 20
integer(pInt), parameter :: nInner = 200_pInt ! inner loop limit 200
real(pReal), parameter :: reltol_Outer = 1.0e-6_pReal ! relative tolerance in outer loop (state)
real(pReal), parameter :: reltol_Inner = 1.0e-6_pReal ! relative tolerance in inner loop (Lp)
real(pReal), parameter :: abstol_Inner = 1.0e-8_pReal ! absolute tolerance in inner loop (Lp)
!