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

corrected determination and evaluation of updateJaco 

introduced check whether FFN1 changed and force recollect in that case
This commit is contained in:
Franz Roters 2008-07-14 14:38:19 +00:00
parent 63d81b92b7
commit 4773cf7d09
4 changed files with 65 additions and 56 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

49
trunk/CPFEM_Taylor.f90
View File

@ -126,7 +126,7 @@
real(pReal), dimension (3,3) :: ffn,ffn1,Kirchhoff_bar real(pReal), dimension (3,3) :: ffn,ffn1,Kirchhoff_bar
real(pReal), dimension (3,3,3,3) :: H_bar real(pReal), dimension (3,3,3,3) :: H_bar
real(pReal), dimension(CPFEM_ngens) :: CPFEM_stress real(pReal), dimension(CPFEM_ngens) :: CPFEM_stress
real(pReal), dimension(CPFEM_ngens,CPFEM_ngens) :: CPFEM_jaco real(pReal), dimension(CPFEM_ngens,CPFEM_ngens) :: CPFEM_jaco, odd_jaco
real(pReal) Temperature,CPFEM_dt,J_inverse real(pReal) Temperature,CPFEM_dt,J_inverse
integer(pInt) CPFEM_mode ! 1: regular computation with aged results& integer(pInt) CPFEM_mode ! 1: regular computation with aged results&
! 2: regular computation& ! 2: regular computation&
@ -152,8 +152,23 @@
'mode',CPFEM_mode 'mode',CPFEM_mode
! !
select case (CPFEM_mode) select case (CPFEM_mode)
case (2,1) ! regular computation (with aging of results) case (2,1) ! regular computation (with aging of results)
if (.not. CPFEM_calc_done) then ! puuh, me needs doing all the work... if (any(abs(ffn1 - CPFEM_ffn1_bar(:,:,CPFEM_in,cp_en)) > relevantStrain)) then
CPFEM_stress_bar(1:CPFEM_ngens,:,:) = CPFEM_odd_stress
odd_jaco = CPFEM_odd_jacobian*math_identity2nd(CPFEM_ngens)
forall (i=1:mesh_NcpElems, j=1:FE_Nips(mesh_element(2,e))) &
CPFEM_jaco_bar(1:CPFEM_ngens,1:CPFEM_ngens,j,i) = odd_jaco
outdatedFFN1 = .true.
if (.not. CPFEM_calc_done .AND.CPFEM_mode == 1) then
CPFEM_Fp_old = CPFEM_Fp_new
constitutive_state_old = constitutive_state_new
endif
!$OMP CRITICAL (write2out)
write(6,*) 'WARNING: FFN1 changed for ip', CPFEM_in, 'of element', cp_en
!$OMP END CRITICAL (write2out)
return
else
if (.not. CPFEM_calc_done) then ! puuh, me needs doing all the work...
if (CPFEM_mode == 1) then ! age results at start of new increment if (CPFEM_mode == 1) then ! age results at start of new increment
CPFEM_Fp_old = CPFEM_Fp_new CPFEM_Fp_old = CPFEM_Fp_new
constitutive_state_old = constitutive_state_new constitutive_state_old = constitutive_state_new
@ -172,24 +187,24 @@
!$OMP END PARALLEL DO !$OMP END PARALLEL DO
call debug_info() ! output of debugging/performance statistics call debug_info() ! output of debugging/performance statistics
CPFEM_calc_done = .true. ! now calc is done CPFEM_calc_done = .true. ! now calc is done
endif endif
! translate from P and dP/dF to CS and dCS/dE ! translate from P and dP/dF to CS and dCS/dE
!$OMP CRITICAL (evilmatmul) !$OMP CRITICAL (evilmatmul)
Kirchhoff_bar = matmul(CPFEM_PK1_bar(:,:,CPFEM_in, cp_en),transpose(CPFEM_ffn1_bar(:,:,CPFEM_in, cp_en))) Kirchhoff_bar = matmul(CPFEM_PK1_bar(:,:,CPFEM_in, cp_en),transpose(CPFEM_ffn1_bar(:,:,CPFEM_in, cp_en)))
!$OMP END CRITICAL (evilmatmul) !$OMP END CRITICAL (evilmatmul)
J_inverse = 1.0_pReal/math_det3x3(CPFEM_ffn1_bar(:,:,CPFEM_in, cp_en)) J_inverse = 1.0_pReal/math_det3x3(CPFEM_ffn1_bar(:,:,CPFEM_in, cp_en))
CPFEM_stress_bar(1:CPFEM_ngens,CPFEM_in,cp_en) = math_Mandel33to6(J_inverse*Kirchhoff_bar) CPFEM_stress_bar(1:CPFEM_ngens,CPFEM_in,cp_en) = math_Mandel33to6(J_inverse*Kirchhoff_bar)
! !
H_bar = 0.0_pReal H_bar = 0.0_pReal
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) &
H_bar(i,j,k,l) = H_bar(i,j,k,l) + & H_bar(i,j,k,l) = H_bar(i,j,k,l) + &
CPFEM_ffn1_bar(j,m,CPFEM_in,cp_en) * & CPFEM_ffn1_bar(j,m,CPFEM_in,cp_en) * &
CPFEM_ffn1_bar(l,n,CPFEM_in,cp_en) * & CPFEM_ffn1_bar(l,n,CPFEM_in,cp_en) * &
CPFEM_dPdF_bar(i,m,k,n,CPFEM_in,cp_en) - & CPFEM_dPdF_bar(i,m,k,n,CPFEM_in,cp_en) - &
math_I3(j,l)*CPFEM_ffn1_bar(i,m,CPFEM_in,cp_en)*CPFEM_PK1_bar(k,m,CPFEM_in,cp_en) + & math_I3(j,l)*CPFEM_ffn1_bar(i,m,CPFEM_in,cp_en)*CPFEM_PK1_bar(k,m,CPFEM_in,cp_en) + &
0.5_pReal*(math_I3(i,k)*Kirchhoff_bar(j,l) + math_I3(j,l)*Kirchhoff_bar(i,k) + & 0.5_pReal*(math_I3(i,k)*Kirchhoff_bar(j,l) + math_I3(j,l)*Kirchhoff_bar(i,k) + &
math_I3(i,l)*Kirchhoff_bar(j,k) + math_I3(j,k)*Kirchhoff_bar(i,l)) math_I3(i,l)*Kirchhoff_bar(j,k) + math_I3(j,k)*Kirchhoff_bar(i,l))
CPFEM_jaco_bar(1:CPFEM_ngens,1:CPFEM_ngens,CPFEM_in,cp_en) = math_Mandel3333to66(J_inverse*H_bar) CPFEM_jaco_bar(1:CPFEM_ngens,1:CPFEM_ngens,CPFEM_in,cp_en) = math_Mandel3333to66(J_inverse*H_bar)
! if (CPFEM_in==8 .and. cp_en==80) then ! if (CPFEM_in==8 .and. cp_en==80) then
! do e=1,80 ! do e=1,80
! do i=1,8 ! do i=1,8
@ -201,7 +216,8 @@
! enddo ! enddo
! enddo ! enddo
! endif ! endif
! !
endif
case (3) ! collect and return odd result case (3) ! collect and return odd result
CPFEM_Temperature(CPFEM_in,cp_en) = Temperature CPFEM_Temperature(CPFEM_in,cp_en) = Temperature
CPFEM_ffn_bar(:,:,CPFEM_in,cp_en) = ffn CPFEM_ffn_bar(:,:,CPFEM_in,cp_en) = ffn
@ -255,10 +271,9 @@
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,ip,jp integer(pInt) cp_en,CPFEM_in,grain
logical updateJaco,error,NRconvergent,failed logical updateJaco,error
real(pReal) CPFEM_dt,volfrac,dTime,shMod,C_kb,resNorm,resMax,subStep,subFrac,temp1,temp2 real(pReal) CPFEM_dt,volfrac
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,dPdF_bar_old real(pReal), dimension(3,3,3,3) :: dPdF,dPdF_bar_old
@ -275,7 +290,7 @@
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_Lp(:,:,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_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_dt,cp_en,CPFEM_in,grain,updateJaco,&
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))

