diff --git a/trunk/CPFEM_Taylor.f90 b/trunk/CPFEM_Taylor.f90 index 9cee7337d..c38d8fc1c 100644 --- a/trunk/CPFEM_Taylor.f90 +++ b/trunk/CPFEM_Taylor.f90 @@ -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 diff --git a/trunk/crystallite.f90 b/trunk/crystallite.f90 index 3f38a7ee9..bffa238aa 100644 --- a/trunk/crystallite.f90 +++ b/trunk/crystallite.f90 @@ -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 diff --git a/trunk/prec.f90 b/trunk/prec.f90 index a74e65936..91f0e9319 100644 --- a/trunk/prec.f90 +++ b/trunk/prec.f90 @@ -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) !