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

clean interface 

still need to get rid of internal converstion to instance and el,ip
arguments
This commit is contained in:
Martin Diehl 2021-02-14 15:29:10 +01:00
parent 5a1ca012f8
commit 4026881e5a
4 changed files with 107 additions and 112 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

6
src/phase.f90
View File

@ -247,9 +247,9 @@ module phase
TDot TDot
end subroutine phase_thermal_getRate end subroutine phase_thermal_getRate
module subroutine plastic_nonlocal_updateCompatibility(orientation,instance,i,e) module subroutine plastic_nonlocal_updateCompatibility(orientation,ph,i,e)
integer, intent(in) :: & integer, intent(in) :: &
instance, & ph, &
i, & i, &
e e
type(rotation), dimension(1,discretization_nIPs,discretization_Nelems), intent(in) :: & type(rotation), dimension(1,discretization_nIPs,discretization_Nelems), intent(in) :: &
@ -616,7 +616,7 @@ subroutine crystallite_orientations(co,ip,el)
if (plasticState(material_phaseAt(1,el))%nonlocal) & if (plasticState(material_phaseAt(1,el))%nonlocal) &
call plastic_nonlocal_updateCompatibility(crystallite_orientation, & call plastic_nonlocal_updateCompatibility(crystallite_orientation, &
phase_plasticInstance(material_phaseAt(1,el)),ip,el) material_phaseAt(1,el),ip,el)
end subroutine crystallite_orientations end subroutine crystallite_orientations

6
src/phase_mechanical.f90
View File

