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

simplified access pattern

This commit is contained in:
Martin Diehl 2021-02-14 17:29:23 +01:00
parent 4026881e5a
commit f46d212e47
3 changed files with 56 additions and 86 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

23
src/phase_mechanical.f90
View File

@ -85,11 +85,8 @@ submodule(phase) mechanical
logical :: broken logical :: broken
end function plastic_dotState end function plastic_dotState
module function plastic_deltaState(co, ip, el, ph, me) result(broken) module function plastic_deltaState(ph, me) result(broken)
integer, intent(in) :: & integer, intent(in) :: &
co, & !< component-ID of integration point
ip, & !< integration point
el, & !< element
ph, & ph, &
me me
logical :: & logical :: &
@ -114,11 +111,9 @@ submodule(phase) mechanical
module subroutine plastic_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, & module subroutine plastic_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, &
S, Fi, co, ip, el) S, Fi, ph,me)
integer, intent(in) :: & integer, intent(in) :: &
co, & !< component-ID of integration point ph,me
ip, & !< integration point
el !< element
real(pReal), intent(in), dimension(3,3) :: & real(pReal), intent(in), dimension(3,3) :: &
S, & !< 2nd Piola-Kirchhoff stress S, & !< 2nd Piola-Kirchhoff stress
Fi !< intermediate deformation gradient Fi !< intermediate deformation gradient
@ -538,7 +533,7 @@ function integrateStress(F,subFp0,subFi0,Delta_t,co,ip,el) result(broken)
Fe, Fi_new, co, ip, el) Fe, Fi_new, co, ip, el)
call plastic_LpAndItsTangents(Lp_constitutive, dLp_dS, dLp_dFi, & call plastic_LpAndItsTangents(Lp_constitutive, dLp_dS, dLp_dFi, &
S, Fi_new, co, ip, el) S, Fi_new, ph,me)
!* update current residuum and check for convergence of loop !* update current residuum and check for convergence of loop
atol_Lp = max(num%rtol_crystalliteStress * max(norm2(Lpguess),norm2(Lp_constitutive)), & ! absolute tolerance from largest acceptable relative error atol_Lp = max(num%rtol_crystalliteStress * max(norm2(Lpguess),norm2(Lp_constitutive)), & ! absolute tolerance from largest acceptable relative error
@ -702,7 +697,7 @@ function integrateStateFPI(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el) resul
plasticState(ph)%state(1:sizeDotState,me) = plasticState(ph)%state(1:sizeDotState,me) & plasticState(ph)%state(1:sizeDotState,me) = plasticState(ph)%state(1:sizeDotState,me) &
- r(1:sizeDotState) - r(1:sizeDotState)
if (converged(r(1:sizeDotState),plasticState(ph)%state(1:sizeDotState,me),plasticState(ph)%atol(1:sizeDotState))) then if (converged(r(1:sizeDotState),plasticState(ph)%state(1:sizeDotState,me),plasticState(ph)%atol(1:sizeDotState))) then
broken = plastic_deltaState(co,ip,el,ph,me) broken = plastic_deltaState(ph,me)
exit iteration exit iteration
endif endif
@ -765,7 +760,7 @@ function integrateStateEuler(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el) res
plasticState(ph)%state(1:sizeDotState,me) = subState0 & plasticState(ph)%state(1:sizeDotState,me) = subState0 &
+ plasticState(ph)%dotState(1:sizeDotState,me) * Delta_t + plasticState(ph)%dotState(1:sizeDotState,me) * Delta_t
broken = plastic_deltaState(co,ip,el,ph,me) broken = plastic_deltaState(ph,me)
if(broken) return if(broken) return
broken = integrateStress(F,subFp0,subFi0,Delta_t,co,ip,el) broken = integrateStress(F,subFp0,subFi0,Delta_t,co,ip,el)
@ -807,7 +802,7 @@ function integrateStateAdaptiveEuler(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip
plasticState(ph)%state(1:sizeDotState,me) = subState0 & plasticState(ph)%state(1:sizeDotState,me) = subState0 &
+ plasticState(ph)%dotstate(1:sizeDotState,me) * Delta_t + plasticState(ph)%dotstate(1:sizeDotState,me) * Delta_t
broken = plastic_deltaState(co,ip,el,ph,me) broken = plastic_deltaState(ph,me)
if(broken) return if(broken) return
broken = integrateStress(F,subFp0,subFi0,Delta_t,co,ip,el) broken = integrateStress(F,subFp0,subFi0,Delta_t,co,ip,el)
@ -956,7 +951,7 @@ function integrateStateRK(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el,A,B,C,D
if(broken) return if(broken) return
broken = plastic_deltaState(co,ip,el,ph,me) broken = plastic_deltaState(ph,me)
if(broken) return if(broken) return
broken = integrateStress(F,subFp0,subFi0,Delta_t,co,ip,el) broken = integrateStress(F,subFp0,subFi0,Delta_t,co,ip,el)
@ -1329,7 +1324,7 @@ module function phase_mechanical_dPdF(dt,co,ip,el) result(dPdF)
call plastic_LpAndItsTangents(devNull,dLpdS,dLpdFi, & call plastic_LpAndItsTangents(devNull,dLpdS,dLpdFi, &
phase_mechanical_S(ph)%data(1:3,1:3,me), & phase_mechanical_S(ph)%data(1:3,1:3,me), &
phase_mechanical_Fi(ph)%data(1:3,1:3,me),co,ip,el) phase_mechanical_Fi(ph)%data(1:3,1:3,me),ph,me)
dLpdS = math_mul3333xx3333(dLpdFi,dFidS) + dLpdS dLpdS = math_mul3333xx3333(dLpdFi,dFidS) + dLpdS
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------

34
src/phase_mechanical_plastic.f90
View File

@ -104,7 +104,7 @@ submodule(phase:mechanical) plastic
end subroutine dislotungsten_LpAndItsTangent end subroutine dislotungsten_LpAndItsTangent
module subroutine nonlocal_LpAndItsTangent(Lp,dLp_dMp, & module subroutine nonlocal_LpAndItsTangent(Lp,dLp_dMp, &
Mp,Temperature,ph,me,ip,el) Mp,Temperature,ph,me)
real(pReal), dimension(3,3), intent(out) :: & real(pReal), dimension(3,3), intent(out) :: &
Lp Lp
real(pReal), dimension(3,3,3,3), intent(out) :: & real(pReal), dimension(3,3,3,3), intent(out) :: &
@ -116,9 +116,7 @@ submodule(phase:mechanical) plastic
Temperature Temperature
integer, intent(in) :: & integer, intent(in) :: &
ph, & ph, &
me, & me
ip, & !< current integration point
el !< current element number
end subroutine nonlocal_LpAndItsTangent end subroutine nonlocal_LpAndItsTangent
@ -209,14 +207,12 @@ submodule(phase:mechanical) plastic
me me
end subroutine plastic_kinehardening_deltaState end subroutine plastic_kinehardening_deltaState
module subroutine plastic_nonlocal_deltaState(Mp,ph,me,ip,el) module subroutine plastic_nonlocal_deltaState(Mp,ph,me)
real(pReal), dimension(3,3), intent(in) :: & real(pReal), dimension(3,3), intent(in) :: &
Mp Mp
integer, intent(in) :: & integer, intent(in) :: &
ph, & ph, &
me, & me
ip, &
el
end subroutine plastic_nonlocal_deltaState end subroutine plastic_nonlocal_deltaState
end interface end interface
@ -244,11 +240,9 @@ end subroutine plastic_init
! Mp in, dLp_dMp out ! Mp in, dLp_dMp out
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module subroutine plastic_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, & module subroutine plastic_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, &
S, Fi, co, ip, el) S, Fi, ph,me)
integer, intent(in) :: & integer, intent(in) :: &
co, & !< component-ID of integration point ph,me
ip, & !< integration point
el !< element
real(pReal), intent(in), dimension(3,3) :: & real(pReal), intent(in), dimension(3,3) :: &
S, & !< 2nd Piola-Kirchhoff stress S, & !< 2nd Piola-Kirchhoff stress
Fi !< intermediate deformation gradient Fi !< intermediate deformation gradient
@ -263,14 +257,13 @@ module subroutine plastic_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, &
real(pReal), dimension(3,3) :: & real(pReal), dimension(3,3) :: &
Mp !< Mandel stress work conjugate with Lp Mp !< Mandel stress work conjugate with Lp
integer :: & integer :: &
i, j, me, ph i, j
Mp = matmul(matmul(transpose(Fi),Fi),S) Mp = matmul(matmul(transpose(Fi),Fi),S)
me = material_phasememberAt(co,ip,el)
ph = material_phaseAt(co,el)
plasticType: select case (phase_plasticity(material_phaseAt(co,el)))
plasticType: select case (phase_plasticity(ph))
case (PLASTICITY_NONE_ID) plasticType case (PLASTICITY_NONE_ID) plasticType
Lp = 0.0_pReal Lp = 0.0_pReal
@ -286,7 +279,7 @@ module subroutine plastic_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, &
call kinehardening_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,me) call kinehardening_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,me)
case (PLASTICITY_NONLOCAL_ID) plasticType case (PLASTICITY_NONLOCAL_ID) plasticType
call nonlocal_LpAndItsTangent(Lp,dLp_dMp,Mp, thermal_T(ph,me),ph,me,ip,el) call nonlocal_LpAndItsTangent(Lp,dLp_dMp,Mp, thermal_T(ph,me),ph,me)
case (PLASTICITY_DISLOTWIN_ID) plasticType case (PLASTICITY_DISLOTWIN_ID) plasticType
call dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp, thermal_T(ph,me),ph,me) call dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp, thermal_T(ph,me),ph,me)
@ -391,12 +384,9 @@ end subroutine plastic_dependentState
!> @brief for constitutive models having an instantaneous change of state !> @brief for constitutive models having an instantaneous change of state
!> will return false if delta state is not needed/supported by the constitutive model !> will return false if delta state is not needed/supported by the constitutive model
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module function plastic_deltaState(co, ip, el, ph, me) result(broken) module function plastic_deltaState(ph, me) result(broken)
integer, intent(in) :: & integer, intent(in) :: &
co, & !< component-ID of integration point
ip, & !< integration point
el, & !< element
ph, & ph, &
me me
logical :: & logical :: &
@ -421,7 +411,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,ph,me,ip,el) call plastic_nonlocal_deltaState(Mp,ph,me)
broken = any(IEEE_is_NaN(plasticState(ph)%deltaState(:,me))) broken = any(IEEE_is_NaN(plasticState(ph)%deltaState(:,me)))
case default case default

