Merge branch 'no-defaults-phenopowerlaw' into 'development'

No defaults phenopowerlaw

See merge request damask/DAMASK!863
This commit is contained in:
Sharan Roongta 2023-12-07 13:19:00 +00:00
commit 78de3e6bde
3 changed files with 79 additions and 46 deletions

View File

@ -21,7 +21,9 @@ module misc
misc_init, &
misc_selfTest, &
misc_optional, &
misc_prefixOptions
misc_prefixOptions, &
misc_ones, &
misc_zeros
contains
@ -136,6 +138,34 @@ pure function misc_prefixOptions(string,prefix) result(prefixed)
end function misc_prefixOptions
!--------------------------------------------------------------------------------------------------
!> @brief 1D array of zeros.
!--------------------------------------------------------------------------------------------------
pure function misc_zeros(N)
integer, intent(in) :: N !< number of zeros
real(pREAL), dimension(N) :: misc_zeros
misc_zeros = 0._pREAL
end function misc_zeros
!--------------------------------------------------------------------------------------------------
!> @brief 1D array of ones.
!--------------------------------------------------------------------------------------------------
pure function misc_ones(N)
integer, intent(in) :: N !< number of ones
real(pREAL), dimension(N) :: misc_ones
misc_ones = 1._pREAL
end function misc_ones
!--------------------------------------------------------------------------------------------------
!> @brief Check correctness of some misc functions.
!--------------------------------------------------------------------------------------------------
@ -143,7 +173,7 @@ subroutine misc_selfTest()
real(pREAL) :: r
character(len=:), allocatable :: str,out
integer :: N
call random_number(r)
if (test_str('DAMASK') /= 'DAMASK') error stop 'optional_str, present'
@ -163,6 +193,12 @@ subroutine misc_selfTest()
out=misc_prefixOptions(str,'p_')
if (out /= '-p_a -1 -p_more 123 -p_flag -') error stop 'misc_prefixOptions'
N = int(r*99._pReal)
if (size(misc_zeros(N)) /= N) error stop 'shape zeros'
if (size(misc_ones(N)) /= N) error stop 'shape ones'
if (any(dNeq(misc_zeros(N),0.0_pReal))) error stop 'value zeros'
if (any(dNeq(misc_ones(N),1.0_pReal))) error stop 'value ones'
contains
function test_str(str_in) result(str_out)

View File

@ -74,7 +74,6 @@ module function plastic_kinehardening_init() result(myPlasticity)
logical, dimension(:), allocatable :: myPlasticity
integer :: &
ph, o, &
N_fam, &
Nmembers, &
sizeState, sizeDeltaState, sizeDotState, &
startIndex, endIndex
@ -138,7 +137,6 @@ module function plastic_kinehardening_init() result(myPlasticity)
N_sl = pl%get_as1dInt('N_sl',defaultVal=emptyIntArray)
prm%sum_N_sl = sum(abs(N_sl))
slipActive: if (prm%sum_N_sl > 0) then
N_fam = size(N_sl)
prm%systems_sl = crystal_labels_slip(N_sl,phase_lattice(ph))
prm%P = crystal_SchmidMatrix_slip(N_sl,phase_lattice(ph),phase_cOverA(ph))
@ -155,15 +153,15 @@ module function plastic_kinehardening_init() result(myPlasticity)
prm%h_sl_sl = crystal_interaction_SlipBySlip(N_sl,pl%get_as1dReal('h_sl-sl'), &
phase_lattice(ph))
xi_0 = math_expand(pl%get_as1dReal('xi_0', requiredSize=N_fam),N_sl)
prm%dot_gamma_0 = math_expand(pl%get_as1dReal('dot_gamma_0', requiredSize=N_fam),N_sl)
prm%n = math_expand(pl%get_as1dReal('n', requiredSize=N_fam),N_sl)
prm%xi_inf = math_expand(pl%get_as1dReal('xi_inf', requiredSize=N_fam),N_sl)
prm%chi_inf = math_expand(pl%get_as1dReal('chi_inf', requiredSize=N_fam),N_sl)
prm%h_0_xi = math_expand(pl%get_as1dReal('h_0_xi', requiredSize=N_fam),N_sl)
prm%h_0_chi = math_expand(pl%get_as1dReal('h_0_chi', requiredSize=N_fam),N_sl)
prm%h_inf_xi = math_expand(pl%get_as1dReal('h_inf_xi', requiredSize=N_fam),N_sl)
prm%h_inf_chi = math_expand(pl%get_as1dReal('h_inf_chi', requiredSize=N_fam),N_sl)
xi_0 = math_expand(pl%get_as1dReal('xi_0', requiredSize=size(N_sl)),N_sl)
prm%dot_gamma_0 = math_expand(pl%get_as1dReal('dot_gamma_0', requiredSize=size(N_sl)),N_sl)
prm%n = math_expand(pl%get_as1dReal('n', requiredSize=size(N_sl)),N_sl)
prm%xi_inf = math_expand(pl%get_as1dReal('xi_inf', requiredSize=size(N_sl)),N_sl)
prm%chi_inf = math_expand(pl%get_as1dReal('chi_inf', requiredSize=size(N_sl)),N_sl)
prm%h_0_xi = math_expand(pl%get_as1dReal('h_0_xi', requiredSize=size(N_sl)),N_sl)
prm%h_0_chi = math_expand(pl%get_as1dReal('h_0_chi', requiredSize=size(N_sl)),N_sl)
prm%h_inf_xi = math_expand(pl%get_as1dReal('h_inf_xi', requiredSize=size(N_sl)),N_sl)
prm%h_inf_chi = math_expand(pl%get_as1dReal('h_inf_chi', requiredSize=size(N_sl)),N_sl)
!--------------------------------------------------------------------------------------------------
! sanity checks

