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

A revision for CPFEM_Taylor.f90.

This commit is contained in:
Denny Tjahjanto 2008-04-10 11:22:17 +00:00
parent d2312e81ff
commit 9d2ce61698
5 changed files with 206 additions and 168 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

228
trunk/CPFEM_GIA8.f90
View File

@ -217,8 +217,9 @@
! return the local stress and the jacobian from storage ! 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_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) 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) write (6,*) 'stress',CPFEM_stress
if (cp_en == 1 .and. CPFEM_in == 1 .and. CPFEM_updateJaco) write (6,*) 'stiffness',CPFEM_jaco ! if (cp_en == 1 .and. CPFEM_in == 1 .and. CPFEM_updateJaco) write (6,*) 'stiffness',CPFEM_jaco
! if (cp_en == 1 .and. CPFEM_in == 1) write (6,*) 'vector',GIA_rVect_new(:,:,1,1)
! !
return return
! !
@ -234,6 +235,7 @@
cp_en) ! Element number cp_en) ! Element number
! !
use prec use prec
use FEsolving, only: theCycle
use debug use debug
use math, only: math_pDecomposition,math_RtoEuler,inDeg,math_I3,math_invert3x3,math_permut,math_invert,math_delta use math, only: math_pDecomposition,math_RtoEuler,inDeg,math_I3,math_invert3x3,math_permut,math_invert,math_delta
use IO, only: IO_error use IO, only: IO_error
@ -243,12 +245,12 @@
implicit none implicit none
! !
character(len=128) msg character(len=128) msg
integer(pInt) cp_en,CPFEM_in,grain,max_cutbacks,i,j,k,l,m,n,iBoun,NRiter,dummy,ii,jj,kk,ll integer(pInt) cp_en,CPFEM_in,grain,max_cutbacks,i,j,k,l,m,n,iBoun,NRiter,dummy,ii,jj,kk,ll,ip,jp
logical updateJaco,error,NRconvergent,failed logical updateJaco,error,NRconvergent,failed
real(pReal) CPFEM_dt,volfrac,dTime,shMod,C_kb,resNorm,resMax,subStep,subFrac,temp1,temp2 real(pReal) CPFEM_dt,volfrac,dTime,shMod,C_kb,resNorm,resMax,subStep,subFrac,temp1,temp2
real(pReal), dimension(3,3) :: F0_bar,F1_bar,dF_bar real(pReal), dimension(3,3) :: F0_bar,F1_bar,dF_bar,PK1_per,F1_per
real(pReal), dimension(3,3) :: U,R real(pReal), dimension(3,3) :: U,R
real(pReal), dimension(3,3,8) :: PK1,Fp0,Fp1,Fe1,F1,F0,dFgrain,dFg_cor real(pReal), dimension(3,3,8) :: PK1,Fp0,Fp1,Fe1,F1,F0
real(pReal), dimension(3,3,12) :: GPK1,GF1,Nye,GRB1 real(pReal), dimension(3,3,12) :: GPK1,GF1,Nye,GRB1
real(pReal), dimension(3,3,3,3,8) :: dPdF real(pReal), dimension(3,3,3,3,8) :: dPdF
real(pReal), dimension(3,3,3,3,12) :: dRdX1 real(pReal), dimension(3,3,3,3,12) :: dRdX1
@ -297,21 +299,20 @@ NRIteration: do
call GIA_RelaxedDeformation(F1,F1_bar,rx) call GIA_RelaxedDeformation(F1,F1_bar,rx)
! !
! -------------- grain loop ----------------- ! -------------- grain loop -----------------
do grain = 1,8 do grain = 1,texture_Ngrains(mesh_element(4,cp_en))
call SingleCrystallite(msg,PK1(:,:,grain),dPdF(:,:,:,:,grain),& call SingleCrystallite(msg,PK1(:,:,grain),dPdF(:,:,:,:,grain),&
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),&
Fp1(:,:,grain),Fe1(:,:,grain),state1(:,grain),& ! output up to here Fp1(:,:,grain),Fe1(:,:,grain),state1(:,grain),& ! output up to here
dTime,cp_en,CPFEM_in,grain,.true.,& dTime,cp_en,CPFEM_in,grain,.true.,&
CPFEM_Temperature(CPFEM_in,cp_en),F1(:,:,grain),F0(:,:,grain),Fp0(:,:,grain),state0(:,grain)) CPFEM_Temperature(CPFEM_in,cp_en),F1(:,:,grain),F0(:,:,grain),Fp0(:,:,grain),state0(:,grain))
if (msg /= 'ok') then ! solution not reached --> exit NRIteration if (msg /= 'ok') then ! solution not reached --> exit NRIteration
write(6,*) 'GIA: grain loop failed to converge within allowable step-size' write(6,*) 'GIA: grain loop failed to converge @ EL:',cp_en,' IP:',CPFEM_in
NRconvergent = .false. NRconvergent = .false.
exit NRiteration exit NRiteration
endif endif
enddo ! grain loop enddo ! grain loop
! !
! calculate the deformation jump and stress jump across the boundaries ! calculate the deformation jump and stress jump across the boundaries
dFgrain = F1 - F0
call GIA_BoundaryJump(GF1,F1) call GIA_BoundaryJump(GF1,F1)
call GIA_BoundaryJump(GPK1,PK1) call GIA_BoundaryJump(GPK1,PK1)
! !
@ -362,7 +363,7 @@ NRIteration: do
enddo enddo
resNorm = sqrt(resNorm) resNorm = sqrt(resNorm)
! !
! write(6,'(x,a,i3,a,i3,a,i3,a,e10.4)')'EL = ',cp_en,':IP = ',CPFEM_in,':Iter = ',NRiter,':RNorm = ',resNorm if (debugger) write(6,'(x,a,i3,a,i3,a,i3,a,e10.4)')'EL:',cp_en,' IP:',CPFEM_in,' Iter:',NRiter,' RNorm:',resNorm
if (NRiter == 1_pInt) resMax = resNorm if (NRiter == 1_pInt) resMax = resNorm
if ((resNorm < resToler*resMax) .or. (resNorm < resAbsol)) then ! resNorm < tolerance ===> convergent if ((resNorm < resToler*resMax) .or. (resNorm < resAbsol)) then ! resNorm < tolerance ===> convergent
NRconvergent = .true. NRconvergent = .true.
@ -402,7 +403,7 @@ NRIteration: do
dvardres = 0.0_pReal dvardres = 0.0_pReal
call math_invert(36,dresdvar,dvardres,dummy,failed) call math_invert(36,dresdvar,dvardres,dummy,failed)
if (failed) then if (failed) then
write(6,*) 'GIA: failed to invert the Jacobian' write(6,*) 'GIA: failed to invert the Jacobian @ EL:',cp_en,' IP:',CPFEM_in
NRconvergent = .false. NRconvergent = .false.
exit NRiteration exit NRiteration
endif endif
@ -430,8 +431,7 @@ NRIteration: do
! !
! return to the general subroutine when convergence is not reached ! return to the general subroutine when convergence is not reached
if (.not. NRconvergent) then if (.not. NRconvergent) then
write(6,'(x,a)') 'GIA: convergence is not reached within allowable step-size' write(6,'(x,a)') 'GIA: convergence is not reached @ EL:',cp_en,' IP:',CPFEM_in
write(6,'(x,a,i3,a,i3)') 'Element: ',cp_en,' @ IP: ',CPFEM_in
call IO_error(600) call IO_error(600)
return return
endif endif
@ -442,19 +442,10 @@ NRIteration: do
constitutive_state_new(:,:,CPFEM_in,cp_en) = state1 constitutive_state_new(:,:,CPFEM_in,cp_en) = state1
! !
! compute the effective stress and consistent tangent ! compute the effective stress and consistent tangent
call GIA_TangentCorrection(dFgrain,dFg_cor)
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) 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) + & CPFEM_PK1_bar(:,:,CPFEM_in,cp_en) = CPFEM_PK1_bar(:,:,CPFEM_in,cp_en) + &
volfrac*PK1(:,:,grain) ! average Cauchy stress volfrac*PK1(:,:,grain) ! average Cauchy stress
if (updateJaco) then ! consistent tangent
do i = 1,3
do j = 1,3
CPFEM_dPdF_bar(:,:,i,j,CPFEM_in,cp_en) = CPFEM_dPdF_bar(:,:,i,j,CPFEM_in,cp_en) + &
volfrac*dPdF(:,:,i,j,grain)*dFg_cor(i,j,grain)
enddo
enddo
endif
! !
! update results plotted in MENTAT ! update results plotted in MENTAT
call math_pDecomposition(Fe1(:,:,grain),U,R,error) ! polar decomposition call math_pDecomposition(Fe1(:,:,grain),U,R,error) ! polar decomposition
@ -470,6 +461,138 @@ NRIteration: do
CPFEM_results(1:3,grain,CPFEM_in,cp_en) = math_RtoEuler(transpose(R))*inDeg ! orientation CPFEM_results(1:3,grain,CPFEM_in,cp_en) = math_RtoEuler(transpose(R))*inDeg ! orientation
CPFEM_results(4 ,grain,CPFEM_in,cp_en) = volfrac ! volume fraction of orientation CPFEM_results(4 ,grain,CPFEM_in,cp_en) = volfrac ! volume fraction of orientation
enddo enddo
!
if (theCycle >= 0_pInt) then
forall (grain=1:texture_Ngrains(mesh_element(4,cp_en))) &
CPFEM_dPdF_bar(:,:,:,:,CPFEM_in,cp_en) = CPFEM_dPdF_bar(:,:,:,:,CPFEM_in,cp_en) + volfrac*dPdF(:,:,:,:,grain)
else
do ip = 1,3
do jp = 1,3
F1_per = F1_bar
F1_per(ip,jp) = F1_per(ip,jp) + 1.0e-5_pReal
forall (iBoun=1:12,i=1:3) var(3_pInt*(iBoun-1_pInt)+i) = rVect(i,iBoun)
NRiter = 1_pInt
!
NRPerturbation: do
forall (iBoun=1:12,i=1:3) rx(i,iBoun) = var(3_pInt*(iBoun-1_pInt)+i) ! relaxation vectors (guess)
call GIA_RelaxedDeformation(F1,F1_bar,rx)
do grain = 1,8
call SingleCrystallite(msg,PK1(:,:,grain),dPdF(:,:,:,:,grain),&
CPFEM_results(5:4+constitutive_Nresults(grain,CPFEM_in,cp_en),grain,CPFEM_in,cp_en),&
Fp1(:,:,grain),Fe1(:,:,grain),state1(:,grain),& ! output up to here
dTime,cp_en,CPFEM_in,grain,.true.,&
CPFEM_Temperature(CPFEM_in,cp_en),F1(:,:,grain),F0(:,:,grain),Fp0(:,:,grain),state0(:,grain))
if (msg /= 'ok') then ! solution not reached --> exit NRIteration
write(6,*) 'GIA: perturbation grain loop failed to converge within allowable step-size'
NRconvergent = .false.
exit NRPerturbation
endif
enddo
call GIA_BoundaryJump(GF1,F1)
call GIA_BoundaryJump(GPK1,PK1)
!
Nye = 0.0_pReal
NyeNorm = 0.0_pReal
do iBoun = 1,12
do i = 1,3
do j = 1,3
do k = 1,3
do l = 1,3
Nye(i,j,iBoun) = Nye(i,j,iBoun) - 0.5_pReal*math_permut(j,k,l)*GIA_bNorm(k,iBoun)*GF1(i,l,iBoun)
enddo
enddo
NyeNorm(iBoun) = NyeNorm(iBoun) + Nye(i,j,iBoun)*Nye(i,j,iBoun)
enddo
enddo
NyeNorm(iBoun) = sqrt(NyeNorm(iBoun))
if (NyeNorm(iBoun) > 1.0e-8_pReal) Nye(:,:,iBoun) = Nye(:,:,iBoun)/NyeNorm(iBoun)
enddo
!
GRB1 = 0.0_pReal
do iBoun = 1,12
do i = 1,3
do j = 1,3
do k = 1,3
do l = 1,3
GRB1(i,j,iBoun) = GRB1(i,j,iBoun) + Nye(i,k,iBoun)*GIA_bNorm(l,iBoun)*math_permut(k,l,j)
enddo
enddo
enddo
enddo
GRB1(:,:,iBoun) = 0.5_pReal*(C_kb + C_kb)*GRB1(:,:,iBoun)
enddo
!
res = 0.0_pReal
resNorm = 0.0_pReal
do iBoun = 1,12
do j = 1,3
do i = 1,3
res(3_pInt*(iBoun-1_pInt)+j) = res(3_pInt*(iBoun-1_pInt)+j) - &
GIA_bNorm(i,iBoun)*(GPK1(i,j,iBoun) - GRB1(i,j,iBoun))
enddo
resNorm = resNorm + res(3_pInt*(iBoun-1_pInt)+j)*res(3_pInt*(iBoun-1_pInt)+j)
enddo
enddo
resNorm = sqrt(resNorm)
!
! if (debugger) write(6,'(x,a,i3,a,i3,a,i3,a,i3,a,e10.4)')'EL = ',cp_en,':IP = ',CPFEM_in,':pert = ',3*(ip-1)+jp,':Iter = ',NRiter,':RNorm = ',resNorm
if (NRiter == 1_pInt) resMax = resNorm
if ((resNorm < resToler*resMax) .or. (resNorm < resAbsol)) then ! resNorm < tolerance ===> convergent
NRconvergent = .true.
exit NRPerturbation
elseif ((NRiter > NRiterMax) .or. (resNorm > resBound*resMax)) then ! resNorm > up. bound ===> substepping
NRconvergent = .false.
exit NRPerturbation
else ! update the residual
dRdX1 = 0.0_pReal
do iBoun = 1,12
if (NyeNorm(iBoun) < 1.0e-8_pReal) NyeNorm(iBoun) = 1.0e-8_pReal
do i = 1,3
do j = 1,3
do k = 1,3
do l = 1,3
temp1 = 0.0_pReal
temp2 = 0.0_pReal
do ii = 1,3
do jj = 1,3
do kk = 1,3
temp1 = temp1 + GIA_bNorm(jj,iBoun)*math_permut(ii,jj,j)*math_delta(i,k)* &
GIA_bNorm(kk,iBoun)*math_permut(ii,kk,l)
do ll = 1,3
temp2 = temp2 + Nye(i,ii,iBoun)*GIA_bNorm(jj,iBoun)*math_permut(ii,jj,j)* &
Nye(k,kk,iBoun)*GIA_bNorm(ll,iBoun)*math_permut(kk,ll,l)
enddo
enddo
enddo
enddo
dRdX1(i,j,k,l,iBoun) = 0.25_pReal*(C_kb + C_kb)*(temp1 - temp2)/NyeNorm(iBoun)
enddo
enddo
enddo
enddo
enddo
call GIA_JacobianMatrix(dresdvar,dPdF,dRdX1)
dvardres = 0.0_pReal
call math_invert(36,dresdvar,dvardres,dummy,failed)
if (failed) then
write(6,*) 'GIA: perturbation failed to invert the Jacobian'
NRconvergent = .false.
exit NRPerturbation
endif
forall (i=1:36,j=1:36) var(i) = var(i) - dvardres(i,j)*res(j)
endif
NRiter = NRiter + 1_pInt
enddo NRPerturbation ! End of N-R iteration blok
!
PK1_per = 0.0_pReal
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)
PK1_per = PK1_per + volfrac*PK1(:,:,grain)
enddo
CPFEM_dPdF_bar(:,:,ip,jp,CPFEM_in,cp_en) = (PK1_per - CPFEM_PK1_bar(:,:,CPFEM_in,cp_en))/1.0e-5_pReal
enddo
enddo
endif
! !
return return
! !
@ -729,69 +852,6 @@ NRIteration: do
! !
END SUBROUTINE END SUBROUTINE
! !
!
!********************************************************************
! Calculate the correction for the effective consisten tangent
!********************************************************************
subroutine GIA_TangentCorrection(&
dF_grain,& ! deformation gradient increment of grains
ddF_corr) !
!
implicit none
!
real(pReal), dimension(3,3,8) :: dF_grain,ddF_corr
integer(pInt) i,j,iBoun,grain
!
ddF_corr = 0.0_pReal
do j = 1,3
!
! first relaxation direction
ddF_corr(1,j,1) = 2.0_pReal*abs(dF_grain(1,j,1))/(abs(dF_grain(1,j,1)) + abs(dF_grain(1,j,2)))
if (abs(dF_grain(1,j,1)) < 1.0e-8_pReal) ddF_corr(1,j,1) = 1.0_pReal
ddF_corr(1,j,2) = 2.0_pReal - ddF_corr(1,j,1)
ddF_corr(1,j,3) = 2.0_pReal*abs(dF_grain(1,j,3))/(abs(dF_grain(1,j,3)) + abs(dF_grain(1,j,4)))
if (abs(dF_grain(1,j,3)) < 1.0e-8_pReal) ddF_corr(1,j,3) = 1.0_pReal
ddF_corr(1,j,4) = 2.0_pReal - ddF_corr(1,j,3)
ddF_corr(1,j,5) = 2.0_pReal*abs(dF_grain(1,j,5))/(abs(dF_grain(1,j,5)) + abs(dF_grain(1,j,6)))
if (abs(dF_grain(1,j,5)) < 1.0e-8_pReal) ddF_corr(1,j,5) = 1.0_pReal
ddF_corr(1,j,6) = 2.0_pReal - ddF_corr(1,j,5)
ddF_corr(1,j,7) = 2.0_pReal*abs(dF_grain(1,j,7))/(abs(dF_grain(1,j,7)) + abs(dF_grain(1,j,8)))
if (abs(dF_grain(1,j,7)) < 1.0e-8_pReal) ddF_corr(1,j,7) = 1.0_pReal
ddF_corr(1,j,8) = 2.0_pReal - ddF_corr(1,j,7)
!
! second relaxation direction
ddF_corr(2,j,1) = 2.0_pReal*abs(dF_grain(2,j,1))/(abs(dF_grain(2,j,1)) + abs(dF_grain(2,j,3)))
if (abs(dF_grain(2,j,1)) < 1.0e-8_pReal) ddF_corr(2,j,1) = 1.0_pReal
ddF_corr(2,j,2) = 2.0_pReal*abs(dF_grain(2,j,2))/(abs(dF_grain(2,j,2)) + abs(dF_grain(2,j,4)))
if (abs(dF_grain(2,j,2)) < 1.0e-8_pReal) ddF_corr(2,j,2) = 1.0_pReal
ddF_corr(2,j,3) = 2.0_pReal - ddF_corr(2,j,1)
ddF_corr(2,j,4) = 2.0_pReal - ddF_corr(2,j,2)
ddF_corr(2,j,5) = 2.0_pReal*abs(dF_grain(2,j,5))/(abs(dF_grain(2,j,5)) + abs(dF_grain(2,j,7)))
if (abs(dF_grain(2,j,5)) < 1.0e-8_pReal) ddF_corr(2,j,5) = 1.0_pReal
ddF_corr(2,j,6) = 2.0_pReal*abs(dF_grain(2,j,6))/(abs(dF_grain(2,j,6)) + abs(dF_grain(2,j,8)))
if (abs(dF_grain(2,j,6)) < 1.0e-8_pReal) ddF_corr(2,j,6) = 1.0_pReal
ddF_corr(2,j,7) = 2.0_pReal - ddF_corr(2,j,5)
ddF_corr(2,j,8) = 2.0_pReal - ddF_corr(2,j,6)
!
! third relaxation direction
ddF_corr(3,j,1) = 2.0_pReal*abs(dF_grain(3,j,1))/(abs(dF_grain(3,j,1)) + abs(dF_grain(3,j,5)))
if (abs(dF_grain(3,j,1)) < 1.0e-8_pReal) ddF_corr(3,j,1) = 1.0_pReal
ddF_corr(3,j,2) = 2.0_pReal*abs(dF_grain(3,j,2))/(abs(dF_grain(3,j,2)) + abs(dF_grain(3,j,6)))
if (abs(dF_grain(3,j,2)) < 1.0e-8_pReal) ddF_corr(3,j,2) = 1.0_pReal
ddF_corr(3,j,3) = 2.0_pReal*abs(dF_grain(3,j,3))/(abs(dF_grain(3,j,3)) + abs(dF_grain(3,j,7)))
if (abs(dF_grain(3,j,3)) < 1.0e-8_pReal) ddF_corr(3,j,3) = 1.0_pReal
ddF_corr(3,j,4) = 2.0_pReal*abs(dF_grain(3,j,4))/(abs(dF_grain(3,j,4)) + abs(dF_grain(3,j,8)))
if (abs(dF_grain(3,j,4)) < 1.0e-8_pReal) ddF_corr(3,j,4) = 1.0_pReal
ddF_corr(3,j,5) = 2.0_pReal - ddF_corr(3,j,1)
ddF_corr(3,j,6) = 2.0_pReal - ddF_corr(3,j,2)
ddF_corr(3,j,7) = 2.0_pReal - ddF_corr(3,j,3)
ddF_corr(3,j,8) = 2.0_pReal - ddF_corr(3,j,4)
!
enddo
!
return
!
END SUBROUTINE
! !
END MODULE END MODULE
!############################################################## !##############################################################

