diff --git a/trunk/crystallite.f90 b/trunk/crystallite.f90 index 6710decef..5e80bfb29 100644 --- a/trunk/crystallite.f90 +++ b/trunk/crystallite.f90 @@ -38,15 +38,16 @@ CONTAINS use debug use constitutive, only: constitutive_Nstatevars,constitutive_Nresults use mesh, only: mesh_element - use math + use math + implicit none ! character(len=*) msg - logical updateJaco,error,cuttedBack,guessNew + logical updateJaco,error,success,guessNew integer(pInt) cp_en,ip,grain,k,l, nCutbacks, maxCutbacks real(pReal) Temperature real(pReal) dt,dt_aim,subFrac,subStep,det - real(pReal), dimension(3,3) :: Lp,Lp_pert,inv + real(pReal), dimension(3,3) :: Lp,Lp_interpolated,Lp_pert,inv 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_current,Fe_new,Fe_pert @@ -55,74 +56,73 @@ CONTAINS 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 + + deltaFg = Fg_new - Fg_old subFrac = 0.0_pReal subStep = 1.0_pReal - nCutbacks = 0_pInt - maxCutbacks = 0_pInt -! - Fg_aim = Fg_old ! make "new", "aim" a synonym for "old" - Fp_new = Fp_old + nCutbacks = 0_pInt + maxCutbacks = 0_pInt + Fg_current = Fg_old ! initialize to start of inc + Fp_current = Fp_old call math_invert3x3(Fp_old,inv,det,error) -!!$OMP CRITICAL (evilmatmul) - Fe_new = math_mul33x33(Fg_old,inv) -!!$OMP END CRITICAL (evilmatmul) - state_bestguess = state_new ! remember potentially available state guess - state_new = state_old -! - cuttedBack = .false. - guessNew = .true. -! + Fe_current = math_mul33x33(Fg_old,inv) + state_current = state_old + success = .false. ! pretend cutback + dt_aim = 0.0_pReal ! prevent initial Lp interpolation + ! begin the cutback loop do while (subStep > subStepMin) ! continue until finished or too much cut backing - if (.not. cuttedBack) then - Fg_current = Fg_aim ! wind forward - Fp_current = Fp_new - Fe_current = Fe_new - state_current = state_new + if (success) then ! wind forward + Fg_current = Fg_aim + Fe_current = Fe_new + Fp_current = Fp_new + state_current = state_new + elseif (dt_aim > 0.0_pReal) then + call math_invert3x3(Fg_aim,inv,det,error) ! inv of Fg_aim + Lp_interpolated = 0.5_pReal*Lp + & + 0.5_pReal*(math_I3 - math_mul33x33(Fp_current,& + math_mul33x33(inv,Fe_current)))/dt_aim ! interpolate Lp and L + if (debugger) then +!$OMP CRITICAL (write2out) + write (6,*) 'Lp interpolation' + write (6,'(a,/,3(3(f12.7,x)/))') 'from',Lp(1:3,:) + write (6,'(a,/,3(3(f12.7,x)/))') 'to',Lp_interpolated(1:3,:) +!$OMP END CRITICAL (write2out) + endif + Lp = Lp_interpolated endif -! + Fg_aim = Fg_current + subStep*deltaFg ! aim for Fg dt_aim = subStep*dt ! aim for dt - - if (guessNew) & - 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 - math_mul33x33(Fp_current,& - math_mul33x33(inv,Fe_current)))/dt_aim ! interpolate Lp and L -!!$OMP END CRITICAL (evilmatmul) - endif -!!$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 + if (debugger) then +!$OMP CRITICAL (write2out) + write (6,*) 'using these values' + write (6,'(a,/,3(4(f9.3,x)/))') 'state current / MPa',state_current/1e6_pReal + write (6,'(a,/,3(4(f9.3,x)/))') 'state new / MPa',state_new/1e6_pReal + write (6,'(a,/,3(3(f12.7,x)/))') 'Fe current',Fe_current(1:3,:) + write (6,'(a,/,3(3(f12.7,x)/))') 'Fp current',Fp_current(1:3,:) + write (6,'(a,/,3(3(f12.7,x)/))') 'Lp (old=new guess)',Lp(1:3,:) + write (6,'(a20,f,x,a2,x,f)') 'integrating from ',subFrac,'to',(subFrac+subStep) +!$OMP END CRITICAL (write2out) + endif + call TimeIntegration(msg,Lp,Fp_new,Fe_new,P,state_new,post_results, & ! def gradients and PK2 at end of time step + maxval(abs(Fg_aim-Fg_new)) < relevantStrain, & ! post results only if asking for final values dt_aim,cp_en,ip,grain,Temperature,Fg_aim,Fp_current,state_current) -!!$OMP END CRITICAL (timeint) -! - if (msg == 'ok') then - cuttedBack = .false. ! no cut back required - guessNew = .false. ! keep the Lp - subFrac = subFrac + subStep - subStep = 1.0_pReal - subFrac ! try one go - nCutbacks = 0_pInt ! rest cutbackcounter - if (debugger) then -!$OMP CRITICAL (write2out) - write (6,*) '--------- one go -----------++##' -!$OMP END CRITICAL (write2out) - endif + if (msg == 'ok') then + subFrac = subFrac + subStep + subStep = min(1.0_pReal-subFrac, subStep*2.0_pReal) ! accelerate + nCutbacks = 0_pInt ! reset cutback counter + success = .true. ! keep current Lp else nCutbacks = nCutbacks + 1 ! record additional cutback - maxCutbacks=max(nCutbacks,maxCutbacks)! remember maximum number of cutbacks - cuttedBack = .true. ! encountered problems --> - guessNew = .true. ! redo plastic Lp guess + maxCutbacks = max(nCutbacks,maxCutbacks)! remember maximum number of cutbacks subStep = subStep / 2.0_pReal ! cut time step in half - state_bestguess = state_current ! current state is then best guess - if (debugger) then -!$OMP CRITICAL (write2out) - write (6,*) '>>>>>>>>>>>>>>>>>>>> cutback <<<<<<<<<<<<<<<<<<<<<<' -!$OMP END CRITICAL (write2out) + success = .false. ! force Lp interpolation + if (debugger) then +!$OMP CRITICAL (write2out) + write (6,*) '>>>>>>>>>>>>>>>>>>>> cutback <<<<<<<<<<<<<<<<<<<<<<' +!$OMP END CRITICAL (write2out) endif endif @@ -133,59 +133,31 @@ CONTAINS !$OMP END CRITICAL (cutback) ! if (msg /= 'ok') return ! solution not reached --> report back -!!$OMP CRITICAL (write2out) -! write(6,*) '*************************' -! write(6,*) '*************************' -! write(6,*) 'updateJaco', updateJaco, cp_en, ip -! write(6,*) '*************************' -! write(6,*) '*************************' -!!$OMP END CRITICAL (write2out) if (updateJaco) then ! consistent tangent using + if (debugger) then +!$OMP CRITICAL (write2out) + write (6,*) 'Jacobian calc' +!$OMP END CRITICAL (write2out) + endif do k=1,3 do l=1,3 Fg_pert = Fg_new ! initialize perturbed Fg 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 -!!$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 - dt_aim,cp_en,ip,grain,Temperature,Fg_pert, Fp_current,state_current) -!!$OMP END CRITICAL (timeint) -!!$OMP CRITICAL (write2out) -! if(cp_en==61 .and. ip==1) then -! write (6,*) k, l -! write (6,*) msg -! write (6,*) Lp_pert -! write (6,*) Fp_pert -! write (6,*) Fe_pert -! write (6,*) P_pert -! write (6,*) state_pert -! write (6,*) post_results -! write (6,*) dt_aim -! write (6,*) cp_en -! write (6,*) ip -! write (6,*) grain -! write (6,*) Temperature -! write (6,*) Fg_pert -! write (6,*) Fp_current -! write (6,*) state_current -! endif -!!$OMP END CRITICAL (write2out) + dt_aim,cp_en,ip,grain,Temperature,Fg_pert,Fp_current,state_current) + 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 -!!$OMP CRITICAL (write2out) -! write(6,*) '*************************' -! write(6,*) cp_en, ip -! write(6,*) 'dPdF_kl', k, l -! write(6,*) dPdF(:,:,k,l) -! write(6,*) '*************************' -!!$OMP END CRITICAL (write2out) enddo enddo endif -! - msg = 'ok' ! a new consistent tangent was computed even if msg was not ok for all components +! + + msg = 'ok' ! a new consistent tangent was computed even if msg was not ok for all components + ! return ! @@ -247,19 +219,29 @@ CONTAINS msg = 'inversion Fp_old' return endif - -!!$OMP CRITICAL (evilmatmul) - A = math_mul33x33(transpose(invFp_old), math_mul33x33(transpose(Fg_new),math_mul33x33(Fg_new,invFp_old))) -!!$OMP END CRITICAL (evilmatmul) + + A = math_mul33x33(transpose(invFp_old), math_mul33x33(transpose(Fg_new),math_mul33x33(Fg_new,invFp_old))) ! - if (all(state == 0.0_pReal)) then - state = state_old ! former state guessed, if none specified - endif + if (all(state == 0.0_pReal)) state = state_old ! former state guessed, if none specified iOuter = 0_pInt ! outer counter ! + if (debugger) then +!$OMP CRITICAL (write2out) + write (6,'(a,/,3(3(f12.7,x)/))') 'Fg to be calculated',Fg_new +!$OMP END CRITICAL (write2out) + endif ! Outer: do ! outer iteration: State iOuter = iOuter+1 + if (debugger) then +!$OMP CRITICAL (write2out) + write (6,'(a,i3)') '---outer ',iOuter + write (6,'(a,/,3(4(f9.3,x)/))') 'state old / MPa',state_old/1e6_pReal + write (6,'(a,/,3(4(f9.3,x)/))') 'state / MPa',state/1e6_pReal + write (6,'(a,/,3(3(f12.7,x)/))') 'Lpguess',Lpguess(1:3,:) +!$OMP END CRITICAL (write2out) + endif + if (iOuter > nOuter) then msg = 'limit Outer iteration' !$OMP CRITICAL (out) @@ -279,8 +261,16 @@ Outer: do ! outer iteration: State ! Inner: do ! inner iteration: Lp iInner = iInner+1 - if (iInner > nInner) then ! too many loops required - Lpguess=Lpguess_old ! do not trust the last update + if (debugger) then +!$OMP CRITICAL (write2out) + write (6,'(a,i3)') 'inner ',iInner + if (wantsConstitutiveResults .and. iOuter == 1 .and. iInner < 3) then + write (6,'(a,/,3(3(f12.7,x)/))') 'Lpguess',Lpguess(1:3,:) + endif +!$OMP END CRITICAL (write2out) + endif + if (iInner > nInner) then ! too many loops required + Lpguess = Lpguess_old ! do not trust the last update but resort to former one msg = 'limit Inner iteration' !$OMP CRITICAL (in) debug_InnerLoopDistribution(nInner) = debug_InnerLoopDistribution(nInner)+1 @@ -290,55 +280,16 @@ Inner: do ! inner iteration: Lp ! B = math_i3 - dt*Lpguess BT = transpose(B) -!!$OMP CRITICAL (evilmatmul) AB = math_mul33x33(A,B) BTA = math_mul33x33(BT,A) Tstar_v = 0.5_pReal*math_mul66x6(C_66,math_mandel33to6(math_mul33x33(BT,AB)-math_I3)) -!!$OMP END CRITICAL (evilmatmul) 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 -!!$OMP CRITICAL (calcLp) call constitutive_LpAndItsTangent(Lp,dLp, & Tstar_v,state,Temperature,grain,ip,cp_en) -!!$OMP END CRITICAL (calcLp) -!!$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,*) 'Tstar_v' -! write(6,*) Tstar_v -! write(6,*) 'A' -! write(6,*) A -! write(6,*) 'B' -! write(6,*) B -! write(6,*) 'AB' -! write(6,*) AB -! write(6,*) 'BTA' -! write(6,*) BTA -! write(6,*) 'C_66' -! write(6,*) C_66 -! write(6,*) '*************************' -! write(6,*) -! call flush(6) -!endif -!!$OMP END CRITICAL (write2out) ! Rinner = Lpguess - Lp ! update current residuum -!!$OMP CRITICAL (write2out) -! write(6,*) -! write(6,*) '*************************' -! write(6,*) iInner, iOuter -! write(6,*) cp_en, ip -! write(6,*) 'Rinner' -! write(6,*) Rinner -! write(6,*) 'Lp' -! write(6,*) Lp -! write(6,*) 'Lpguess' -! write(6,*) Lpguess -!!$OMP END CRITICAL (write2out) ! if (.not.(any(Rinner/=Rinner)) .and. & ! exclude any NaN in residuum ( (maxval(abs(Rinner)) < abstol_Inner) .or. & ! below abs tol .or. @@ -370,9 +321,7 @@ Inner: do ! inner iteration: Lp dTdLp(3*(i-1)+j,3*(k-1)+l) = dTdLp(3*(i-1)+j,3*(k-1)+l) + & C(i,j,l,n)*AB(k,n)+C(i,j,m,l)*BTA(m,k) dTdLp = -0.5_pReal*dt*dTdLp -!!$OMP CRITICAL (evilmatmul) - dRdLp = eye2 - math_mul99x99(dLp,dTdLp) ! calc dR/dLp -!!$OMP END CRITICAL (evilmatmul) + dRdLp = eye2 - math_mul99x99(dLp,dTdLp) ! calc dR/dLp invdRdLp = 0.0_pReal call math_invert(9,dRdLp,invdRdLp,dummy,failed) ! invert dR/dLp --> dLp/dR if (failed) then @@ -381,9 +330,10 @@ Inner: do ! inner iteration: Lp !$OMP CRITICAL (write2out) write (6,*) msg write (6,'(a,/,9(9(e9.3,x)/))') 'dRdLp', dRdLp(1:9,:) - write (6,*) 'state',state + write (6,'(a,/,3(4(f9.3,x)/))') 'state / MPa',state/1e6_pReal write (6,'(a,/,3(3(f12.7,x)/))') 'Lpguess',Lpguess(1:3,:) - write (6,*) 'Tstar',Tstar_v + write (6,'(a,/,3(3(e12.7,x)/))') 'Lp',Lp(1:3,:) + write (6,'(a,/,6(f9.3,x))') 'Tstar / MPa',Tstar_v/1e6_pReal !$OMP END CRITICAL (write2out) endif return @@ -392,7 +342,6 @@ Inner: do ! inner iteration: Lp 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 ok endif ! @@ -405,28 +354,24 @@ Inner: do ! inner iteration: Lp !$OMP CRITICAL (in) debug_InnerLoopDistribution(iInner) = debug_InnerLoopDistribution(iInner)+1 !$OMP END CRITICAL (in) -!!$OMP CRITICAL (stateupdate) ROuter = state - state_old - & dt*constitutive_dotState(Tstar_v,state,Temperature,& grain,ip,cp_en) ! residuum from evolution of microstructure -!!$OMP END CRITICAL (stateupdate) - state = state - ROuter ! update of microstructure -! if (iOuter==nOuter) then -!!$OMP CRITICAL (write2out) -! write (6,*) 'WARNING: Outer loop has not really converged' -!!$OMP END CRITICAL (write2out) -! exit Outer -! endif + state = state - ROuter ! update of microstructure + + if (iOuter==nOuter) then +!!$OMP CRITICAL (write2out) + write (6,*) 'Terminated outer loop at el,ip,grain',cp_en,ip,grain +!!$OMP END CRITICAL (write2out) + exit Outer + endif if (maxval(abs(Router/state),state /= 0.0_pReal) < reltol_Outer) exit Outer enddo Outer ! !$OMP CRITICAL (out) debug_OuterLoopDistribution(iOuter) = debug_OuterLoopDistribution(iOuter)+1 !$OMP END CRITICAL (out) - -!!$OMP CRITICAL (evilmatmul) invFp_new = math_mul33x33(invFp_old,B) -!!$OMP END CRITICAL (evilmatmul) call math_invert3x3(invFp_new,Fp_new,det,failed) if (failed) then msg = 'inversion Fp_new^-1' @@ -440,20 +385,9 @@ Inner: do ! inner iteration: Lp ! Fp_new = Fp_new*det**(1.0_pReal/3.0_pReal) ! regularize Fp by det = det(InvFp_new) !! forall (i=1:3) Tstar_v(i) = Tstar_v(i)+p_hydro ! add hydrostatic component back -!!$OMP CRITICAL (evilmatmul) - Fe_new = math_mul33x33(Fg_new,invFp_new) ! calc resulting Fe + Fe_new = math_mul33x33(Fg_new,invFp_new) ! calc resulting Fe P = math_mul33x33(Fe_new,math_mul33x33(Tstar,transpose(invFp_new))) ! first PK stress -!!$OMP END CRITICAL (evilmatmul) -! -!!$OMP CRITICAL (write2out) -! write(6,*) -! write(6,*) '*************************' -! write(6,*) cp_en, ip -! write(6,*) iInner, iOuter -! write(6,*) 'P' -! write(6,*) P -! write(6,*) '*************************' -!!$OMP END CRITICAL (write2out) + return ! END SUBROUTINE