diff --git a/code/IO.f90 b/code/IO.f90 index c33b63210..2e0e6a2cc 100644 --- a/code/IO.f90 +++ b/code/IO.f90 @@ -1048,6 +1048,10 @@ endfunction case (277) msg = 'Non-positive relevant mismatch in RGC' +!* Error message when selected perturbation method is not defined + case (299) + msg = 'Chosen prturbation method does not exist' + case (300) msg = 'This material can only be used with elements with three direct stress components' case (500) diff --git a/code/crystallite.f90 b/code/crystallite.f90 index e237e8b28..9cb706d51 100644 --- a/code/crystallite.f90 +++ b/code/crystallite.f90 @@ -227,7 +227,10 @@ subroutine crystallite_stressAndItsTangent(updateJaco) use prec, only: pInt, & pReal use numerics, only: subStepMinCryst, & + subStepSizeCryst, & + stepIncreaseCryst, & pert_Fg, & + pert_method, & nState, & nCryst use debug, only: debugger, & @@ -293,11 +296,12 @@ subroutine crystallite_stressAndItsTangent(updateJaco) temperatureConverged, & ! flag indicating if temperature converged stateConverged, & ! flag indicating if state converged converged ! flag indicating if iteration converged - - real(pReal), dimension(9,9) :: dPdF99 + real(pReal), dimension(9,9) :: dPdF99 + real(pReal), dimension(3,3,3,3) :: dPdF_pos,dPdF_neg ! ------ initialize to starting condition ------ - + centralDifference = .true. + !$OMP CRITICAL (write2out) ! write (6,*) ! write (6,*) 'Crystallite request from Materialpoint' @@ -323,7 +327,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco) crystallite_subTstar0_v(:,g,i,e) = crystallite_partionedTstar0_v(:,g,i,e) ! ...2nd PK stress crystallite_subFrac(g,i,e) = 0.0_pReal - crystallite_subStep(g,i,e) = 2.0_pReal + crystallite_subStep(g,i,e) = 1.0_pReal/subStepSizeCryst ! <> crystallite_onTrack(g,i,e) = .true. crystallite_converged(g,i,e) = .false. ! pretend failed step of twice the required size endif @@ -361,7 +365,8 @@ subroutine crystallite_stressAndItsTangent(updateJaco) !$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) + crystallite_subStep(g,i,e) = min(1.0_pReal-crystallite_subFrac(g,i,e), & + stepIncreaseCryst*crystallite_subStep(g,i,e)) ! <> 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 @@ -376,7 +381,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco) !$OMPEND CRITICAL (distributionCrystallite) endif else - crystallite_subStep(g,i,e) = 0.5_pReal * crystallite_subStep(g,i,e) ! cut step in half and restore... + crystallite_subStep(g,i,e) = subStepSizeCryst*crystallite_subStep(g,i,e) ! cut step in half and restore... crystallite_Temperature(g,i,e) = crystallite_subTemperature0(g,i,e) ! ...temperature crystallite_Fp(:,:,g,i,e) = crystallite_subFp0(:,:,g,i,e) ! ...plastic def grad crystallite_invFp(:,:,g,i,e) = math_inv3x3(crystallite_Fp(:,:,g,i,e)) @@ -625,8 +630,11 @@ subroutine crystallite_stressAndItsTangent(updateJaco) write (6,'(a,/,16(6(e12.4,x)/),2(f12.4,x))') 'state of 1 1 1',myState/1e6 !$OMPEND CRITICAL (write2out) endif - do k = 1,3 ! perturbation... - do l = 1,3 ! ...components + +! begin perturbation of components of F + if (pert_method == 1_pInt .or. pert_method == 3_pInt) then ! <<< when forward or central difference is desired >>> + do k = 1,3 ! perturbation... + do l = 1,3 ! ...components to the positive direction crystallite_subF(:,:,g,i,e) = myF ! initialize perturbed F to match converged crystallite_subF(k,l,g,i,e) = crystallite_subF(k,l,g,i,e) + pert_Fg ! perturb single component if (debugger) then @@ -663,8 +671,10 @@ subroutine crystallite_stressAndItsTangent(updateJaco) !$OMPEND CRITICAL (write2out) endif enddo - if (converged) & ! converged state warrants stiffness update - crystallite_dPdF(:,:,k,l,g,i,e) = (crystallite_P(:,:,g,i,e) - myP)/pert_Fg ! tangent dP_ij/dFg_kl + if (converged) then ! converged state warrants stiffness update + dPdF_pos(:,:,k,l) = (crystallite_P(:,:,g,i,e) - myP)/pert_Fg ! tangent dP_ij/dFg_kl + if (pert_method == 1_pInt) crystallite_dPdF(:,:,k,l,g,i,e) = dPdF_pos(:,:,k,l) + endif constitutive_state(g,i,e)%p = myState ! restore unperturbed, converged state, ... constitutive_dotState(g,i,e)%p = myDotState ! ... dotState, ... crystallite_Temperature(g,i,e) = myTemperature ! ... temperature, ... @@ -679,7 +689,68 @@ subroutine crystallite_stressAndItsTangent(updateJaco) debug_StiffnessstateLoopDistribution(NiterationState) + 1 !$OMPEND CRITICAL (out) enddo - enddo + enddo + endif + if (pert_method == 2_pInt .or. pert_method == 3_pInt) then ! <<< when backward or central difference is desired >>> + do k = 1,3 ! perturbation... + do l = 1,3 ! ...components to the negative direction + crystallite_subF(:,:,g,i,e) = myF ! initialize perturbed F to match converged + crystallite_subF(k,l,g,i,e) = crystallite_subF(k,l,g,i,e) - pert_Fg ! perturb single component + if (debugger) then + !$OMP CRITICAL (write2out) + write (6,*) '=============' + 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) + !$OMPEND CRITICAL (write2out) + endif + onTrack = .true. + converged = .false. + NiterationState = 0_pInt + do while(.not. converged .and. onTrack .and. NiterationState < nState) ! keep cycling until done (potentially non-converged) + NiterationState = NiterationState + 1_pInt + onTrack = crystallite_integrateStress(g,i,e) ! stress of perturbed situation (overwrites _P,_Tstar_v,_Fp,_Lp,_Fe) + if (onTrack) then + call constitutive_collectDotState(crystallite_Tstar_v(:,g,i,e), crystallite_subTstar0_v(:,g,i,e), & + crystallite_Fe, crystallite_Fp, crystallite_Temperature(g,i,e), & + crystallite_subdt(g,i,e), g, i, e) + + stateConverged = crystallite_updateState(g,i,e) ! update state + temperatureConverged = crystallite_updateTemperature(g,i,e) ! update temperature + converged = stateConverged .and. temperatureConverged + endif + if (debugger) then + !$OMP CRITICAL (write2out) + write (6,*) '-------------' + write (6,'(a,x,l,x,l)') 'ontrack + converged:',onTrack,converged + write (6,'(a,/,3(3(f12.4,x)/))') 'pertP of 1 1 1',crystallite_P(1:3,:,g,i,e)/1e6 + 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 + !$OMPEND CRITICAL (write2out) + endif + enddo + if (converged) then ! converged state warrants stiffness update + dPdF_neg(:,:,k,l) = (myP - crystallite_P(:,:,g,i,e))/pert_Fg ! tangent dP_ij/dFg_kl + if (pert_method == 2_pInt) crystallite_dPdF(:,:,k,l,g,i,e) = dPdF_neg(:,:,k,l) + endif + constitutive_state(g,i,e)%p = myState ! restore unperturbed, converged state, ... + constitutive_dotState(g,i,e)%p = myDotState ! ... dotState, ... + crystallite_Temperature(g,i,e) = myTemperature ! ... temperature, ... + crystallite_Fp(:,:,g,i,e) = myFp ! ... and kinematics + crystallite_invFp(:,:,g,i,e) = myInvFp + crystallite_Fe(:,:,g,i,e) = myFe + crystallite_Lp(:,:,g,i,e) = myLp + crystallite_Tstar_v(:,g,i,e) = myTstar_v + crystallite_P(:,:,g,i,e) = myP + !$OMP CRITICAL (out) + debug_StiffnessStateLoopDistribution(NiterationState) = & + debug_StiffnessstateLoopDistribution(NiterationState) + 1 + !$OMPEND CRITICAL (out) + enddo + enddo + endif + if (pert_method == 3_pInt) crystallite_dPdF(:,:,:,:,g,i,e) = 0.5_pReal*(dPdF_neg + dPdF_pos) else ! grain did not converged crystallite_dPdF(:,:,:,:,g,i,e) = crystallite_fallbackdPdF(:,:,:,:,g,i,e) ! use (elastic) fallback endif ! grain convergence diff --git a/code/homogenization.f90 b/code/homogenization.f90 index 635a52b21..74d27717f 100644 --- a/code/homogenization.f90 +++ b/code/homogenization.f90 @@ -188,6 +188,8 @@ subroutine materialpoint_stressAndItsTangent(& use prec, only: pInt, & pReal use numerics, only: subStepMinHomog, & + subStepSizeHomog, & + stepIncreaseHomog, & nHomog, & nMPstate use FEsolving, only: FEsolving_execElem, & @@ -256,7 +258,7 @@ subroutine materialpoint_stressAndItsTangent(& materialpoint_subF0(:,:,i,e) = materialpoint_F0(:,:,i,e) ! ...def grad materialpoint_subFrac(i,e) = 0.0_pReal - materialpoint_subStep(i,e) = 8.0_pReal + materialpoint_subStep(i,e) = 1.0_pReal/subStepSizeHomog ! <> materialpoint_converged(i,e) = .false. ! pretend failed step of twice the required size materialpoint_requested(i,e) = .true. ! everybody requires calculation enddo @@ -289,7 +291,8 @@ subroutine materialpoint_stressAndItsTangent(& ! calculate new subStep and new subFrac materialpoint_subFrac(i,e) = materialpoint_subFrac(i,e) + materialpoint_subStep(i,e) - 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) + materialpoint_subStep(i,e) = min(1.0_pReal-materialpoint_subFrac(i,e), & + stepIncreaseHomog*materialpoint_subStep(i,e)) ! <> ! still stepping needed if (materialpoint_subStep(i,e) > subStepMinHomog) then @@ -314,7 +317,8 @@ subroutine materialpoint_stressAndItsTangent(& ! materialpoint didn't converge, so we need a cutback here else - materialpoint_subStep(i,e) = 0.125_pReal * materialpoint_subStep(i,e) ! crystallite had severe trouble, so do a significant cutback + materialpoint_subStep(i,e) = subStepSizeHomog * materialpoint_subStep(i,e) ! crystallite had severe trouble, so do a significant cutback + ! <> if (debugger) then !$OMP CRITICAL (write2out) @@ -436,6 +440,7 @@ elementLoop: do e = FEsolving_execElem(1),FEsolving_execElem(2) ! iterate call homogenization_averageTemperature(i,e) else terminallyIll = .true. + write(6,'(a48,i4,i4,/)') 'homogenization terminally-ill ',i,e exit elementLoop endif enddo diff --git a/code/numerics.config b/code/numerics.config index 12891256c..e5fcbc966 100644 --- a/code/numerics.config +++ b/code/numerics.config @@ -4,22 +4,33 @@ 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 (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 +pert_Fg 1.0e-7 # deformation gradient perturbation for grain tangent +pert_method 1 # perturbation method (1 = forward, 2 = backward or 3 = central) + +## crystallite numerical parameters ## 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 +subStepSizeCryst 0.25 # size of substep when cutback introduced in crystallite (value between 0 and 1) +stepIncreaseCryst 1.5 # increase of next substep size when previous substep converged in crystallite (value higher than 1) nState 50 # state loop limit nStress 200 # stress loop limit 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!) +## homogenization numerical parameters ## +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 +subStepSizeHomog 0.25 # size of substep when cutback introduced in homogenization (value between 0 and 1) +stepIncreaseHomog 1.5 # increase of next substep size when previous substep converged in homogenization (value higher than 1) +nMPstate 10 # materialpoint state loop limit + +## RGC scheme numerical parameters ## aTol_RGC 1.0e+4 # absolute tolerance of RGC residuum (in Pa) rTol_RGC 1.0e-3 # relative ... aMax_RGC 1.0e+10 # absolute upper-limit of RGC residuum (in Pa) -rMax_RGC 1.0e+2 # relative ... +rMax_RGC 1.0e+3 # relative ... perturbPenalty_RGC 1.0e-7 # perturbation for computing penalty tangent relevantMismatch_RGC 1.0e-5 # minimum threshold of mismatch + fixed_seed 1234 # put any number larger than zero, integer, if you want to have a pseudo random distribution diff --git a/code/numerics.f90 b/code/numerics.f90 index 0e9da150f..ef4773553 100644 --- a/code/numerics.f90 +++ b/code/numerics.f90 @@ -13,11 +13,16 @@ integer(pInt) iJacoStiffness, & ! freque nMPstate, & ! materialpoint state loop limit nCryst, & ! crystallite loop limit (only for debugging info, loop limit is determined by "subStepMinCryst") nState, & ! state loop limit - nStress ! stress loop limit + nStress, & ! stress loop limit + pert_method ! method used in perturbation technique for tangent real(pReal) relevantStrain, & ! strain increment considered significant pert_Fg, & ! strain perturbation for FEM Jacobi subStepMinCryst, & ! minimum (relative) size of sub-step allowed during cutback in crystallite subStepMinHomog, & ! minimum (relative) size of sub-step allowed during cutback in homogenization + subStepSizeCryst, & ! size of first substep when cutback in crystallite + subStepSizeHomog, & ! size of first substep when cutback in homogenization + stepIncreaseCryst, & ! increase of next substep size when previous substep converged in crystallite + stepIncreaseHomog, & ! increase of next substep size when previous substep converged 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 @@ -33,7 +38,7 @@ real(pReal) relevantStrain, & ! strain !* Random seeding parameters: added <<>> integer(pInt) fixedSeed ! fixed seeding for pseudo-random number generator - + CONTAINS !******************************************* @@ -75,12 +80,17 @@ subroutine numerics_init() relevantStrain = 1.0e-7_pReal iJacoStiffness = 1_pInt iJacoLpresiduum = 1_pInt - pert_Fg = 1.0e-6_pReal + pert_Fg = 1.0e-7_pReal + pert_method = 1 nHomog = 20_pInt subStepMinHomog = 1.0e-3_pReal + subStepSizeHomog = 0.25 + stepIncreaseHomog = 1.5 nMPstate = 10_pInt nCryst = 20_pInt subStepMinCryst = 1.0e-3_pReal + subStepsizeCryst = 0.25 + stepIncreaseCryst = 1.5 nState = 10_pInt nStress = 40_pInt rTol_crystalliteState = 1.0e-6_pReal @@ -121,6 +131,8 @@ subroutine numerics_init() iJacoLpresiduum = IO_intValue(line,positions,2) case ('pert_fg') pert_Fg = IO_floatValue(line,positions,2) + case ('pert_method') + pert_method = IO_intValue(line,positions,2) case ('nhomog') nHomog = IO_intValue(line,positions,2) case ('nmpstate') @@ -133,8 +145,16 @@ subroutine numerics_init() nStress = IO_intValue(line,positions,2) case ('substepmincryst') subStepMinCryst = IO_floatValue(line,positions,2) + case ('substepsizecryst') + subStepSizeCryst = IO_floatValue(line,positions,2) + case ('stepincreasecryst') + stepIncreaseCryst = IO_floatValue(line,positions,2) case ('substepminhomog') subStepMinHomog = IO_floatValue(line,positions,2) + case ('substepsizehomog') + subStepSizeHomog = IO_floatValue(line,positions,2) + case ('stepincreasehomog') + stepIncreaseHomog = IO_floatValue(line,positions,2) case ('rtol_crystallitestate') rTol_crystalliteState = IO_floatValue(line,positions,2) case ('rtol_crystallitetemperature') @@ -178,17 +198,25 @@ subroutine numerics_init() write(6,'(a24,x,i8)') 'iJacoStiffness: ',iJacoStiffness 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)') 'pert_method: ',pert_method write(6,'(a24,x,i8)') 'nCryst: ',nCryst write(6,'(a24,x,e8.1)') 'subStepMinCryst: ',subStepMinCryst + write(6,'(a24,x,e8.1)') 'subStepSizeCryst: ',subStepSizeCryst + write(6,'(a24,x,e8.1)') 'stepIncreaseCryst: ',stepIncreaseCryst write(6,'(a24,x,i8)') 'nState: ',nState write(6,'(a24,x,i8)') 'nStress: ',nStress 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 write(6,'(a24,x,e8.1)') 'aTol_crystalliteStress: ',aTol_crystalliteStress + write(6,*) + + write(6,'(a24,x,i8)') 'nHomog: ',nHomog + write(6,'(a24,x,e8.1)') 'subStepMinHomog: ',subStepMinHomog + write(6,'(a24,x,e8.1)') 'subStepSizeHomog: ',subStepSizeHomog + write(6,'(a24,x,e8.1)') 'stepIncreaseHomog: ',stepIncreaseHomog + write(6,'(a24,x,i8)') 'nMPstate: ',nMPstate + write(6,*) !* RGC parameters: added <<>> write(6,'(a24,x,e8.1)') 'aTol_RGC: ',absTol_RGC @@ -197,6 +225,7 @@ subroutine numerics_init() write(6,'(a24,x,e8.1)') 'rMax_RGC: ',relMax_RGC write(6,'(a24,x,e8.1)') 'perturbPenalty_RGC: ',pPert_RGC write(6,'(a24,x,e8.1)') 'relevantMismatch_RGC: ',xSmoo_RGC + write(6,*) !* Random seeding parameters: added <<>> write(6,'(a24,x,i8)') 'fixed_seed: ',fixedSeed @@ -207,13 +236,19 @@ subroutine numerics_init() if (iJacoStiffness < 1_pInt) call IO_error(261) if (iJacoLpresiduum < 1_pInt) call IO_error(262) if (pert_Fg <= 0.0_pReal) call IO_error(263) + if (pert_method <= 0_pInt .or. pert_method >= 4_pInt) & + call IO_error(299) if (nHomog < 1_pInt) call IO_error(264) if (nMPstate < 1_pInt) call IO_error(279) !! missing in IO !! 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 (subStepMinCryst <= 0.0_pReal) call IO_error(268) + if (subStepSizeCryst <= 0.0_pReal) call IO_error(268) + if (stepIncreaseCryst <= 0.0_pReal) call IO_error(268) if (subStepMinHomog <= 0.0_pReal) call IO_error(268) + if (subStepSizeHomog <= 0.0_pReal) call IO_error(268) + if (stepIncreaseHomog <= 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)