From c766ba2e3a1e3fdd60f147679512e6ad4b4bd6fd Mon Sep 17 00:00:00 2001 From: Christoph Kords Date: Mon, 26 Oct 2009 16:43:43 +0000 Subject: [PATCH] minimum substepping size in homogenization and crystallite is now controlled by two independent parameters subStepMinHomog and subStepMinCryst respectively. deleted flush statements where they are not needed, since they seriously slow down computation. --- code/crystallite.f90 | 26 ++++---------------------- code/homogenization.f90 | 13 +++++-------- code/numerics.config | 8 +++++--- code/numerics.f90 | 22 ++++++++++++++-------- 4 files changed, 28 insertions(+), 41 deletions(-) diff --git a/code/crystallite.f90 b/code/crystallite.f90 index 8411d58f5..e237e8b28 100644 --- a/code/crystallite.f90 +++ b/code/crystallite.f90 @@ -226,7 +226,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco) !*** variables and functions from other modules ***! use prec, only: pInt, & pReal - use numerics, only: subStepMin, & + use numerics, only: subStepMinCryst, & pert_Fg, & nState, & nCryst @@ -337,7 +337,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco) NiterationCrystallite = 0_pInt - do while (any(crystallite_subStep(:,:,FEsolving_execELem(1):FEsolving_execElem(2)) > subStepMin)) ! cutback loop for crystallites + do while (any(crystallite_subStep(:,:,FEsolving_execELem(1):FEsolving_execElem(2)) > subStepMinCryst)) ! cutback loop for crystallites if (any(.not. crystallite_converged(:,:,FEsolving_execELem(1):FEsolving_execElem(2)) .and. & ! any non-converged grain .not. crystallite_localConstitution(:,:,FEsolving_execELem(1):FEsolving_execElem(2))) ) & ! has non-local constitution? @@ -358,12 +358,11 @@ subroutine crystallite_stressAndItsTangent(updateJaco) crystallite_subFrac(g,i,e),' to current crystallite_subfrac ', & crystallite_subFrac(g,i,e)+crystallite_subStep(g,i,e),' in crystallite_stressAndItsTangent' write(6,*) - call flush(6) !$OMPEND CRITICAL (write2out) endif crystallite_subFrac(g,i,e) = crystallite_subFrac(g,i,e) + crystallite_subStep(g,i,e) crystallite_subStep(g,i,e) = min(1.0_pReal-crystallite_subFrac(g,i,e), 1.0_pReal * crystallite_subStep(g,i,e)) ! keep cut back step size (no acceleration) - if (crystallite_subStep(g,i,e) > subStepMin) then + if (crystallite_subStep(g,i,e) > subStepMinCryst) then crystallite_subTemperature0(g,i,e) = crystallite_Temperature(g,i,e) ! wind forward... crystallite_subF0(:,:,g,i,e) = crystallite_subF(:,:,g,i,e) ! ...def grad crystallite_subFp0(:,:,g,i,e) = crystallite_Fp(:,:,g,i,e) ! ...plastic def grad @@ -389,12 +388,11 @@ subroutine crystallite_stressAndItsTangent(updateJaco) write(6,'(a78,f10.8)') 'cutback step in crystallite_stressAndItsTangent with new crystallite_subStep: ',& crystallite_subStep(g,i,e) write(6,*) - call flush(6) !$OMPEND CRITICAL (write2out) endif endif - crystallite_onTrack(g,i,e) = crystallite_subStep(g,i,e) > subStepMin ! still on track or already done (beyond repair) + crystallite_onTrack(g,i,e) = crystallite_subStep(g,i,e) > subStepMinCryst ! still on track or already done (beyond repair) if (crystallite_onTrack(g,i,e)) then ! specify task (according to substep) crystallite_subF(:,:,g,i,e) = crystallite_subF0(:,:,g,i,e) + & crystallite_subStep(g,i,e) * & @@ -457,7 +455,6 @@ subroutine crystallite_stressAndItsTangent(updateJaco) if (debugger) then !$OMP CRITICAL (write2out) write(6,*) count(crystallite_onTrack(1,:,:)),'IPs onTrack after preguess for state' - call flush(6) !$OMPEND CRITICAL (write2out) endif @@ -492,7 +489,6 @@ subroutine crystallite_stressAndItsTangent(updateJaco) if (debugger) then !$OMP CRITICAL (write2out) write(6,*) count(crystallite_onTrack(1,:,:)),'IPs onTrack after stress integration' - call flush(6) !$OMPEND CRITICAL (write2out) endif @@ -561,7 +557,6 @@ subroutine crystallite_stressAndItsTangent(updateJaco) write(6,*) count(crystallite_converged(1,:,:)),'IPs converged' write(6,*) count(crystallite_todo(1,:,:)),'IPs todo' write(6,*) - call flush(6) !$OMPEND CRITICAL (write2out) endif @@ -628,7 +623,6 @@ subroutine crystallite_stressAndItsTangent(updateJaco) write (6,'(a,/,3(3(f12.8,x)/))') ' Fp of 1 1 1',myFp(1:3,:) write (6,'(a,/,3(3(f12.8,x)/))') ' Lp of 1 1 1',myLp(1:3,:) write (6,'(a,/,16(6(e12.4,x)/),2(f12.4,x))') 'state of 1 1 1',myState/1e6 - call flush(6) !$OMPEND CRITICAL (write2out) endif do k = 1,3 ! perturbation... @@ -641,7 +635,6 @@ subroutine crystallite_stressAndItsTangent(updateJaco) write (6,'(i1,x,i1)') k,l write (6,*) '=============' write (6,'(a,/,3(3(f12.6,x)/))') 'pertF of 1 1 1',crystallite_subF(1:3,:,g,i,e) - call flush(6) !$OMPEND CRITICAL (write2out) endif onTrack = .true. @@ -667,7 +660,6 @@ subroutine crystallite_stressAndItsTangent(updateJaco) write (6,'(a,/,3(3(f12.4,x)/))') 'DP of 1 1 1',(crystallite_P(1:3,:,g,i,e)-myP(1:3,:))/1e6 write (6,'(a,/,16(6(e12.4,x)/),/,2(f12.4,x))') 'state of 1 1 1',constitutive_state(g,i,e)%p/1e6 write (6,'(a,/,16(6(e12.4,x)/),/,2(f12.4,x))') 'Dstate of 1 1 1',(constitutive_state(g,i,e)%p-myState)/1e6 - call flush(6) !$OMPEND CRITICAL (write2out) endif enddo @@ -762,7 +754,6 @@ endsubroutine if (debugger) then !$OMP CRITICAL (write2out) write(6,*) '::: updateState encountered NaN',g,i,e - call flush(6) !$OMPEND CRITICAL (write2out) endif return @@ -788,7 +779,6 @@ endsubroutine write(6,*) write(6,'(a,/,12(f12.5,x))') 'resid tolerance',abs(residuum/rTol_crystalliteState/constitutive_state(g,i,e)%p(1:mySize)) write(6,*) - call flush(6) !$OMPEND CRITICAL (write2out) endif return @@ -846,7 +836,6 @@ endsubroutine crystallite_updateTemperature = .false. ! indicate update failed !$OMP CRITICAL (write2out) write(6,*) '::: updateTemperature encountered NaN',g,i,e - call flush(6) !$OMPEND CRITICAL (write2out) return endif @@ -974,7 +963,6 @@ endsubroutine write(6,*) '::: integrateStress failed on invFp_current inversion',g,i,e write(6,*) write(6,'(a11,3(i3,x),/,3(3(f12.7,x)/))') 'invFp_new at ',g,i,e,invFp_new - call flush(6) !$OMPEND CRITICAL (write2out) endif return @@ -1004,7 +992,6 @@ LpLoop: do !$OMP CRITICAL (write2out) write(6,*) '::: integrateStress reached loop limit',g,i,e write(6,*) - call flush(6) !$OMPEND CRITICAL (write2out) endif return @@ -1033,7 +1020,6 @@ LpLoop: do write(6,*) write(6,'(a19,3(i3,x),/,3(3(f20.7,x)/))') 'Lp_constitutive at ',g,i,e,Lp_constitutive write(6,'(a11,3(i3,x),/,3(3(f20.7,x)/))') 'Lpguess at ',g,i,e,Lpguess - call flush(6) !$OMPEND CRITICAL (write2out) endif @@ -1055,7 +1041,6 @@ LpLoop: do if (debugger) then !$OMP CRITICAL (write2out) write(6,*) '::: integrateStress encountered NaN at iteration', NiterationStress,'at',g,i,e - call flush(6) !$OMPEND CRITICAL (write2out) endif return @@ -1095,7 +1080,6 @@ LpLoop: do write(6,'(a20,3(i3,x),/,9(9(f15.3,x)/))') 'dLpdT_constitutive at ',g,i,e,dLpdT_constitutive write(6,'(a19,3(i3,x),/,3(3(f20.7,x)/))') 'Lp_constitutive at ',g,i,e,Lp_constitutive write(6,'(a11,3(i3,x),/,3(3(f20.7,x)/))') 'Lpguess at ',g,i,e,Lpguess - call flush(6) !$OMPEND CRITICAL (write2out) endif return @@ -1126,7 +1110,6 @@ LpLoop: do write(6,*) '::: integrateStress failed on invFp_new inversion at iteration', NiterationStress write(6,*) write(6,'(a11,3(i3,x),/,3(3(f12.7,x)/))') 'invFp_new at ',g,i,e,invFp_new - call flush(6) !$OMPEND CRITICAL (write2out) endif return @@ -1155,7 +1138,6 @@ LpLoop: do write(6,'(a,/,3(3(f12.7,x)/))') 'P / MPa',crystallite_P(:,:,g,i,e)/1e6 write(6,'(a,/,3(3(f12.7,x)/))') 'Lp',crystallite_Lp(:,:,g,i,e) write(6,'(a,/,3(3(f12.7,x)/))') 'Fp',crystallite_Fp(:,:,g,i,e) - call flush(6) !$OMP CRITICAL (write2out) endif diff --git a/code/homogenization.f90 b/code/homogenization.f90 index af2f95e22..b02d2ba1d 100644 --- a/code/homogenization.f90 +++ b/code/homogenization.f90 @@ -187,7 +187,7 @@ subroutine materialpoint_stressAndItsTangent(& use prec, only: pInt, & pReal - use numerics, only: subStepMin, & + use numerics, only: subStepMinHomog, & nHomog, & nMPstate use FEsolving, only: FEsolving_execElem, & @@ -236,7 +236,6 @@ subroutine materialpoint_stressAndItsTangent(& write (6,'(a,/,3(3(f12.7,x)/))') 'F of 1 1',materialpoint_F(1:3,:,1,1) write (6,'(a,/,3(3(f12.7,x)/))') 'Fp0 of 1 1 1',crystallite_Fp0(1:3,:,1,1,1) write (6,'(a,/,3(3(f12.7,x)/))') 'Lp0 of 1 1 1',crystallite_Lp0(1:3,:,1,1,1) - call flush(6) !$OMP PARALLEL DO do e = FEsolving_execElem(1),FEsolving_execElem(2) ! iterate over elements to be processed @@ -268,7 +267,7 @@ subroutine materialpoint_stressAndItsTangent(& ! ------ cutback loop ------ - do while (any(materialpoint_subStep(:,FEsolving_execELem(1):FEsolving_execElem(2)) > subStepMin)) ! cutback loop for material points + do while (any(materialpoint_subStep(:,FEsolving_execELem(1):FEsolving_execElem(2)) > subStepMinHomog)) ! cutback loop for material points ! write(6,'(a,/,125(8(f8.5,x),/))') 'mp_subSteps',materialpoint_subStep(:,FEsolving_execELem(1):FEsolving_execElem(2)) !$OMP PARALLEL DO @@ -285,7 +284,6 @@ subroutine materialpoint_stressAndItsTangent(& write(6,'(a21,f10.8,a34,f10.8,a37,/)') 'winding forward from ', & materialpoint_subFrac(i,e), ' to current materialpoint_subFrac ', & materialpoint_subFrac(i,e)+materialpoint_subStep(i,e),' in materialpoint_stressAndItsTangent' - call flush(6) !$OMPEND CRITICAL (write2out) endif @@ -294,7 +292,7 @@ subroutine materialpoint_stressAndItsTangent(& materialpoint_subStep(i,e) = min(1.0_pReal-materialpoint_subFrac(i,e), 1.0_pReal * materialpoint_subStep(i,e)) ! keep cut back time step (no acceleration) ! still stepping needed - if (materialpoint_subStep(i,e) > subStepMin) then + if (materialpoint_subStep(i,e) > subStepMinHomog) then ! wind forward grain starting point of... crystallite_partionedTemperature0(1:myNgrains,i,e) = crystallite_Temperature(1:myNgrains,i,e) ! ...temperatures @@ -322,7 +320,6 @@ subroutine materialpoint_stressAndItsTangent(& !$OMP CRITICAL (write2out) write(6,'(a82,f10.8,/)') 'cutback step in materialpoint_stressAndItsTangent with new materialpoint_subStep: ',& materialpoint_subStep(i,e) - call flush(6) !$OMPEND CRITICAL (write2out) endif @@ -337,7 +334,7 @@ subroutine materialpoint_stressAndItsTangent(& endif - materialpoint_requested(i,e) = materialpoint_subStep(i,e) > subStepMin + materialpoint_requested(i,e) = materialpoint_subStep(i,e) > subStepMinHomog if (materialpoint_requested(i,e)) then materialpoint_subF(:,:,i,e) = materialpoint_subF0(:,:,i,e) + & materialpoint_subStep(i,e) * (materialpoint_F(:,:,i,e) - materialpoint_F0(:,:,i,e)) @@ -349,7 +346,7 @@ subroutine materialpoint_stressAndItsTangent(& !$OMP END PARALLEL DO !* Checks for cutback/substepping loops: added <<>> - ! write (6,'(a,/,8(L,x))') 'MP exceeds substep min',materialpoint_subStep(:,FEsolving_execELem(1):FEsolving_execElem(2)) > subStepMin + ! write (6,'(a,/,8(L,x))') 'MP exceeds substep min',materialpoint_subStep(:,FEsolving_execELem(1):FEsolving_execElem(2)) > subStepMinHomog ! write (6,'(a,/,8(L,x))') 'MP requested',materialpoint_requested(:,FEsolving_execELem(1):FEsolving_execElem(2)) ! write (6,'(a,/,8(f6.4,x))') 'MP subFrac',materialpoint_subFrac(:,FEsolving_execELem(1):FEsolving_execElem(2)) ! write (6,'(a,/,8(f6.4,x))') 'MP subStep',materialpoint_subStep(:,FEsolving_execELem(1):FEsolving_execElem(2)) diff --git a/code/numerics.config b/code/numerics.config index ec807373e..3cb7af291 100644 --- a/code/numerics.config +++ b/code/numerics.config @@ -5,11 +5,13 @@ relevantStrain 1.0e-9 # strain increment considered significant iJacoStiffness 1 # frequency of stiffness update iJacoLpresiduum 1 # frequency of Jacobian update of residuum in Lp pert_Fg 1.0e-7 # strain perturbation for FEM Jacobi -nHomog 25 # homogenization loop limit -nCryst 20 # crystallite loop limit (only for debugging info, real loop limit is "subStepMin") +nHomog 25 # homogenization loop limit (only for debugging info, loop limit is determined by "subStepMinHomog") +subStepMinHomog 1.0e-3 # minimum (relative) size of sub-step allowed during cutback in homogenization +nMPstate 10 # materialpoint state loop limit +nCryst 20 # crystallite loop limit (only for debugging info, loop limit is determined by "subStepMinCryst") +subStepMinCryst 1.0e-3 # minimum (relative) size of sub-step allowed during cutback in crystallite nState 50 # state loop limit nStress 200 # stress loop limit -subStepMin 1.0e-3 # minimum (relative) size of sub-step allowed during cutback in crystallite rTol_crystalliteState 1.0e-6 # relative tolerance in crystallite state loop (abs tol provided by constitutive law) rTol_crystalliteStress 1.0e-6 # relative tolerance in crystallite stress loop (Lp residuum) aTol_crystalliteStress 1.0e-8 # absolute tolerance in crystallite stress loop (Lp residuum!) diff --git a/code/numerics.f90 b/code/numerics.f90 index 6dd5f8b2c..c550dc4cd 100644 --- a/code/numerics.f90 +++ b/code/numerics.f90 @@ -9,14 +9,15 @@ implicit none character(len=64), parameter :: numerics_configFile = 'numerics.config' ! name of configuration file integer(pInt) iJacoStiffness, & ! frequency of stiffness update iJacoLpresiduum, & ! frequency of Jacobian update of residuum in Lp - nHomog, & ! homogenization loop limit + nHomog, & ! homogenization loop limit (only for debugging info, loop limit is determined by "subStepMinHomog") nMPstate, & ! materialpoint state loop limit - nCryst, & ! crystallite loop limit (only for debugging info, real loop limit is "subStepMin") + nCryst, & ! crystallite loop limit (only for debugging info, loop limit is determined by "subStepMinCryst") nState, & ! state loop limit nStress ! stress loop limit real(pReal) relevantStrain, & ! strain increment considered significant pert_Fg, & ! strain perturbation for FEM Jacobi - subStepMin, & ! minimum (relative) size of sub-step allowed during cutback in crystallite + subStepMinCryst, & ! minimum (relative) size of sub-step allowed during cutback in crystallite + subStepMinHomog, & ! minimum (relative) size of sub-step allowed during cutback in homogenization rTol_crystalliteState, & ! relative tolerance in crystallite state loop rTol_crystalliteTemperature, & ! relative tolerance in crystallite temperature loop rTol_crystalliteStress, & ! relative tolerance in crystallite stress loop @@ -76,11 +77,12 @@ subroutine numerics_init() iJacoLpresiduum = 1_pInt pert_Fg = 1.0e-6_pReal nHomog = 20_pInt + subStepMinHomog = 1.0e-3_pReal nMPstate = 10_pInt nCryst = 20_pInt + subStepMinCryst = 1.0e-3_pReal nState = 10_pInt nStress = 40_pInt - subStepMin = 1.0e-3_pReal rTol_crystalliteState = 1.0e-6_pReal rTol_crystalliteTemperature = 1.0e-6_pReal rTol_crystalliteStress = 1.0e-6_pReal @@ -129,8 +131,10 @@ subroutine numerics_init() nState = IO_intValue(line,positions,2) case ('nstress') nStress = IO_intValue(line,positions,2) - case ('substepmin') - subStepMin = IO_floatValue(line,positions,2) + case ('substepmincryst') + subStepMinCryst = IO_floatValue(line,positions,2) + case ('substepminhomog') + subStepMinHomog = IO_floatValue(line,positions,2) case ('rtol_crystallitestate') rTol_crystalliteState = IO_floatValue(line,positions,2) case ('rtol_crystallitetemperature') @@ -175,11 +179,12 @@ subroutine numerics_init() write(6,'(a24,x,i8)') 'iJacoLpresiduum: ',iJacoLpresiduum write(6,'(a24,x,e8.1)') 'pert_Fg: ',pert_Fg write(6,'(a24,x,i8)') 'nHomog: ',nHomog + write(6,'(a24,x,e8.1)') 'subStepMinHomog: ',subStepMinHomog write(6,'(a24,x,i8)') 'nMPstate: ',nMPstate write(6,'(a24,x,i8)') 'nCryst: ',nCryst + write(6,'(a24,x,e8.1)') 'subStepMinCryst: ',subStepMinCryst write(6,'(a24,x,i8)') 'nState: ',nState write(6,'(a24,x,i8)') 'nStress: ',nStress - write(6,'(a24,x,e8.1)') 'subStepMin: ',subStepMin write(6,'(a24,x,e8.1)') 'rTol_crystalliteState: ',rTol_crystalliteState write(6,'(a24,x,e8.1)') 'rTol_crystalliteTemp: ',rTol_crystalliteTemperature write(6,'(a24,x,e8.1)') 'rTol_crystalliteStress: ',rTol_crystalliteStress @@ -207,7 +212,8 @@ subroutine numerics_init() if (nCryst < 1_pInt) call IO_error(265) if (nState < 1_pInt) call IO_error(266) if (nStress < 1_pInt) call IO_error(267) - if (subStepMin <= 0.0_pReal) call IO_error(268) + if (subStepMinCryst <= 0.0_pReal) call IO_error(268) + if (subStepMinHomog <= 0.0_pReal) call IO_error(268) if (rTol_crystalliteState <= 0.0_pReal) call IO_error(269) if (rTol_crystalliteTemperature <= 0.0_pReal) call IO_error(276) !! oops !! if (rTol_crystalliteStress <= 0.0_pReal) call IO_error(270)