223031012
/
DAMASK_EICMD
mirror of https://github.com/achalhp/DAMASK_EICMD.git
1
0
Fork 0
Code Issues Projects Releases Wiki Activity

more systematic naming

This commit is contained in:
Martin Diehl 2023-12-29 18:02:52 +01:00
parent b07c2a3191
commit 1f1baa8b48
No known key found for this signature in database
GPG Key ID: 1FD50837275A0A9B
1 changed files with 214 additions and 224 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

438
src/phase_mechanical_plastic_nonlocal.f90
View File

@ -6,18 +6,18 @@
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
submodule(phase:plastic) nonlocal submodule(phase:plastic) nonlocal
use geometry_plastic_nonlocal, only: & use geometry_plastic_nonlocal, only: &
nIPneighbors => geometry_plastic_nonlocal_nIPneighbors, & nCellNeighbors => geometry_plastic_nonlocal_nIPneighbors, &
IPneighborhood => geometry_plastic_nonlocal_IPneighborhood, & IPneighborhood => geometry_plastic_nonlocal_IPneighborhood, &
IPvolume => geometry_plastic_nonlocal_IPvolume0, & IPvolume0 => geometry_plastic_nonlocal_IPvolume0, &
IParea => geometry_plastic_nonlocal_IParea0, & IParea0 => geometry_plastic_nonlocal_IParea0, &
IPareaNormal => geometry_plastic_nonlocal_IPareaNormal0, & IPareaNormal0 => geometry_plastic_nonlocal_IPareaNormal0, &
geometry_plastic_nonlocal_disable geometry_plastic_nonlocal_disable
type :: tGeometry type :: tGeometry
real(pREAL), dimension(:), allocatable :: V_0 real(pREAL), dimension(:), allocatable :: v_0
real(pREAL), dimension(:,:), allocatable :: a_0, x_0
real(pREAL), dimension(:,:,:), allocatable :: n_0
integer, dimension(:,:,:), allocatable :: IPneighborhood integer, dimension(:,:,:), allocatable :: IPneighborhood
real(pREAL), dimension(:,:), allocatable :: IParea, IPcoordinates
real(pREAL), dimension(:,:,:), allocatable :: IPareaNormal
end type tGeometry end type tGeometry
type(tGeometry), dimension(:), allocatable :: geom type(tGeometry), dimension(:), allocatable :: geom
@ -133,18 +133,18 @@ submodule(phase:plastic) nonlocal
type :: tNonlocalState type :: tNonlocalState
real(pREAL), pointer, dimension(:,:) :: & real(pREAL), pointer, dimension(:,:) :: &
rho, & ! < all dislocations rho, & ! < all dislocations
rhoSgl, & rho_sgl, &
rhoSglMobile, & ! iRhoU rho_sgl_mob, & ! iRhoU
rho_sgl_mob_edg_pos, & rho_sgl_mob_edg_pos, &
rho_sgl_mob_edg_neg, & rho_sgl_mob_edg_neg, &
rho_sgl_mob_scr_pos, & rho_sgl_mob_scr_pos, &
rho_sgl_mob_scr_neg, & rho_sgl_mob_scr_neg, &
rhoSglImmobile, & rho_sgl_imm, &
rho_sgl_imm_edg_pos, & rho_sgl_imm_edg_pos, &
rho_sgl_imm_edg_neg, & rho_sgl_imm_edg_neg, &
rho_sgl_imm_scr_pos, & rho_sgl_imm_scr_pos, &
rho_sgl_imm_scr_neg, & rho_sgl_imm_scr_neg, &
rhoDip, & rho_dip, &
rho_dip_edg, & rho_dip_edg, &
rho_dip_scr, & rho_dip_scr, &
gamma, & gamma, &
@ -383,12 +383,12 @@ module function plastic_nonlocal_init() result(myPlasticity)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! allocate state arrays ! allocate state arrays
Nmembers = count(material_ID_phase == ph) Nmembers = count(material_ID_phase == ph)
sizeDotState = size([ 'rhoSglEdgePosMobile ','rhoSglEdgeNegMobile ', & sizeDotState = size(['rho_sgl_mob_edg_pos', 'rho_sgl_mob_edg_neg', &
'rhoSglScrewPosMobile ','rhoSglScrewNegMobile ', & 'rho_sgl_mob_scr_pos', 'rho_sgl_mob_scr_neg', &
'rhoSglEdgePosImmobile ','rhoSglEdgeNegImmobile ', & 'rho_sgl_imm_edg_pos', 'rho_sgl_imm_edg_neg', &
'rhoSglScrewPosImmobile','rhoSglScrewNegImmobile', & 'rho_sgl_imm_scr_pos', 'rho_sgl_imm_scr_neg', &
'rhoDipEdge ','rhoDipScrew ', & 'rho_dip_edg ', 'rho_dip_scr ', &
'gamma ' ]) * prm%sum_N_sl !< "basic" microstructural state variables that are independent from other state variables 'gamma ' ]) * prm%sum_N_sl !< "basic" microstructural state variables that are independent from other state variables
sizeState = sizeDotState & sizeState = sizeDotState &
+ size([ 'velocityEdgePos ','velocityEdgeNeg ', & + size([ 'velocityEdgePos ','velocityEdgeNeg ', &
'velocityScrewPos ','velocityScrewNeg ']) * prm%sum_N_sl !< other dependent state variables that are not updated by microstructure 'velocityScrewPos ','velocityScrewNeg ']) * prm%sum_N_sl !< other dependent state variables that are not updated by microstructure
@ -396,11 +396,11 @@ module function plastic_nonlocal_init() result(myPlasticity)
call phase_allocateState(plasticState(ph),Nmembers,sizeState,sizeDotState,sizeDeltaState,0) ! ToDo: state structure does not follow convention call phase_allocateState(plasticState(ph),Nmembers,sizeState,sizeDotState,sizeDeltaState,0) ! ToDo: state structure does not follow convention
allocate(geom(ph)%V_0(Nmembers)) allocate(geom(ph)%v_0(Nmembers))
allocate(geom(ph)%IPneighborhood(3,nIPneighbors,Nmembers)) allocate(geom(ph)%a_0(nCellNeighbors,Nmembers))
allocate(geom(ph)%IPareaNormal(3,nIPneighbors,Nmembers)) allocate(geom(ph)%x_0(3,Nmembers))
allocate(geom(ph)%IParea(nIPneighbors,Nmembers)) allocate(geom(ph)%n_0(3,nCellNeighbors,Nmembers))
allocate(geom(ph)%IPcoordinates(3,Nmembers)) allocate(geom(ph)%IPneighborhood(3,nCellNeighbors,Nmembers))
call storeGeometry(ph) call storeGeometry(ph)
if (plasticState(ph)%nonlocal .and. .not. allocated(IPneighborhood)) & if (plasticState(ph)%nonlocal .and. .not. allocated(IPneighborhood)) &
@ -412,13 +412,13 @@ module function plastic_nonlocal_init() result(myPlasticity)
del%rho => plasticState(ph)%deltaState (0*prm%sum_N_sl+1:10*prm%sum_N_sl,:) del%rho => plasticState(ph)%deltaState (0*prm%sum_N_sl+1:10*prm%sum_N_sl,:)
plasticState(ph)%atol(1:10*prm%sum_N_sl) = prm%atol_rho plasticState(ph)%atol(1:10*prm%sum_N_sl) = prm%atol_rho
stt%rhoSgl => plasticState(ph)%state (0*prm%sum_N_sl+1: 8*prm%sum_N_sl,:) stt%rho_sgl => plasticState(ph)%state (0*prm%sum_N_sl+1: 8*prm%sum_N_sl,:)
dot%rhoSgl => plasticState(ph)%dotState (0*prm%sum_N_sl+1: 8*prm%sum_N_sl,:) dot%rho_sgl => plasticState(ph)%dotState (0*prm%sum_N_sl+1: 8*prm%sum_N_sl,:)
del%rhoSgl => plasticState(ph)%deltaState (0*prm%sum_N_sl+1: 8*prm%sum_N_sl,:) del%rho_sgl => plasticState(ph)%deltaState (0*prm%sum_N_sl+1: 8*prm%sum_N_sl,:)
stt%rhoSglMobile => plasticState(ph)%state (0*prm%sum_N_sl+1: 4*prm%sum_N_sl,:) stt%rho_sgl_mob => plasticState(ph)%state (0*prm%sum_N_sl+1: 4*prm%sum_N_sl,:)
dot%rhoSglMobile => plasticState(ph)%dotState (0*prm%sum_N_sl+1: 4*prm%sum_N_sl,:) dot%rho_sgl_mob => plasticState(ph)%dotState (0*prm%sum_N_sl+1: 4*prm%sum_N_sl,:)
del%rhoSglMobile => plasticState(ph)%deltaState (0*prm%sum_N_sl+1: 4*prm%sum_N_sl,:) del%rho_sgl_mob => plasticState(ph)%deltaState (0*prm%sum_N_sl+1: 4*prm%sum_N_sl,:)
stt%rho_sgl_mob_edg_pos => plasticState(ph)%state (0*prm%sum_N_sl+1: 1*prm%sum_N_sl,:) stt%rho_sgl_mob_edg_pos => plasticState(ph)%state (0*prm%sum_N_sl+1: 1*prm%sum_N_sl,:)
dot%rho_sgl_mob_edg_pos => plasticState(ph)%dotState (0*prm%sum_N_sl+1: 1*prm%sum_N_sl,:) dot%rho_sgl_mob_edg_pos => plasticState(ph)%dotState (0*prm%sum_N_sl+1: 1*prm%sum_N_sl,:)
@ -436,9 +436,9 @@ module function plastic_nonlocal_init() result(myPlasticity)
dot%rho_sgl_mob_scr_neg => plasticState(ph)%dotState (3*prm%sum_N_sl+1: 4*prm%sum_N_sl,:) dot%rho_sgl_mob_scr_neg => plasticState(ph)%dotState (3*prm%sum_N_sl+1: 4*prm%sum_N_sl,:)
del%rho_sgl_mob_scr_neg => plasticState(ph)%deltaState (3*prm%sum_N_sl+1: 4*prm%sum_N_sl,:) del%rho_sgl_mob_scr_neg => plasticState(ph)%deltaState (3*prm%sum_N_sl+1: 4*prm%sum_N_sl,:)
stt%rhoSglImmobile => plasticState(ph)%state (4*prm%sum_N_sl+1: 8*prm%sum_N_sl,:) stt%rho_sgl_imm => plasticState(ph)%state (4*prm%sum_N_sl+1: 8*prm%sum_N_sl,:)
dot%rhoSglImmobile => plasticState(ph)%dotState (4*prm%sum_N_sl+1: 8*prm%sum_N_sl,:) dot%rho_sgl_imm => plasticState(ph)%dotState (4*prm%sum_N_sl+1: 8*prm%sum_N_sl,:)
del%rhoSglImmobile => plasticState(ph)%deltaState (4*prm%sum_N_sl+1: 8*prm%sum_N_sl,:) del%rho_sgl_imm => plasticState(ph)%deltaState (4*prm%sum_N_sl+1: 8*prm%sum_N_sl,:)
stt%rho_sgl_imm_edg_pos => plasticState(ph)%state (4*prm%sum_N_sl+1: 5*prm%sum_N_sl,:) stt%rho_sgl_imm_edg_pos => plasticState(ph)%state (4*prm%sum_N_sl+1: 5*prm%sum_N_sl,:)
dot%rho_sgl_imm_edg_pos => plasticState(ph)%dotState (4*prm%sum_N_sl+1: 5*prm%sum_N_sl,:) dot%rho_sgl_imm_edg_pos => plasticState(ph)%dotState (4*prm%sum_N_sl+1: 5*prm%sum_N_sl,:)
@ -456,9 +456,9 @@ module function plastic_nonlocal_init() result(myPlasticity)
dot%rho_sgl_imm_scr_neg => plasticState(ph)%dotState (7*prm%sum_N_sl+1: 8*prm%sum_N_sl,:) dot%rho_sgl_imm_scr_neg => plasticState(ph)%dotState (7*prm%sum_N_sl+1: 8*prm%sum_N_sl,:)
del%rho_sgl_imm_scr_neg => plasticState(ph)%deltaState (7*prm%sum_N_sl+1: 8*prm%sum_N_sl,:) del%rho_sgl_imm_scr_neg => plasticState(ph)%deltaState (7*prm%sum_N_sl+1: 8*prm%sum_N_sl,:)
stt%rhoDip => plasticState(ph)%state (8*prm%sum_N_sl+1:10*prm%sum_N_sl,:) stt%rho_dip => plasticState(ph)%state (8*prm%sum_N_sl+1:10*prm%sum_N_sl,:)
dot%rhoDip => plasticState(ph)%dotState (8*prm%sum_N_sl+1:10*prm%sum_N_sl,:) dot%rho_dip => plasticState(ph)%dotState (8*prm%sum_N_sl+1:10*prm%sum_N_sl,:)
del%rhoDip => plasticState(ph)%deltaState (8*prm%sum_N_sl+1:10*prm%sum_N_sl,:) del%rho_dip => plasticState(ph)%deltaState (8*prm%sum_N_sl+1:10*prm%sum_N_sl,:)
stt%rho_dip_edg => plasticState(ph)%state (8*prm%sum_N_sl+1: 9*prm%sum_N_sl,:) stt%rho_dip_edg => plasticState(ph)%state (8*prm%sum_N_sl+1: 9*prm%sum_N_sl,:)
dot%rho_dip_edg => plasticState(ph)%dotState (8*prm%sum_N_sl+1: 9*prm%sum_N_sl,:) dot%rho_dip_edg => plasticState(ph)%dotState (8*prm%sum_N_sl+1: 9*prm%sum_N_sl,:)
@ -486,7 +486,7 @@ module function plastic_nonlocal_init() result(myPlasticity)
allocate(dst%tau_back(prm%sum_N_sl,Nmembers),source=0.0_pREAL) allocate(dst%tau_back(prm%sum_N_sl,Nmembers),source=0.0_pREAL)
allocate(dst%rho_forest(prm%sum_N_sl,Nmembers),source=0.0_pREAL) allocate(dst%rho_forest(prm%sum_N_sl,Nmembers),source=0.0_pREAL)
allocate(dst%max_dipole_height(2*prm%sum_N_sl,Nmembers),source=0.0_pREAL) ! edge and screw allocate(dst%max_dipole_height(2*prm%sum_N_sl,Nmembers),source=0.0_pREAL) ! edge and screw
allocate(dst%compatibility(2,maxval(param%sum_N_sl),maxval(param%sum_N_sl),nIPneighbors,Nmembers),source=0.0_pREAL) allocate(dst%compatibility(2,maxval(param%sum_N_sl),maxval(param%sum_N_sl),nCellNeighbors,Nmembers),source=0.0_pREAL)
end associate end associate
if (Nmembers > 0) call stateInit(ini,ph,Nmembers) if (Nmembers > 0) call stateInit(ini,ph,Nmembers)
@ -497,7 +497,7 @@ module function plastic_nonlocal_init() result(myPlasticity)
end do end do
allocate(compatibility(2,maxval(param%sum_N_sl),maxval(param%sum_N_sl),nIPneighbors,& allocate(compatibility(2,maxval(param%sum_N_sl),maxval(param%sum_N_sl),nCellNeighbors,&
discretization_nIPs,discretization_Nelems), source=0.0_pREAL) discretization_nIPs,discretization_Nelems), source=0.0_pREAL)
! BEGIN DEPRECATED---------------------------------------------------------------------------------- ! BEGIN DEPRECATED----------------------------------------------------------------------------------
@ -541,9 +541,9 @@ module subroutine nonlocal_dependentState(ph, en)
en en
integer :: & integer :: &
no, & !< neighbor offset en_nbr, & ! neighbor phase entry
neighbor_el, & ! element number of neighboring material point el_nbr, & ! element number of neighboring material point
neighbor_ip, & ! integration point of neighboring material point ip_nbr, & ! integration point of neighboring material point
c, & ! index of dilsocation character (edge, screw) c, & ! index of dilsocation character (edge, screw)
s, & ! slip system index s, & ! slip system index
dir, & dir, &
@ -558,7 +558,7 @@ module subroutine nonlocal_dependentState(ph, en)
real(pREAL), dimension(2) :: & real(pREAL), dimension(2) :: &
rhoExcessGradient, & rhoExcessGradient, &
rhoExcessGradient_over_rho, & rhoExcessGradient_over_rho, &
rhoTotal rho_sum
real(pREAL), dimension(3) :: & real(pREAL), dimension(3) :: &
rhoExcessDifferences, & rhoExcessDifferences, &
normal_latticeConf normal_latticeConf
@ -567,25 +567,23 @@ module subroutine nonlocal_dependentState(ph, en)
invFp, & !< inverse of plastic deformation gradient invFp, & !< inverse of plastic deformation gradient
connections, & connections, &
invConnections invConnections
real(pREAL), dimension(3,nIPneighbors) :: & real(pREAL), dimension(3,nCellNeighbors) :: &
connection_latticeConf connection_latticeConf
real(pREAL), dimension(2,param(ph)%sum_N_sl) :: &
rhoExcess
real(pREAL), dimension(param(ph)%sum_N_sl) :: & real(pREAL), dimension(param(ph)%sum_N_sl) :: &
rho_edg_delta, & rho_edg_delta, &
rho_scr_delta rho_scr_delta
real(pREAL), dimension(param(ph)%sum_N_sl,10) :: & real(pREAL), dimension(param(ph)%sum_N_sl,10) :: &
rho, & rho, &
rho0, & rho_0, &
rho_neighbor0 rho_0_nbr
real(pREAL), dimension(param(ph)%sum_N_sl,param(ph)%sum_N_sl) :: & real(pREAL), dimension(param(ph)%sum_N_sl,param(ph)%sum_N_sl) :: &
myInteractionMatrix ! corrected slip interaction matrix myInteractionMatrix ! corrected slip interaction matrix
real(pREAL), dimension(param(ph)%sum_N_sl,nIPneighbors) :: & real(pREAL), dimension(param(ph)%sum_N_sl,nCellNeighbors) :: &
rho_edg_delta_neighbor, & rho_edg_delta_nbr, &
rho_scr_delta_neighbor rho_scr_delta_nbr
real(pREAL), dimension(2,maxval(param%sum_N_sl),nIPneighbors) :: & real(pREAL), dimension(2,maxval(param%sum_N_sl),nCellNeighbors) :: &
neighbor_rhoExcess, & ! excess density at neighboring material point rho_delta_nbr, & ! excess density at neighboring material point
neighbor_rhoTotal ! total density at neighboring material point rho_sum_nbr ! total density at neighboring material point
real(pREAL), dimension(3,param(ph)%sum_N_sl,2) :: & real(pREAL), dimension(3,param(ph)%sum_N_sl,2) :: &
m ! direction of dislocation motion m ! direction of dislocation motion
@ -621,60 +619,57 @@ module subroutine nonlocal_dependentState(ph, en)
! ToDo: MD: this is most likely only correct for F_i = I ! ToDo: MD: this is most likely only correct for F_i = I
!################################################################################################# !#################################################################################################
rho0 = getRho0(ph,en) rho_0 = getRho0(ph,en)
if (plasticState(ph)%nonlocal .and. prm%shortRangeStressCorrection) then if (plasticState(ph)%nonlocal .and. prm%shortRangeStressCorrection) then
invFp = math_inv33(phase_mechanical_Fp(ph)%data(1:3,1:3,en)) invFp = math_inv33(phase_mechanical_Fp(ph)%data(1:3,1:3,en))
invFe = math_inv33(phase_mechanical_Fe(ph)%data(1:3,1:3,en)) invFe = math_inv33(phase_mechanical_Fe(ph)%data(1:3,1:3,en))
rho_edg_delta = rho0(:,mob_edg_pos) - rho0(:,mob_edg_neg) rho_edg_delta = rho_0(:,mob_edg_pos) - rho_0(:,mob_edg_neg)
rho_scr_delta = rho0(:,mob_scr_pos) - rho0(:,mob_scr_neg) rho_scr_delta = rho_0(:,mob_scr_pos) - rho_0(:,mob_scr_neg)
rhoExcess(1,:) = rho_edg_delta FVsize = geom(ph)%v_0(en)**(1.0_pREAL/3.0_pREAL)
rhoExcess(2,:) = rho_scr_delta
FVsize = geom(ph)%V_0(en)**(1.0_pREAL/3.0_pREAL)
!* loop through my neighborhood and get the connection vectors (in lattice frame) and the excess densities !* loop through my neighborhood and get the connection vectors (in lattice frame) and the excess densities
nRealNeighbors = 0.0_pREAL nRealNeighbors = 0.0_pREAL
neighbor_rhoTotal = 0.0_pREAL rho_sum_nbr = 0.0_pREAL
do n = 1,nIPneighbors do n = 1,nCellNeighbors
neighbor_el = geom(ph)%IPneighborhood(1,n,en) el_nbr = geom(ph)%IPneighborhood(1,n,en)
neighbor_ip = geom(ph)%IPneighborhood(2,n,en) ip_nbr = geom(ph)%IPneighborhood(2,n,en)
if (neighbor_el > 0 .and. neighbor_ip > 0) then if (el_nbr > 0 .and. ip_nbr > 0) then
if (material_ID_phase(1,(neighbor_el-1)*discretization_nIPs + neighbor_ip) == ph) then if (material_ID_phase(1,(el_nbr-1)*discretization_nIPs + ip_nbr) == ph) then
no = material_entry_phase(1,(neighbor_el-1)*discretization_nIPs + neighbor_ip) en_nbr = material_entry_phase(1,(el_nbr-1)*discretization_nIPs + ip_nbr)
nRealNeighbors = nRealNeighbors + 1.0_pREAL nRealNeighbors = nRealNeighbors + 1.0_pREAL
rho_neighbor0 = getRho0(ph,no) rho_0_nbr = getRho0(ph,en_nbr)
rho_edg_delta_neighbor(:,n) = rho_neighbor0(:,mob_edg_pos) - rho_neighbor0(:,mob_edg_neg) rho_edg_delta_nbr(:,n) = rho_0_nbr(:,mob_edg_pos) - rho_0_nbr(:,mob_edg_neg)
rho_scr_delta_neighbor(:,n) = rho_neighbor0(:,mob_scr_pos) - rho_neighbor0(:,mob_scr_neg) rho_scr_delta_nbr(:,n) = rho_0_nbr(:,mob_scr_pos) - rho_0_nbr(:,mob_scr_neg)
neighbor_rhoTotal(1,:,n) = sum(abs(rho_neighbor0(:,edg)),2) rho_sum_nbr(1,:,n) = sum(abs(rho_0_nbr(:,edg)),2)
neighbor_rhoTotal(2,:,n) = sum(abs(rho_neighbor0(:,scr)),2) rho_sum_nbr(2,:,n) = sum(abs(rho_0_nbr(:,scr)),2)
connection_latticeConf(1:3,n) = matmul(invFe, geom(ph)%IPcoordinates(1:3,no) & connection_latticeConf(1:3,n) = matmul(invFe, geom(ph)%x_0(1:3,en_nbr) &
- geom(ph)%IPcoordinates(1:3,en)) - geom(ph)%x_0(1:3,en))
normal_latticeConf = matmul(transpose(invFp), geom(ph)%IPareaNormal(1:3,n,en)) normal_latticeConf = matmul(transpose(invFp), geom(ph)%n_0(1:3,n,en))
if (math_inner(normal_latticeConf,connection_latticeConf(1:3,n)) < 0.0_pREAL) & ! neighboring connection points in opposite direction to face normal: must be periodic image if (math_inner(normal_latticeConf,connection_latticeConf(1:3,n)) < 0.0_pREAL) & ! neighboring connection points in opposite direction to face normal: must be periodic image
connection_latticeConf(1:3,n) = normal_latticeConf * geom(ph)%V_0(en)/geom(ph)%IParea(n,en) ! instead take the surface normal scaled with the diameter of the cell connection_latticeConf(1:3,n) = normal_latticeConf * geom(ph)%v_0(en)/geom(ph)%a_0(n,en) ! instead take the surface normal scaled with the diameter of the cell
else else
! local neighbor or different lattice structure or different constitution instance -> use central values instead ! local neighbor or different lattice structure or different constitution instance -> use central values instead
connection_latticeConf(1:3,n) = 0.0_pREAL connection_latticeConf(1:3,n) = 0.0_pREAL
rho_edg_delta_neighbor(:,n) = rho_edg_delta rho_edg_delta_nbr(:,n) = rho_edg_delta
rho_scr_delta_neighbor(:,n) = rho_scr_delta rho_scr_delta_nbr(:,n) = rho_scr_delta
end if end if
else else
! free surface -> use central values instead ! free surface -> use central values instead
connection_latticeConf(1:3,n) = 0.0_pREAL connection_latticeConf(1:3,n) = 0.0_pREAL
rho_edg_delta_neighbor(:,n) = rho_edg_delta rho_edg_delta_nbr(:,n) = rho_edg_delta
rho_scr_delta_neighbor(:,n) = rho_scr_delta rho_scr_delta_nbr(:,n) = rho_scr_delta
end if end if
end do end do
neighbor_rhoExcess(1,:,:) = rho_edg_delta_neighbor rho_delta_nbr(1,:,:) = rho_edg_delta_nbr
neighbor_rhoExcess(2,:,:) = rho_scr_delta_neighbor rho_delta_nbr(2,:,:) = rho_scr_delta_nbr
!* loop through the slip systems and calculate the dislocation gradient by !* loop through the slip systems and calculate the dislocation gradient by
!* 1. interpolation of the excess density in the neighorhood !* 1. interpolation of the excess density in the neighorhood
@ -691,8 +686,8 @@ module subroutine nonlocal_dependentState(ph, en)
neighbors(2) = 2 * dir neighbors(2) = 2 * dir
connections(dir,1:3) = connection_latticeConf(1:3,neighbors(1)) & connections(dir,1:3) = connection_latticeConf(1:3,neighbors(1)) &
- connection_latticeConf(1:3,neighbors(2)) - connection_latticeConf(1:3,neighbors(2))
rhoExcessDifferences(dir) = neighbor_rhoExcess(c,s,neighbors(1)) & rhoExcessDifferences(dir) = rho_delta_nbr(c,s,neighbors(1)) &
- neighbor_rhoExcess(c,s,neighbors(2)) - rho_delta_nbr(c,s,neighbors(2))
end do end do
invConnections = math_inv33(connections) invConnections = math_inv33(connections)
if (all(dEq0(invConnections))) call IO_error(-1,ext_msg='back stress calculation: inversion error') if (all(dEq0(invConnections))) call IO_error(-1,ext_msg='back stress calculation: inversion error')
@ -705,11 +700,11 @@ module subroutine nonlocal_dependentState(ph, en)
rhoExcessGradient(2) = rhoExcessGradient(2) + sum(rho(s,imm_scr)) / FVsize rhoExcessGradient(2) = rhoExcessGradient(2) + sum(rho(s,imm_scr)) / FVsize
! ... normalized with the total density ... ! ... normalized with the total density ...
rhoTotal(1) = (sum(abs(rho(s,edg))) + sum(neighbor_rhoTotal(1,s,:))) / (1.0_pREAL + nRealNeighbors) rho_sum(1) = (sum(abs(rho(s,edg))) + sum(rho_sum_nbr(1,s,:))) / (1.0_pREAL + nRealNeighbors)
rhoTotal(2) = (sum(abs(rho(s,scr))) + sum(neighbor_rhoTotal(2,s,:))) / (1.0_pREAL + nRealNeighbors) rho_sum(2) = (sum(abs(rho(s,scr))) + sum(rho_sum_nbr(2,s,:))) / (1.0_pREAL + nRealNeighbors)
rhoExcessGradient_over_rho = 0.0_pREAL rhoExcessGradient_over_rho = 0.0_pREAL
where(rhoTotal > 0.0_pREAL) rhoExcessGradient_over_rho = rhoExcessGradient / rhoTotal where(rho_sum > 0.0_pREAL) rhoExcessGradient_over_rho = rhoExcessGradient / rho_sum
! ... gives the local stress correction when multiplied with a factor ! ... gives the local stress correction when multiplied with a factor
dst%tau_back(s,en) = - mu * prm%b_sl(s) / (2.0_pREAL * PI) & dst%tau_back(s,en) = - mu * prm%b_sl(s) / (2.0_pREAL * PI) &
@ -743,7 +738,7 @@ module subroutine nonlocal_LpAndItsTangent(Lp,dLp_dMp, &
t, & !< dislocation type t, & !< dislocation type
s !< index of my current slip system s !< index of my current slip system
real(pREAL), dimension(param(ph)%sum_N_sl,8) :: & real(pREAL), dimension(param(ph)%sum_N_sl,8) :: &
rhoSgl !< single dislocation densities (including blocked) rho_sgl !< single dislocation densities (including blocked)
real(pREAL), dimension(param(ph)%sum_N_sl,10) :: & real(pREAL), dimension(param(ph)%sum_N_sl,10) :: &
rho rho
real(pREAL), dimension(param(ph)%sum_N_sl,4) :: & real(pREAL), dimension(param(ph)%sum_N_sl,4) :: &
@ -766,7 +761,7 @@ module subroutine nonlocal_LpAndItsTangent(Lp,dLp_dMp, &
!*** shortcut to state variables !*** shortcut to state variables
rho = getRho(ph,en) rho = getRho(ph,en)
rhoSgl = rho(:,sgl) rho_sgl = rho(:,sgl)
do s = 1,prm%sum_N_sl do s = 1,prm%sum_N_sl
tau(s) = math_tensordot(Mp, prm%P_sl(1:3,1:3,s)) tau(s) = math_tensordot(Mp, prm%P_sl(1:3,1:3,s))
@ -799,20 +794,20 @@ module subroutine nonlocal_LpAndItsTangent(Lp,dLp_dMp, &
stt%v(:,en) = pack(v,.true.) stt%v(:,en) = pack(v,.true.)
!*** Bauschinger effect !*** Bauschinger effect
forall (s = 1:prm%sum_N_sl, t = 5:8, rhoSgl(s,t) * v(s,t-4) < 0.0_pREAL) & forall (s = 1:prm%sum_N_sl, t = 5:8, rho_sgl(s,t) * v(s,t-4) < 0.0_pREAL) &
rhoSgl(s,t-4) = rhoSgl(s,t-4) + abs(rhoSgl(s,t)) rho_sgl(s,t-4) = rho_sgl(s,t-4) + abs(rho_sgl(s,t))
dot_gamma = sum(rhoSgl(:,1:4) * v, 2) * prm%b_sl dot_gamma = sum(rho_sgl(:,1:4) * v, 2) * prm%b_sl
do s = 1,prm%sum_N_sl do s = 1,prm%sum_N_sl
Lp = Lp + dot_gamma(s) * prm%P_sl(1:3,1:3,s) Lp = Lp + dot_gamma(s) * prm%P_sl(1:3,1:3,s)
forall (i=1:3,j=1:3,k=1:3,l=1:3) & forall (i=1:3,j=1:3,k=1:3,l=1:3) &
dLp_dMp(i,j,k,l) = dLp_dMp(i,j,k,l) & dLp_dMp(i,j,k,l) = dLp_dMp(i,j,k,l) &
+ prm%P_sl(i,j,s) * prm%P_sl(k,l,s) & + prm%P_sl(i,j,s) * prm%P_sl(k,l,s) &
* sum(rhoSgl(s,1:4) * dv_dtau(s,1:4)) * prm%b_sl(s) & * sum(rho_sgl(s,1:4) * dv_dtau(s,1:4)) * prm%b_sl(s) &
+ prm%P_sl(i,j,s) & + prm%P_sl(i,j,s) &
* (+ prm%P_nS_pos(k,l,s) * rhoSgl(s,3) * dv_dtauNS(s,3) & * (+ prm%P_nS_pos(k,l,s) * rho_sgl(s,3) * dv_dtauNS(s,3) &
- prm%P_nS_neg(k,l,s) * rhoSgl(s,4) * dv_dtauNS(s,4)) * prm%b_sl(s) - prm%P_nS_neg(k,l,s) * rho_sgl(s,4) * dv_dtauNS(s,4)) * prm%b_sl(s)
end do end do
end associate end associate
@ -848,7 +843,7 @@ module subroutine plastic_nonlocal_deltaState(Mp,ph,en)
real(pREAL), dimension(param(ph)%sum_N_sl) :: & real(pREAL), dimension(param(ph)%sum_N_sl) :: &
tau ! current resolved shear stress tau ! current resolved shear stress
real(pREAL), dimension(param(ph)%sum_N_sl,2) :: & real(pREAL), dimension(param(ph)%sum_N_sl,2) :: &
rhoDip, & ! current dipole dislocation densities (screw and edge dipoles) rho_dip, & ! current dipole dislocation densities (screw and edge dipoles)
dUpper, & ! current maximum stable dipole distance for edges and screws dUpper, & ! current maximum stable dipole distance for edges and screws
dUpperOld, & ! old maximum stable dipole distance for edges and screws dUpperOld, & ! old maximum stable dipole distance for edges and screws
deltaDUpper ! change in maximum stable dipole distance for edges and screws deltaDUpper ! change in maximum stable dipole distance for edges and screws
@ -864,7 +859,7 @@ module subroutine plastic_nonlocal_deltaState(Mp,ph,en)
dUpperOld = reshape(dst%max_dipole_height(:,en),[prm%sum_N_sl,2]) dUpperOld = reshape(dst%max_dipole_height(:,en),[prm%sum_N_sl,2])
rho = getRho(ph,en) rho = getRho(ph,en)
rhoDip = rho(:,dip) rho_dip = rho(:,dip)
!**************************************************************************** !****************************************************************************
!*** dislocation remobilization (bauschinger effect) !*** dislocation remobilization (bauschinger effect)
@ -904,7 +899,7 @@ module subroutine plastic_nonlocal_deltaState(Mp,ph,en)
deltaRhoDipole2SingleStress = 0.0_pREAL deltaRhoDipole2SingleStress = 0.0_pREAL
forall (c=1:2, s=1:prm%sum_N_sl, deltaDUpper(s,c) < 0.0_pREAL .and. & forall (c=1:2, s=1:prm%sum_N_sl, deltaDUpper(s,c) < 0.0_pREAL .and. &
dNeq0(dUpperOld(s,c) - prm%minDipoleHeight(s,c))) & dNeq0(dUpperOld(s,c) - prm%minDipoleHeight(s,c))) &
deltaRhoDipole2SingleStress(s,8+c) = rhoDip(s,c) * deltaDUpper(s,c) & deltaRhoDipole2SingleStress(s,8+c) = rho_dip(s,c) * deltaDUpper(s,c) &
/ (dUpperOld(s,c) - prm%minDipoleHeight(s,c)) / (dUpperOld(s,c) - prm%minDipoleHeight(s,c))
forall (t=1:4) deltaRhoDipole2SingleStress(:,t) = -0.5_pREAL * deltaRhoDipole2SingleStress(:,(t-1)/2+9) forall (t=1:4) deltaRhoDipole2SingleStress(:,t) = -0.5_pREAL * deltaRhoDipole2SingleStress(:,(t-1)/2+9)
@ -938,24 +933,21 @@ module subroutine nonlocal_dotState(Mp,timestep, &
s !< index of my current slip system s !< index of my current slip system
real(pREAL), dimension(param(ph)%sum_N_sl,10) :: & real(pREAL), dimension(param(ph)%sum_N_sl,10) :: &
rho, & rho, &
rho0, & !< dislocation density at beginning of time step
rhoDot, & !< density evolution rhoDot, & !< density evolution
rhoDotMultiplication, & !< density evolution by multiplication rhoDotMultiplication, & !< density evolution by multiplication
rhoDotSingle2DipoleGlide, & !< density evolution by dipole formation (by glide) rhoDotSingle2DipoleGlide, & !< density evolution by dipole formation (by glide)
rhoDotAthermalAnnihilation, & !< density evolution by athermal annihilation rhoDotAthermalAnnihilation, & !< density evolution by athermal annihilation
rhoDotThermalAnnihilation !< density evolution by thermal annihilation rhoDotThermalAnnihilation !< density evolution by thermal annihilation
real(pREAL), dimension(param(ph)%sum_N_sl,8) :: & real(pREAL), dimension(param(ph)%sum_N_sl,8) :: &
rhoSgl, & !< current single dislocation densities (positive/negative screw and edge without dipoles) rho_sgl !< current single dislocation densities (positive/negative screw and edge without dipoles)
my_rhoSgl0 !< single dislocation densities of central ip (positive/negative screw and edge without dipoles)
real(pREAL), dimension(param(ph)%sum_N_sl,4) :: & real(pREAL), dimension(param(ph)%sum_N_sl,4) :: &
v, & !< current dislocation glide velocity v, & !< current dislocation glide velocity
v0, &
dot_gamma !< shear rates dot_gamma !< shear rates
real(pREAL), dimension(param(ph)%sum_N_sl) :: & real(pREAL), dimension(param(ph)%sum_N_sl) :: &
tau, & !< current resolved shear stress tau, & !< current resolved shear stress
v_climb !< climb velocity of edge dipoles v_climb !< climb velocity of edge dipoles
real(pREAL), dimension(param(ph)%sum_N_sl,2) :: & real(pREAL), dimension(param(ph)%sum_N_sl,2) :: &
rhoDip, & !< current dipole dislocation densities (screw and edge dipoles) rho_dip, & !< current dipole dislocation densities (screw and edge dipoles)
dLower, & !< minimum stable dipole distance for edges and screws dLower, & !< minimum stable dipole distance for edges and screws
dUpper !< current maximum stable dipole distance for edges and screws dUpper !< current maximum stable dipole distance for edges and screws
real(pREAL) :: & real(pREAL) :: &
@ -978,14 +970,12 @@ module subroutine nonlocal_dotState(Mp,timestep, &
tau = 0.0_pREAL tau = 0.0_pREAL
dot_gamma = 0.0_pREAL dot_gamma = 0.0_pREAL
rho = getRho(ph,en) rho = getRho(ph,en)
rhoSgl = rho(:,sgl) rho_sgl = rho(:,sgl)
rhoDip = rho(:,dip) rho_dip = rho(:,dip)
rho0 = getRho0(ph,en)
my_rhoSgl0 = rho0(:,sgl)
v = reshape(stt%v(:,en),[prm%sum_N_sl,4]) v = reshape(stt%v(:,en),[prm%sum_N_sl,4])
dot_gamma = rhoSgl(:,1:4) * v * spread(prm%b_sl,2,4) dot_gamma = rho_sgl(:,1:4) * v * spread(prm%b_sl,2,4)
! limits for stable dipole height ! limits for stable dipole height
@ -1024,8 +1014,6 @@ module subroutine nonlocal_dotState(Mp,timestep, &
* sqrt(dst%rho_forest(:,en)) / prm%i_sl / prm%b_sl, 2, 4) ! eq. 3.26 * sqrt(dst%rho_forest(:,en)) / prm%i_sl / prm%b_sl, 2, 4) ! eq. 3.26
end if isBCC end if isBCC
v0 = reshape(st0%v(:,en),[prm%sum_N_sl,4])
!**************************************************************************** !****************************************************************************
! dipole formation and annihilation ! dipole formation and annihilation
@ -1033,20 +1021,20 @@ module subroutine nonlocal_dotState(Mp,timestep, &
! formation by glide ! formation by glide
do c = 1,2 do c = 1,2
rhoDotSingle2DipoleGlide(:,2*c-1) = -2.0_pREAL * dUpper(:,c) / prm%b_sl & rhoDotSingle2DipoleGlide(:,2*c-1) = -2.0_pREAL * dUpper(:,c) / prm%b_sl &
* ( rhoSgl(:,2*c-1) * abs(dot_gamma(:,2*c)) & ! negative mobile --> positive mobile * ( rho_sgl(:,2*c-1) * abs(dot_gamma(:,2*c)) & ! negative mobile --> positive mobile
+ rhoSgl(:,2*c) * abs(dot_gamma(:,2*c-1)) & ! positive mobile --> negative mobile + rho_sgl(:,2*c) * abs(dot_gamma(:,2*c-1)) & ! positive mobile --> negative mobile
+ abs(rhoSgl(:,2*c+4)) * abs(dot_gamma(:,2*c-1))) ! positive mobile --> negative immobile + abs(rho_sgl(:,2*c+4)) * abs(dot_gamma(:,2*c-1))) ! positive mobile --> negative immobile
rhoDotSingle2DipoleGlide(:,2*c) = -2.0_pREAL * dUpper(:,c) / prm%b_sl & rhoDotSingle2DipoleGlide(:,2*c) = -2.0_pREAL * dUpper(:,c) / prm%b_sl &
* ( rhoSgl(:,2*c-1) * abs(dot_gamma(:,2*c)) & ! negative mobile --> positive mobile * ( rho_sgl(:,2*c-1) * abs(dot_gamma(:,2*c)) & ! negative mobile --> positive mobile
+ rhoSgl(:,2*c) * abs(dot_gamma(:,2*c-1)) & ! positive mobile --> negative mobile + rho_sgl(:,2*c) * abs(dot_gamma(:,2*c-1)) & ! positive mobile --> negative mobile
+ abs(rhoSgl(:,2*c+3)) * abs(dot_gamma(:,2*c))) ! negative mobile --> positive immobile + abs(rho_sgl(:,2*c+3)) * abs(dot_gamma(:,2*c))) ! negative mobile --> positive immobile
rhoDotSingle2DipoleGlide(:,2*c+3) = -2.0_pREAL * dUpper(:,c) / prm%b_sl & rhoDotSingle2DipoleGlide(:,2*c+3) = -2.0_pREAL * dUpper(:,c) / prm%b_sl &
* rhoSgl(:,2*c+3) * abs(dot_gamma(:,2*c)) ! negative mobile --> positive immobile * rho_sgl(:,2*c+3) * abs(dot_gamma(:,2*c)) ! negative mobile --> positive immobile
rhoDotSingle2DipoleGlide(:,2*c+4) = -2.0_pREAL * dUpper(:,c) / prm%b_sl & rhoDotSingle2DipoleGlide(:,2*c+4) = -2.0_pREAL * dUpper(:,c) / prm%b_sl &
* rhoSgl(:,2*c+4) * abs(dot_gamma(:,2*c-1)) ! positive mobile --> negative immobile * rho_sgl(:,2*c+4) * abs(dot_gamma(:,2*c-1)) ! positive mobile --> negative immobile
rhoDotSingle2DipoleGlide(:,c+8) = abs(rhoDotSingle2DipoleGlide(:,2*c+3)) & rhoDotSingle2DipoleGlide(:,c+8) = abs(rhoDotSingle2DipoleGlide(:,2*c+3)) &
+ abs(rhoDotSingle2DipoleGlide(:,2*c+4)) & + abs(rhoDotSingle2DipoleGlide(:,2*c+4)) &
@ -1059,9 +1047,9 @@ module subroutine nonlocal_dotState(Mp,timestep, &
rhoDotAthermalAnnihilation = 0.0_pREAL rhoDotAthermalAnnihilation = 0.0_pREAL
forall (c=1:2) & forall (c=1:2) &
rhoDotAthermalAnnihilation(:,c+8) = -2.0_pREAL * dLower(:,c) / prm%b_sl & rhoDotAthermalAnnihilation(:,c+8) = -2.0_pREAL * dLower(:,c) / prm%b_sl &
* ( 2.0_pREAL * (rhoSgl(:,2*c-1) * abs(dot_gamma(:,2*c)) + rhoSgl(:,2*c) * abs(dot_gamma(:,2*c-1))) & ! was single hitting single * ( 2.0_pREAL * (rho_sgl(:,2*c-1) * abs(dot_gamma(:,2*c)) + rho_sgl(:,2*c) * abs(dot_gamma(:,2*c-1))) & ! was single hitting single
+ 2.0_pREAL * (abs(rhoSgl(:,2*c+3)) * abs(dot_gamma(:,2*c)) + abs(rhoSgl(:,2*c+4)) * abs(dot_gamma(:,2*c-1))) & ! was single hitting immobile single or was immobile single hit by single + 2.0_pREAL * (abs(rho_sgl(:,2*c+3)) * abs(dot_gamma(:,2*c)) + abs(rho_sgl(:,2*c+4)) * abs(dot_gamma(:,2*c-1))) & ! was single hitting immobile single or was immobile single hit by single
+ rhoDip(:,c) * (abs(dot_gamma(:,2*c-1)) + abs(dot_gamma(:,2*c)))) ! single knocks dipole constituent + rho_dip(:,c) * (abs(dot_gamma(:,2*c-1)) + abs(dot_gamma(:,2*c)))) ! single knocks dipole constituent
! annihilated screw dipoles leave edge jogs behind on the colinear system ! annihilated screw dipoles leave edge jogs behind on the colinear system
if (phase_lattice(ph) == 'cF') & if (phase_lattice(ph) == 'cF') &
@ -1076,8 +1064,8 @@ module subroutine nonlocal_dotState(Mp,timestep, &
v_climb = D_SD * mu * prm%V_at & v_climb = D_SD * mu * prm%V_at &
/ (PI * (1.0_pREAL-nu) * (dUpper(:,1) + dLower(:,1)) * K_B * Temperature) ! eq. 3.54 / (PI * (1.0_pREAL-nu) * (dUpper(:,1) + dLower(:,1)) * K_B * Temperature) ! eq. 3.54
forall (s = 1:prm%sum_N_sl, dUpper(s,1) > dLower(s,1)) & forall (s = 1:prm%sum_N_sl, dUpper(s,1) > dLower(s,1)) &
rhoDotThermalAnnihilation(s,9) = max(- 4.0_pREAL * rhoDip(s,1) * v_climb(s) / (dUpper(s,1) - dLower(s,1)), & rhoDotThermalAnnihilation(s,9) = max(- 4.0_pREAL * rho_dip(s,1) * v_climb(s) / (dUpper(s,1) - dLower(s,1)), &
- rhoDip(s,1) / timestep - rhoDotAthermalAnnihilation(s,9) & - rho_dip(s,1) / timestep - rhoDotAthermalAnnihilation(s,9) &
- rhoDotSingle2DipoleGlide(s,9)) ! make sure that we do not annihilate more dipoles than we have - rhoDotSingle2DipoleGlide(s,9)) ! make sure that we do not annihilate more dipoles than we have
rhoDot = rhoDotFlux(timestep, ph,en) & rhoDot = rhoDotFlux(timestep, ph,en) &
@ -1118,38 +1106,37 @@ function rhoDotFlux(timestep,ph,en)
ns, & !< short notation for the total number of active slip systems ns, & !< short notation for the total number of active slip systems
c, & !< character of dislocation c, & !< character of dislocation
n, & !< index of my current neighbor n, & !< index of my current neighbor
neighbor_el, & !< element number of my neighbor el_nbr, & !< element number of my neighbor
neighbor_ip, & !< integration point of my neighbor ip_nbr, & !< integration point of my neighbor
neighbor_n, & !< neighbor index pointing to en when looking from my neighbor n_nbr, & !< neighbor index pointing to en when looking from my neighbor
opposite_neighbor, & !< index of my opposite neighbor opposite_neighbor, & !< index of my opposite neighbor
opposite_ip, & !< ip of my opposite neighbor opposite_ip, & !< ip of my opposite neighbor
opposite_el, & !< element index of my opposite neighbor opposite_el, & !< element index of my opposite neighbor
opposite_n, & !< neighbor index pointing to en when looking from my opposite neighbor opposite_n, & !< neighbor index pointing to en when looking from my opposite neighbor
t, & !< type of dislocation t, & !< type of dislocation
no,& !< neighbor offset shortcut en_nbr,& !< neighbor phase entry
np,& !< neighbor phase shortcut ph_nbr,& !< neighbor phase ID
topp, & !< type of dislocation with opposite sign to t topp, & !< type of dislocation with opposite sign to t
s !< index of my current slip system s !< index of my current slip system
real(pREAL), dimension(param(ph)%sum_N_sl,10) :: & real(pREAL), dimension(param(ph)%sum_N_sl,10) :: &
rho, & rho, &
rho0, & !< dislocation density at beginning of time step rho_0, & !< dislocation density at beginning of time step
rhoDotFlux !< density evolution by flux rhoDotFlux !< density evolution by flux
real(pREAL), dimension(param(ph)%sum_N_sl,8) :: & real(pREAL), dimension(param(ph)%sum_N_sl,8) :: &
rhoSgl, & !< current single dislocation densities (positive/negative screw and edge without dipoles) rho_0_sgl_nbr, & !< current single dislocation densities of neighboring ip (positive/negative screw and edge without dipoles)
neighbor_rhoSgl0, & !< current single dislocation densities of neighboring ip (positive/negative screw and edge without dipoles) rho_0_sgl !< single dislocation densities of central ip (positive/negative screw and edge without dipoles)
my_rhoSgl0 !< single dislocation densities of central ip (positive/negative screw and edge without dipoles)
real(pREAL), dimension(param(ph)%sum_N_sl,4) :: & real(pREAL), dimension(param(ph)%sum_N_sl,4) :: &
v, & !< current dislocation glide velocity v, & !< current dislocation glide velocity
v0, & v_0, &
neighbor_v0, & !< dislocation glide velocity of enighboring ip v_0_nbr, & !< dislocation glide velocity of enighboring ip
dot_gamma !< shear rates dot_gamma !< shear rates
real(pREAL), dimension(3,param(ph)%sum_N_sl,4) :: & real(pREAL), dimension(3,param(ph)%sum_N_sl,4) :: &
m !< direction of dislocation motion m !< direction of dislocation motion
real(pREAL), dimension(3,3) :: & real(pREAL), dimension(3,3) :: &
my_F, & !< my total deformation gradient my_F, & !< my total deformation gradient
neighbor_F, & !< total deformation gradient of my neighbor F_nbr, & !< total deformation gradient of my neighbor
my_Fe, & !< my elastic deformation gradient my_Fe, & !< my elastic deformation gradient
neighbor_Fe, & !< elastic deformation gradient of my neighbor F_e_nbr, & !< elastic deformation gradient of my neighbor
Favg !< average total deformation gradient of en and my neighbor Favg !< average total deformation gradient of en and my neighbor
real(pREAL), dimension(3) :: & real(pREAL), dimension(3) :: &
normal_neighbor2me, & !< interface normal pointing from my neighbor to en in neighbor's lattice configuration normal_neighbor2me, & !< interface normal pointing from my neighbor to en in neighbor's lattice configuration
@ -1157,7 +1144,7 @@ function rhoDotFlux(timestep,ph,en)
normal_me2neighbor, & !< interface normal pointing from en to my neighbor in my lattice configuration normal_me2neighbor, & !< interface normal pointing from en to my neighbor in my lattice configuration
normal_me2neighbor_defConf !< interface normal pointing from en to my neighbor in shared deformed configuration normal_me2neighbor_defConf !< interface normal pointing from en to my neighbor in shared deformed configuration
real(pREAL) :: & real(pREAL) :: &
area, & !< area of the current interface a, & !< area of the current interface
transmissivity, & !< overall transmissivity of dislocation flux to neighboring material point transmissivity, & !< overall transmissivity of dislocation flux to neighboring material point
lineLength !< dislocation line length leaving the current interface lineLength !< dislocation line length leaving the current interface
@ -1166,19 +1153,19 @@ function rhoDotFlux(timestep,ph,en)
dst => dependentState(ph), & dst => dependentState(ph), &
stt => state(ph), & stt => state(ph), &
st0 => state0(ph)) st0 => state0(ph))
ns = prm%sum_N_sl ns = prm%sum_N_sl
dot_gamma = 0.0_pREAL dot_gamma = 0.0_pREAL
rho = getRho(ph,en) rho = getRho(ph,en)
rhoSgl = rho(:,sgl) rho_0 = getRho0(ph,en)
rho0 = getRho0(ph,en) rho_0_sgl = rho_0(:,sgl)
my_rhoSgl0 = rho0(:,sgl)
v = reshape(stt%v(:,en),[prm%sum_N_sl,4]) !ToDo: MD: I think we should use state0 here v = reshape(stt%v(:,en),[prm%sum_N_sl,4]) !ToDo: MD: I think we should use state0 here
dot_gamma = rhoSgl(:,1:4) * v * spread(prm%b_sl,2,4) dot_gamma = rho(:,mob) * v * spread(prm%b_sl,2,4)
v0 = reshape(st0%v(:,en),[prm%sum_N_sl,4]) v_0 = reshape(st0%v(:,en),[prm%sum_N_sl,4])
!**************************************************************************** !****************************************************************************
!*** calculate dislocation fluxes (only for nonlocal plasticity) !*** calculate dislocation fluxes (only for nonlocal plasticity)
@ -1187,8 +1174,8 @@ function rhoDotFlux(timestep,ph,en)
!*** check CFL (Courant-Friedrichs-Lewy) condition for flux !*** check CFL (Courant-Friedrichs-Lewy) condition for flux
if (any( abs(dot_gamma) > 0.0_pREAL & ! any active slip system ... if (any( abs(dot_gamma) > 0.0_pREAL & ! any active slip system ...
.and. prm%C_CFL * abs(v0) * timestep & .and. prm%C_CFL * abs(v_0) * timestep &
> geom(ph)%V_0(en)/ maxval(geom(ph)%IParea(:,en)))) then ! ...with velocity above critical value (we use the reference volume and area for simplicity here) > geom(ph)%v_0(en)/ maxval(geom(ph)%a_0(:,en)))) then ! ...with velocity above critical value (we use the reference volume and area for simplicity here)
rhoDotFlux = IEEE_value(1.0_pREAL,IEEE_quiet_NaN) ! enforce cutback rhoDotFlux = IEEE_value(1.0_pREAL,IEEE_quiet_NaN) ! enforce cutback
return return
end if end if
@ -1205,28 +1192,28 @@ function rhoDotFlux(timestep,ph,en)
my_F = phase_mechanical_F(ph)%data(1:3,1:3,en) my_F = phase_mechanical_F(ph)%data(1:3,1:3,en)
my_Fe = matmul(my_F, math_inv33(phase_mechanical_Fp(ph)%data(1:3,1:3,en))) my_Fe = matmul(my_F, math_inv33(phase_mechanical_Fp(ph)%data(1:3,1:3,en)))
neighbors: do n = 1,nIPneighbors neighbors: do n = 1,nCellNeighbors
neighbor_el = geom(ph)%IPneighborhood(1,n,en) el_nbr = geom(ph)%IPneighborhood(1,n,en)
neighbor_ip = geom(ph)%IPneighborhood(2,n,en) ip_nbr = geom(ph)%IPneighborhood(2,n,en)
neighbor_n = geom(ph)%IPneighborhood(3,n,en) n_nbr = geom(ph)%IPneighborhood(3,n,en)
np = material_ID_phase(1,(neighbor_el-1)*discretization_nIPs + neighbor_ip) ph_nbr = material_ID_phase(1,(el_nbr-1)*discretization_nIPs + ip_nbr)
no = material_entry_phase(1,(neighbor_el-1)*discretization_nIPs + neighbor_ip) en_nbr = material_entry_phase(1,(el_nbr-1)*discretization_nIPs + ip_nbr)
opposite_neighbor = n + mod(n,2) - mod(n+1,2) opposite_neighbor = n + mod(n,2) - mod(n+1,2)
opposite_el = geom(ph)%IPneighborhood(1,opposite_neighbor,en) opposite_el = geom(ph)%IPneighborhood(1,opposite_neighbor,en)
opposite_ip = geom(ph)%IPneighborhood(2,opposite_neighbor,en) opposite_ip = geom(ph)%IPneighborhood(2,opposite_neighbor,en)
opposite_n = geom(ph)%IPneighborhood(3,opposite_neighbor,en) opposite_n = geom(ph)%IPneighborhood(3,opposite_neighbor,en)
if (neighbor_n > 0) then ! if neighbor exists, average deformation gradient if (n_nbr > 0) then ! if neighbor exists, average deformation gradient
neighbor_F = phase_mechanical_F(np)%data(1:3,1:3,no) F_nbr = phase_mechanical_F(ph_nbr)%data(1:3,1:3,en_nbr)
neighbor_Fe = matmul(neighbor_F, math_inv33(phase_mechanical_Fp(np)%data(1:3,1:3,no))) F_e_nbr = matmul(F_nbr, math_inv33(phase_mechanical_Fp(ph_nbr)%data(1:3,1:3,en_nbr)))
Favg = 0.5_pREAL * (my_F + neighbor_F) Favg = 0.5_pREAL * (my_F + F_nbr)
else ! if no neighbor, take my value as average else ! if no neighbor, take my value as average
Favg = my_F Favg = my_F
end if end if
neighbor_v0 = 0.0_pREAL ! needed for check of sign change in flux density below v_0_nbr = 0.0_pREAL ! needed for check of sign change in flux density below
!* FLUX FROM MY NEIGHBOR TO ME !* FLUX FROM MY NEIGHBOR TO ME
!* This is only considered, if I have a neighbor of nonlocal plasticity !* This is only considered, if I have a neighbor of nonlocal plasticity
@ -1236,38 +1223,38 @@ function rhoDotFlux(timestep,ph,en)
!* my neighbor's interface. !* my neighbor's interface.
!* The entering flux from my neighbor will be distributed on my slip systems according to the !* The entering flux from my neighbor will be distributed on my slip systems according to the
!* compatibility !* compatibility
if (neighbor_n > 0) then if (n_nbr > 0) then
if (mechanical_plasticity_type(np) == MECHANICAL_PLASTICITY_NONLOCAL .and. & if (mechanical_plasticity_type(ph_nbr) == MECHANICAL_PLASTICITY_NONLOCAL .and. &
any(dependentState(ph)%compatibility(:,:,:,n,en) > 0.0_pREAL)) then any(dependentState(ph)%compatibility(:,:,:,n,en) > 0.0_pREAL)) then
forall (s = 1:ns, t = 1:4) forall (s = 1:ns, t = 1:4)
neighbor_v0(s,t) = plasticState(np)%state0(iV (s,t,np),no) v_0_nbr(s,t) = plasticState(ph_nbr)%state0(iV (s,t,ph_nbr),en_nbr)
neighbor_rhoSgl0(s,t) = max(plasticState(np)%state0(iRhoU(s,t,np),no),0.0_pREAL) rho_0_sgl_nbr(s,t) = max(plasticState(ph_nbr)%state0(iRhoU(s,t,ph_nbr),en_nbr),0.0_pREAL)
endforall endforall
where (neighbor_rhoSgl0 * IPvolume(neighbor_ip,neighbor_el) ** 0.667_pREAL < prm%rho_min & where (rho_0_sgl_nbr * IPvolume0(ip_nbr,el_nbr) ** 0.667_pREAL < prm%rho_min &
.or. neighbor_rhoSgl0 < prm%rho_significant) & .or. rho_0_sgl_nbr < prm%rho_significant) &
neighbor_rhoSgl0 = 0.0_pREAL rho_0_sgl_nbr = 0.0_pREAL
normal_neighbor2me_defConf = math_det33(Favg) * matmul(math_inv33(transpose(Favg)), & normal_neighbor2me_defConf = math_det33(Favg) * matmul(math_inv33(transpose(Favg)), &
IPareaNormal(1:3,neighbor_n,neighbor_ip,neighbor_el)) ! normal of the interface in (average) deformed configuration (pointing neighbor => en) IPareaNormal0(1:3,n_nbr,ip_nbr,el_nbr)) ! normal of the interface in (average) deformed configuration (pointing neighbor => en)
normal_neighbor2me = matmul(transpose(neighbor_Fe), normal_neighbor2me_defConf) & normal_neighbor2me = matmul(transpose(F_e_nbr), normal_neighbor2me_defConf) &
/ math_det33(neighbor_Fe) ! interface normal in the lattice configuration of my neighbor / math_det33(F_e_nbr) ! interface normal in the lattice configuration of my neighbor
area = IParea(neighbor_n,neighbor_ip,neighbor_el) * norm2(normal_neighbor2me) a = IParea0(n_nbr,ip_nbr,el_nbr) * norm2(normal_neighbor2me)
normal_neighbor2me = normal_neighbor2me / norm2(normal_neighbor2me) ! normalize the surface normal to unit length normal_neighbor2me = normal_neighbor2me / norm2(normal_neighbor2me) ! normalize the surface normal to unit length
do s = 1,ns do s = 1,ns
do t = 1,4 do t = 1,4
c = (t + 1) / 2 c = (t + 1) / 2
topp = t + mod(t,2) - mod(t+1,2) topp = t + mod(t,2) - mod(t+1,2)
if (neighbor_v0(s,t) * math_inner(m(1:3,s,t), normal_neighbor2me) > 0.0_pREAL & ! flux from my neighbor to en == entering flux for en if (v_0_nbr(s,t) * math_inner(m(1:3,s,t), normal_neighbor2me) > 0.0_pREAL & ! flux from my neighbor to en == entering flux for en
.and. v0(s,t) * neighbor_v0(s,t) >= 0.0_pREAL ) then ! ... only if no sign change in flux density .and. v_0(s,t) * v_0_nbr(s,t) >= 0.0_pREAL ) then ! ... only if no sign change in flux density
lineLength = neighbor_rhoSgl0(s,t) * neighbor_v0(s,t) & lineLength = rho_0_sgl_nbr(s,t) * v_0_nbr(s,t) &
* math_inner(m(1:3,s,t), normal_neighbor2me) * area ! positive line length that wants to enter through this interface * math_inner(m(1:3,s,t), normal_neighbor2me) * a ! positive line length that wants to enter through this interface
where (dependentState(ph)%compatibility(c,:,s,n,en) > 0.0_pREAL) & where (dependentState(ph)%compatibility(c,:,s,n,en) > 0.0_pREAL) &
rhoDotFlux(:,t) = rhoDotFlux(1:ns,t) & rhoDotFlux(:,t) = rhoDotFlux(1:ns,t) &
+ lineLength/geom(ph)%V_0(en)*dependentState(ph)%compatibility(c,:,s,n,en)**2 ! transferring to equally signed mobile dislocation type + lineLength/geom(ph)%v_0(en)*dependentState(ph)%compatibility(c,:,s,n,en)**2 ! transferring to equally signed mobile dislocation type
where (dependentState(ph)%compatibility(c,:,s,n,en) < 0.0_pREAL) & where (dependentState(ph)%compatibility(c,:,s,n,en) < 0.0_pREAL) &
rhoDotFlux(:,topp) = rhoDotFlux(:,topp) & rhoDotFlux(:,topp) = rhoDotFlux(:,topp) &
+ lineLength/geom(ph)%V_0(en)*dependentState(ph)%compatibility(c,:,s,n,en)**2 ! transferring to opposite signed mobile dislocation type + lineLength/geom(ph)%v_0(en)*dependentState(ph)%compatibility(c,:,s,n,en)**2 ! transferring to opposite signed mobile dislocation type
end if end if
end do end do
@ -1283,29 +1270,29 @@ function rhoDotFlux(timestep,ph,en)
!* In case of reduced transmissivity, part of the leaving flux is stored as dead dislocation density. !* In case of reduced transmissivity, part of the leaving flux is stored as dead dislocation density.
!* That means for an interface of zero transmissivity the leaving flux is fully converted to dead dislocations. !* That means for an interface of zero transmissivity the leaving flux is fully converted to dead dislocations.
if (opposite_n > 0) then if (opposite_n > 0) then
if (mechanical_plasticity_type(np) == MECHANICAL_PLASTICITY_NONLOCAL) then if (mechanical_plasticity_type(ph_nbr) == MECHANICAL_PLASTICITY_NONLOCAL) then
normal_me2neighbor_defConf = math_det33(Favg) & normal_me2neighbor_defConf = math_det33(Favg) &
* matmul(math_inv33(transpose(Favg)),geom(ph)%IPareaNormal(1:3,n,en)) ! normal of the interface in (average) deformed configuration (pointing en => neighbor) * matmul(math_inv33(transpose(Favg)),geom(ph)%n_0(1:3,n,en)) ! normal of the interface in (average) deformed configuration (pointing en => neighbor)
normal_me2neighbor = matmul(transpose(my_Fe), normal_me2neighbor_defConf) & normal_me2neighbor = matmul(transpose(my_Fe), normal_me2neighbor_defConf) &
/ math_det33(my_Fe) ! interface normal in my lattice configuration / math_det33(my_Fe) ! interface normal in my lattice configuration
area = geom(ph)%IParea(n,en) * norm2(normal_me2neighbor) a = geom(ph)%a_0(n,en) * norm2(normal_me2neighbor)
normal_me2neighbor = normal_me2neighbor / norm2(normal_me2neighbor) ! normalize the surface normal to unit length normal_me2neighbor = normal_me2neighbor / norm2(normal_me2neighbor) ! normalize the surface normal to unit length
do s = 1,ns do s = 1,ns
do t = 1,4 do t = 1,4
c = (t + 1) / 2 c = (t + 1) / 2
if (v0(s,t) * math_inner(m(1:3,s,t), normal_me2neighbor) > 0.0_pREAL ) then ! flux from en to my neighbor == leaving flux for en (might also be a pure flux from my mobile density to dead density if interface not at all transmissive) if (v_0(s,t) * math_inner(m(1:3,s,t), normal_me2neighbor) > 0.0_pREAL ) then ! flux from en to my neighbor == leaving flux for en (might also be a pure flux from my mobile density to dead density if interface not at all transmissive)
if (v0(s,t) * neighbor_v0(s,t) >= 0.0_pREAL) then ! no sign change in flux density if (v_0(s,t) * v_0_nbr(s,t) >= 0.0_pREAL) then ! no sign change in flux density
transmissivity = sum(dependentState(ph)%compatibility(c,:,s,n,en)**2) ! overall transmissivity from this slip system to my neighbor transmissivity = sum(dependentState(ph)%compatibility(c,:,s,n,en)**2) ! overall transmissivity from this slip system to my neighbor
else ! sign change in flux density means sign change in stress which does not allow for dislocations to arive at the neighbor else ! sign change in flux density means sign change in stress which does not allow for dislocations to arive at the neighbor
transmissivity = 0.0_pREAL transmissivity = 0.0_pREAL
end if end if
lineLength = my_rhoSgl0(s,t) * v0(s,t) & lineLength = rho_0_sgl(s,t) * v_0(s,t) &
* math_inner(m(1:3,s,t), normal_me2neighbor) * area ! positive line length of mobiles that wants to leave through this interface * math_inner(m(1:3,s,t), normal_me2neighbor) * a ! positive line length of mobiles that wants to leave through this interface
rhoDotFlux(s,t) = rhoDotFlux(s,t) - lineLength / geom(ph)%V_0(en) ! subtract dislocation flux from current type rhoDotFlux(s,t) = rhoDotFlux(s,t) - lineLength / geom(ph)%v_0(en) ! subtract dislocation flux from current type
rhoDotFlux(s,t+4) = rhoDotFlux(s,t+4) & rhoDotFlux(s,t+4) = rhoDotFlux(s,t+4) &
+ lineLength / geom(ph)%V_0(en) * (1.0_pREAL - transmissivity) & + lineLength / geom(ph)%v_0(en) * (1.0_pREAL - transmissivity) &
* sign(1.0_pREAL, v0(s,t)) ! dislocation flux that is not able to leave through interface (because of low transmissivity) will remain as immobile single density at the material point * sign(1.0_pREAL, v_0(s,t)) ! dislocation flux that is not able to leave through interface (because of low transmissivity) will remain as immobile single density at the material point
end if end if
end do end do
end do end do
@ -1338,14 +1325,14 @@ module subroutine plastic_nonlocal_updateCompatibility(orientation,ce)
en, & en, &
ip, & ip, &
el, & el, &
neighbor_e, & ! element index of my neighbor el_nbr, & ! element index of my neighbor
neighbor_i, & ! integration point index of my neighbor ip_nbr, & ! integration point index of my neighbor
neighbor_me, & en_nbr, &
neighbor_phase, & ph_nbr, &
ns, & ! number of active slip systems ns, & ! number of active slip systems
s1, & ! slip system index (en) s1, & ! slip system index (en)
s2 ! slip system index (my neighbor) s2 ! slip system index (my neighbor)
real(pREAL), dimension(2,param(material_ID_phase(1,ce))%sum_N_sl,param(material_ID_phase(1,ce))%sum_N_sl,nIPneighbors) :: & real(pREAL), dimension(2,param(material_ID_phase(1,ce))%sum_N_sl,param(material_ID_phase(1,ce))%sum_N_sl,nCellNeighbors) :: &
my_compatibility ! my_compatibility for current element and ip my_compatibility ! my_compatibility for current element and ip
real(pREAL) :: & real(pREAL) :: &
my_compatibilitySum, & my_compatibilitySum, &
@ -1368,24 +1355,24 @@ module subroutine plastic_nonlocal_updateCompatibility(orientation,ce)
my_compatibility = 0.0_pREAL my_compatibility = 0.0_pREAL
forall(s1 = 1:ns) my_compatibility(:,s1,s1,:) = 1.0_pREAL forall(s1 = 1:ns) my_compatibility(:,s1,s1,:) = 1.0_pREAL
neighbors: do n = 1,nIPneighbors neighbors: do n = 1,nCellNeighbors
neighbor_e = IPneighborhood(1,n,ip,el) el_nbr = IPneighborhood(1,n,ip,el)
neighbor_i = IPneighborhood(2,n,ip,el) ip_nbr = IPneighborhood(2,n,ip,el)
neighbor_me = material_entry_phase(1,(neighbor_e-1)*discretization_nIPs + neighbor_i) en_nbr = material_entry_phase(1,(el_nbr-1)*discretization_nIPs + ip_nbr)
neighbor_phase = material_ID_phase(1,(neighbor_e-1)*discretization_nIPs + neighbor_i) ph_nbr = material_ID_phase(1,(el_nbr-1)*discretization_nIPs + ip_nbr)
if (neighbor_e <= 0 .or. neighbor_i <= 0) then if (el_nbr <= 0 .or. ip_nbr <= 0) then
!* FREE SURFACE !* FREE SURFACE
forall(s1 = 1:ns) my_compatibility(:,s1,s1,n) = sqrt(prm%chi_surface) forall(s1 = 1:ns) my_compatibility(:,s1,s1,n) = sqrt(prm%chi_surface)
elseif (neighbor_phase /= ph) then elseif (ph_nbr /= ph) then
!* PHASE BOUNDARY !* PHASE BOUNDARY
if (plasticState(neighbor_phase)%nonlocal .and. plasticState(ph)%nonlocal) & if (plasticState(ph_nbr)%nonlocal .and. plasticState(ph)%nonlocal) &
forall(s1 = 1:ns) my_compatibility(:,s1,s1,n) = 0.0_pREAL forall(s1 = 1:ns) my_compatibility(:,s1,s1,n) = 0.0_pREAL
elseif (prm%chi_GB >= 0.0_pREAL) then elseif (prm%chi_GB >= 0.0_pREAL) then
!* GRAIN BOUNDARY !* GRAIN BOUNDARY
if (any(dNeq(phase_O_0(ph)%data(en)%asQuaternion(), & if (any(dNeq(phase_O_0(ph)%data(en)%asQuaternion(), &
phase_O_0(neighbor_phase)%data(neighbor_me)%asQuaternion())) .and. & phase_O_0(ph_nbr)%data(en_nbr)%asQuaternion())) .and. &
plasticState(neighbor_phase)%nonlocal) & plasticState(ph_nbr)%nonlocal) &
forall(s1 = 1:ns) my_compatibility(:,s1,s1,n) = sqrt(prm%chi_GB) forall(s1 = 1:ns) my_compatibility(:,s1,s1,n) = sqrt(prm%chi_GB)
else else
!* GRAIN BOUNDARY ? !* GRAIN BOUNDARY ?
@ -1397,7 +1384,7 @@ module subroutine plastic_nonlocal_updateCompatibility(orientation,ce)
!* the number of compatible slip systems is minimized with the sum of the original compatibility values exceeding one. !* the number of compatible slip systems is minimized with the sum of the original compatibility values exceeding one.
!* Finally the smallest compatibility value is decreased until the sum is exactly equal to one. !* Finally the smallest compatibility value is decreased until the sum is exactly equal to one.
!* All values below the threshold are set to zero. !* All values below the threshold are set to zero.
mis = orientation(ph)%data(en)%misorientation(orientation(neighbor_phase)%data(neighbor_me)) mis = orientation(ph)%data(en)%misorientation(orientation(ph_nbr)%data(en_nbr))
mySlipSystems: do s1 = 1,ns mySlipSystems: do s1 = 1,ns
neighborSlipSystems: do s2 = 1,ns neighborSlipSystems: do s2 = 1,ns
my_compatibility(1,s2,s1,n) = math_inner(prm%slip_normal(1:3,s1), & my_compatibility(1,s2,s1,n) = math_inner(prm%slip_normal(1:3,s1), &
@ -1543,8 +1530,8 @@ subroutine stateInit(ini,phase,Nentries)
associate(stt => state(phase)) associate(stt => state(phase))
if (ini%random_rho_u > 0.0_pREAL) then ! randomly distribute dislocation segments on random slip system and of random type in the volume if (ini%random_rho_u > 0.0_pREAL) then ! randomly distribute dislocation segments on random slip system and of random type in the volume
totalVolume = sum(geom(phase)%V_0) totalVolume = sum(geom(phase)%v_0)
minimumIPVolume = minval(geom(phase)%V_0) minimumIPVolume = minval(geom(phase)%v_0)
densityBinning = ini%random_rho_u_binning / minimumIpVolume ** (2.0_pREAL / 3.0_pREAL) densityBinning = ini%random_rho_u_binning / minimumIpVolume ** (2.0_pREAL / 3.0_pREAL)
! fill random material points with dislocation segments until the desired overall density is reached ! fill random material points with dislocation segments until the desired overall density is reached
@ -1553,8 +1540,8 @@ subroutine stateInit(ini,phase,Nentries)
call random_number(rnd) call random_number(rnd)
e = nint(rnd(1)*real(Nentries,pREAL) + 0.5_pREAL) e = nint(rnd(1)*real(Nentries,pREAL) + 0.5_pREAL)
s = nint(rnd(2)*real(sum(ini%N_sl),pREAL)*4.0_pREAL + 0.5_pREAL) s = nint(rnd(2)*real(sum(ini%N_sl),pREAL)*4.0_pREAL + 0.5_pREAL)
meanDensity = meanDensity + densityBinning * geom(phase)%V_0(e) / totalVolume meanDensity = meanDensity + densityBinning * geom(phase)%v_0(e) / totalVolume
stt%rhoSglMobile(s,e) = densityBinning stt%rho_sgl_mob(s,e) = densityBinning
end do end do
else ! homogeneous distribution with noise else ! homogeneous distribution with noise
do f = 1,size(ini%N_sl,1) do f = 1,size(ini%N_sl,1)
@ -1684,7 +1671,7 @@ pure function getRho(ph,en) result(rho)
rho(:,mob) = max(rho(:,mob),0.0_pREAL) rho(:,mob) = max(rho(:,mob),0.0_pREAL)
rho(:,dip) = max(rho(:,dip),0.0_pREAL) rho(:,dip) = max(rho(:,dip),0.0_pREAL)
where(abs(rho) < max(prm%rho_min/geom(ph)%V_0(en)**(2.0_pREAL/3.0_pREAL),prm%rho_significant)) & where(abs(rho) < max(prm%rho_min/geom(ph)%v_0(en)**(2.0_pREAL/3.0_pREAL),prm%rho_significant)) &
rho = 0.0_pREAL rho = 0.0_pREAL
end associate end associate
@ -1696,54 +1683,57 @@ end function getRho
!> @brief returns copy of current dislocation densities from state !> @brief returns copy of current dislocation densities from state
!> @details raw values is rectified !> @details raw values is rectified
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
pure function getRho0(ph,en) result(rho0) pure function getRho0(ph,en) result(rho_0)
integer, intent(in) :: ph, en integer, intent(in) :: ph, en
real(pREAL), dimension(param(ph)%sum_N_sl,10) :: rho0 real(pREAL), dimension(param(ph)%sum_N_sl,10) :: rho_0
associate(prm => param(ph)) associate(prm => param(ph))
rho0 = reshape(state0(ph)%rho(:,en),[prm%sum_N_sl,10]) rho_0 = reshape(state0(ph)%rho(:,en),[prm%sum_N_sl,10])
! ensure positive densities (not for imm, they have a sign) ! ensure positive densities (not for imm, they have a sign)
rho0(:,mob) = max(rho0(:,mob),0.0_pREAL) rho_0(:,mob) = max(rho_0(:,mob),0.0_pREAL)
rho0(:,dip) = max(rho0(:,dip),0.0_pREAL) rho_0(:,dip) = max(rho_0(:,dip),0.0_pREAL)
where (abs(rho0) < max(prm%rho_min/geom(ph)%V_0(en)**(2.0_pREAL/3.0_pREAL),prm%rho_significant)) & where (abs(rho_0) < max(prm%rho_min/geom(ph)%v_0(en)**(2.0_pREAL/3.0_pREAL),prm%rho_significant)) &
rho0 = 0.0_pREAL rho_0 = 0.0_pREAL
end associate end associate
end function getRho0 end function getRho0
!--------------------------------------------------------------------------------------------------
!--------------------------------------------------------------------------------------------------
subroutine storeGeometry(ph) subroutine storeGeometry(ph)
integer, intent(in) :: ph integer, intent(in) :: ph
integer :: ce, co, nCell integer :: ce, co, nCell
real(pREAL), dimension(:), allocatable :: V real(pREAL), dimension(:), allocatable :: v_0
real(pREAL), dimension(:,:), allocatable :: a_0, x_0
real(pREAL), dimension(:,:,:), allocatable :: n_0
integer, dimension(:,:,:), allocatable :: neighborhood integer, dimension(:,:,:), allocatable :: neighborhood
real(pREAL), dimension(:,:), allocatable :: area, coords
real(pREAL), dimension(:,:,:), allocatable :: areaNormal
nCell = product(shape(IPvolume))
V = reshape(IPvolume,[nCell]) nCell = product(shape(IPVolume0))
neighborhood = reshape(IPneighborhood,[3,nIPneighbors,nCell])
area = reshape(IParea,[nIPneighbors,nCell]) v_0 = reshape(IPVolume0,[nCell])
areaNormal = reshape(IPareaNormal,[3,nIPneighbors,nCell]) a_0 = reshape(IPArea0,[nCellNeighbors,nCell])
coords = reshape(discretization_IPcoords,[3,nCell]) x_0 = reshape(discretization_IPcoords,[3,nCell])
n_0 = reshape(IPAreaNormal0,[3,nCellNeighbors,nCell])
neighborhood = reshape(IPneighborhood,[3,nCellNeighbors,nCell])
do ce = 1, size(material_entry_homogenization,1) do ce = 1, size(material_entry_homogenization,1)
do co = 1, homogenization_maxNconstituents do co = 1, homogenization_maxNconstituents
if (material_ID_phase(co,ce) == ph) then if (material_ID_phase(co,ce) == ph) then
geom(ph)%V_0(material_entry_phase(co,ce)) = V(ce) geom(ph)%v_0(material_entry_phase(co,ce)) = v_0(ce)
geom(ph)%a_0(:,material_entry_phase(co,ce)) = a_0(:,ce)
geom(ph)%x_0(:,material_entry_phase(co,ce)) = x_0(:,ce)
geom(ph)%n_0(:,:,material_entry_phase(co,ce)) = n_0(:,:,ce)
geom(ph)%IPneighborhood(:,:,material_entry_phase(co,ce)) = neighborhood(:,:,ce) geom(ph)%IPneighborhood(:,:,material_entry_phase(co,ce)) = neighborhood(:,:,ce)
geom(ph)%IParea(:,material_entry_phase(co,ce)) = area(:,ce)
geom(ph)%IPareaNormal(:,:,material_entry_phase(co,ce)) = areaNormal(:,:,ce)
geom(ph)%IPcoordinates(:,material_entry_phase(co,ce)) = coords(:,ce)
end if end if
end do end do
end do end do