2
trunk/FEsolving.f90
View File

@ -9,6 +9,6 @@
integer(pInt) cycleCounter integer(pInt) cycleCounter
integer(pInt) theInc,theCycle,theLovl integer(pInt) theInc,theCycle,theLovl
real(pReal) theTime real(pReal) theTime
logical :: lastIncConverged = .false.,outdatedByNewInc = .false. logical :: lastIncConverged = .false.,outdatedByNewInc = .false., outdatedFFN1 = .false.
END MODULE FEsolving END MODULE FEsolving

51
trunk/crystallite.f90
View File

@ -87,12 +87,17 @@ CONTAINS
Fg_aim = Fg_current + subStep*deltaFg ! aim for Fg Fg_aim = Fg_current + subStep*deltaFg ! aim for Fg
dt_aim = subStep*dt ! aim for dt dt_aim = subStep*dt ! aim for dt
if (guessNew) & if (guessNew) &
state_new = state_bestguess ! try best available guess for state state_new = state_bestguess ! try best available guess for state
if (dt_aim > 0.0_pReal) then
call math_invert3x3(Fg_aim,inv,det,error) ! inv of Fg_aim
!$OMP CRITICAL (evilmatmul)
Lp = 0.5_pReal*Lp + 0.5_pReal*(math_I3 - matmul(Fp_current,matmul(inv,Fe_current)))/dt_aim ! interpolate Lp and L
!$OMP END CRITICAL (evilmatmul)
endif
!!$OMP CRITICAL (timeint) !!$OMP CRITICAL (timeint)
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,0_pInt) dt_aim,cp_en,ip,grain,Temperature,Fg_old,Fg_aim,state_current)
!!$OMP END CRITICAL (timeint) !!$OMP END CRITICAL (timeint)
! !
if (msg == 'ok') then if (msg == 'ok') then
@ -144,7 +149,7 @@ CONTAINS
state_pert = state_new ! initial guess from end of time step state_pert = state_new ! initial guess from end of time step
!!$OMP CRITICAL (timeint) !!$OMP CRITICAL (timeint)
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 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_pInt) dt_aim,cp_en,ip,grain,Temperature,Fg_old,Fg_pert,state_current)
!!$OMP END CRITICAL (timeint) !!$OMP END CRITICAL (timeint)
!!$OMP CRITICAL (write2out) !!$OMP CRITICAL (write2out)
! if(cp_en==61 .and. ip==1) then ! if(cp_en==61 .and. ip==1) then
@ -208,8 +213,7 @@ CONTAINS
Temperature,& ! temperature Temperature,& ! temperature
Fg_new,& ! new total def gradient Fg_new,& ! new total def gradient
Fp_old,& ! former plastic def gradient Fp_old,& ! former plastic def gradient
state_old,& ! former microstructure state_old) ! former microstructure
flag)
use prec use prec
use debug use debug
use mesh, only: mesh_element use mesh, only: mesh_element
@ -224,7 +228,7 @@ CONTAINS
character(len=*) msg character(len=*) msg
logical failed,wantsConstitutiveResults logical failed,wantsConstitutiveResults
integer(pInt) cp_en, ip, grain integer(pInt) cp_en, ip, grain
integer(pInt) iOuter,iInner,dummy, i,j,k,l,m,n,flag integer(pInt) iOuter,iInner,dummy, i,j,k,l,m,n
real(pReal) dt, Temperature, det, p_hydro, leapfrog,maxleap real(pReal) dt, Temperature, det, p_hydro, leapfrog,maxleap
real(pReal), dimension(6) :: Tstar_v real(pReal), dimension(6) :: Tstar_v
real(pReal), dimension(9,9) :: dLp,dTdLp,dRdLp,invdRdLp,eye2 real(pReal), dimension(9,9) :: dLp,dTdLp,dRdLp,invdRdLp,eye2
@ -244,27 +248,10 @@ CONTAINS
msg = 'inversion Fp_old' msg = 'inversion Fp_old'
return return
endif endif
!$OMP CRITICAL (evilmatmul) !$OMP CRITICAL (evilmatmul)
A = matmul(transpose(invFp_old), matmul(transpose(Fg_new),matmul(Fg_new,invFp_old))) A = matmul(transpose(invFp_old), matmul(transpose(Fg_new),matmul(Fg_new,invFp_old)))
!$OMP END CRITICAL (evilmatmul) !$OMP END CRITICAL (evilmatmul)
!!$OMP CRITICAL (write2out)
!if(cp_en==61 .and. ip==1) then
! write(6,*)
! write(6,*) '*************************'
! write(6,*) iInner, iOuter
! write(6,*) cp_en, ip
! write(6,*) 'invFp_old'
! write(6,*) invFp_old
! write(6,*) 'Fg_new'
! write(6,*) Fg_new
! write(6,*) 'A'
! write(6,*) A
! write(6,*) '*************************'
! write(6,*)
! call flush(6)
!endif
!!$OMP END CRITICAL (write2out)
! !
if (all(state == 0.0_pReal)) then if (all(state == 0.0_pReal)) then
state = state_old ! former state guessed, if none specified state = state_old ! former state guessed, if none specified
@ -425,12 +412,12 @@ Inner: do ! inner iteration: Lp
grain,ip,cp_en) ! residuum from evolution of microstructure grain,ip,cp_en) ! residuum from evolution of microstructure
!!$OMP END CRITICAL (stateupdate) !!$OMP END CRITICAL (stateupdate)
state = state - ROuter ! update of microstructure state = state - ROuter ! update of microstructure
if (iOuter==nOuter) then ! if (iOuter==nOuter) then
!$OMP CRITICAL (write2out) !!$OMP CRITICAL (write2out)
write (6,*) 'WARNING: Outer loop has not really converged' ! write (6,*) 'WARNING: Outer loop has not really converged'
!$OMP END CRITICAL (write2out) !!$OMP END CRITICAL (write2out)
exit Outer ! exit Outer
endif ! endif
if (maxval(abs(Router/state),state /= 0.0_pReal) < reltol_Outer) exit Outer if (maxval(abs(Router/state),state /= 0.0_pReal) < reltol_Outer) exit Outer
enddo Outer enddo Outer
! !