View File

@ -77,7 +77,6 @@ module function plastic_phenopowerlaw_init() result(myPlasticity)
logical, dimension(:), allocatable :: myPlasticity
integer :: &
ph, i, &
N_fam, &
Nmembers, &
sizeState, sizeDotState, &
startIndex, endIndex
@ -86,8 +85,7 @@ module function plastic_phenopowerlaw_init() result(myPlasticity)
N_tw !< number of twin-systems for a given twin family
real(pREAL), dimension(:), allocatable :: &
xi_0_sl, & !< initial critical shear stress for slip
xi_0_tw, & !< initial critical shear stress for twin
ones
xi_0_tw !< initial critical shear stress for twin
real(pREAL), dimension(:,:), allocatable :: &
a_nS !< non-Schmid coefficients
character(len=:), allocatable :: &
@ -142,17 +140,18 @@ module function plastic_phenopowerlaw_init() result(myPlasticity)
!--------------------------------------------------------------------------------------------------
! slip related parameters
slipActive: if (prm%sum_N_sl > 0) then
N_fam = size(N_sl)
ones = [(1.0_pREAL,i=1,N_fam)]
prm%dot_gamma_0_sl = math_expand(pl%get_as1dReal('dot_gamma_0_sl',requiredSize=N_fam,defaultVal= ones), N_sl)
prm%n_sl = math_expand(pl%get_as1dReal('n_sl', requiredSize=N_fam,defaultVal= ones), N_sl)
prm%a_sl = math_expand(pl%get_as1dReal('a_sl', requiredSize=N_fam,defaultVal= ones), N_sl)
prm%h_0_sl_sl = math_expand(pl%get_as1dReal('h_0_sl-sl', requiredSize=N_fam,defaultVal=0.0_pREAL*ones), N_sl)
prm%c_1 = math_expand(pl%get_as1dReal('c_1', requiredSize=N_fam,defaultVal=0.0_pREAL*ones), N_sl)
prm%c_2 = math_expand(pl%get_as1dReal('c_2', requiredSize=N_fam,defaultVal= ones), N_sl)
xi_0_sl = math_expand(pl%get_as1dReal('xi_0_sl', requiredSize=N_fam), N_sl)
prm%xi_inf_sl = math_expand(pl%get_as1dReal('xi_inf_sl', requiredSize=N_fam), N_sl)
prm%f_sat_sl_tw = math_expand(pl%get_as1dReal('f_sat_sl-tw',requiredSize=N_fam,defaultVal=0.0_pREAL*ones), N_sl)
prm%dot_gamma_0_sl = math_expand(pl%get_as1dReal('dot_gamma_0_sl',requiredSize=size(N_sl)), N_sl)
prm%n_sl = math_expand(pl%get_as1dReal('n_sl', requiredSize=size(N_sl)), N_sl)
prm%a_sl = math_expand(pl%get_as1dReal('a_sl', requiredSize=size(N_sl)), N_sl)
prm%h_0_sl_sl = math_expand(pl%get_as1dReal('h_0_sl-sl', requiredSize=size(N_sl)), N_sl)
xi_0_sl = math_expand(pl%get_as1dReal('xi_0_sl', requiredSize=size(N_sl)), N_sl)
prm%xi_inf_sl = math_expand(pl%get_as1dReal('xi_inf_sl', requiredSize=size(N_sl)), N_sl)
prm%c_1 = math_expand(pl%get_as1dReal('c_1', requiredSize=size(N_sl), &
defaultVal=misc_zeros(size(N_sl))), N_sl)
prm%c_2 = math_expand(pl%get_as1dReal('c_2', requiredSize=size(N_sl), &
defaultVal=misc_ones(size(N_sl))), N_sl)
prm%f_sat_sl_tw = math_expand(pl%get_as1dReal('f_sat_sl-tw', requiredSize=size(N_sl), &
defaultVal=misc_zeros(size(N_sl))), N_sl)
prm%h_sl_sl = crystal_interaction_SlipBySlip(N_sl,pl%get_as1dReal('h_sl-sl'),phase_lattice(ph))
@ -179,14 +178,14 @@ module function plastic_phenopowerlaw_init() result(myPlasticity)
else slipActive
xi_0_sl = emptyRealArray
allocate(prm%dot_gamma_0_sl, & !< reference shear strain rate for slip
allocate(prm%dot_gamma_0_sl, &
prm%a_sl, &
prm%n_sl, & !< stress exponent for slip
prm%xi_inf_sl, & !< maximum critical shear stress for slip
prm%f_sat_sl_tw, & !< push-up factor for slip saturation due to twinning
prm%n_sl, &
prm%xi_inf_sl, &
prm%f_sat_sl_tw, &
prm%c_1, &
prm%c_2, &
prm%h_0_sl_sl, & !< reference hardening slip - slip
prm%h_0_sl_sl, &
source=emptyRealArray)
allocate(prm%h_sl_sl(0,0))
end if slipActive
@ -194,14 +193,14 @@ module function plastic_phenopowerlaw_init() result(myPlasticity)
!--------------------------------------------------------------------------------------------------
! twin related parameters
twinActive: if (prm%sum_N_tw > 0) then
N_fam = size(N_tw)
ones = [(1.0_pREAL,i=1,N_fam)]
prm%dot_gamma_0_tw = math_expand(pl%get_as1dReal('dot_gamma_0_tw', requiredSize=N_fam,defaultVal=ones), N_tw)
prm%n_tw = math_expand(pl%get_as1dReal('n_tw', requiredSize=N_fam,defaultVal= ones), N_tw)
prm%c_3 = math_expand(pl%get_as1dReal('c_3', requiredSize=N_fam,defaultVal=0.0_pREAL*ones), N_tw)
prm%c_4 = math_expand(pl%get_as1dReal('c_4', requiredSize=N_fam,defaultVal=0.0_pREAL*ones), N_tw)
prm%h_0_tw_tw = math_expand(pl%get_as1dReal('h_0_tw-tw', requiredSize=N_fam,defaultVal=0.0_pReal*ones), N_tw)
xi_0_tw = math_expand(pl%get_as1dReal('xi_0_tw', requiredSize=N_fam), N_tw)
prm%dot_gamma_0_tw = math_expand(pl%get_as1dReal('dot_gamma_0_tw', requiredSize=size(N_tw)), N_tw)
prm%n_tw = math_expand(pl%get_as1dReal('n_tw', requiredSize=size(N_tw)), N_tw)
prm%h_0_tw_tw = math_expand(pl%get_as1dReal('h_0_tw-tw', requiredSize=size(N_tw)), N_tw)
xi_0_tw = math_expand(pl%get_as1dReal('xi_0_tw', requiredSize=size(N_tw)), N_tw)
prm%c_3 = math_expand(pl%get_as1dReal('c_3', requiredSize=size(N_tw), &
defaultVal=misc_ones(size(N_tw))), N_tw)
prm%c_4 = math_expand(pl%get_as1dReal('c_4', requiredSize=size(N_tw), &
defaultVal=misc_zeros(size(N_tw))), N_tw)
prm%gamma_char = crystal_characteristicShear_twin(N_tw,phase_lattice(ph),phase_cOverA(ph))
prm%h_tw_tw = crystal_interaction_TwinByTwin(N_tw,pl%get_as1dReal('h_tw-tw'),phase_lattice(ph))
@ -216,13 +215,13 @@ module function plastic_phenopowerlaw_init() result(myPlasticity)
else twinActive
xi_0_tw = emptyRealArray
allocate(prm%dot_gamma_0_tw, & !< reference shear strain rate for twin
prm%n_tw, & !< stress exponent for twin
allocate(prm%dot_gamma_0_tw, &
prm%n_tw, &
prm%c_3, &
prm%c_4, &
prm%gamma_char, & !< characteristic shear for twins
prm%h_0_tw_sl, & !< reference hardening twin - slip
prm%h_0_tw_tw, & !< reference hardening twin - twin
prm%gamma_char, &
prm%h_0_tw_sl, &
prm%h_0_tw_tw, &
source=emptyRealArray)
allocate(prm%h_tw_tw(0,0))
end if twinActive
@ -230,7 +229,7 @@ module function plastic_phenopowerlaw_init() result(myPlasticity)
!--------------------------------------------------------------------------------------------------
! slip-twin related parameters
slipAndTwinActive: if (prm%sum_N_sl > 0 .and. prm%sum_N_tw > 0) then
prm%h_0_tw_sl = math_expand(pl%get_as1dReal('h_0_tw-sl',requiredSize=N_fam,defaultVal=0.0_pReal*ones), N_tw)
prm%h_0_tw_sl = math_expand(pl%get_as1dReal('h_0_tw-sl',requiredSize=size(N_tw)), N_tw)
prm%h_sl_tw = crystal_interaction_SlipByTwin(N_sl,N_tw,pl%get_as1dReal('h_sl-tw'),phase_lattice(ph))
prm%h_tw_sl = crystal_interaction_TwinBySlip(N_tw,N_sl,pl%get_as1dReal('h_tw-sl'),phase_lattice(ph))
else slipAndTwinActive