atol=0.0 is fine

crystallite takes max(atol, rtol*X), so atol=0.0 means that convergence
is based on rtol only
This commit is contained in:
Martin Diehl 2020-03-16 00:27:52 +01:00
parent e3bbd32b1e
commit 23c6510faa
12 changed files with 39 additions and 53 deletions

@ -1 +1 @@
Subproject commit ffb01b351ab907b61885560f452f7f0970a1a744
Subproject commit cdc06984fa68a9873d2de752dfe99e720e18a4f1

View File

@ -3,23 +3,16 @@
elasticity hooke
plasticity dislotwin
#(output) edge_density
#(output) dipole_density
#(output) shearrate_slip
#(output) accumulated_shear_slip
#(output) mfp_slip
#(output) resolved_stress_slip
#(output) threshold_stress_slip
#(output) twin_fraction
#(output) shearrate_twin
#(output) accumulated_shear_twin
#(output) mfp_twin
#(output) resolved_stress_twin
#(output) threshold_stress_twin
#(output) shearrate_shearband
#(output) resolved_stress_shearband
#(output) sb_eigenvalues
#(output) sb_eigenvectors
(output) rho_mob
(output) rho_dip
(output) gamma_sl
(output) lambda_sl
(output) tau_pass
(output) f_tw
(output) lambda_tw
(output) tau_hat_tw
(output) f_tr
### Material parameters ###
lattice_structure fcc
@ -45,7 +38,6 @@ D0 4.0e-5 # Vacancy diffusion prefactor [m**2/s]
Qsd 4.5e-19 # Activation energy for climb [J]
Catomicvolume 1.0 # Adj. parameter controlling the atomic volume [in b^3]
Cedgedipmindistance 1.0 # Adj. parameter controlling the minimum dipole distance [in b]
atol_rho 1.0
interactionSlipSlip 0.122 0.122 0.625 0.07 0.137 0.122 # Interaction coefficients (Kubin et al. 2008)
### Shearband parameters ###
@ -68,6 +60,5 @@ Cmfptwin 1.0 # Adj. parameter controlling twin mean free
Cthresholdtwin 1.0 # Adj. parameter controlling twin threshold stress
interactionSlipTwin 0.0 1.0 1.0 # Dislocation-Twin interaction coefficients
interactionTwinTwin 0.0 1.0 # Twin-Twin interaction coefficients
atol_twinFrac 1.0e-7
SFE_0K -0.0396 # stacking fault energy at zero K; TWIP steel: -0.0526; Cu: -0.0396
dSFE_dT 0.0002 # temperature dependance of stacking fault energy

View File

@ -41,8 +41,6 @@ interaction_twinslip 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
####################################################
# open for discussion
####################################################
atol_shear 1e-6
atol_twinfrac 1e-6
n_twin 20
n_slip 20

View File

