introduced a flexibility in cut-backing scheme in homogenization.f90 and in crystallite.f90:

(1) subStepSizeHomog and subStepSizeCryst := size of substep when cut-back is applied (initially was hard-coded).
(2) stepIncreaseHomog and stepIncreaseCryst := step increase when calculation for substep converge (was also hardcoded).

introduced a possibility to choose different finite difference scheme, i.e., forward-, backward- and central-difference, for computing grain numerical tangent. note that central-difference scheme will slow down the computation significantly. please use it only if necessary.

parameters to set these new features have been included in numerics.f90 and numerics.config, whereas corresponding error messages have been introduced in the IO.f90
This commit is contained in:
Denny Tjahjanto 2009-11-10 13:36:27 +00:00
parent 67f87486b1
commit cb88019aa6
5 changed files with 151 additions and 25 deletions

View File

@ -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)

View File

@ -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 ! <<added flexibility in cutback size>>
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)) ! <<introduce possibility for acceleration>>
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

View File

@ -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 ! <<added to adopt flexibility in cutback size>>
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)) ! <<introduce flexibility for step increase/acceleration>>
! 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
! <<modified to add more flexibility in 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

View File

@ -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

View File

@ -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 <<<updated 27.08.2009>>>
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 <<<updated 31.07.2009>>>
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 <<<updated 27.08.2009>>>
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)