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

store IP neighborhood for ph/en access

This commit is contained in:
Martin Diehl 2022-02-04 07:33:12 +01:00
parent 8b12814eec
commit 900ef0a2c9
1 changed files with 40 additions and 23 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

63
src/phase_mechanical_plastic_nonlocal.f90
View File

@ -15,6 +15,9 @@ submodule(phase:plastic) nonlocal
type :: tGeometry type :: tGeometry
real(pReal), dimension(:), allocatable :: V_0 real(pReal), dimension(:), allocatable :: V_0
integer, dimension(:,:,:), allocatable :: IPneighborhood
real(pReal), dimension(:,:), allocatable :: IParea
real(pReal), dimension(:,:,:), allocatable :: IPareaNormal
end type tGeometry end type tGeometry
type(tGeometry), dimension(:), allocatable :: geom type(tGeometry), dimension(:), allocatable :: geom
@ -407,6 +410,9 @@ 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)%IPareaNormal(3,nIPneighbors,Nmembers))
allocate(geom(ph)%IParea(nIPneighbors,Nmembers))
call storeGeometry(ph) call storeGeometry(ph)
if(plasticState(ph)%nonlocal .and. .not. allocated(IPneighborhood)) & if(plasticState(ph)%nonlocal .and. .not. allocated(IPneighborhood)) &
@ -652,8 +658,8 @@ module subroutine nonlocal_dependentState(ph, en, ip, el)
nRealNeighbors = 0.0_pReal nRealNeighbors = 0.0_pReal
neighbor_rhoTotal = 0.0_pReal neighbor_rhoTotal = 0.0_pReal
do n = 1,nIPneighbors do n = 1,nIPneighbors
neighbor_el = IPneighborhood(1,n,ip,el) neighbor_el = geom(ph)%IPneighborhood(1,n,en)
neighbor_ip = IPneighborhood(2,n,ip,el) neighbor_ip = geom(ph)%IPneighborhood(2,n,en)
no = material_phasememberAt(1,neighbor_ip,neighbor_el) no = material_phasememberAt(1,neighbor_ip,neighbor_el)
if (neighbor_el > 0 .and. neighbor_ip > 0) then if (neighbor_el > 0 .and. neighbor_ip > 0) then
if (material_phaseAt(1,neighbor_el) == ph) then if (material_phaseAt(1,neighbor_el) == ph) then
@ -669,9 +675,9 @@ module subroutine nonlocal_dependentState(ph, en, ip, el)
connection_latticeConf(1:3,n) = matmul(invFe, discretization_IPcoords(1:3,neighbor_el+neighbor_ip-1) & connection_latticeConf(1:3,n) = matmul(invFe, discretization_IPcoords(1:3,neighbor_el+neighbor_ip-1) &
- discretization_IPcoords(1:3,el+neighbor_ip-1)) - discretization_IPcoords(1:3,el+neighbor_ip-1))
normal_latticeConf = matmul(transpose(invFp), IPareaNormal(1:3,n,ip,el)) normal_latticeConf = matmul(transpose(invFp), geom(ph)%IPareaNormal(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)/IParea(n,ip,el) ! 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)%IParea(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
@ -1101,7 +1107,7 @@ module subroutine nonlocal_dotState(Mp,timestep, &
- rhoDip(s,1) / timestep - rhoDotAthermalAnnihilation(s,9) & - rhoDip(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,ip,el) & rhoDot = rhoDotFlux(timestep, ph,en) &
+ rhoDotMultiplication & + rhoDotMultiplication &
+ rhoDotSingle2DipoleGlide & + rhoDotSingle2DipoleGlide &
+ rhoDotAthermalAnnihilation & + rhoDotAthermalAnnihilation &
@ -1131,17 +1137,15 @@ end subroutine nonlocal_dotState
!> @brief calculates the rate of change of microstructure !> @brief calculates the rate of change of microstructure
!--------------------------------------------------------------------------------------------------- !---------------------------------------------------------------------------------------------------
#if __INTEL_COMPILER >= 2020 #if __INTEL_COMPILER >= 2020
non_recursive function rhoDotFlux(timestep,ph,en,ip,el) non_recursive function rhoDotFlux(timestep,ph,en)
#else #else
function rhoDotFlux(timestep,ph,en,ip,el) function rhoDotFlux(timestep,ph,en)
#endif #endif
real(pReal), intent(in) :: & real(pReal), intent(in) :: &
timestep !< substepped crystallite time increment timestep !< substepped crystallite time increment
integer, intent(in) :: & integer, intent(in) :: &
ph, & ph, &
en, & en
ip, & !< current integration point
el !< current element number
integer :: & integer :: &
neighbor_ph, & !< phase of my neighbor's plasticity neighbor_ph, & !< phase of my neighbor's plasticity
@ -1217,14 +1221,14 @@ function rhoDotFlux(timestep,ph,en,ip,el)
!*** 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(v0) * timestep &
> geom(ph)%V_0(en)/ maxval(IParea(:,ip,el)))) then ! ...with velocity above critical value (we use the reference volume and area for simplicity here) > 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)
#ifdef DEBUG #ifdef DEBUG
if (debugConstitutive%extensive) then if (debugConstitutive%extensive) then
print'(a,i5,a,i2)', '<< CONST >> CFL condition not fullfilled at ph ',ph,' en ',en print'(a,i5,a,i2)', '<< CONST >> CFL condition not fullfilled at ph ',ph,' en ',en
print'(a,e10.3,a,e10.3)', '<< CONST >> velocity is at ', & print'(a,e10.3,a,e10.3)', '<< CONST >> velocity is at ', &
maxval(abs(v0), abs(dot_gamma) > 0.0_pReal & maxval(abs(v0), abs(dot_gamma) > 0.0_pReal &
.and. prm%C_CFL * abs(v0) * timestep & .and. prm%C_CFL * abs(v0) * timestep &
> geom(ph)%V_0(en) / maxval(IParea(:,ip,el))), & > geom(ph)%V_0(en) / maxval(geom(ph)%IParea(:,en))), &
' at a timestep of ',timestep ' at a timestep of ',timestep
print*, '<< CONST >> enforcing cutback !!!' print*, '<< CONST >> enforcing cutback !!!'
end if end if
@ -1247,16 +1251,16 @@ function rhoDotFlux(timestep,ph,en,ip,el)
neighbors: do n = 1,nIPneighbors neighbors: do n = 1,nIPneighbors
neighbor_el = IPneighborhood(1,n,ip,el) neighbor_el = geom(ph)%IPneighborhood(1,n,en)
neighbor_ip = IPneighborhood(2,n,ip,el) neighbor_ip = geom(ph)%IPneighborhood(2,n,en)
neighbor_n = IPneighborhood(3,n,ip,el) neighbor_n = geom(ph)%IPneighborhood(3,n,en)
np = material_phaseAt(1,neighbor_el) np = material_phaseAt(1,neighbor_el)
no = material_phasememberAt(1,neighbor_ip,neighbor_el) no = material_phasememberAt(1,neighbor_ip,neighbor_el)
opposite_neighbor = n + mod(n,2) - mod(n+1,2) opposite_neighbor = n + mod(n,2) - mod(n+1,2)
opposite_el = IPneighborhood(1,opposite_neighbor,ip,el) opposite_el = geom(ph)%IPneighborhood(1,opposite_neighbor,en)
opposite_ip = IPneighborhood(2,opposite_neighbor,ip,el) opposite_ip = geom(ph)%IPneighborhood(2,opposite_neighbor,en)
opposite_n = IPneighborhood(3,opposite_neighbor,ip,el) opposite_n = geom(ph)%IPneighborhood(3,opposite_neighbor,en)
if (neighbor_n > 0) then ! if neighbor exists, average deformation gradient if (neighbor_n > 0) then ! if neighbor exists, average deformation gradient
neighbor_ph = material_phaseAt(1,neighbor_el) neighbor_ph = material_phaseAt(1,neighbor_el)
@ -1327,10 +1331,10 @@ function rhoDotFlux(timestep,ph,en,ip,el)
if (phase_plasticity(material_phaseAt(1,opposite_el)) == PLASTIC_NONLOCAL_ID) then if (phase_plasticity(material_phaseAt(1,opposite_el)) == PLASTIC_NONLOCAL_ID) then
normal_me2neighbor_defConf = math_det33(Favg) & normal_me2neighbor_defConf = math_det33(Favg) &
* matmul(math_inv33(transpose(Favg)),IPareaNormal(1:3,n,ip,el)) ! normal of the interface in (average) deformed configuration (pointing en => neighbor) * matmul(math_inv33(transpose(Favg)),geom(ph)%IPareaNormal(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 = IParea(n,ip,el) * norm2(normal_me2neighbor) area = geom(ph)%IParea(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
@ -1758,14 +1762,27 @@ subroutine storeGeometry(ph)
integer, intent(in) :: ph integer, intent(in) :: ph
integer :: ce, co integer :: ce, co, nCell
real(pReal), dimension(:), allocatable :: V real(pReal), dimension(:), allocatable :: V
integer, dimension(:,:,:), allocatable :: neighborhood
real(pReal), dimension(:,:), allocatable :: area
real(pReal), dimension(:,:,:), allocatable :: areaNormal
nCell = product(shape(IPvolume))
V = reshape(IPvolume,[nCell])
neighborhood = reshape(IPneighborhood,[3,nIPneighbors,nCell])
area = reshape(IParea,[nIPneighbors,nCell])
areaNormal = reshape(IPareaNormal,[3,nIPneighbors,nCell])
V = reshape(IPvolume,[product(shape(IPvolume))])
do ce = 1, size(material_homogenizationEntry,1) do ce = 1, size(material_homogenizationEntry,1)
do co = 1, homogenization_maxNconstituents do co = 1, homogenization_maxNconstituents
if (material_phaseID(co,ce) == ph) geom(ph)%V_0(material_phaseEntry(co,ce)) = V(ce) if (material_phaseID(co,ce) == ph) then
geom(ph)%V_0(material_phaseEntry(co,ce)) = V(ce)
geom(ph)%IPneighborhood(:,:,material_phaseEntry(co,ce)) = neighborhood(:,:,ce)
geom(ph)%IParea(:,material_phaseEntry(co,ce)) = area(:,ce)
geom(ph)%IPareaNormal(:,:,material_phaseEntry(co,ce)) = areaNormal(:,:,ce)
end if
end do end do
end do end do