Merge branch '298-forest-projections' into 'development'

Resolve "Forest projections"

Closes #298

See merge request damask/DAMASK!768
This commit is contained in:
Franz Roters 2023-06-27 15:33:57 +00:00
commit 5f82e975c8
5 changed files with 152 additions and 145 deletions

@ -1 +1 @@
Subproject commit 3ba790ed2d3d8d8cf66fa0ad7be7abe8f77c0d54
Subproject commit 6b0e09f1f31a2892c8e48c8ca688e8dfc05cd70a

View File

@ -12,6 +12,7 @@ output: [rho_dip, rho_mob]
N_sl: [12, 12]
f_edge: [1.0, 1.0]
b_sl: [2.49e-10, 2.49e-10]
rho_mob_0: [2.81e+12, 2.8e+12]
rho_dip_0: [1.0, 1.0] # not given

View File

@ -15,6 +15,8 @@ output: [rho_mob, rho_dip, gamma_sl, Lambda_sl, tau_pass, f_tw, Lambda_tw, f_tr]
# Glide
N_sl: [12]
f_edge: [1.0]
b_sl: [2.56e-10] # a/sqrt(2)
Q_sl: [3.5e-19]
p_sl: [0.325]

View File

@ -83,13 +83,14 @@ module function plastic_dislotungsten_init() result(myPlasticity)
logical, dimension(:), allocatable :: myPlasticity
integer :: &
ph, &
ph, i, &
Nmembers, &
sizeState, sizeDotState, &
startIndex, endIndex
integer, dimension(:), allocatable :: &
N_sl
real(pREAL),dimension(:), allocatable :: &
f_edge, & !< edge character fraction of total dislocation density
rho_mob_0, & !< initial dislocation density
rho_dip_0, & !< initial dipole density
a !< non-Schmid coefficients
@ -123,14 +124,15 @@ module function plastic_dislotungsten_init() result(myPlasticity)
do ph = 1, phases%length
if (.not. myPlasticity(ph)) cycle
associate(prm => param(ph), stt => state(ph), dst => dependentState(ph), &
associate(prm => param(ph), &
stt => state(ph), dst => dependentState(ph), &
idx_dot => indexDotState(ph))
phase => phases%get_dict(ph)
mech => phase%get_dict('mechanical')
pl => mech%get_dict('plastic')
print'(/,1x,a,i0,a)', 'phase ',ph,': '//phases%key(ph)
print'(/,1x,a,1x,i0,a)', 'phase',ph,': '//phases%key(ph)
refs = config_listReferences(pl,indent=3)
if (len(refs) > 0) print'(/,1x,a)', refs
@ -159,48 +161,36 @@ module function plastic_dislotungsten_init() result(myPlasticity)
prm%P_nS_neg = prm%P_sl
end if
prm%h_sl_sl = lattice_interaction_SlipBySlip(N_sl,pl%get_as1dReal('h_sl-sl'), &
phase_lattice(ph))
prm%forestProjection = lattice_forestProjection_edge(N_sl,phase_lattice(ph),&
phase_cOverA(ph))
prm%forestProjection = transpose(prm%forestProjection)
rho_mob_0 = pl%get_as1dReal('rho_mob_0', requiredSize=size(N_sl))
rho_dip_0 = pl%get_as1dReal('rho_dip_0', requiredSize=size(N_sl))
prm%b_sl = pl%get_as1dReal('b_sl', requiredSize=size(N_sl))
prm%Q_s = pl%get_as1dReal('Q_s', requiredSize=size(N_sl))
prm%i_sl = pl%get_as1dReal('i_sl', requiredSize=size(N_sl))
prm%tau_Peierls = pl%get_as1dReal('tau_Peierls', requiredSize=size(N_sl))
prm%p = pl%get_as1dReal('p_sl', requiredSize=size(N_sl))
prm%q = pl%get_as1dReal('q_sl', requiredSize=size(N_sl))
prm%h = pl%get_as1dReal('h', requiredSize=size(N_sl))
prm%w = pl%get_as1dReal('w', requiredSize=size(N_sl))
prm%omega = pl%get_as1dReal('omega', requiredSize=size(N_sl))
prm%B = pl%get_as1dReal('B', requiredSize=size(N_sl))
prm%dipoleformation = .not. pl%get_asBool('no_dipole_formation', defaultVal=.false.)
prm%D = pl%get_asReal('D')
prm%D_0 = pl%get_asReal('D_0')
prm%Q_cl = pl%get_asReal('Q_cl')
prm%f_at = pl%get_asReal('f_at') * prm%b_sl**3
prm%dipoleformation = .not. pl%get_asBool('no_dipole_formation', defaultVal = .false.)
f_edge = math_expand(pl%get_as1dReal('f_edge', requiredSize=size(N_sl), &
defaultVal=[(0.5_pREAL, i=1,size(N_sl))]),N_sl)
rho_mob_0 = math_expand(pl%get_as1dReal('rho_mob_0', requiredSize=size(N_sl)),N_sl)
rho_dip_0 = math_expand(pl%get_as1dReal('rho_dip_0', requiredSize=size(N_sl)),N_sl)
prm%b_sl = math_expand(pl%get_as1dReal('b_sl', requiredSize=size(N_sl)),N_sl)
prm%Q_s = math_expand(pl%get_as1dReal('Q_s', requiredSize=size(N_sl)),N_sl)
prm%i_sl = math_expand(pl%get_as1dReal('i_sl', requiredSize=size(N_sl)),N_sl)
prm%tau_Peierls = math_expand(pl%get_as1dReal('tau_Peierls', requiredSize=size(N_sl)),N_sl)
prm%p = math_expand(pl%get_as1dReal('p_sl', requiredSize=size(N_sl)),N_sl)
prm%q = math_expand(pl%get_as1dReal('q_sl', requiredSize=size(N_sl)),N_sl)
prm%h = math_expand(pl%get_as1dReal('h', requiredSize=size(N_sl)),N_sl)
prm%w = math_expand(pl%get_as1dReal('w', requiredSize=size(N_sl)),N_sl)
prm%omega = math_expand(pl%get_as1dReal('omega', requiredSize=size(N_sl)),N_sl)
prm%B = math_expand(pl%get_as1dReal('B', requiredSize=size(N_sl)),N_sl)
prm%d_caron = prm%b_sl * pl%get_asReal('D_a')
prm%f_at = prm%b_sl**3*pl%get_asReal('f_at')
! expand: family => system
rho_mob_0 = math_expand(rho_mob_0, N_sl)
rho_dip_0 = math_expand(rho_dip_0, N_sl)
prm%q = math_expand(prm%q, N_sl)
prm%p = math_expand(prm%p, N_sl)
prm%Q_s = math_expand(prm%Q_s, N_sl)
prm%b_sl = math_expand(prm%b_sl, N_sl)
prm%h = math_expand(prm%h, N_sl)
prm%w = math_expand(prm%w, N_sl)
prm%omega = math_expand(prm%omega, N_sl)
prm%tau_Peierls = math_expand(prm%tau_Peierls, N_sl)
prm%B = math_expand(prm%B, N_sl)
prm%i_sl = math_expand(prm%i_sl, N_sl)
prm%f_at = math_expand(prm%f_at, N_sl)
prm%d_caron = pl%get_asReal('D_a') * prm%b_sl
prm%h_sl_sl = lattice_interaction_SlipBySlip(N_sl,pl%get_as1dReal('h_sl-sl'), &
phase_lattice(ph))
prm%forestProjection = spread( f_edge,1,prm%sum_N_sl) &
* lattice_forestProjection_edge (N_sl,phase_lattice(ph),phase_cOverA(ph)) &
+ spread(1.0_pREAL-f_edge,1,prm%sum_N_sl) &
* lattice_forestProjection_screw(N_sl,phase_lattice(ph),phase_cOverA(ph))
! sanity checks
if ( prm%D_0 < 0.0_pREAL) extmsg = trim(extmsg)//' D_0'
@ -215,9 +205,20 @@ module function plastic_dislotungsten_init() result(myPlasticity)
if (any(prm%f_at <= 0.0_pREAL)) extmsg = trim(extmsg)//' f_at or b_sl'
else slipActive
rho_mob_0 = emptyRealArray; rho_dip_0 = emptyRealArray
allocate(prm%b_sl,prm%d_caron,prm%i_sl,prm%f_at,prm%tau_Peierls, &
prm%Q_s,prm%p,prm%q,prm%B,prm%h,prm%w,prm%omega, &
rho_mob_0 = emptyRealArray
rho_dip_0 = emptyRealArray
allocate(prm%b_sl, &
prm%d_caron, &
prm%i_sl, &
prm%f_at, &
prm%tau_Peierls, &
prm%Q_s, &
prm%p, &
prm%q, &
prm%B, &
prm%h, &
prm%w, &
prm%omega, &
source = emptyRealArray)
allocate(prm%forestProjection(0,0))
allocate(prm%h_sl_sl (0,0))

View File

@ -34,15 +34,15 @@ submodule(phase:plastic) dislotwin
Gamma_sf, & !< stacking fault energy
Delta_G !< free energy difference between austensite and martensite
real(pREAL), allocatable, dimension(:) :: &
b_sl, & !< absolute length of Burgers vector [m] for each slip system
b_tw, & !< absolute length of Burgers vector [m] for each twin system
b_tr, & !< absolute length of Burgers vector [m] for each transformation system
Q_sl,& !< activation energy for glide [J] for each slip system
v_0, & !< dislocation velocity prefactor [m/s] for each slip system
dot_N_0_tw, & !< twin nucleation rate [1/m³s] for each twin system
t_tw, & !< twin thickness [m] for each twin system
b_sl, & !< magnitude of Burgers vector (m) for each slip system
b_tw, & !< magnitude of Burgers vector (m) for each twin system
b_tr, & !< magnitude of Burgers vector (m) for each transformation system
Q_sl,& !< activation energy for glide (J) for each slip system
v_0, & !< dislocation velocity prefactor (m/s) for each slip system
dot_N_0_tw, & !< twin nucleation rate (1/m³s) for each twin system
t_tw, & !< twin thickness (m) for each twin system
i_sl, & !< Adj. parameter for distance between 2 forest dislocations for each slip system
t_tr, & !< martensite lamellar thickness [m] for each trans system
t_tr, & !< martensite lamellar thickness (m) for each trans system
p, & !< p-exponent in glide velocity
q, & !< q-exponent in glide velocity
r, & !< exponent in twin nucleation rate
@ -71,7 +71,7 @@ submodule(phase:plastic) dislotwin
N_tw, &
N_tr
integer, allocatable, dimension(:,:) :: &
fcc_twinNucleationSlipPair ! ToDo: Better name? Is also use for trans
fcc_twinNucleationSlipPair ! ToDo: Better name? Is also used for trans
character(len=:), allocatable :: &
lattice_tr, &
isotropic_bound
@ -138,6 +138,7 @@ module function plastic_dislotwin_init() result(myPlasticity)
N_sl
real(pREAL) :: a_cF
real(pREAL), allocatable, dimension(:) :: &
f_edge, & !< edge character fraction of total dislocation density
rho_mob_0, & !< initial unipolar dislocation density per slip system
rho_dip_0 !< initial dipole dislocation density per slip system
character(len=:), allocatable :: &
@ -176,14 +177,15 @@ module function plastic_dislotwin_init() result(myPlasticity)
do ph = 1, phases%length
if (.not. myPlasticity(ph)) cycle
associate(prm => param(ph), stt => state(ph), dst => dependentState(ph), &
associate(prm => param(ph), &
stt => state(ph), dst => dependentState(ph), &
idx_dot => indexDotState(ph))
phase => phases%get_dict(ph)
mech => phase%get_dict('mechanical')
pl => mech%get_dict('plastic')
print'(/,1x,a,i0,a)', 'phase ',ph,': '//phases%key(ph)
print'(/,1x,a,1x,i0,a)', 'phase',ph,': '//phases%key(ph)
refs = config_listReferences(pl,indent=3)
if (len(refs) > 0) print'(/,1x,a)', refs
@ -202,49 +204,43 @@ module function plastic_dislotwin_init() result(myPlasticity)
slipActive: if (prm%sum_N_sl > 0) then
prm%systems_sl = lattice_labels_slip(N_sl,phase_lattice(ph))
prm%P_sl = lattice_SchmidMatrix_slip(N_sl,phase_lattice(ph),phase_cOverA(ph))
prm%h_sl_sl = lattice_interaction_SlipBySlip(N_sl,pl%get_as1dReal('h_sl-sl'),phase_lattice(ph))
prm%forestProjection = lattice_forestProjection_edge(N_sl,phase_lattice(ph),phase_cOverA(ph))
prm%forestProjection = transpose(prm%forestProjection)
prm%n0_sl = lattice_slip_normal(N_sl,phase_lattice(ph),phase_cOverA(ph))
prm%fccTwinTransNucleation = phase_lattice(ph) == 'cF' .and. (N_sl(1) == 12)
if (prm%fccTwinTransNucleation) prm%fcc_twinNucleationSlipPair = lattice_CF_TWINNUCLEATIONSLIPPAIR
rho_mob_0 = pl%get_as1dReal('rho_mob_0', requiredSize=size(N_sl))
rho_dip_0 = pl%get_as1dReal('rho_dip_0', requiredSize=size(N_sl))
prm%v_0 = pl%get_as1dReal('v_0', requiredSize=size(N_sl))
prm%b_sl = pl%get_as1dReal('b_sl', requiredSize=size(N_sl))
prm%Q_sl = pl%get_as1dReal('Q_sl', requiredSize=size(N_sl))
prm%i_sl = pl%get_as1dReal('i_sl', requiredSize=size(N_sl))
prm%p = pl%get_as1dReal('p_sl', requiredSize=size(N_sl))
prm%q = pl%get_as1dReal('q_sl', requiredSize=size(N_sl))
prm%tau_0 = pl%get_as1dReal('tau_0', requiredSize=size(N_sl))
prm%B = pl%get_as1dReal('B', requiredSize=size(N_sl), &
defaultVal=[(0.0_pREAL, i=1,size(N_sl))])
prm%Q_cl = pl%get_asReal('Q_cl')
prm%extendedDislocations = pl%get_asBool('extend_dislocations',defaultVal=.false.)
prm%omitDipoles = pl%get_asBool('omit_dipoles', defaultVal=.false.)
prm%Q_cl = pl%get_asReal('Q_cl')
f_edge = math_expand(pl%get_as1dReal('f_edge', requiredSize=size(N_sl), &
defaultVal=[(0.5_pREAL,i=1,size(N_sl))]),N_sl)
rho_mob_0 = math_expand(pl%get_as1dReal('rho_mob_0', requiredSize=size(N_sl)),N_sl)
rho_dip_0 = math_expand(pl%get_as1dReal('rho_dip_0', requiredSize=size(N_sl)),N_sl)
prm%v_0 = math_expand(pl%get_as1dReal('v_0', requiredSize=size(N_sl)),N_sl)
prm%b_sl = math_expand(pl%get_as1dReal('b_sl', requiredSize=size(N_sl)),N_sl)
prm%Q_sl = math_expand(pl%get_as1dReal('Q_sl', requiredSize=size(N_sl)),N_sl)
prm%i_sl = math_expand(pl%get_as1dReal('i_sl', requiredSize=size(N_sl)),N_sl)
prm%p = math_expand(pl%get_as1dReal('p_sl', requiredSize=size(N_sl)),N_sl)
prm%q = math_expand(pl%get_as1dReal('q_sl', requiredSize=size(N_sl)),N_sl)
prm%tau_0 = math_expand(pl%get_as1dReal('tau_0', requiredSize=size(N_sl)),N_sl)
prm%B = math_expand(pl%get_as1dReal('B', requiredSize=size(N_sl), &
defaultVal=[(0.0_pREAL,i=1,size(N_sl))]),N_sl)
prm%d_caron = prm%b_sl * pl%get_asReal('D_a')
prm%h_sl_sl = lattice_interaction_SlipBySlip(N_sl,pl%get_as1dReal('h_sl-sl'),phase_lattice(ph))
prm%forestProjection = spread( f_edge,1,prm%sum_N_sl) &
* lattice_forestProjection_edge (N_sl,phase_lattice(ph),phase_cOverA(ph)) &
+ spread(1.0_pREAL-f_edge,1,prm%sum_N_sl) &
* lattice_forestProjection_screw(N_sl,phase_lattice(ph),phase_cOverA(ph))
prm%fccTwinTransNucleation = phase_lattice(ph) == 'cF' .and. N_sl(1) == 12
if (prm%fccTwinTransNucleation) prm%fcc_twinNucleationSlipPair = lattice_CF_TWINNUCLEATIONSLIPPAIR
! multiplication factor according to crystal structure (nearest neighbors bcc vs fcc/hex)
! details: Argon & Moffat, Acta Metallurgica, Vol. 29, pg 293 to 299, 1981
prm%omega = pl%get_asReal('omega', defaultVal=1000.0_pREAL) &
* merge(12.0_pREAL,8.0_pREAL,any(phase_lattice(ph) == ['cF','hP']))
! expand: family => system
rho_mob_0 = math_expand(rho_mob_0, N_sl)
rho_dip_0 = math_expand(rho_dip_0, N_sl)
prm%v_0 = math_expand(prm%v_0, N_sl)
prm%b_sl = math_expand(prm%b_sl, N_sl)
prm%Q_sl = math_expand(prm%Q_sl, N_sl)
prm%i_sl = math_expand(prm%i_sl, N_sl)
prm%p = math_expand(prm%p, N_sl)
prm%q = math_expand(prm%q, N_sl)
prm%tau_0 = math_expand(prm%tau_0, N_sl)
prm%B = math_expand(prm%B, N_sl)
prm%d_caron = pl%get_asReal('D_a') * prm%b_sl
! sanity checks
if ( prm%Q_cl <= 0.0_pREAL) extmsg = trim(extmsg)//' Q_cl'
if (any(rho_mob_0 < 0.0_pREAL)) extmsg = trim(extmsg)//' rho_mob_0'
@ -258,9 +254,19 @@ module function plastic_dislotwin_init() result(myPlasticity)
if (any(prm%p<=0.0_pREAL .or. prm%p>1.0_pREAL)) extmsg = trim(extmsg)//' p_sl'
if (any(prm%q< 1.0_pREAL .or. prm%q>2.0_pREAL)) extmsg = trim(extmsg)//' q_sl'
else slipActive
rho_mob_0 = emptyRealArray; rho_dip_0 = emptyRealArray
allocate(prm%b_sl,prm%Q_sl,prm%v_0,prm%i_sl,prm%p,prm%q,prm%B,source=emptyRealArray)
allocate(prm%forestProjection(0,0),prm%h_sl_sl(0,0))
rho_mob_0 = emptyRealArray
rho_dip_0 = emptyRealArray
allocate(prm%v_0, &
prm%b_sl, &
prm%Q_sl, &
prm%i_sl, &
prm%p, &
prm%q, &
prm%tau_0, &
prm%B, &
source=emptyRealArray)
allocate(prm%forestProjection(0,0), &
prm%h_sl_sl(0,0)) ! PE: What about P_sl and systems_sl?
end if slipActive
!--------------------------------------------------------------------------------------------------
@ -270,22 +276,17 @@ module function plastic_dislotwin_init() result(myPlasticity)
twinActive: if (prm%sum_N_tw > 0) then
prm%systems_tw = lattice_labels_twin(prm%N_tw,phase_lattice(ph))
prm%P_tw = lattice_SchmidMatrix_twin(prm%N_tw,phase_lattice(ph),phase_cOverA(ph))
prm%h_tw_tw = lattice_interaction_TwinByTwin(prm%N_tw,pl%get_as1dReal('h_tw-tw'), &
phase_lattice(ph))
prm%b_tw = pl%get_as1dReal('b_tw', requiredSize=size(prm%N_tw))
prm%t_tw = pl%get_as1dReal('t_tw', requiredSize=size(prm%N_tw))
prm%r = pl%get_as1dReal('p_tw', requiredSize=size(prm%N_tw))
prm%gamma_char_tw = lattice_characteristicShear_Twin(prm%N_tw,phase_lattice(ph),phase_cOverA(ph))
prm%L_tw = pl%get_asReal('L_tw')
prm%i_tw = pl%get_asReal('i_tw')
prm%gamma_char_tw = lattice_characteristicShear_Twin(prm%N_tw,phase_lattice(ph),phase_cOverA(ph))
prm%b_tw = math_expand(pl%get_as1dReal('b_tw', requiredSize=size(prm%N_tw)),prm%N_tw)
prm%t_tw = math_expand(pl%get_as1dReal('t_tw', requiredSize=size(prm%N_tw)),prm%N_tw)
prm%r = math_expand(pl%get_as1dReal('p_tw', requiredSize=size(prm%N_tw)),prm%N_tw)
! expand: family => system
prm%b_tw = math_expand(prm%b_tw,prm%N_tw)
prm%t_tw = math_expand(prm%t_tw,prm%N_tw)
prm%r = math_expand(prm%r,prm%N_tw)
prm%h_tw_tw = lattice_interaction_TwinByTwin(prm%N_tw,pl%get_as1dReal('h_tw-tw'), &
phase_lattice(ph))
! sanity checks
if (.not. prm%fccTwinTransNucleation) extmsg = trim(extmsg)//' TWIP for non-fcc'
@ -295,7 +296,11 @@ module function plastic_dislotwin_init() result(myPlasticity)
if (any(prm%t_tw < 0.0_pREAL)) extmsg = trim(extmsg)//' t_tw'
if (any(prm%r < 0.0_pREAL)) extmsg = trim(extmsg)//' p_tw'
else twinActive
allocate(prm%gamma_char_tw,prm%b_tw,prm%t_tw,prm%r,source=emptyRealArray)
allocate(prm%gamma_char_tw, &
prm%b_tw, &
prm%t_tw, &
prm%r, &
source=emptyRealArray)
allocate(prm%h_tw_tw(0,0))
end if twinActive
@ -304,26 +309,24 @@ module function plastic_dislotwin_init() result(myPlasticity)
prm%N_tr = pl%get_as1dInt('N_tr', defaultVal=emptyIntArray)
prm%sum_N_tr = sum(abs(prm%N_tr))
transActive: if (prm%sum_N_tr > 0) then
prm%b_tr = pl%get_as1dReal('b_tr')
prm%b_tr = math_expand(prm%b_tr,prm%N_tr)
prm%P_tr = lattice_SchmidMatrix_trans(prm%N_tr,'hP',prm%cOverA_hP)
prm%i_tr = pl%get_asReal('i_tr')
prm%Delta_G = polynomial(pl,'Delta_G','T')
prm%i_tr = pl%get_asReal('i_tr')
prm%L_tr = pl%get_asReal('L_tr')
prm%cOverA_hP = pl%get_asReal('c/a_hP')
prm%V_mol = pl%get_asReal('V_mol')
prm%b_tr = math_expand(pl%get_as1dReal('b_tr'),prm%N_tr)
prm%t_tr = math_expand(pl%get_as1dReal('t_tr'),prm%N_tr)
prm%s = math_expand(pl%get_as1dReal('p_tr'),prm%N_tr)
a_cF = prm%b_tr(1)*sqrt(6.0_pREAL) ! b_tr is Shockley partial
prm%h = 5.0_pREAL * a_cF/sqrt(3.0_pREAL)
prm%cOverA_hP = pl%get_asReal('c/a_hP')
prm%rho = 4.0_pREAL/(sqrt(3.0_pREAL)*a_cF**2)/N_A
prm%V_mol = pl%get_asReal('V_mol')
prm%h_tr_tr = lattice_interaction_TransByTrans(prm%N_tr,pl%get_as1dReal('h_tr-tr'),&
phase_lattice(ph))
prm%P_tr = lattice_SchmidMatrix_trans(prm%N_tr,'hP',prm%cOverA_hP)
prm%t_tr = pl%get_as1dReal('t_tr')
prm%t_tr = math_expand(prm%t_tr,prm%N_tr)
prm%s = pl%get_as1dReal('p_tr')
prm%s = math_expand(prm%s,prm%N_tr)
! sanity checks
if (.not. prm%fccTwinTransNucleation) extmsg = trim(extmsg)//' TRIP for non-fcc'