68
trunk/CPFEM_Taylor.f90
View File

@ -199,8 +199,8 @@
! return the local stress and the jacobian from storage ! 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_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) 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) write (6,*) 'stress',CPFEM_stress
if (cp_en == 1 .and. CPFEM_in == 1 .and. CPFEM_updateJaco) write (6,*) 'stiffness',CPFEM_jaco ! if (cp_en == 1 .and. CPFEM_in == 1 .and. CPFEM_updateJaco) write (6,*) 'stiffness',CPFEM_jaco
! !
return return
! !
@ -216,6 +216,7 @@
cp_en) ! Element number cp_en) ! Element number
! !
use prec use prec
use FEsolving, only: theCycle
use debug use debug
use math, only: math_pDecomposition,math_RtoEuler,inDeg,math_I3,math_invert3x3,math_permut,math_invert,math_delta use math, only: math_pDecomposition,math_RtoEuler,inDeg,math_I3,math_invert3x3,math_permut,math_invert,math_delta
use IO, only: IO_error use IO, only: IO_error
@ -225,39 +226,34 @@
implicit none implicit none
! !
character(len=128) msg character(len=128) msg
integer(pInt) cp_en,CPFEM_in,grain,max_cutbacks,i,j,k,l,m,n integer(pInt) cp_en,CPFEM_in,grain,max_cutbacks,i,j,k,l,m,n,iBoun,NRiter,dummy,ii,jj,kk,ll,ip,jp
logical updateJaco,error logical updateJaco,error,NRconvergent,failed
real(pReal) CPFEM_dt,volfrac real(pReal) CPFEM_dt,volfrac,dTime,shMod,C_kb,resNorm,resMax,subStep,subFrac,temp1,temp2
real(pReal), dimension(3,3) :: F0_bar,F1_bar,dF_bar real(pReal), dimension(3,3) :: PK1_pert,F1_pert
real(pReal), dimension(3,3) :: U,R real(pReal), dimension(3,3) :: U,R,Fe1
real(pReal), dimension(3,3,8) :: PK1,Fp0,Fp1,Fe1 real(pReal), dimension(3,3) :: PK1
real(pReal), dimension(3,3,3,3,8) :: dPdF real(pReal), dimension(3,3,3,3) :: dPdF
real(pReal), dimension(constitutive_maxNstatevars,8) :: state0,state1
! !
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) CPFEM_dPdF_bar(:,:,:,:,CPFEM_in,cp_en) = 0.0_pReal ! zero out average consistent tangent
!
F0_bar = CPFEM_ffn_bar(:,:,CPFEM_in,cp_en) ! effective deformation gradient at t_n
F1_bar = CPFEM_ffn1_bar(:,:,CPFEM_in,cp_en) ! effective deformation gradient at t_n+1
state0 = constitutive_state_old(:,:,CPFEM_in,cp_en) ! state variables at t_n
Fp0 = CPFEM_Fp_old(:,:,:,CPFEM_in,cp_en) ! grain plastic def. gradient at t_n
!
! -------------- grain loop -----------------
do grain = 1,texture_Ngrains(mesh_element(4,cp_en))
call SingleCrystallite(msg,PK1(:,:,grain),dPdF(:,:,:,:,grain),&
CPFEM_results(5:4+constitutive_Nresults(grain,CPFEM_in,cp_en),grain,CPFEM_in,cp_en),&
Fp1(:,:,grain),Fe1(:,:,grain),state1(:,grain),& ! output up to here
CPFEM_dt,cp_en,CPFEM_in,grain,.true.,&
CPFEM_Temperature(CPFEM_in,cp_en),F1_bar,F0_bar,Fp0(:,:,grain),state0(:,grain))
if (msg /= 'ok') then ! solution not reached
call IO_error(600)
return
endif
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)
call SingleCrystallite(msg,PK1,dPdF,&
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_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))
!
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
!
! update results plotted in MENTAT ! update results plotted in MENTAT
call math_pDecomposition(Fe1(:,:,grain),U,R,error) ! polar decomposition call math_pDecomposition(Fe1,U,R,error) ! polar decomposition
if (error) then if (error) then
write(6,*) Fe1(:,:,grain) write(6,*) Fe1
write(6,*) 'polar decomposition' write(6,*) 'polar decomposition'
write(6,*) 'Grain: ',grain write(6,*) 'Grain: ',grain
write(6,*) 'Integration point: ',CPFEM_in write(6,*) 'Integration point: ',CPFEM_in
@ -265,21 +261,9 @@
call IO_error(650) call IO_error(650)
return return
endif endif
!
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(:,:,grain) ! average Cauchy stress
if (updateJaco) then ! consistent tangent
CPFEM_dPdF_bar(:,:,:,:,CPFEM_in,cp_en) = CPFEM_dPdF_bar(:,:,:,:,CPFEM_in,cp_en) + &
volfrac*dPdF(:,:,:,:,grain)
endif
CPFEM_results(1:3,grain,CPFEM_in,cp_en) = math_RtoEuler(transpose(R))*inDeg ! orientation CPFEM_results(1:3,grain,CPFEM_in,cp_en) = math_RtoEuler(transpose(R))*inDeg ! orientation
CPFEM_results(4 ,grain,CPFEM_in,cp_en) = volfrac ! volume fraction of orientation CPFEM_results(4 ,grain,CPFEM_in,cp_en) = volfrac ! volume fraction of orientation
enddo ! grain loop enddo ! grain
!
! updates all variables, deformation gradients, and vectors
CPFEM_Fp_new(:,:,:,CPFEM_in,cp_en) = Fp1
constitutive_state_new(:,:,CPFEM_in,cp_en) = state1
! !
return return
! !