85
src/phase_mechanical_plastic_nonlocal.f90
View File

@ -512,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) if (Nconstituents > 0) call stateInit(ini,p,Nconstituents,i)
plasticState(p)%state0 = plasticState(p)%state plasticState(p)%state0 = plasticState(p)%state
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -576,7 +576,6 @@ module subroutine nonlocal_dependentState(ph, me, ip, el)
el el
integer :: & integer :: &
ph, &
no, & !< neighbor offset no, & !< neighbor offset
neighbor_el, & ! element number of neighboring material point neighbor_el, & ! element number of neighboring material point
neighbor_ip, & ! integration point of neighboring material point neighbor_ip, & ! integration point of neighboring material point
@ -626,7 +625,7 @@ module subroutine nonlocal_dependentState(ph, me, ip, el)
associate(prm => param(ins(ph)),dst => microstructure(ins(ph)), stt => state(ins(ph))) associate(prm => param(ins(ph)),dst => microstructure(ins(ph)), stt => state(ins(ph)))
rho = getRho(ph,me,ip,el) rho = getRho(ph,me)
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))
@ -654,7 +653,7 @@ module subroutine nonlocal_dependentState(ph, 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(ph,me,ip,el) rho0 = getRho0(ph,me)
if (.not. phase_localPlasticity(material_phaseAt(1,el)) .and. prm%shortRangeStressCorrection) then if (.not. phase_localPlasticity(material_phaseAt(1,el)) .and. prm%shortRangeStressCorrection) then
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))
@ -665,7 +664,7 @@ module subroutine nonlocal_dependentState(ph, me, ip, el)
rhoExcess(1,:) = rho_edg_delta rhoExcess(1,:) = rho_edg_delta
rhoExcess(2,:) = rho_scr_delta rhoExcess(2,:) = rho_scr_delta
FVsize = IPvolume(ip,el) ** (1.0_pReal/3.0_pReal) FVsize = geom(ph)%V_0(me) ** (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
@ -680,7 +679,7 @@ module subroutine nonlocal_dependentState(ph, me, ip, el)
if (neighbor_instance == ins(ph)) then if (neighbor_instance == ins(ph)) then
nRealNeighbors = nRealNeighbors + 1.0_pReal nRealNeighbors = nRealNeighbors + 1.0_pReal
rho_neighbor0 = getRho0(ph,no,neighbor_ip,neighbor_el) rho_neighbor0 = getRho0(ph,no)
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)
@ -772,16 +771,14 @@ end subroutine nonlocal_dependentState
!> @brief calculates plastic velocity gradient and its tangent !> @brief calculates plastic velocity gradient and its tangent
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module subroutine nonlocal_LpAndItsTangent(Lp,dLp_dMp, & module subroutine nonlocal_LpAndItsTangent(Lp,dLp_dMp, &
Mp,Temperature,ph,me,ip,el) Mp,Temperature,ph,me)
real(pReal), dimension(3,3), intent(out) :: & real(pReal), dimension(3,3), intent(out) :: &
Lp !< plastic velocity gradient Lp !< plastic velocity gradient
real(pReal), dimension(3,3,3,3), intent(out) :: & real(pReal), dimension(3,3,3,3), intent(out) :: &
dLp_dMp dLp_dMp
integer, intent(in) :: & integer, intent(in) :: &
ph, & ph, &
me, & me
ip, & !< current integration point
el !< current element number
real(pReal), intent(in) :: & real(pReal), intent(in) :: &
Temperature !< temperature Temperature !< temperature
@ -814,7 +811,7 @@ module subroutine nonlocal_LpAndItsTangent(Lp,dLp_dMp, &
ns = prm%sum_N_sl ns = prm%sum_N_sl
!*** shortcut to state variables !*** shortcut to state variables
rho = getRho(ph,me,ip,el) rho = getRho(ph,me)
rhoSgl = rho(:,sgl) rhoSgl = rho(:,sgl)
do s = 1,ns do s = 1,ns
@ -880,18 +877,15 @@ end subroutine nonlocal_LpAndItsTangent
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief (instantaneous) incremental change of microstructure !> @brief (instantaneous) incremental change of microstructure
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module subroutine plastic_nonlocal_deltaState(Mp,ph,me,ip,el) module subroutine plastic_nonlocal_deltaState(Mp,ph,me)
real(pReal), dimension(3,3), intent(in) :: & real(pReal), dimension(3,3), intent(in) :: &
Mp !< MandelStress Mp !< MandelStress
integer, intent(in) :: & integer, intent(in) :: &
ph, & ph, &
me, & !< offset me
ip, &
el
integer :: & integer :: &
ph, & !< phase
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
t, & ! type of dislocation t, & ! type of dislocation
@ -918,7 +912,7 @@ module subroutine plastic_nonlocal_deltaState(Mp,ph,me,ip,el)
forall (s = 1:ns, t = 1:4) v(s,t) = plasticState(ph)%state(iV(s,t,ins(ph)),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,ins(ph)),me) forall (s = 1:ns, c = 1:2) dUpperOld(s,c) = plasticState(ph)%state(iD(s,c,ins(ph)),me)
rho = getRho(ph,me,ip,el) rho = getRho(ph,me)
rhoDip = rho(:,dip) rhoDip = rho(:,dip)
!**************************************************************************** !****************************************************************************
@ -968,15 +962,6 @@ module subroutine plastic_nonlocal_deltaState(Mp,ph,me,ip,el)
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])
#ifdef DEBUG
if (debugConstitutive%extensive &
.and. ((debugConstitutive%element == el .and. debugConstitutive%ip == ip)&
.or. .not. debugConstitutive%selective)) then
print'(a,/,8(12x,12(e12.5,1x),/))', '<< CONST >> dislocation remobilization', deltaRhoRemobilization(:,1:8)
print'(a,/,10(12x,12(e12.5,1x),/),/)', '<< CONST >> dipole dissociation by stress increase', deltaRhoDipole2SingleStress
endif
#endif
end associate end associate
end subroutine plastic_nonlocal_deltaState end subroutine plastic_nonlocal_deltaState
@ -1043,10 +1028,10 @@ module subroutine nonlocal_dotState(Mp, Temperature,timestep, &
tau = 0.0_pReal tau = 0.0_pReal
gdot = 0.0_pReal gdot = 0.0_pReal
rho = getRho(ph,me,ip,el) rho = getRho(ph,me)
rhoSgl = rho(:,sgl) rhoSgl = rho(:,sgl)
rhoDip = rho(:,dip) rhoDip = rho(:,dip)
rho0 = getRho0(ph,me,ip,el) rho0 = getRho0(ph,me)
my_rhoSgl0 = rho0(:,sgl) my_rhoSgl0 = rho0(:,sgl)
forall (s = 1:ns, t = 1:4) v(s,t) = plasticState(ph)%state(iV(s,t,ins(ph)),me) forall (s = 1:ns, t = 1:4) v(s,t) = plasticState(ph)%state(iV(s,t,ins(ph)),me)
@ -1194,7 +1179,6 @@ function rhoDotFlux(timestep,ph,me,ip,el)
el !< current element number el !< current element number
integer :: & integer :: &
ph, &
neighbor_instance, & !< instance of my neighbor's plasticity neighbor_instance, & !< instance of my neighbor's plasticity
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
@ -1251,9 +1235,9 @@ function rhoDotFlux(timestep,ph,me,ip,el)
gdot = 0.0_pReal gdot = 0.0_pReal
rho = getRho(ph,me,ip,el) rho = getRho(ph,me)
rhoSgl = rho(:,sgl) rhoSgl = rho(:,sgl)
rho0 = getRho0(ph,me,ip,el) rho0 = getRho0(ph,me)
my_rhoSgl0 = rho0(:,sgl) my_rhoSgl0 = rho0(:,sgl)
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 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
@ -1432,7 +1416,6 @@ module subroutine plastic_nonlocal_updateCompatibility(orientation,ph,i,e)
n, & ! neighbor index n, & ! neighbor index
neighbor_e, & ! element index of my neighbor neighbor_e, & ! element index of my neighbor
neighbor_i, & ! integration point index of my neighbor neighbor_i, & ! integration point index of my neighbor
ph, &
neighbor_phase, & neighbor_phase, &
ns, & ! number of active slip systems ns, & ! number of active slip systems
s1, & ! slip system index (me) s1, & ! slip system index (me)
@ -1604,7 +1587,7 @@ end subroutine plastic_nonlocal_results
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief populates the initial dislocation density !> @brief populates the initial dislocation density
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine stateInit(ini,phase,Nconstituents) subroutine stateInit(ini,phase,Nconstituents,i)
type(tInitialParameters) :: & type(tInitialParameters) :: &
ini ini
@ -1631,7 +1614,7 @@ subroutine stateInit(ini,phase,Nconstituents)
volume volume
associate(stt => state(ins(phase))) associate(stt => state(i))
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
@ -1794,21 +1777,22 @@ 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(ph,me,ip,el) pure function getRho(ph,me)
integer, intent(in) :: ph, me,ip,el integer, intent(in) :: ph, me
real(pReal), dimension(param(ins(ph))%sum_N_sl,10) :: getRho real(pReal), dimension(param(ins(ph))%sum_N_sl,10) :: getRho
associate(prm => param(ins(ph))) associate(prm => param(ins(ph)))
getRho = reshape(state(ins(ph))%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)
getRho(:,dip) = max(getRho(:,dip),0.0_pReal) getRho(:,dip) = max(getRho(:,dip),0.0_pReal)
where(abs(getRho) < max(prm%rho_min/IPvolume(ip,el)**(2.0_pReal/3.0_pReal),prm%rho_significant)) & where(abs(getRho) < max(prm%rho_min/geom(ph)%V_0(me)**(2.0_pReal/3.0_pReal),prm%rho_significant)) &
getRho = 0.0_pReal getRho = 0.0_pReal
end associate end associate
@ -1819,21 +1803,22 @@ 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,me,ip,el) pure function getRho0(ph,me)
integer, intent(in) :: ph, me,ip,el integer, intent(in) :: ph, me
real(pReal), dimension(param(ins(ph))%sum_N_sl,10) :: getRho0 real(pReal), dimension(param(ins(ph))%sum_N_sl,10) :: getRho0
associate(prm => param(ins(ph))) associate(prm => param(ins(ph)))
getRho0 = reshape(state0(ins(ph))%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)
getRho0(:,dip) = max(getRho0(:,dip),0.0_pReal) getRho0(:,dip) = max(getRho0(:,dip),0.0_pReal)
where(abs(getRho0) < max(prm%rho_min/IPvolume(ip,el)**(2.0_pReal/3.0_pReal),prm%rho_significant)) & where(abs(getRho0) < max(prm%rho_min/geom(ph)%V_0(me)**(2.0_pReal/3.0_pReal),prm%rho_significant)) &
getRho0 = 0.0_pReal getRho0 = 0.0_pReal
end associate end associate