@ -155,8 +155,8 @@ submodule(phase) mechanical
character(len=*), intent(in) :: group character(len=*), intent(in) :: group
end subroutine plastic_dislotungsten_results end subroutine plastic_dislotungsten_results
module subroutine plastic_nonlocal_results(instance,group) module subroutine plastic_nonlocal_results(ph,group)
integer, intent(in) :: instance integer, intent(in) :: ph
character(len=*), intent(in) :: group character(len=*), intent(in) :: group
end subroutine plastic_nonlocal_results end subroutine plastic_nonlocal_results
@ -418,7 +418,7 @@ module subroutine mechanical_results(group,ph)
call plastic_dislotungsten_results(ph,group//'plastic/') call plastic_dislotungsten_results(ph,group//'plastic/')
case(PLASTICITY_NONLOCAL_ID) case(PLASTICITY_NONLOCAL_ID)
call plastic_nonlocal_results(phase_plasticInstance(ph),group//'plastic/') call plastic_nonlocal_results(ph,group//'plastic/')
end select end select

12
src/phase_mechanical_plastic.f90
View File

@ -193,9 +193,9 @@ submodule(phase:mechanical) plastic
me me
end subroutine dislotungsten_dependentState end subroutine dislotungsten_dependentState
module subroutine nonlocal_dependentState(instance, me, ip, el) module subroutine nonlocal_dependentState(ph, me, ip, el)
integer, intent(in) :: & integer, intent(in) :: &
instance, & ph, &
me, & me, &
ip, & !< current integration point ip, & !< current integration point
el !< current element number el !< current element number
@ -209,11 +209,11 @@ submodule(phase:mechanical) plastic
me me
end subroutine plastic_kinehardening_deltaState end subroutine plastic_kinehardening_deltaState
module subroutine plastic_nonlocal_deltaState(Mp,instance,me,ip,el) module subroutine plastic_nonlocal_deltaState(Mp,ph,me,ip,el)
real(pReal), dimension(3,3), intent(in) :: & real(pReal), dimension(3,3), intent(in) :: &
Mp Mp
integer, intent(in) :: & integer, intent(in) :: &
instance, & ph, &
me, & me, &
ip, & ip, &
el el
@ -380,7 +380,7 @@ module subroutine plastic_dependentState(co, ip, el)
call dislotungsten_dependentState(ph,me) call dislotungsten_dependentState(ph,me)
case (PLASTICITY_NONLOCAL_ID) plasticType case (PLASTICITY_NONLOCAL_ID) plasticType
call nonlocal_dependentState(instance,me,ip,el) call nonlocal_dependentState(ph,me,ip,el)
end select plasticType end select plasticType
@ -421,7 +421,7 @@ module function plastic_deltaState(co, ip, el, ph, me) result(broken)
broken = any(IEEE_is_NaN(plasticState(ph)%deltaState(:,me))) broken = any(IEEE_is_NaN(plasticState(ph)%deltaState(:,me)))
case (PLASTICITY_NONLOCAL_ID) plasticType case (PLASTICITY_NONLOCAL_ID) plasticType
call plastic_nonlocal_deltaState(Mp,instance,me,ip,el) call plastic_nonlocal_deltaState(Mp,ph,me,ip,el)
broken = any(IEEE_is_NaN(plasticState(ph)%deltaState(:,me))) broken = any(IEEE_is_NaN(plasticState(ph)%deltaState(:,me)))
case default case default

195
src/phase_mechanical_plastic_nonlocal.f90
View File

@ -12,6 +12,8 @@ submodule(phase:plastic) nonlocal
IParea => geometry_plastic_nonlocal_IParea0, & IParea => geometry_plastic_nonlocal_IParea0, &
IPareaNormal => geometry_plastic_nonlocal_IPareaNormal0, & IPareaNormal => geometry_plastic_nonlocal_IPareaNormal0, &
geometry_plastic_nonlocal_disable geometry_plastic_nonlocal_disable
use phase, &
ins => phase_plasticInstance
type :: tGeometry type :: tGeometry
real(pReal), dimension(:), allocatable :: V_0 real(pReal), dimension(:), allocatable :: V_0
@ -160,7 +162,7 @@ submodule(phase:plastic) nonlocal
state, & state, &
state0 state0
type(tParameters), dimension(:), allocatable :: param !< containers of constitutive parameters (len Ninstances) type(tParameters), dimension(:), allocatable :: param !< containers of constitutive parameters
type(tNonlocalMicrostructure), dimension(:), allocatable :: microstructure type(tNonlocalMicrostructure), dimension(:), allocatable :: microstructure
@ -510,7 +512,7 @@ module function plastic_nonlocal_init() result(myPlasticity)
allocate(dst%tau_back(prm%sum_N_sl,Nconstituents),source=0.0_pReal) allocate(dst%tau_back(prm%sum_N_sl,Nconstituents),source=0.0_pReal)
end associate end associate
if (Nconstituents > 0) call stateInit(ini,p,Nconstituents,i) if (Nconstituents > 0) call stateInit(ini,p,Nconstituents)
plasticState(p)%state0 = plasticState(p)%state plasticState(p)%state0 = plasticState(p)%state
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -565,10 +567,10 @@ end function plastic_nonlocal_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief calculates quantities characterizing the microstructure !> @brief calculates quantities characterizing the microstructure
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module subroutine nonlocal_dependentState(instance, me, ip, el) module subroutine nonlocal_dependentState(ph, me, ip, el)
integer, intent(in) :: & integer, intent(in) :: &
instance, & ph, &
me, & me, &
ip, & ip, &
el el
@ -602,29 +604,29 @@ module subroutine nonlocal_dependentState(instance, me, ip, el)
invConnections invConnections
real(pReal), dimension(3,nIPneighbors) :: & real(pReal), dimension(3,nIPneighbors) :: &
connection_latticeConf connection_latticeConf
real(pReal), dimension(2,param(instance)%sum_N_sl) :: & real(pReal), dimension(2,param(ins(ph))%sum_N_sl) :: &
rhoExcess rhoExcess
real(pReal), dimension(param(instance)%sum_N_sl) :: & real(pReal), dimension(param(ins(ph))%sum_N_sl) :: &
rho_edg_delta, & rho_edg_delta, &
rho_scr_delta rho_scr_delta
real(pReal), dimension(param(instance)%sum_N_sl,10) :: & real(pReal), dimension(param(ins(ph))%sum_N_sl,10) :: &
rho, & rho, &
rho0, & rho0, &
rho_neighbor0 rho_neighbor0
real(pReal), dimension(param(instance)%sum_N_sl,param(instance)%sum_N_sl) :: & real(pReal), dimension(param(ins(ph))%sum_N_sl,param(ins(ph))%sum_N_sl) :: &
myInteractionMatrix ! corrected slip interaction matrix myInteractionMatrix ! corrected slip interaction matrix
real(pReal), dimension(param(instance)%sum_N_sl,nIPneighbors) :: & real(pReal), dimension(param(ins(ph))%sum_N_sl,nIPneighbors) :: &
rho_edg_delta_neighbor, & rho_edg_delta_neighbor, &
rho_scr_delta_neighbor rho_scr_delta_neighbor
real(pReal), dimension(2,maxval(param%sum_N_sl),nIPneighbors) :: & real(pReal), dimension(2,maxval(param%sum_N_sl),nIPneighbors) :: &
neighbor_rhoExcess, & ! excess density at neighboring material point neighbor_rhoExcess, & ! excess density at neighboring material point
neighbor_rhoTotal ! total density at neighboring material point neighbor_rhoTotal ! total density at neighboring material point
real(pReal), dimension(3,param(instance)%sum_N_sl,2) :: & real(pReal), dimension(3,param(ins(ph))%sum_N_sl,2) :: &
m ! direction of dislocation motion m ! direction of dislocation motion
associate(prm => param(instance),dst => microstructure(instance), stt => state(instance)) associate(prm => param(ins(ph)),dst => microstructure(ins(ph)), stt => state(ins(ph)))
rho = getRho(instance,me,ip,el) rho = getRho(ph,me,ip,el)
stt%rho_forest(:,me) = matmul(prm%forestProjection_Edge, sum(abs(rho(:,edg)),2)) & stt%rho_forest(:,me) = matmul(prm%forestProjection_Edge, sum(abs(rho(:,edg)),2)) &
+ matmul(prm%forestProjection_Screw,sum(abs(rho(:,scr)),2)) + matmul(prm%forestProjection_Screw,sum(abs(rho(:,scr)),2))
@ -652,9 +654,8 @@ module subroutine nonlocal_dependentState(instance, me, ip, el)
! 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(instance,me,ip,el) rho0 = getRho0(ph,me,ip,el)
if (.not. phase_localPlasticity(material_phaseAt(1,el)) .and. prm%shortRangeStressCorrection) then if (.not. phase_localPlasticity(material_phaseAt(1,el)) .and. prm%shortRangeStressCorrection) then
ph = material_phaseAt(1,el)
invFp = math_inv33(phase_mechanical_Fp(ph)%data(1:3,1:3,me)) invFp = math_inv33(phase_mechanical_Fp(ph)%data(1:3,1:3,me))
invFe = math_inv33(phase_mechanical_Fe(ph)%data(1:3,1:3,me)) invFe = math_inv33(phase_mechanical_Fe(ph)%data(1:3,1:3,me))
@ -675,11 +676,11 @@ module subroutine nonlocal_dependentState(instance, me, ip, el)
neighbor_ip = IPneighborhood(2,n,ip,el) neighbor_ip = IPneighborhood(2,n,ip,el)
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
neighbor_instance = phase_plasticInstance(material_phaseAt(1,neighbor_el)) neighbor_instance = ins(material_phaseAt(1,neighbor_el))
if (neighbor_instance == instance) then if (neighbor_instance == ins(ph)) then
nRealNeighbors = nRealNeighbors + 1.0_pReal nRealNeighbors = nRealNeighbors + 1.0_pReal
rho_neighbor0 = getRho0(instance,no,neighbor_ip,neighbor_el) rho_neighbor0 = getRho0(ph,no,neighbor_ip,neighbor_el)
rho_edg_delta_neighbor(:,n) = rho_neighbor0(:,mob_edg_pos) - rho_neighbor0(:,mob_edg_neg) rho_edg_delta_neighbor(:,n) = rho_neighbor0(:,mob_edg_pos) - rho_neighbor0(:,mob_edg_neg)
rho_scr_delta_neighbor(:,n) = rho_neighbor0(:,mob_scr_pos) - rho_neighbor0(:,mob_scr_neg) rho_scr_delta_neighbor(:,n) = rho_neighbor0(:,mob_scr_pos) - rho_neighbor0(:,mob_scr_neg)
@ -795,25 +796,25 @@ module subroutine nonlocal_LpAndItsTangent(Lp,dLp_dMp, &
l, & l, &
t, & !< dislocation type t, & !< dislocation type
s !< index of my current slip system s !< index of my current slip system
real(pReal), dimension(param(phase_plasticInstance(ph))%sum_N_sl,8) :: & real(pReal), dimension(param(ins(ph))%sum_N_sl,8) :: &
rhoSgl !< single dislocation densities (including blocked) rhoSgl !< single dislocation densities (including blocked)
real(pReal), dimension(param(phase_plasticInstance(ph))%sum_N_sl,10) :: & real(pReal), dimension(param(ins(ph))%sum_N_sl,10) :: &
rho rho
real(pReal), dimension(param(phase_plasticInstance(ph))%sum_N_sl,4) :: & real(pReal), dimension(param(ins(ph))%sum_N_sl,4) :: &
v, & !< velocity v, & !< velocity
tauNS, & !< resolved shear stress including non Schmid and backstress terms tauNS, & !< resolved shear stress including non Schmid and backstress terms
dv_dtau, & !< velocity derivative with respect to the shear stress dv_dtau, & !< velocity derivative with respect to the shear stress
dv_dtauNS !< velocity derivative with respect to the shear stress dv_dtauNS !< velocity derivative with respect to the shear stress
real(pReal), dimension(param(phase_plasticInstance(ph))%sum_N_sl) :: & real(pReal), dimension(param(ins(ph))%sum_N_sl) :: &
tau, & !< resolved shear stress including backstress terms tau, & !< resolved shear stress including backstress terms
gdotTotal !< shear rate gdotTotal !< shear rate
associate(prm => param(phase_plasticInstance(ph)),dst=>microstructure(phase_plasticInstance(ph)),& associate(prm => param(ins(ph)),dst=>microstructure(ins(ph)),&
stt=>state(phase_plasticInstance(ph))) stt=>state(ins(ph)))
ns = prm%sum_N_sl ns = prm%sum_N_sl
!*** shortcut to state variables !*** shortcut to state variables
rho = getRho(phase_plasticInstance(ph),me,ip,el) rho = getRho(ph,me,ip,el)
rhoSgl = rho(:,sgl) rhoSgl = rho(:,sgl)
do s = 1,ns do s = 1,ns
@ -833,7 +834,7 @@ module subroutine nonlocal_LpAndItsTangent(Lp,dLp_dMp, &
! edges ! edges
call kinetics(v(:,1), dv_dtau(:,1), dv_dtauNS(:,1), & call kinetics(v(:,1), dv_dtau(:,1), dv_dtauNS(:,1), &
tau, tauNS(:,1), dst%tau_pass(:,me),1,Temperature, phase_plasticInstance(ph)) tau, tauNS(:,1), dst%tau_pass(:,me),1,Temperature, ph)
v(:,2) = v(:,1) v(:,2) = v(:,1)
dv_dtau(:,2) = dv_dtau(:,1) dv_dtau(:,2) = dv_dtau(:,1)
dv_dtauNS(:,2) = dv_dtauNS(:,1) dv_dtauNS(:,2) = dv_dtauNS(:,1)
@ -846,7 +847,7 @@ module subroutine nonlocal_LpAndItsTangent(Lp,dLp_dMp, &
else else
do t = 3,4 do t = 3,4
call kinetics(v(:,t), dv_dtau(:,t), dv_dtauNS(:,t), & call kinetics(v(:,t), dv_dtau(:,t), dv_dtauNS(:,t), &
tau, tauNS(:,t), dst%tau_pass(:,me),2,Temperature, phase_plasticInstance(ph)) tau, tauNS(:,t), dst%tau_pass(:,me),2,Temperature, ph)
enddo enddo
endif endif
@ -879,12 +880,12 @@ end subroutine nonlocal_LpAndItsTangent
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief (instantaneous) incremental change of microstructure !> @brief (instantaneous) incremental change of microstructure
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module subroutine plastic_nonlocal_deltaState(Mp,instance,me,ip,el) module subroutine plastic_nonlocal_deltaState(Mp,ph,me,ip,el)
real(pReal), dimension(3,3), intent(in) :: & real(pReal), dimension(3,3), intent(in) :: &
Mp !< MandelStress Mp !< MandelStress
integer, intent(in) :: & integer, intent(in) :: &
instance, & ! current instance of this plasticity ph, &
me, & !< offset me, & !< offset
ip, & ip, &
el el
@ -895,31 +896,29 @@ module subroutine plastic_nonlocal_deltaState(Mp,instance,me,ip,el)
c, & ! character of dislocation c, & ! character of dislocation
t, & ! type of dislocation t, & ! type of dislocation
s ! index of my current slip system s ! index of my current slip system
real(pReal), dimension(param(instance)%sum_N_sl,10) :: & real(pReal), dimension(param(ins(ph))%sum_N_sl,10) :: &
deltaRhoRemobilization, & ! density increment by remobilization deltaRhoRemobilization, & ! density increment by remobilization
deltaRhoDipole2SingleStress ! density increment by dipole dissociation (by stress change) deltaRhoDipole2SingleStress ! density increment by dipole dissociation (by stress change)
real(pReal), dimension(param(instance)%sum_N_sl,10) :: & real(pReal), dimension(param(ins(ph))%sum_N_sl,10) :: &
rho ! current dislocation densities rho ! current dislocation densities
real(pReal), dimension(param(instance)%sum_N_sl,4) :: & real(pReal), dimension(param(ins(ph))%sum_N_sl,4) :: &
v ! dislocation glide velocity v ! dislocation glide velocity
real(pReal), dimension(param(instance)%sum_N_sl) :: & real(pReal), dimension(param(ins(ph))%sum_N_sl) :: &
tau ! current resolved shear stress tau ! current resolved shear stress
real(pReal), dimension(param(instance)%sum_N_sl,2) :: & real(pReal), dimension(param(ins(ph))%sum_N_sl,2) :: &
rhoDip, & ! current dipole dislocation densities (screw and edge dipoles) rhoDip, & ! 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
ph = material_phaseAt(1,el) associate(prm => param(ins(ph)),dst => microstructure(ins(ph)),del => deltaState(ins(ph)))
ns = prm%sum_N_sl
associate(prm => param(instance),dst => microstructure(instance),del => deltaState(instance))
ns = prm%sum_N_sl
!*** shortcut to state variables !*** shortcut to state variables
forall (s = 1:ns, t = 1:4) v(s,t) = plasticState(ph)%state(iV(s,t,instance),me) forall (s = 1:ns, t = 1:4) v(s,t) = plasticState(ph)%state(iV(s,t,ins(ph)),me)
forall (s = 1:ns, c = 1:2) dUpperOld(s,c) = plasticState(ph)%state(iD(s,c,instance),me) forall (s = 1:ns, c = 1:2) dUpperOld(s,c) = plasticState(ph)%state(iD(s,c,ins(ph)),me)
rho = getRho(instance,me,ip,el) rho = getRho(ph,me,ip,el)
rhoDip = rho(:,dip) rhoDip = rho(:,dip)
!**************************************************************************** !****************************************************************************
@ -964,7 +963,7 @@ module subroutine plastic_nonlocal_deltaState(Mp,instance,me,ip,el)
/ (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)
forall (s = 1:ns, c = 1:2) plasticState(ph)%state(iD(s,c,instance),me) = dUpper(s,c) forall (s = 1:ns, c = 1:2) plasticState(ph)%state(iD(s,c,ins(ph)),me) = dUpper(s,c)
plasticState(ph)%deltaState(:,me) = 0.0_pReal plasticState(ph)%deltaState(:,me) = 0.0_pReal
del%rho(:,me) = reshape(deltaRhoRemobilization + deltaRhoDipole2SingleStress, [10*ns]) del%rho(:,me) = reshape(deltaRhoRemobilization + deltaRhoDipole2SingleStress, [10*ns])
@ -1005,7 +1004,7 @@ module subroutine nonlocal_dotState(Mp, Temperature,timestep, &
c, & !< character of dislocation c, & !< character of dislocation
t, & !< type of dislocation t, & !< type of dislocation
s !< index of my current slip system s !< index of my current slip system
real(pReal), dimension(param(phase_plasticInstance(ph))%sum_N_sl,10) :: & real(pReal), dimension(param(ins(ph))%sum_N_sl,10) :: &
rho, & rho, &
rho0, & !< dislocation density at beginning of time step rho0, & !< dislocation density at beginning of time step
rhoDot, & !< density evolution rhoDot, & !< density evolution
@ -1013,17 +1012,17 @@ module subroutine nonlocal_dotState(Mp, Temperature,timestep, &
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(phase_plasticInstance(ph))%sum_N_sl,8) :: & real(pReal), dimension(param(ins(ph))%sum_N_sl,8) :: &
rhoSgl, & !< current single dislocation densities (positive/negative screw and edge without dipoles) rhoSgl, & !< 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) my_rhoSgl0 !< single dislocation densities of central ip (positive/negative screw and edge without dipoles)
real(pReal), dimension(param(phase_plasticInstance(ph))%sum_N_sl,4) :: & real(pReal), dimension(param(ins(ph))%sum_N_sl,4) :: &
v, & !< current dislocation glide velocity v, & !< current dislocation glide velocity
v0, & v0, &
gdot !< shear rates gdot !< shear rates
real(pReal), dimension(param(phase_plasticInstance(ph))%sum_N_sl) :: & real(pReal), dimension(param(ins(ph))%sum_N_sl) :: &
tau, & !< current resolved shear stress tau, & !< current resolved shear stress
vClimb !< climb velocity of edge dipoles vClimb !< climb velocity of edge dipoles
real(pReal), dimension(param(phase_plasticInstance(ph))%sum_N_sl,2) :: & real(pReal), dimension(param(ins(ph))%sum_N_sl,2) :: &
rhoDip, & !< current dipole dislocation densities (screw and edge dipoles) rhoDip, & !< 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
@ -1035,22 +1034,22 @@ module subroutine nonlocal_dotState(Mp, Temperature,timestep, &
return return
endif endif
associate(prm => param(phase_plasticInstance(ph)), & associate(prm => param(ins(ph)), &
dst => microstructure(phase_plasticInstance(ph)), & dst => microstructure(ins(ph)), &
dot => dotState(phase_plasticInstance(ph)), & dot => dotState(ins(ph)), &
stt => state(phase_plasticInstance(ph))) stt => state(ins(ph)))
ns = prm%sum_N_sl ns = prm%sum_N_sl
tau = 0.0_pReal tau = 0.0_pReal
gdot = 0.0_pReal gdot = 0.0_pReal
rho = getRho(phase_plasticInstance(ph),me,ip,el) rho = getRho(ph,me,ip,el)
rhoSgl = rho(:,sgl) rhoSgl = rho(:,sgl)
rhoDip = rho(:,dip) rhoDip = rho(:,dip)
rho0 = getRho0(phase_plasticInstance(ph),me,ip,el) rho0 = getRho0(ph,me,ip,el)
my_rhoSgl0 = rho0(:,sgl) my_rhoSgl0 = rho0(:,sgl)
forall (s = 1:ns, t = 1:4) v(s,t) = plasticState(ph)%state(iV(s,t,phase_plasticInstance(ph)),me) forall (s = 1:ns, t = 1:4) v(s,t) = plasticState(ph)%state(iV(s,t,ins(ph)),me)
gdot = rhoSgl(:,1:4) * v * spread(prm%b_sl,2,4) gdot = rhoSgl(:,1:4) * v * spread(prm%b_sl,2,4)
#ifdef DEBUG #ifdef DEBUG
@ -1099,7 +1098,7 @@ module subroutine nonlocal_dotState(Mp, Temperature,timestep, &
* sqrt(stt%rho_forest(:,me)) / prm%i_sl / prm%b_sl, 2, 4) * sqrt(stt%rho_forest(:,me)) / prm%i_sl / prm%b_sl, 2, 4)
endif isBCC endif isBCC
forall (s = 1:ns, t = 1:4) v0(s,t) = plasticState(ph)%state0(iV(s,t,phase_plasticInstance(ph)),me) forall (s = 1:ns, t = 1:4) v0(s,t) = plasticState(ph)%state0(iV(s,t,ins(ph)),me)
!**************************************************************************** !****************************************************************************
@ -1155,7 +1154,7 @@ module subroutine nonlocal_dotState(Mp, Temperature,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, phase_plasticInstance(ph),me,ip,el) & rhoDot = rhoDotFlux(timestep, ph,me,ip,el) &
+ rhoDotMultiplication & + rhoDotMultiplication &
+ rhoDotSingle2DipoleGlide & + rhoDotSingle2DipoleGlide &
+ rhoDotAthermalAnnihilation & + rhoDotAthermalAnnihilation &
@ -1184,12 +1183,12 @@ end subroutine nonlocal_dotState
!--------------------------------------------------------------------------------------------------- !---------------------------------------------------------------------------------------------------
!> @brief calculates the rate of change of microstructure !> @brief calculates the rate of change of microstructure
!--------------------------------------------------------------------------------------------------- !---------------------------------------------------------------------------------------------------
function rhoDotFlux(timestep,instance,me,ip,el) function rhoDotFlux(timestep,ph,me,ip,el)
real(pReal), intent(in) :: & real(pReal), intent(in) :: &
timestep !< substepped crystallite time increment timestep !< substepped crystallite time increment
integer, intent(in) :: & integer, intent(in) :: &
instance, & ph, &
me, & me, &
ip, & !< current integration point ip, & !< current integration point
el !< current element number el !< current element number
@ -1212,20 +1211,20 @@ function rhoDotFlux(timestep,instance,me,ip,el)
np,& !< neighbor phase shortcut np,& !< neighbor phase shortcut
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(instance)%sum_N_sl,10) :: & real(pReal), dimension(param(ins(ph))%sum_N_sl,10) :: &
rho, & rho, &
rho0, & !< dislocation density at beginning of time step rho0, & !< dislocation density at beginning of time step
rhoDotFlux !< density evolution by flux rhoDotFlux !< density evolution by flux
real(pReal), dimension(param(instance)%sum_N_sl,8) :: & real(pReal), dimension(param(ins(ph))%sum_N_sl,8) :: &
rhoSgl, & !< current single dislocation densities (positive/negative screw and edge without dipoles) rhoSgl, & !< current single dislocation densities (positive/negative screw and edge without dipoles)
neighbor_rhoSgl0, & !< 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)
my_rhoSgl0 !< 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(instance)%sum_N_sl,4) :: & real(pReal), dimension(param(ins(ph))%sum_N_sl,4) :: &
v, & !< current dislocation glide velocity v, & !< current dislocation glide velocity
v0, & v0, &
neighbor_v0, & !< dislocation glide velocity of enighboring ip neighbor_v0, & !< dislocation glide velocity of enighboring ip
gdot !< shear rates gdot !< shear rates
real(pReal), dimension(3,param(instance)%sum_N_sl,4) :: & real(pReal), dimension(3,param(ins(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
@ -1243,26 +1242,25 @@ function rhoDotFlux(timestep,instance,me,ip,el)
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
ph = material_phaseAt(1,el)
associate(prm => param(instance), & associate(prm => param(ins(ph)), &
dst => microstructure(instance), & dst => microstructure(ins(ph)), &
dot => dotState(instance), & dot => dotState(ins(ph)), &
stt => state(instance)) stt => state(ins(ph)))
ns = prm%sum_N_sl ns = prm%sum_N_sl
gdot = 0.0_pReal gdot = 0.0_pReal
rho = getRho(instance,me,ip,el) rho = getRho(ph,me,ip,el)
rhoSgl = rho(:,sgl) rhoSgl = rho(:,sgl)
rho0 = getRho0(instance,me,ip,el) rho0 = getRho0(ph,me,ip,el)
my_rhoSgl0 = rho0(:,sgl) my_rhoSgl0 = rho0(:,sgl)
forall (s = 1:ns, t = 1:4) v(s,t) = plasticState(ph)%state(iV(s,t,instance),me) !ToDo: MD: I think we should use state0 here forall (s = 1:ns, t = 1:4) v(s,t) = plasticState(ph)%state(iV(s,t,ins(ph)),me) !ToDo: MD: I think we should use state0 here
gdot = rhoSgl(:,1:4) * v * spread(prm%b_sl,2,4) gdot = rhoSgl(:,1:4) * v * spread(prm%b_sl,2,4)
forall (s = 1:ns, t = 1:4) v0(s,t) = plasticState(ph)%state0(iV(s,t,instance),me) forall (s = 1:ns, t = 1:4) v0(s,t) = plasticState(ph)%state0(iV(s,t,ins(ph)),me)
!**************************************************************************** !****************************************************************************
!*** calculate dislocation fluxes (only for nonlocal plasticity) !*** calculate dislocation fluxes (only for nonlocal plasticity)
@ -1314,7 +1312,7 @@ function rhoDotFlux(timestep,instance,me,ip,el)
opposite_n = IPneighborhood(3,opposite_neighbor,ip,el) opposite_n = IPneighborhood(3,opposite_neighbor,ip,el)
if (neighbor_n > 0) then ! if neighbor exists, average deformation gradient if (neighbor_n > 0) then ! if neighbor exists, average deformation gradient
neighbor_instance = phase_plasticInstance(material_phaseAt(1,neighbor_el)) neighbor_instance = ins(material_phaseAt(1,neighbor_el))
neighbor_F = phase_mechanical_F(np)%data(1:3,1:3,no) neighbor_F = phase_mechanical_F(np)%data(1:3,1:3,no)
neighbor_Fe = matmul(neighbor_F, math_inv33(phase_mechanical_Fp(np)%data(1:3,1:3,no))) neighbor_Fe = matmul(neighbor_F, math_inv33(phase_mechanical_Fp(np)%data(1:3,1:3,no)))
Favg = 0.5_pReal * (my_F + neighbor_F) Favg = 0.5_pReal * (my_F + neighbor_F)
@ -1421,12 +1419,12 @@ end function rhoDotFlux
! plane normals and signed cosine of the angle between the slip directions. Only the largest values ! plane normals and signed cosine of the angle between the slip directions. Only the largest values
! that sum up to a total of 1 are considered, all others are set to zero. ! that sum up to a total of 1 are considered, all others are set to zero.
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module subroutine plastic_nonlocal_updateCompatibility(orientation,instance,i,e) module subroutine plastic_nonlocal_updateCompatibility(orientation,ph,i,e)
type(rotation), dimension(1,discretization_nIPs,discretization_Nelems), intent(in) :: & type(rotation), dimension(1,discretization_nIPs,discretization_Nelems), intent(in) :: &
orientation ! crystal orientation orientation ! crystal orientation
integer, intent(in) :: & integer, intent(in) :: &
instance, & ph, &
i, & i, &
e e
@ -1439,19 +1437,17 @@ module subroutine plastic_nonlocal_updateCompatibility(orientation,instance,i,e)
ns, & ! number of active slip systems ns, & ! number of active slip systems
s1, & ! slip system index (me) s1, & ! slip system index (me)
s2 ! slip system index (my neighbor) s2 ! slip system index (my neighbor)
real(pReal), dimension(2,param(instance)%sum_N_sl,param(instance)%sum_N_sl,nIPneighbors) :: & real(pReal), dimension(2,param(ins(ph))%sum_N_sl,param(ins(ph))%sum_N_sl,nIPneighbors) :: &
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, &
thresholdValue, & thresholdValue, &
nThresholdValues nThresholdValues
logical, dimension(param(instance)%sum_N_sl) :: & logical, dimension(param(ins(ph))%sum_N_sl) :: &
belowThreshold belowThreshold
type(rotation) :: mis type(rotation) :: mis
ph = material_phaseAt(1,e) associate(prm => param(ins(ph)))
associate(prm => param(instance))
ns = prm%sum_N_sl ns = prm%sum_N_sl
!*** start out fully compatible !*** start out fully compatible
@ -1537,14 +1533,14 @@ end subroutine plastic_nonlocal_updateCompatibility
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief writes results to HDF5 output file !> @brief writes results to HDF5 output file
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module subroutine plastic_nonlocal_results(instance,group) module subroutine plastic_nonlocal_results(ph,group)
integer, intent(in) :: instance integer, intent(in) :: ph
character(len=*),intent(in) :: group character(len=*),intent(in) :: group
integer :: o integer :: o
associate(prm => param(instance),dst => microstructure(instance),stt=>state(instance)) associate(prm => param(ins(ph)),dst => microstructure(ins(ph)),stt=>state(ins(ph)))
outputsLoop: do o = 1,size(prm%output) outputsLoop: do o = 1,size(prm%output)
select case(trim(prm%output(o))) select case(trim(prm%output(o)))
case('rho_u_ed_pos') case('rho_u_ed_pos')
@ -1608,14 +1604,13 @@ end subroutine plastic_nonlocal_results
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief populates the initial dislocation density !> @brief populates the initial dislocation density
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine stateInit(ini,phase,Nconstituents,instance) subroutine stateInit(ini,phase,Nconstituents)
type(tInitialParameters) :: & type(tInitialParameters) :: &
ini ini
integer,intent(in) :: & integer,intent(in) :: &
phase, & phase, &
Nconstituents, & Nconstituents
instance
integer :: & integer :: &
e, & e, &
i, & i, &
@ -1636,7 +1631,7 @@ subroutine stateInit(ini,phase,Nconstituents,instance)
volume volume
associate(stt => state(instance)) associate(stt => state(ins(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
do e = 1,discretization_Nelems do e = 1,discretization_Nelems
@ -1684,18 +1679,18 @@ end subroutine stateInit
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief calculates kinetics !> @brief calculates kinetics
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
pure subroutine kinetics(v, dv_dtau, dv_dtauNS, tau, tauNS, tauThreshold, c, Temperature, instance) pure subroutine kinetics(v, dv_dtau, dv_dtauNS, tau, tauNS, tauThreshold, c, Temperature, ph)
integer, intent(in) :: & integer, intent(in) :: &
c, & !< dislocation character (1:edge, 2:screw) c, & !< dislocation character (1:edge, 2:screw)
instance ph
real(pReal), intent(in) :: & real(pReal), intent(in) :: &
Temperature !< temperature Temperature !< temperature
real(pReal), dimension(param(instance)%sum_N_sl), intent(in) :: & real(pReal), dimension(param(ins(ph))%sum_N_sl), intent(in) :: &
tau, & !< resolved external shear stress (without non Schmid effects) tau, & !< resolved external shear stress (without non Schmid effects)
tauNS, & !< resolved external shear stress (including non Schmid effects) tauNS, & !< resolved external shear stress (including non Schmid effects)
tauThreshold !< threshold shear stress tauThreshold !< threshold shear stress
real(pReal), dimension(param(instance)%sum_N_sl), intent(out) :: & real(pReal), dimension(param(ins(ph))%sum_N_sl), intent(out) :: &
v, & !< velocity v, & !< velocity
dv_dtau, & !< velocity derivative with respect to resolved shear stress (without non Schmid contributions) dv_dtau, & !< velocity derivative with respect to resolved shear stress (without non Schmid contributions)
dv_dtauNS !< velocity derivative with respect to resolved shear stress (including non Schmid contributions) dv_dtauNS !< velocity derivative with respect to resolved shear stress (including non Schmid contributions)
@ -1726,7 +1721,7 @@ pure subroutine kinetics(v, dv_dtau, dv_dtauNS, tau, tauNS, tauThreshold, c, Tem
criticalStress_S, & !< maximum obstacle strength criticalStress_S, & !< maximum obstacle strength
mobility !< dislocation mobility mobility !< dislocation mobility
associate(prm => param(instance)) associate(prm => param(ins(ph)))
ns = prm%sum_N_sl ns = prm%sum_N_sl
v = 0.0_pReal v = 0.0_pReal
dv_dtau = 0.0_pReal dv_dtau = 0.0_pReal
@ -1799,14 +1794,14 @@ end subroutine kinetics
!> @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 getRho(instance,me,ip,el) pure function getRho(ph,me,ip,el)
integer, intent(in) :: instance, me,ip,el integer, intent(in) :: ph, me,ip,el
real(pReal), dimension(param(instance)%sum_N_sl,10) :: getRho real(pReal), dimension(param(ins(ph))%sum_N_sl,10) :: getRho
associate(prm => param(instance)) associate(prm => param(ins(ph)))
getRho = reshape(state(instance)%rho(:,me),[prm%sum_N_sl,10]) getRho = reshape(state(ins(ph))%rho(:,me),[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)
getRho(:,mob) = max(getRho(:,mob),0.0_pReal) getRho(:,mob) = max(getRho(:,mob),0.0_pReal)
@ -1824,14 +1819,14 @@ 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(instance,me,ip,el) pure function getRho0(ph,me,ip,el)
integer, intent(in) :: instance, me,ip,el integer, intent(in) :: ph, me,ip,el
real(pReal), dimension(param(instance)%sum_N_sl,10) :: getRho0 real(pReal), dimension(param(ins(ph))%sum_N_sl,10) :: getRho0
associate(prm => param(instance)) associate(prm => param(ins(ph)))
getRho0 = reshape(state0(instance)%rho(:,me),[prm%sum_N_sl,10]) getRho0 = reshape(state0(ins(ph))%rho(:,me),[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)
getRho0(:,mob) = max(getRho0(:,mob),0.0_pReal) getRho0(:,mob) = max(getRho0(:,mob),0.0_pReal)