@ -218,15 +218,15 @@ module subroutine plastic_disloUCLA_init
stt%rho_mob=>plasticState(p)%state(startIndex:endIndex,:)
stt%rho_mob= spread(prm%rho_mob_0,2,NipcMyPhase)
dot%rho_mob=>plasticState(p)%dotState(startIndex:endIndex,:)
plasticState(p)%atol(startIndex:endIndex) = config%getFloat('atol_rho')
if (any(plasticState(p)%atol(startIndex:endIndex) <= 0.0_pReal)) extmsg = trim(extmsg)//' atol_rho'
plasticState(p)%atol(startIndex:endIndex) = config%getFloat('atol_rho',defaultVal=1.0_pReal)
if (any(plasticState(p)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_rho'
startIndex = endIndex + 1
endIndex = endIndex + prm%sum_N_sl
stt%rho_dip=>plasticState(p)%state(startIndex:endIndex,:)
stt%rho_dip= spread(prm%rho_dip_0,2,NipcMyPhase)
dot%rho_dip=>plasticState(p)%dotState(startIndex:endIndex,:)
plasticState(p)%atol(startIndex:endIndex) = config%getFloat('atol_rho')
plasticState(p)%atol(startIndex:endIndex) = config%getFloat('atol_rho',defaultVal=1.0_pReal)
startIndex = endIndex + 1
endIndex = endIndex + prm%sum_N_sl
@ -401,7 +401,7 @@ module subroutine plastic_disloUCLA_results(instance,group)
if(prm%sum_N_sl>0) call results_writeDataset(group,stt%rho_dip,'rho_dip',&
'dislocation dipole density''1/m²')
case('shear_rate_slip') ! should be gamma
if(prm%sum_N_sl>0) call results_writeDataset(group,stt%gamma_sl,'dot_gamma_sl',& ! this is not dot!!
if(prm%sum_N_sl>0) call results_writeDataset(group,stt%gamma_sl,'dot_gamma_sl',& ! this is not dot!!
'plastic shear','1')
case('mfp_slip') !ToDo: should be Lambda
if(prm%sum_N_sl>0) call results_writeDataset(group,dst%Lambda_sl,'Lambda_sl',&

View File

@ -173,9 +173,9 @@ module subroutine plastic_dislotwin_init
dst => dependentState(phase_plasticityInstance(p)), &
config => config_phase(p))
prm%aTol_rho = config%getFloat('atol_rho', defaultVal=0.0_pReal)
prm%aTol_f_tw = config%getFloat('atol_twinfrac', defaultVal=0.0_pReal)
prm%aTol_f_tr = config%getFloat('atol_transfrac', defaultVal=0.0_pReal)
prm%aTol_rho = config%getFloat('atol_rho', defaultVal=1.0_pReal)
prm%aTol_f_tw = config%getFloat('atol_twinfrac', defaultVal=1.0e-7_pReal)
prm%aTol_f_tr = config%getFloat('atol_transfrac', defaultVal=1.0e-6_pReal)
prm%output = config%getStrings('(output)', defaultVal=emptyStringArray)
@ -399,13 +399,13 @@ module subroutine plastic_dislotwin_init
if (any(prm%atomicVolume <= 0.0_pReal)) &
call IO_error(211,el=p,ext_msg='cAtomicVolume ('//PLASTICITY_DISLOTWIN_LABEL//')')
if (prm%sum_N_tw > 0) then
if (prm%aTol_rho <= 0.0_pReal) &
if (prm%aTol_rho < 0.0_pReal) &
call IO_error(211,el=p,ext_msg='aTol_rho ('//PLASTICITY_DISLOTWIN_label//')')
if (prm%aTol_f_tw <= 0.0_pReal) &
if (prm%aTol_f_tw < 0.0_pReal) &
call IO_error(211,el=p,ext_msg='aTol_f_tw ('//PLASTICITY_DISLOTWIN_label//')')
endif
if (prm%sum_N_tr > 0) then
if (prm%aTol_f_tr <= 0.0_pReal) &
if (prm%aTol_f_tr < 0.0_pReal) &
call IO_error(211,el=p,ext_msg='aTol_f_tr ('//PLASTICITY_DISLOTWIN_label//')')
endif

View File

@ -125,12 +125,12 @@ module subroutine plastic_isotropic_init
stt%xi = prm%xi_0
dot%xi => plasticState(p)%dotState(1,:)
plasticState(p)%atol(1) = config%getFloat('atol_flowstress',defaultVal=1.0_pReal)
if (plasticState(p)%atol(1) <= 0.0_pReal) extmsg = trim(extmsg)//' atol_xi'
if (plasticState(p)%atol(1) < 0.0_pReal) extmsg = trim(extmsg)//' atol_xi'
stt%gamma => plasticState(p)%state (2,:)
dot%gamma => plasticState(p)%dotState(2,:)
plasticState(p)%atol(2) = config%getFloat('atol_shear',defaultVal=1.0e-6_pReal)
if (plasticState(p)%atol(2) <= 0.0_pReal) extmsg = trim(extmsg)//' atol_gamma'
if (plasticState(p)%atol(2) < 0.0_pReal) extmsg = trim(extmsg)//' atol_gamma'
! global alias
plasticState(p)%slipRate => plasticState(p)%dotState(2:2,:)

View File

@ -170,7 +170,7 @@ module subroutine plastic_kinehardening_init
stt%crss = spread(prm%crss0, 2, NipcMyPhase)
dot%crss => plasticState(p)%dotState(startIndex:endIndex,:)
plasticState(p)%atol(startIndex:endIndex) = config%getFloat('atol_xi',defaultVal=1.0_pReal)
if(any(plasticState(p)%atol(startIndex:endIndex) <= 0.0_pReal)) extmsg = trim(extmsg)//' atol_xi'
if(any(plasticState(p)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_xi'
startIndex = endIndex + 1
endIndex = endIndex + prm%sum_N_sl
@ -183,7 +183,7 @@ module subroutine plastic_kinehardening_init
stt%accshear => plasticState(p)%state (startIndex:endIndex,:)
dot%accshear => plasticState(p)%dotState(startIndex:endIndex,:)
plasticState(p)%atol(startIndex:endIndex) = config%getFloat('atol_gamma',defaultVal=1.0e-6_pReal)
if(any(plasticState(p)%atol(startIndex:endIndex) <= 0.0_pReal)) extmsg = trim(extmsg)//' atol_gamma'
if(any(plasticState(p)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_gamma'
! global alias
plasticState(p)%slipRate => plasticState(p)%dotState(startIndex:endIndex,:)

View File

@ -207,8 +207,8 @@ module subroutine plastic_phenopowerlaw_init
!--------------------------------------------------------------------------------------------------
! allocate state arrays
NipcMyPhase = count(material_phaseAt == p) * discretization_nIP
sizeDotState = size(['tau_slip ','gamma_slip']) * prm%sum_N_sl &
+ size(['tau_twin ','gamma_twin']) * prm%sum_N_tw
sizeDotState = size(['xi_sl ','gamma_sl']) * prm%sum_N_sl &
+ size(['xi_tw ','gamma_tw']) * prm%sum_N_tw
sizeState = sizeDotState
call material_allocatePlasticState(p,NipcMyPhase,sizeState,sizeDotState,0)
@ -220,31 +220,32 @@ module subroutine plastic_phenopowerlaw_init
stt%xi_slip => plasticState(p)%state (startIndex:endIndex,:)
stt%xi_slip = spread(prm%xi_slip_0, 2, NipcMyPhase)
dot%xi_slip => plasticState(p)%dotState(startIndex:endIndex,:)
plasticState(p)%atol(startIndex:endIndex) = config%getFloat('atol_resistance',defaultVal=1.0_pReal)
if(any(plasticState(p)%atol(startIndex:endIndex)<=0.0_pReal)) extmsg = trim(extmsg)//' atol_xi'
plasticState(p)%atol(startIndex:endIndex) = config%getFloat('atol_xi',defaultVal=1.0_pReal)
if(any(plasticState(p)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_xi'
startIndex = endIndex + 1
endIndex = endIndex + prm%sum_N_tw
stt%xi_twin => plasticState(p)%state (startIndex:endIndex,:)
stt%xi_twin = spread(prm%xi_twin_0, 2, NipcMyPhase)
dot%xi_twin => plasticState(p)%dotState(startIndex:endIndex,:)
plasticState(p)%atol(startIndex:endIndex) = config%getFloat('atol_resistance',defaultVal=1.0_pReal)
plasticState(p)%atol(startIndex:endIndex) = config%getFloat('atol_xi',defaultVal=1.0_pReal)
if(any(plasticState(p)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_xi'
startIndex = endIndex + 1
endIndex = endIndex + prm%sum_N_sl
stt%gamma_slip => plasticState(p)%state (startIndex:endIndex,:)
dot%gamma_slip => plasticState(p)%dotState(startIndex:endIndex,:)
plasticState(p)%atol(startIndex:endIndex) = config%getFloat('atol_shear',defaultVal=1.0e-6_pReal)
if(any(plasticState(p)%atol(startIndex:endIndex)<=0.0_pReal)) extmsg = trim(extmsg)//' atol_gamma_slip'
plasticState(p)%atol(startIndex:endIndex) = config%getFloat('atol_gamma',defaultVal=1.0e-6_pReal)
if(any(plasticState(p)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_gamma'
! global alias
plasticState(p)%slipRate => plasticState(p)%dotState(startIndex:endIndex,:)
plasticState(p)%slipRate => plasticState(p)%dotState(startIndex:endIndex,:)
startIndex = endIndex + 1
endIndex = endIndex + prm%sum_N_tw
stt%gamma_twin => plasticState(p)%state (startIndex:endIndex,:)
dot%gamma_twin => plasticState(p)%dotState(startIndex:endIndex,:)
plasticState(p)%atol(startIndex:endIndex) = config%getFloat('atol_gamma',defaultVal=1.0e-6_pReal)
if(any(plasticState(p)%atol(startIndex:endIndex)<=0.0_pReal)) extmsg = trim(extmsg)//' atol_gamma'
if(any(plasticState(p)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_gamma'
plasticState(p)%state0 = plasticState(p)%state ! ToDo: this could be done centrally

View File

@ -110,8 +110,7 @@ subroutine source_damage_anisoBrittle_init
NofMyPhase = count(material_phaseAt==p) * discretization_nIP
call material_allocateSourceState(p,sourceOffset,NofMyPhase,1,1,0)
sourceState(p)%p(sourceOffset)%atol = config%getFloat('anisobrittle_atol',defaultVal=1.0e-3_pReal)
if(any(sourceState(p)%p(sourceOffset)%atol <= 0.0_pReal)) &
extmsg = trim(extmsg)//' anisobrittle_atol'
if(any(sourceState(p)%p(sourceOffset)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' anisobrittle_atol'
end associate

View File

@ -97,8 +97,7 @@ subroutine source_damage_anisoDuctile_init
NofMyPhase=count(material_phaseAt==p) * discretization_nIP
call material_allocateSourceState(p,sourceOffset,NofMyPhase,1,1,0)
sourceState(p)%p(sourceOffset)%atol = config%getFloat('anisoductile_atol',defaultVal=1.0e-3_pReal)
if(any(sourceState(p)%p(sourceOffset)%atol <=0.0_pReal)) &
extmsg = trim(extmsg)//' anisoductile_atol'
if(any(sourceState(p)%p(sourceOffset)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' anisoductile_atol'
end associate

View File

@ -85,8 +85,7 @@ subroutine source_damage_isoBrittle_init
NofMyPhase = count(material_phaseAt==p) * discretization_nIP
call material_allocateSourceState(p,sourceOffset,NofMyPhase,1,1,1)
sourceState(p)%p(sourceOffset)%atol = config%getFloat('isobrittle_atol',defaultVal=1.0e-3_pReal)
if(any(sourceState(p)%p(sourceOffset)%atol <= 0.0_pReal)) &
extmsg = trim(extmsg)//' isobrittle_atol'
if(any(sourceState(p)%p(sourceOffset)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' isobrittle_atol'
end associate

View File

@ -84,8 +84,7 @@ subroutine source_damage_isoDuctile_init
NofMyPhase=count(material_phaseAt==p) * discretization_nIP
call material_allocateSourceState(p,sourceOffset,NofMyPhase,1,1,0)
sourceState(p)%p(sourceOffset)%atol = config%getFloat('isoductile_atol',defaultVal=1.0e-3_pReal)
if(any(sourceState(p)%p(sourceOffset)%atol <= 0.0_pReal)) &
extmsg = trim(extmsg)//' isoductile_atol'
if(any(sourceState(p)%p(sourceOffset)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' isoductile_atol'
end associate