polishing

This commit is contained in:
Martin Diehl 2020-03-14 19:11:26 +01:00
parent 33aaa94865
commit eb08f9f0b2
3 changed files with 53 additions and 53 deletions

View File

@ -19,7 +19,7 @@ module subroutine plastic_none_init
p, &
NipcMyPhase
write(6,'(/,a)') ' <<<+- plastic_'//PLASTICITY_NONE_label//' init -+>>>'; flush(6)
write(6,'(/,a)') ' <<<+- plastic_'//PLASTICITY_NONE_LABEL//' init -+>>>'; flush(6)
Ninstance = count(phase_plasticity == PLASTICITY_NONE_ID)
if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0) &

View File

@ -53,8 +53,8 @@ submodule(constitutive) plastic_nonlocal
atomicVolume, & !< atomic volume
Dsd0, & !< prefactor for self-diffusion coefficient
selfDiffusionEnergy, & !< activation enthalpy for diffusion
aTolRho, & !< absolute tolerance for dislocation density in state integration
aTolShear, & !< absolute tolerance for accumulated shear in state integration
atol_rho, & !< absolute tolerance for dislocation density in state integration
atol_gamma, & !< absolute tolerance for accumulated shear in state integration
significantRho, & !< density considered significant
significantN, & !< number of dislocations considered significant
doublekinkwidth, & !< width of a doubkle kink in multiples of the burgers vector length b
@ -182,7 +182,7 @@ module subroutine plastic_nonlocal_init
structure
integer :: NofMyPhase
write(6,'(/,a)') ' <<<+- constitutive_'//PLASTICITY_NONLOCAL_label//' init -+>>>'; flush(6)
write(6,'(/,a)') ' <<<+- constitutive_'//PLASTICITY_NONLOCAL_LABEL//' init -+>>>'; flush(6)
write(6,'(/,a)') ' Reuber et al., Acta Materialia 71:333348, 2014'
write(6,'(a)') ' https://doi.org/10.1016/j.actamat.2014.03.012'
@ -214,8 +214,8 @@ module subroutine plastic_nonlocal_init
dst => microstructure(phase_plasticityInstance(p)), &
config => config_phase(p))
prm%aTolRho = config%getFloat('atol_rho', defaultVal=0.0_pReal)
prm%aTolShear = config%getFloat('atol_shear', defaultVal=0.0_pReal)
prm%atol_rho = config%getFloat('atol_rho', defaultVal=0.0_pReal)
prm%atol_gamma = config%getFloat('atol_shear', defaultVal=0.0_pReal)
structure = config%getString('lattice_structure')
@ -351,8 +351,8 @@ module subroutine plastic_nonlocal_init
if (prm%significantN < 0.0_pReal) extmsg = trim(extmsg)//' significantN'
if (prm%significantrho < 0.0_pReal) extmsg = trim(extmsg)//' significantrho'
if (prm%atolshear <= 0.0_pReal) extmsg = trim(extmsg)//' atolshear'
if (prm%atolrho <= 0.0_pReal) extmsg = trim(extmsg)//' atolrho'
if (prm%atol_rho <= 0.0_pReal) extmsg = trim(extmsg)//' atol_rho'
if (prm%atol_gamma <= 0.0_pReal) extmsg = trim(extmsg)//' atol_gamma'
if (prm%CFLfactor < 0.0_pReal) extmsg = trim(extmsg)//' CFLfactor'
if (prm%p <= 0.0_pReal .or. prm%p > 1.0_pReal) extmsg = trim(extmsg)//' p'
@ -405,7 +405,7 @@ module subroutine plastic_nonlocal_init
stt%rho => plasticState(p)%state (0*prm%totalNslip+1:10*prm%totalNslip,:)
dot%rho => plasticState(p)%dotState (0*prm%totalNslip+1:10*prm%totalNslip,:)
del%rho => plasticState(p)%deltaState (0*prm%totalNslip+1:10*prm%totalNslip,:)
plasticState(p)%aTolState(1:10*prm%totalNslip) = prm%aTolRho
plasticState(p)%atolState(1:10*prm%totalNslip) = prm%atol_rho
stt%rhoSgl => plasticState(p)%state (0*prm%totalNslip+1: 8*prm%totalNslip,:)
dot%rhoSgl => plasticState(p)%dotState (0*prm%totalNslip+1: 8*prm%totalNslip,:)
@ -466,7 +466,7 @@ module subroutine plastic_nonlocal_init
stt%gamma => plasticState(p)%state (10*prm%totalNslip + 1:11*prm%totalNslip ,1:NofMyPhase)
dot%gamma => plasticState(p)%dotState (10*prm%totalNslip + 1:11*prm%totalNslip ,1:NofMyPhase)
del%gamma => plasticState(p)%deltaState (10*prm%totalNslip + 1:11*prm%totalNslip ,1:NofMyPhase)
plasticState(p)%aTolState(10*prm%totalNslip + 1:11*prm%totalNslip ) = prm%aTolShear
plasticState(p)%atolState(10*prm%totalNslip + 1:11*prm%totalNslip ) = prm%atol_gamma
plasticState(p)%slipRate => plasticState(p)%dotState (10*prm%totalNslip + 1:11*prm%totalNslip ,1:NofMyPhase)
stt%rho_forest => plasticState(p)%state (11*prm%totalNslip + 1:12*prm%totalNslip ,1:NofMyPhase)
@ -537,7 +537,7 @@ module subroutine plastic_nonlocal_init
enddo
enddo
if (iD(param(phase_plasticityInstance(p))%totalNslip,2,phase_plasticityInstance(p)) /= plasticState(p)%sizeState) &
call IO_error(0, ext_msg = 'state indices not properly set ('//PLASTICITY_NONLOCAL_label//')')
call IO_error(0, ext_msg = 'state indices not properly set ('//PLASTICITY_NONLOCAL_LABEL//')')
endif myPhase2
@ -1665,8 +1665,8 @@ module subroutine plastic_nonlocal_dotState(Mp, F, Fp, Temperature, &
#endif
if ( any(rho(:,mob) + rhoDot(1:ns,1:4) * timestep < -prm%aTolRho) &
.or. any(rho(:,dip) + rhoDot(1:ns,9:10) * timestep < -prm%aTolRho)) then
if ( any(rho(:,mob) + rhoDot(1:ns,1:4) * timestep < -prm%atol_rho) &
.or. any(rho(:,dip) + rhoDot(1:ns,9:10) * timestep < -prm%atol_rho)) then
#ifdef DEBUG
if (iand(debug_level(debug_constitutive),debug_levelExtensive) /= 0) then
write(6,'(a,i5,a,i2)') '<< CONST >> evolution rate leads to negative density at el ',el,' ip ',ip

View File

@ -91,9 +91,9 @@ module crystallite
subStepSizeLp, & !< size of first substep when cutback in Lp calculation
subStepSizeLi, & !< size of first substep when cutback in Li calculation
stepIncreaseCryst, & !< increase of next substep size when previous substep converged
rTol_crystalliteState, & !< relative tolerance in state loop
rTol_crystalliteStress, & !< relative tolerance in stress loop
aTol_crystalliteStress !< absolute tolerance in stress loop
rtol_crystalliteState, & !< relative tolerance in state loop
rtol_crystalliteStress, & !< relative tolerance in stress loop
atol_crystalliteStress !< absolute tolerance in stress loop
end type tNumerics
type(tNumerics) :: num ! numerics parameters. Better name?
@ -171,9 +171,9 @@ subroutine crystallite_init
num%subStepSizeLp = config_numerics%getFloat('substepsizelp', defaultVal=0.5_pReal)
num%subStepSizeLi = config_numerics%getFloat('substepsizeli', defaultVal=0.5_pReal)
num%rTol_crystalliteState = config_numerics%getFloat('rtol_crystallitestate', defaultVal=1.0e-6_pReal)
num%rTol_crystalliteStress = config_numerics%getFloat('rtol_crystallitestress',defaultVal=1.0e-6_pReal)
num%aTol_crystalliteStress = config_numerics%getFloat('atol_crystallitestress',defaultVal=1.0e-8_pReal)
num%rtol_crystalliteState = config_numerics%getFloat('rtol_crystallitestate', defaultVal=1.0e-6_pReal)
num%rtol_crystalliteStress = config_numerics%getFloat('rtol_crystallitestress',defaultVal=1.0e-6_pReal)
num%atol_crystalliteStress = config_numerics%getFloat('atol_crystallitestress',defaultVal=1.0e-8_pReal)
num%iJacoLpresiduum = config_numerics%getInt ('ijacolpresiduum', defaultVal=1)
@ -187,9 +187,9 @@ subroutine crystallite_init
if(num%subStepSizeLp <= 0.0_pReal) call IO_error(301,ext_msg='subStepSizeLp')
if(num%subStepSizeLi <= 0.0_pReal) call IO_error(301,ext_msg='subStepSizeLi')
if(num%rTol_crystalliteState <= 0.0_pReal) call IO_error(301,ext_msg='rTol_crystalliteState')
if(num%rTol_crystalliteStress <= 0.0_pReal) call IO_error(301,ext_msg='rTol_crystalliteStress')
if(num%aTol_crystalliteStress <= 0.0_pReal) call IO_error(301,ext_msg='aTol_crystalliteStress')
if(num%rtol_crystalliteState <= 0.0_pReal) call IO_error(301,ext_msg='rtol_crystalliteState')
if(num%rtol_crystalliteStress <= 0.0_pReal) call IO_error(301,ext_msg='rtol_crystalliteStress')
if(num%atol_crystalliteStress <= 0.0_pReal) call IO_error(301,ext_msg='atol_crystalliteStress')
if(num%iJacoLpresiduum < 1) call IO_error(301,ext_msg='iJacoLpresiduum')
@ -651,7 +651,7 @@ subroutine crystallite_results
integer :: p,o
real(pReal), allocatable, dimension(:,:,:) :: selected_tensors
type(rotation), allocatable, dimension(:) :: selected_rotations
character(len=pStringLen) :: group,lattice_label
character(len=pStringLen) :: group,structureLabel
do p=1,size(config_name_phase)
group = trim('current/constituent')//'/'//trim(config_name_phase(p))//'/generic'
@ -694,22 +694,22 @@ subroutine crystallite_results
'2nd Piola-Kirchoff stress','Pa')
case('orientation')
select case(lattice_structure(p))
case(LATTICE_iso_ID)
lattice_label = 'iso'
case(LATTICE_fcc_ID)
lattice_label = 'fcc'
case(LATTICE_bcc_ID)
lattice_label = 'bcc'
case(LATTICE_bct_ID)
lattice_label = 'bct'
case(LATTICE_hex_ID)
lattice_label = 'hex'
case(LATTICE_ort_ID)
lattice_label = 'ort'
case(lattice_ISO_ID)
structureLabel = 'iso'
case(lattice_FCC_ID)
structureLabel = 'fcc'
case(lattice_BCC_ID)
structureLabel = 'bcc'
case(lattice_BCT_ID)
structureLabel = 'bct'
case(lattice_HEX_ID)
structureLabel = 'hex'
case(lattice_ORT_ID)
structureLabel = 'ort'
end select
selected_rotations = select_rotations(crystallite_orientation,p)
call results_writeDataset(group,selected_rotations,'orientation',&
'crystal orientation as quaternion',lattice_label)
'crystal orientation as quaternion',structureLabel)
end select
enddo
enddo
@ -824,8 +824,8 @@ logical function integrateStress(ipc,ip,el,timeFraction)
real(pReal) steplengthLp, &
steplengthLi, &
dt, & ! time increment
aTolLp, &
aTolLi, &
atol_Lp, &
atol_Li, &
devNull
integer NiterationStressLp, & ! number of stress integrations
NiterationStressLi, & ! number of inner stress integrations
@ -891,13 +891,13 @@ logical function integrateStress(ipc,ip,el,timeFraction)
S, Fi_new, ipc, ip, el)
!* update current residuum and check for convergence of loop
aTolLp = max(num%rTol_crystalliteStress * max(norm2(Lpguess),norm2(Lp_constitutive)), & ! absolute tolerance from largest acceptable relative error
num%aTol_crystalliteStress) ! minimum lower cutoff
atol_Lp = max(num%rtol_crystalliteStress * max(norm2(Lpguess),norm2(Lp_constitutive)), & ! absolute tolerance from largest acceptable relative error
num%atol_crystalliteStress) ! minimum lower cutoff
residuumLp = Lpguess - Lp_constitutive
if (any(IEEE_is_NaN(residuumLp))) then
return ! error
elseif (norm2(residuumLp) < aTolLp) then ! converged if below absolute tolerance
elseif (norm2(residuumLp) < atol_Lp) then ! converged if below absolute tolerance
exit LpLoop
elseif (NiterationStressLp == 1 .or. norm2(residuumLp) < norm2(residuumLp_old)) then ! not converged, but improved norm of residuum (always proceed in first iteration)...
residuumLp_old = residuumLp ! ...remember old values and...
@ -934,12 +934,12 @@ logical function integrateStress(ipc,ip,el,timeFraction)
S, Fi_new, ipc, ip, el)
!* update current residuum and check for convergence of loop
aTolLi = max(num%rTol_crystalliteStress * max(norm2(Liguess),norm2(Li_constitutive)), & ! absolute tolerance from largest acceptable relative error
num%aTol_crystalliteStress) ! minimum lower cutoff
atol_Li = max(num%rtol_crystalliteStress * max(norm2(Liguess),norm2(Li_constitutive)), & ! absolute tolerance from largest acceptable relative error
num%atol_crystalliteStress) ! minimum lower cutoff
residuumLi = Liguess - Li_constitutive
if (any(IEEE_is_NaN(residuumLi))) then
return ! error
elseif (norm2(residuumLi) < aTolLi) then ! converged if below absolute tolerance
elseif (norm2(residuumLi) < atol_Li) then ! converged if below absolute tolerance
exit LiLoop
elseif (NiterationStressLi == 1 .or. norm2(residuumLi) < norm2(residuumLi_old)) then ! not converged, but improved norm of residuum (always proceed in first iteration)...
residuumLi_old = residuumLi ! ...remember old values and...
@ -1089,7 +1089,7 @@ subroutine integrateStateFPI
crystallite_converged(g,i,e) = converged(residuum_plastic(1:sizeDotState), &
plasticState(p)%state(1:sizeDotState,c), &
plasticState(p)%aTolState(1:sizeDotState))
plasticState(p)%atolState(1:sizeDotState))
do s = 1, phase_Nsources(p)
@ -1113,7 +1113,7 @@ subroutine integrateStateFPI
crystallite_converged(g,i,e) = &
crystallite_converged(g,i,e) .and. converged(residuum_source(1:sizeDotState), &
sourceState(p)%p(s)%state(1:sizeDotState,c), &
sourceState(p)%p(s)%aTolState(1:sizeDotState))
sourceState(p)%p(s)%atolState(1:sizeDotState))
enddo
endif
enddo; enddo; enddo
@ -1269,7 +1269,7 @@ subroutine integrateStateAdaptiveEuler
crystallite_converged(g,i,e) = converged(residuum_plastic(1:sizeDotState,g,i,e), &
plasticState(p)%state(1:sizeDotState,c), &
plasticState(p)%aTolState(1:sizeDotState))
plasticState(p)%atolState(1:sizeDotState))
do s = 1, phase_Nsources(p)
sizeDotState = sourceState(p)%p(s)%sizeDotState
@ -1280,7 +1280,7 @@ subroutine integrateStateAdaptiveEuler
crystallite_converged(g,i,e) = &
crystallite_converged(g,i,e) .and. converged(residuum_source(1:sizeDotState,s,g,i,e), &
sourceState(p)%p(s)%state(1:sizeDotState,c), &
sourceState(p)%p(s)%aTolState(1:sizeDotState))
sourceState(p)%p(s)%atolState(1:sizeDotState))
enddo
endif
@ -1497,7 +1497,7 @@ subroutine integrateStateRKCK45
crystallite_todo(g,i,e) = converged(residuum_plastic(1:sizeDotState,g,i,e), &
plasticState(p)%state(1:sizeDotState,cc), &
plasticState(p)%aTolState(1:sizeDotState))
plasticState(p)%atolState(1:sizeDotState))
do s = 1, phase_Nsources(p)
sizeDotState = sourceState(p)%p(s)%sizeDotState
@ -1505,7 +1505,7 @@ subroutine integrateStateRKCK45
crystallite_todo(g,i,e) = &
crystallite_todo(g,i,e) .and. converged(residuum_source(1:sizeDotState,s,g,i,e), &
sourceState(p)%p(s)%state(1:sizeDotState,cc), &
sourceState(p)%p(s)%aTolState(1:sizeDotState))
sourceState(p)%p(s)%atolState(1:sizeDotState))
enddo
endif
enddo; enddo; enddo
@ -1558,16 +1558,16 @@ end subroutine setConvergenceFlag
!--------------------------------------------------------------------------------------------------
!> @brief determines whether a point is converged
!--------------------------------------------------------------------------------------------------
logical pure function converged(residuum,state,aTol)
logical pure function converged(residuum,state,atol)
real(pReal), intent(in), dimension(:) ::&
residuum, state, aTol
residuum, state, atol
real(pReal) :: &
rTol
rTol = num%rTol_crystalliteState
converged = all(abs(residuum) <= max(aTol, rTol*abs(state)))
converged = all(abs(residuum) <= max(atol, rtol*abs(state)))
end function converged