19
trunk/mpie_cpfem_marc2007r1.f90
View File

@ -4,7 +4,7 @@
! written by F. Roters, P. Eisenlohr, L. Hantcherli, W.A. Counts ! written by F. Roters, P. Eisenlohr, L. Hantcherli, W.A. Counts
! MPI fuer Eisenforschung, Duesseldorf ! MPI fuer Eisenforschung, Duesseldorf
! !
! last modified: 08.11.2007 ! last modified: 09.07.2008
!******************************************************************** !********************************************************************
! Usage: ! Usage:
! - choose material as hypela2 ! - choose material as hypela2
@ -170,10 +170,13 @@
if (inc == 0) then if (inc == 0) then
cycleCounter = 0 cycleCounter = 4
else else
if (theCycle > ncycle) cycleCounter = 0 ! reset counter for each cutback if (theCycle > ncycle .or. theInc /= inc) cycleCounter = 0 ! reset counter for each cutback or new inc
if (theCycle /= ncycle .or. theLovl /= lovl) cycleCounter = cycleCounter+1 ! ping pong if (theCycle /= ncycle .or. theLovl /= lovl) then
cycleCounter = cycleCounter+1 ! ping pong
outdatedFFN1 = .false.
endif
endif endif
if (cptim > theTime .or. theInc /= inc) then ! reached convergence if (cptim > theTime .or. theInc /= inc) then ! reached convergence
lastIncConverged = .true. lastIncConverged = .true.
@ -192,14 +195,18 @@
if (computationMode == 2 .and. outdatedByNewInc) then if (computationMode == 2 .and. outdatedByNewInc) then
computationMode = 1 ! compute and age former results computationMode = 1 ! compute and age former results
outdatedByNewInc = .false. outdatedByNewInc = .false.
endif endif
if (computationMode == 2 .and. outdatedFFN1) then
computationMode = 4 ! return odd results to force new vyvle
endif
theTime = cptim ! record current starting time theTime = cptim ! record current starting time
theInc = inc ! record current increment number theInc = inc ! record current increment number
theCycle = ncycle ! record current cycle count theCycle = ncycle ! record current cycle count
theLovl = lovl ! record current lovl theLovl = lovl ! record current lovl
call CPFEM_general(computationMode,ffn,ffn1,t(1),timinc,n(1),nn,s,mod(theCycle,2_pInt*ijaco)==0,d,ngens) call CPFEM_general(computationMode,ffn,ffn1,t(1),timinc,n(1),nn,s,mod(cycleCounter-4,4_pInt*ijaco)==0,d,ngens)
! Mandel: 11, 22, 33, SQRT(2)*12, SQRT(2)*23, SQRT(2)*13 ! Mandel: 11, 22, 33, SQRT(2)*12, SQRT(2)*23, SQRT(2)*13