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

Make use of average consistent tangent dPdF_bar to fill eventual holes in the individual dPdF of one orientation (during consistent tangent calculation)

This commit is contained in:
Luc Hantcherli 2008-04-11 13:04:25 +00:00
parent 664a29c183
commit fa2d6b9b6d
2 changed files with 53 additions and 24 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

21
trunk/CPFEM_Taylor.f90
View File

@ -232,13 +232,16 @@
real(pReal), dimension(3,3) :: PK1_pert,F1_pert real(pReal), dimension(3,3) :: PK1_pert,F1_pert
real(pReal), dimension(3,3) :: U,R,Fe1 real(pReal), dimension(3,3) :: U,R,Fe1
real(pReal), dimension(3,3) :: PK1 real(pReal), dimension(3,3) :: PK1
real(pReal), dimension(3,3,3,3) :: dPdF 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 CPFEM_PK1_bar(:,:,CPFEM_in,cp_en) = 0.0_pReal ! zero out average first PK stress
if (updateJaco) CPFEM_dPdF_bar(:,:,:,:,CPFEM_in,cp_en) = 0.0_pReal ! zero out average consistent tangent 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
do grain = 1,texture_Ngrains(mesh_element(4,cp_en)) do grain = 1,texture_Ngrains(mesh_element(4,cp_en))
volfrac = constitutive_matVolFrac(grain,CPFEM_in,cp_en)*constitutive_texVolFrac(grain,CPFEM_in,cp_en) dPdF = dPdF_bar_old ! preguess consistent tangent of grain with avg
call SingleCrystallite(msg,PK1,dPdF,& call SingleCrystallite(msg,PK1,dPdF,&
CPFEM_results(5:4+constitutive_Nresults(grain,CPFEM_in,cp_en),grain,CPFEM_in,cp_en),& CPFEM_results(5:4+constitutive_Nresults(grain,CPFEM_in,cp_en),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_Fp_new(:,:,grain,CPFEM_in,cp_en),Fe1,constitutive_state_new(:,grain,CPFEM_in,cp_en),& ! output up to here
@ -246,9 +249,13 @@
CPFEM_Temperature(CPFEM_in,cp_en),& CPFEM_Temperature(CPFEM_in,cp_en),&
CPFEM_ffn1_bar(:,:,CPFEM_in,cp_en),CPFEM_ffn_bar(:,:,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)) 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)
CPFEM_PK1_bar(:,:,CPFEM_in,cp_en) = CPFEM_PK1_bar(:,:,CPFEM_in,cp_en) + volfrac*PK1 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 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 ! update results plotted in MENTAT
call math_pDecomposition(Fe1,U,R,error) ! polar decomposition call math_pDecomposition(Fe1,U,R,error) ! polar decomposition

56
trunk/crystallite.f90
View File

@ -53,7 +53,8 @@ CONTAINS
real(pReal), dimension(3,3) :: Fe_old,Fe_current,Fe_new,Fe_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) :: Tstar,tau,P,P_pert
real(pReal), dimension(3,3,3,3) :: dPdF real(pReal), dimension(3,3,3,3) :: dPdF
real(pReal), dimension(constitutive_Nstatevars(grain,ip,cp_en)) :: state_old,state_current,state_new,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 real(pReal), dimension(constitutive_Nresults(grain,ip,cp_en)) :: post_results
! !
deltaFg = Fg_new-Fg_old deltaFg = Fg_new-Fg_old
@ -64,7 +65,8 @@ CONTAINS
Fg_aim = Fg_old ! make "new", "aim" a synonym for "old" Fg_aim = Fg_old ! make "new", "aim" a synonym for "old"
Fp_new = Fp_old Fp_new = Fp_old
call math_invert3x3(Fp_old,inv,det,error) 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 state_new = state_old
! !
cuttedBack = .false. cuttedBack = .false.
@ -76,7 +78,7 @@ CONTAINS
Fg_current = Fg_aim ! wind forward Fg_current = Fg_aim ! wind forward
Fp_current = Fp_new Fp_current = Fp_new
Fe_current = Fe_new Fe_current = Fe_new
state_current = state_new state_current = state_new
endif endif
! !
Fg_aim = Fg_current + subStep*deltaFg ! aim for Fg Fg_aim = Fg_current + subStep*deltaFg ! aim for Fg
@ -85,10 +87,11 @@ CONTAINS
if (guessNew) then ! calculate new Lp guess when cutted back if (guessNew) then ! calculate new Lp guess when cutted back
if (dt_aim /= 0.0_pReal) then if (dt_aim /= 0.0_pReal) then
call math_invert3x3(Fg_aim,inv,det,error) call math_invert3x3(Fg_aim,inv,det,error)
Lp = (math_I3-matmul(Fp_current,matmul(inv,Fe_current)))/dt ! fully plastic initial guess Lp = (math_I3-matmul(Fp_current,matmul(inv,Fe_current)))/dt ! fully plastic initial guess
else else
Lp = 0.0_pReal ! fully elastic guess Lp = 0.0_pReal ! fully elastic guess
endif endif
state_new = state_bestguess ! try best available guess for state
endif endif
call TimeIntegration(msg,Lp,Fp_new,Fe_new,P,state_new,post_results,.true., & ! def gradients and PK2 at end of time step 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)
@ -102,7 +105,8 @@ CONTAINS
nCutbacks = nCutbacks + 1 ! record additional cutback nCutbacks = nCutbacks + 1 ! record additional cutback
cuttedBack = .true. ! encountered problems --> cuttedBack = .true. ! encountered problems -->
guessNew = .true. ! redo plastic Lp guess guessNew = .true. ! redo plastic Lp guess
subStep = subStep / 2.0_pReal ! cut time step in half subStep = subStep / 2.0_pReal ! cut time step in half
state_bestguess = state_current ! current state is then best guess
endif endif
enddo ! potential substepping enddo ! potential substepping
! !
@ -118,12 +122,8 @@ CONTAINS
state_pert = state_new ! initial guess from end of time step 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 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) dt_aim,cp_en,ip,grain,Temperature,Fg_pert,Fp_current,state_current)
if (msg /= 'ok') then if (msg == 'ok') &
write(6,*)'killed myself at component',k,l dPdF(:,:,k,l) = (P_pert-P)/pert_Fg ! constructing tangent dP_ij/dFg_kl only if valid forward difference
msg = 'consistent tangent --> '//msg
return
endif
dPdF(:,:,k,l) = (P_pert-P)/pert_Fg ! constructing the tangent dP_ij/dFg_kl from forward differences
enddo enddo
enddo enddo
endif endif
@ -216,7 +216,15 @@ Inner: do ! inner iteration: Lp
msg = 'limit Inner iteration' msg = 'limit Inner iteration'
debug_InnerLoopDistribution(nInner) = debug_InnerLoopDistribution(nInner)+1 debug_InnerLoopDistribution(nInner) = debug_InnerLoopDistribution(nInner)+1
return return
endif 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
B = math_i3 - dt*Lpguess B = math_i3 - dt*Lpguess
BT = transpose(B) BT = transpose(B)
AB = matmul(A,B) AB = matmul(A,B)
@ -226,8 +234,19 @@ Inner: do ! inner iteration: Lp
p_hydro=(Tstar_v(1)+Tstar_v(2)+Tstar_v(3))/3.0_pReal 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, & call constitutive_LpAndItsTangent(Lp,dLp, &
Tstar_v,state,Temperature,grain,ip,cp_en) Tstar_v,state,Temperature,grain,ip,cp_en)
Rinner = Lpguess - Lp ! update current residuum ! 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
Rinner = Lpguess - Lp ! update current residuum
! if (debugger) then
! write (6,'(a,/,3(3(e9.3,x)/))') 'Rinner',Rinner(1:3,:)
! endif
if ((maxval(abs(Rinner)) < abstol_Inner) .or. & if ((maxval(abs(Rinner)) < abstol_Inner) .or. &
(any(abs(dt*Lpguess) > relevantStrain) .and. & (any(abs(dt*Lpguess) > relevantStrain) .and. &
maxval(abs(Rinner/Lpguess),abs(dt*Lpguess) > relevantStrain) < reltol_Inner)) & maxval(abs(Rinner/Lpguess),abs(dt*Lpguess) > relevantStrain) < reltol_Inner)) &
@ -236,9 +255,11 @@ Inner: do ! inner iteration: Lp
! check for acceleration/deceleration in Newton--Raphson correction ! check for acceleration/deceleration in Newton--Raphson correction
if (leapfrog > 1.0_pReal .and. & if (leapfrog > 1.0_pReal .and. &
(sum(Rinner*Rinner) > sum(Rinner_old*Rinner_old) .or. & ! worse residuum (sum(Rinner*Rinner) > sum(Rinner_old*Rinner_old) .or. & ! worse residuum
sum(Rinner*Rinner_old) < 0.0_pReal)) then ! residuum changed sign (overshoot) sum(Rinner*Rinner_old) < 0.0_pReal) .or. & ! residuum changed sign (overshoot)
any(Rinner /= Rinner) ) then ! NaN
maxleap = 0.5_pReal * leapfrog ! limit next acceleration 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 else ! better residuum
dTdLp = 0.0_pReal ! calc dT/dLp dTdLp = 0.0_pReal ! calc dT/dLp
forall (i=1:3,j=1:3,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,m=1:3,n=1:3) &
@ -262,7 +283,8 @@ Inner: do ! inner iteration: Lp
! !
Rinner_old = Rinner ! remember current residuum Rinner_old = Rinner ! remember current residuum
Lpguess_old = Lpguess ! remember current Lp guess Lpguess_old = Lpguess ! remember current Lp guess
if (iInner > 1 .and. leapfrog < maxleap) leapfrog = 2.0_pReal * leapfrog ! accelerate if (iInner > 1 .and. leapfrog < maxleap) leapfrog = 2.0_pReal * leapfrog ! accelerate
! if (debugger) write(6,*) 'leapfrogging at',leapfrog
endif endif
! !
Lpguess = Lpguess_old ! start from current guess Lpguess = Lpguess_old ! start from current guess