11
trunk/IO.f90
View File

@ -40,6 +40,7 @@
character(256) path character(256) path
inquire(6, name=path) ! determine outputfile inquire(6, name=path) ! determine outputfile
write(6,*)'trying to open mattex at',path(1:scan(path,pathSep,back=.true.))//relPath
open(unit,status='old',err=100,file=path(1:scan(path,pathSep,back=.true.))//relPath) open(unit,status='old',err=100,file=path(1:scan(path,pathSep,back=.true.))//relPath)
IO_open_file = .true. IO_open_file = .true.
return return
@ -568,25 +569,17 @@
select case (ID) select case (ID)
case (100) case (100)
msg='Unable to open input file.' msg='Unable to open input file.'
case (110) case (110)
msg='No materials specified via State Variable 2.' msg='No materials specified via State Variable 2.'
case (120) case (120)
msg='No textures specified via State Variable 3.' msg='No textures specified via State Variable 3.'
case (200) case (200)
msg='Error reading from material+texture file' msg='Error reading from material+texture file'
case (300) case (300)
msg='This material can only be used with & msg='This material can only be used with elements with three direct stress components'
&elements with three direct stress components'
case (400) case (400)
msg='Unknown alloy number specified' msg='Unknown alloy number specified'
case (500) case (500)
msg='Unknown lattice type specified' msg='Unknown lattice type specified'
case (600) case (600)

23
trunk/crystallite.f90
View File

@ -104,9 +104,9 @@ CONTAINS
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
endif endif
enddo enddo ! potential substepping
! !
debug_cutbackDistribution(nCutbacks+1) = debug_cutbackDistribution(nCutbacks+1)+1 debug_cutbackDistribution(min(nCutback,nCutbacks)+1) = debug_cutbackDistribution(min(nCutback,nCutbacks)+1)+1
! !
if (msg /= 'ok') return ! solution not reached --> report back if (msg /= 'ok') return ! solution not reached --> report back
if (updateJaco) then ! consistent tangent using if (updateJaco) then ! consistent tangent using
@ -119,6 +119,7 @@ CONTAINS
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') then
write(6,*)'killed myself at component',k,l
msg = 'consistent tangent --> '//msg msg = 'consistent tangent --> '//msg
return return
endif endif
@ -249,13 +250,13 @@ Inner: do ! inner iteration: Lp
call math_invert(9,dRdLp,invdRdLp,dummy,failed) ! invert dR/dLp --> dLp/dR call math_invert(9,dRdLp,invdRdLp,dummy,failed) ! invert dR/dLp --> dLp/dR
if (failed) then if (failed) then
msg = 'inversion dR/dLp' msg = 'inversion dR/dLp'
if (debugger) then ! if (debugger) then
write (6,*) msg ! write (6,*) msg
write (6,'(a,/,9(9(e9.3,x)/))') 'dRdLp', dRdLp(1:9,:) ! write (6,'(a,/,9(9(e9.3,x)/))') 'dRdLp', dRdLp(1:9,:)
write (6,*) 'state',state ! write (6,*) 'state',state
write (6,'(a,/,3(3(f12.7,x)/))') 'Lpguess',Lpguess(1:3,:) ! write (6,'(a,/,3(3(f12.7,x)/))') 'Lpguess',Lpguess(1:3,:)
write (6,*) 'Tstar',Tstar_v ! write (6,*) 'Tstar',Tstar_v
endif ! endif
return return
endif endif
! !
@ -280,7 +281,6 @@ Inner: do ! inner iteration: Lp
! !
debug_OuterLoopDistribution(iOuter) = debug_OuterLoopDistribution(iOuter)+1 debug_OuterLoopDistribution(iOuter) = debug_OuterLoopDistribution(iOuter)+1
invFp_new = matmul(invFp_old,B) invFp_new = matmul(invFp_old,B)
call math_invert3x3(invFp_new,Fp_new,det,failed) call math_invert3x3(invFp_new,Fp_new,det,failed)
if (failed) then if (failed) then
@ -301,7 +301,8 @@ Inner: do ! inner iteration: Lp
return return
! !
END SUBROUTINE END SUBROUTINE
!! !
!
END MODULE END MODULE
!############################################################## !##############################################################

14
trunk/prec.f90
View File

@ -14,18 +14,18 @@
! *** Numerical parameters *** ! *** Numerical parameters ***
integer(pInt), parameter :: ijaco = 1_pInt ! frequency of FEM Jacobi update integer(pInt), parameter :: ijaco = 1_pInt ! frequency of FEM Jacobi update
integer(pInt), parameter :: nCutback = 10_pInt ! cutbacks in time-step integration integer(pInt), parameter :: nCutback = 20_pInt ! max cutbacks accounted for in debug distribution
integer(pInt), parameter :: nReg = 1_pInt ! regularization attempts for Jacobi inversion 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 real(pReal), parameter :: pert_Fg = 1.0e-6_pReal ! strain perturbation for FEM Jacobi
integer(pInt), parameter :: nOuter = 10_pInt ! outer loop limit integer(pInt), parameter :: nOuter = 20_pInt ! outer loop limit
integer(pInt), parameter :: nInner = 1000_pInt ! inner loop limit integer(pInt), parameter :: nInner = 200_pInt ! inner loop limit
real(pReal), parameter :: reltol_Outer = 1.0e-6_pReal ! relative tolerance in outer loop (state) real(pReal), parameter :: reltol_Outer = 1.0e-4_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 :: 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) real(pReal), parameter :: abstol_Inner = 1.0e-8_pReal ! absolute tolerance in inner loop (Lp)
! !
real(pReal), parameter :: resToler = 1.0e-6_pReal ! relative tolerance of residual in GIA iteration real(pReal), parameter :: resToler = 1.0e-4_pReal ! relative tolerance of residual in GIA iteration
real(pReal), parameter :: resAbsol = 1.0e+0_pReal ! absolute tolerance of residual in GIA iteration (corresponds to 1 Pa) real(pReal), parameter :: resAbsol = 1.0e+2_pReal ! absolute tolerance of residual in GIA iteration (corresponds to ~1 Pa)
real(pReal), parameter :: resBound = 1.0e+2_pReal ! relative maximum value (upper bound) for GIA residual real(pReal), parameter :: resBound = 1.0e+1_pReal ! relative maximum value (upper bound) for GIA residual
integer(pInt), parameter :: NRiterMax = 24_pInt ! maximum number of GIA iteration integer(pInt), parameter :: NRiterMax = 24_pInt ! maximum number of GIA iteration
real(pReal), parameter :: subStepMin = 1.0e-3_pReal ! minimum (relative) size of sub-step allowed during cutback real(pReal), parameter :: subStepMin = 1.0e-3_pReal ! minimum (relative) size of sub-step allowed during cutback