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

explicit dotState for nonlocal 

all flux related quantities are calculated based on the converged
quantities
This commit is contained in:
Martin Diehl 2020-02-07 12:23:22 +01:00
parent 4f4c6c5949
commit f854dc27e9
3 changed files with 226 additions and 227 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

194
src/constitutive.f90
View File

@ -25,44 +25,44 @@ module constitutive
use kinematics_cleavage_opening
use kinematics_slipplane_opening
use kinematics_thermal_expansion
implicit none
private
integer, public, protected :: &
constitutive_plasticity_maxSizeDotState, &
constitutive_source_maxSizeDotState
interface
module subroutine plastic_none_init
end subroutine plastic_none_init
module subroutine plastic_isotropic_init
end subroutine plastic_isotropic_init
module subroutine plastic_phenopowerlaw_init
end subroutine plastic_phenopowerlaw_init
module subroutine plastic_kinehardening_init
end subroutine plastic_kinehardening_init
module subroutine plastic_dislotwin_init
end subroutine plastic_dislotwin_init
module subroutine plastic_disloUCLA_init
end subroutine plastic_disloUCLA_init
module subroutine plastic_nonlocal_init
end subroutine plastic_nonlocal_init
module subroutine plastic_isotropic_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,of)
real(pReal), dimension(3,3), intent(out) :: &
Lp !< plastic velocity gradient
real(pReal), dimension(3,3,3,3), intent(out) :: &
dLp_dMp !< derivative of Lp with respect to the Mandel stress
real(pReal), dimension(3,3), intent(in) :: &
Mp !< Mandel stress
integer, intent(in) :: &
@ -75,20 +75,20 @@ module constitutive
Lp !< plastic velocity gradient
real(pReal), dimension(3,3,3,3), intent(out) :: &
dLp_dMp !< derivative of Lp with respect to the Mandel stress
real(pReal), dimension(3,3), intent(in) :: &
Mp !< Mandel stress
integer, intent(in) :: &
instance, &
of
end subroutine plastic_phenopowerlaw_LpAndItsTangent
pure module subroutine plastic_kinehardening_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,of)
real(pReal), dimension(3,3), intent(out) :: &
Lp !< plastic velocity gradient
real(pReal), dimension(3,3,3,3), intent(out) :: &
dLp_dMp !< derivative of Lp with respect to the Mandel stress
real(pReal), dimension(3,3), intent(in) :: &
Mp !< Mandel stress
integer, intent(in) :: &
@ -101,7 +101,7 @@ module constitutive
Lp !< plastic velocity gradient
real(pReal), dimension(3,3,3,3), intent(out) :: &
dLp_dMp !< derivative of Lp with respect to the Mandel stress
real(pReal), dimension(3,3), intent(in) :: &
Mp !< Mandel stress
real(pReal), intent(in) :: &
@ -110,13 +110,13 @@ module constitutive
instance, &
of
end subroutine plastic_dislotwin_LpAndItsTangent
pure module subroutine plastic_disloUCLA_LpAndItsTangent(Lp,dLp_dMp,Mp,T,instance,of)
real(pReal), dimension(3,3), intent(out) :: &
Lp !< plastic velocity gradient
real(pReal), dimension(3,3,3,3), intent(out) :: &
dLp_dMp !< derivative of Lp with respect to the Mandel stress
real(pReal), dimension(3,3), intent(in) :: &
Mp !< Mandel stress
real(pReal), intent(in) :: &
@ -125,14 +125,14 @@ module constitutive
instance, &
of
end subroutine plastic_disloUCLA_LpAndItsTangent
module subroutine plastic_nonlocal_LpAndItsTangent(Lp, dLp_dMp, &
Mp, Temperature, volume, ip, el)
real(pReal), dimension(3,3), intent(out) :: &
Lp !< plastic velocity gradient
real(pReal), dimension(3,3,3,3), intent(out) :: &
dLp_dMp !< derivative of Lp with respect to the Mandel stress
real(pReal), dimension(3,3), intent(in) :: &
Mp !< Mandel stress
real(pReal), intent(in) :: &
@ -149,15 +149,15 @@ module constitutive
Li !< inleastic velocity gradient
real(pReal), dimension(3,3,3,3), intent(out) :: &
dLi_dMi !< derivative of Li with respect to Mandel stress
real(pReal), dimension(3,3), intent(in) :: &
Mi !< Mandel stress
Mi !< Mandel stress
integer, intent(in) :: &
instance, &
of
end subroutine plastic_isotropic_LiAndItsTangent
module subroutine plastic_isotropic_dotState(Mp,instance,of)
real(pReal), dimension(3,3), intent(in) :: &
Mp !< Mandel stress
@ -165,7 +165,7 @@ module constitutive
instance, &
of
end subroutine plastic_isotropic_dotState
module subroutine plastic_phenopowerlaw_dotState(Mp,instance,of)
real(pReal), dimension(3,3), intent(in) :: &
Mp !< Mandel stress
@ -173,7 +173,7 @@ module constitutive
instance, &
of
end subroutine plastic_phenopowerlaw_dotState
module subroutine plastic_kinehardening_dotState(Mp,instance,of)
real(pReal), dimension(3,3), intent(in) :: &
Mp !< Mandel stress
@ -181,7 +181,7 @@ module constitutive
instance, &
of
end subroutine plastic_kinehardening_dotState
module subroutine plastic_dislotwin_dotState(Mp,T,instance,of)
real(pReal), dimension(3,3), intent(in) :: &
Mp !< Mandel stress
@ -201,8 +201,8 @@ module constitutive
instance, &
of
end subroutine plastic_disloUCLA_dotState
module subroutine plastic_nonlocal_dotState(Mp, Fe, Fp, Temperature, &
module subroutine plastic_nonlocal_dotState(Mp, F, Fp, Temperature, &
timestep,ip,el)
integer, intent(in) :: &
ip, & !< current integration point
@ -213,11 +213,11 @@ module constitutive
real(pReal), dimension(3,3), intent(in) ::&
Mp !< MandelStress
real(pReal), dimension(3,3,homogenization_maxNgrains,discretization_nIP,discretization_nElem), intent(in) :: &
Fe, & !< elastic deformation gradient
F, & !< deformation gradient
Fp !< plastic deformation gradient
end subroutine plastic_nonlocal_dotState
module subroutine plastic_dislotwin_dependentState(T,instance,of)
integer, intent(in) :: &
instance, &
@ -225,13 +225,13 @@ module constitutive
real(pReal), intent(in) :: &
T
end subroutine plastic_dislotwin_dependentState
module subroutine plastic_disloUCLA_dependentState(instance,of)
integer, intent(in) :: &
instance, &
of
end subroutine plastic_disloUCLA_dependentState
module subroutine plastic_nonlocal_dependentState(Fe, Fp, ip, el)
integer, intent(in) :: &
ip, &
@ -249,7 +249,7 @@ module constitutive
instance, &
of
end subroutine plastic_kinehardening_deltaState
module subroutine plastic_nonlocal_deltaState(Mp,ip,el)
integer, intent(in) :: &
ip, &
@ -267,48 +267,48 @@ module constitutive
ip, & !< integration point
el !< element
end function plastic_dislotwin_homogenizedC
module subroutine plastic_nonlocal_updateCompatibility(orientation,i,e)
module subroutine plastic_nonlocal_updateCompatibility(orientation,i,e)
integer, intent(in) :: &
i, &
e
type(rotation), dimension(1,discretization_nIP,discretization_nElem), intent(in) :: &
orientation !< crystal orientation
end subroutine plastic_nonlocal_updateCompatibility
module subroutine plastic_isotropic_results(instance,group)
integer, intent(in) :: instance
character(len=*), intent(in) :: group
end subroutine plastic_isotropic_results
module subroutine plastic_phenopowerlaw_results(instance,group)
integer, intent(in) :: instance
character(len=*), intent(in) :: group
end subroutine plastic_phenopowerlaw_results
module subroutine plastic_kinehardening_results(instance,group)
integer, intent(in) :: instance
character(len=*), intent(in) :: group
end subroutine plastic_kinehardening_results
module subroutine plastic_dislotwin_results(instance,group)
integer, intent(in) :: instance
character(len=*), intent(in) :: group
end subroutine plastic_dislotwin_results
module subroutine plastic_disloUCLA_results(instance,group)
integer, intent(in) :: instance
character(len=*), intent(in) :: group
end subroutine plastic_disloUCLA_results
module subroutine plastic_nonlocal_results(instance,group)
integer, intent(in) :: instance
character(len=*), intent(in) :: group
end subroutine plastic_nonlocal_results
end interface
public :: &
plastic_nonlocal_updateCompatibility, &
constitutive_init, &
@ -321,7 +321,7 @@ module constitutive
constitutive_collectDotState, &
constitutive_collectDeltaState, &
constitutive_results
contains
@ -355,18 +355,18 @@ subroutine constitutive_init
if (any(phase_source == SOURCE_damage_isoDuctile_ID)) call source_damage_isoDuctile_init
if (any(phase_source == SOURCE_damage_anisoBrittle_ID)) call source_damage_anisoBrittle_init
if (any(phase_source == SOURCE_damage_anisoDuctile_ID)) call source_damage_anisoDuctile_init
!--------------------------------------------------------------------------------------------------
! initialize kinematic mechanisms
if (any(phase_kinematics == KINEMATICS_cleavage_opening_ID)) call kinematics_cleavage_opening_init
if (any(phase_kinematics == KINEMATICS_slipplane_opening_ID)) call kinematics_slipplane_opening_init
if (any(phase_kinematics == KINEMATICS_thermal_expansion_ID)) call kinematics_thermal_expansion_init
write(6,'(/,a)') ' <<<+- constitutive init -+>>>'; flush(6)
constitutive_plasticity_maxSizeDotState = 0
constitutive_source_maxSizeDotState = 0
PhaseLoop2:do ph = 1,material_Nphase
!--------------------------------------------------------------------------------------------------
! partition and inititalize state
@ -398,7 +398,7 @@ function constitutive_homogenizedC(ipc,ip,el)
ipc, & !< component-ID of integration point
ip, & !< integration point
el !< element
plasticityType: select case (phase_plasticity(material_phaseAt(ipc,el)))
case (PLASTICITY_DISLOTWIN_ID) plasticityType
constitutive_homogenizedC = plastic_dislotwin_homogenizedC(ipc,ip,el)
@ -425,10 +425,10 @@ subroutine constitutive_microstructure(Fe, Fp, ipc, ip, el)
ho, & !< homogenization
tme, & !< thermal member position
instance, of
ho = material_homogenizationAt(el)
tme = thermalMapping(ho)%p(ip,el)
plasticityType: select case (phase_plasticity(material_phaseAt(ipc,el)))
case (PLASTICITY_DISLOTWIN_ID) plasticityType
of = material_phasememberAt(ipc,ip,el)
@ -451,7 +451,7 @@ end subroutine constitutive_microstructure
! Mp in, dLp_dMp out
!--------------------------------------------------------------------------------------------------
subroutine constitutive_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, &
S, Fi, ipc, ip, el)
S, Fi, ipc, ip, el)
integer, intent(in) :: &
ipc, & !< component-ID of integration point
ip, & !< integration point
@ -473,49 +473,49 @@ subroutine constitutive_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, &
tme !< thermal member position
integer :: &
i, j, instance, of
ho = material_homogenizationAt(el)
tme = thermalMapping(ho)%p(ip,el)
Mp = matmul(matmul(transpose(Fi),Fi),S)
plasticityType: select case (phase_plasticity(material_phaseAt(ipc,el)))
case (PLASTICITY_NONE_ID) plasticityType
Lp = 0.0_pReal
dLp_dMp = 0.0_pReal
case (PLASTICITY_ISOTROPIC_ID) plasticityType
of = material_phasememberAt(ipc,ip,el)
instance = phase_plasticityInstance(material_phaseAt(ipc,el))
call plastic_isotropic_LpAndItsTangent (Lp,dLp_dMp,Mp,instance,of)
case (PLASTICITY_PHENOPOWERLAW_ID) plasticityType
of = material_phasememberAt(ipc,ip,el)
instance = phase_plasticityInstance(material_phaseAt(ipc,el))
call plastic_phenopowerlaw_LpAndItsTangent (Lp,dLp_dMp,Mp,instance,of)
case (PLASTICITY_KINEHARDENING_ID) plasticityType
of = material_phasememberAt(ipc,ip,el)
instance = phase_plasticityInstance(material_phaseAt(ipc,el))
call plastic_kinehardening_LpAndItsTangent (Lp,dLp_dMp, Mp,instance,of)
case (PLASTICITY_NONLOCAL_ID) plasticityType
call plastic_nonlocal_LpAndItsTangent (Lp,dLp_dMp,Mp, &
temperature(ho)%p(tme),geometry_plastic_nonlocal_IPvolume0(ip,el),ip,el)
case (PLASTICITY_DISLOTWIN_ID) plasticityType
of = material_phasememberAt(ipc,ip,el)
instance = phase_plasticityInstance(material_phaseAt(ipc,el))
call plastic_dislotwin_LpAndItsTangent (Lp,dLp_dMp,Mp,temperature(ho)%p(tme),instance,of)
case (PLASTICITY_DISLOUCLA_ID) plasticityType
of = material_phasememberAt(ipc,ip,el)
instance = phase_plasticityInstance(material_phaseAt(ipc,el))
call plastic_disloucla_LpAndItsTangent (Lp,dLp_dMp,Mp,temperature(ho)%p(tme),instance,of)
end select plasticityType
do i=1,3; do j=1,3
dLp_dFi(i,j,1:3,1:3) = matmul(matmul(Fi,S),transpose(dLp_dMp(i,j,1:3,1:3))) + &
matmul(matmul(Fi,dLp_dMp(i,j,1:3,1:3)),S)
@ -545,7 +545,7 @@ subroutine constitutive_LiAndItsTangents(Li, dLi_dS, dLi_dFi, &
real(pReal), intent(out), dimension(3,3,3,3) :: &
dLi_dS, & !< derivative of Li with respect to S
dLi_dFi
real(pReal), dimension(3,3) :: &
my_Li, & !< intermediate velocity gradient
FiInv, &
@ -557,11 +557,11 @@ subroutine constitutive_LiAndItsTangents(Li, dLi_dS, dLi_dFi, &
integer :: &
k, i, j, &
instance, of
Li = 0.0_pReal
dLi_dS = 0.0_pReal
dLi_dFi = 0.0_pReal
plasticityType: select case (phase_plasticity(material_phaseAt(ipc,el)))
case (PLASTICITY_isotropic_ID) plasticityType
of = material_phasememberAt(ipc,ip,el)
@ -571,10 +571,10 @@ subroutine constitutive_LiAndItsTangents(Li, dLi_dS, dLi_dFi, &
my_Li = 0.0_pReal
my_dLi_dS = 0.0_pReal
end select plasticityType
Li = Li + my_Li
dLi_dS = dLi_dS + my_dLi_dS
KinematicsLoop: do k = 1, phase_Nkinematics(material_phaseAt(ipc,el))
kinematicsType: select case (phase_kinematics(k,material_phaseAt(ipc,el)))
case (KINEMATICS_cleavage_opening_ID) kinematicsType
@ -590,12 +590,12 @@ subroutine constitutive_LiAndItsTangents(Li, dLi_dS, dLi_dFi, &
Li = Li + my_Li
dLi_dS = dLi_dS + my_dLi_dS
enddo KinematicsLoop
FiInv = math_inv33(Fi)
detFi = math_det33(Fi)
Li = matmul(matmul(Fi,Li),FiInv)*detFi !< push forward to intermediate configuration
temp_33 = matmul(FiInv,Li)
do i = 1,3; do j = 1,3
dLi_dS(1:3,1:3,i,j) = matmul(matmul(Fi,dLi_dS(1:3,1:3,i,j)),FiInv)*detFi
dLi_dFi(1:3,1:3,i,j) = dLi_dFi(1:3,1:3,i,j) + Li*FiInv(j,i)
@ -621,10 +621,10 @@ pure function constitutive_initialFi(ipc, ip, el)
integer :: &
phase, &
homog, offset
constitutive_initialFi = math_I3
phase = material_phaseAt(ipc,el)
KinematicsLoop: do k = 1, phase_Nkinematics(phase) !< Warning: small initial strain assumption
kinematicsType: select case (phase_kinematics(k,phase))
case (KINEMATICS_thermal_expansion_ID) kinematicsType
@ -640,7 +640,7 @@ end function constitutive_initialFi
!--------------------------------------------------------------------------------------------------
!> @brief returns the 2nd Piola-Kirchhoff stress tensor and its tangent with respect to
!> the elastic/intermediate deformation gradients depending on the selected elastic law
!> the elastic/intermediate deformation gradients depending on the selected elastic law
!! (so far no case switch because only Hooke is implemented)
!--------------------------------------------------------------------------------------------------
subroutine constitutive_SandItsTangents(S, dS_dFe, dS_dFi, Fe, Fi, ipc, ip, el)
@ -690,17 +690,17 @@ subroutine constitutive_hooke_SandItsTangents(S, dS_dFe, dS_dFi, &
d !< counter in degradation loop
integer :: &
i, j
ho = material_homogenizationAt(el)
C = math_66toSym3333(constitutive_homogenizedC(ipc,ip,el))
DegradationLoop: do d = 1, phase_NstiffnessDegradations(material_phaseAt(ipc,el))
degradationType: select case(phase_stiffnessDegradation(d,material_phaseAt(ipc,el)))
case (STIFFNESS_DEGRADATION_damage_ID) degradationType
C = C * damage(ho)%p(damageMapping(ho)%p(ip,el))**2
end select degradationType
enddo DegradationLoop
E = 0.5_pReal*(matmul(transpose(Fe),Fe)-math_I3) !< Green-Lagrange strain in unloaded configuration
S = math_mul3333xx33(C,matmul(matmul(transpose(Fi),E),Fi)) !< 2PK stress in lattice configuration in work conjugate with GL strain pulled back to lattice configuration
@ -715,7 +715,7 @@ end subroutine constitutive_hooke_SandItsTangents
!--------------------------------------------------------------------------------------------------
!> @brief contains the constitutive equation for calculating the rate of change of microstructure
!--------------------------------------------------------------------------------------------------
subroutine constitutive_collectDotState(S, FeArray, Fi, FpArray, subdt, ipc, ip, el)
subroutine constitutive_collectDotState(S, FArray, Fi, FpArray, subdt, ipc, ip, el)
integer, intent(in) :: &
ipc, & !< component-ID of integration point
@ -724,7 +724,7 @@ subroutine constitutive_collectDotState(S, FeArray, Fi, FpArray, subdt, ipc, ip,
real(pReal), intent(in) :: &
subdt !< timestep
real(pReal), intent(in), dimension(3,3,homogenization_maxNgrains,discretization_nIP,discretization_nElem) :: &
FeArray, & !< elastic deformation gradient
FArray, & !< elastic deformation gradient
FpArray !< plastic deformation gradient
real(pReal), intent(in), dimension(3,3) :: &
Fi !< intermediate deformation gradient
@ -771,7 +771,7 @@ subroutine constitutive_collectDotState(S, FeArray, Fi, FpArray, subdt, ipc, ip,
call plastic_disloucla_dotState (Mp,temperature(ho)%p(tme),instance,of)
case (PLASTICITY_NONLOCAL_ID) plasticityType
call plastic_nonlocal_dotState (Mp,FeArray,FpArray,temperature(ho)%p(tme), &
call plastic_nonlocal_dotState (Mp,FArray,FpArray,temperature(ho)%p(tme), &
subdt,ip,el)
end select plasticityType
@ -858,32 +858,32 @@ subroutine constitutive_results
do p=1,size(config_name_phase)
group = trim('current/constituent')//'/'//trim(config_name_phase(p))
call results_closeGroup(results_addGroup(group))
group = trim(group)//'/plastic'
call results_closeGroup(results_addGroup(group))
call results_closeGroup(results_addGroup(group))
select case(phase_plasticity(p))
case(PLASTICITY_ISOTROPIC_ID)
call plastic_isotropic_results(phase_plasticityInstance(p),group)
call plastic_isotropic_results(phase_plasticityInstance(p),group)
case(PLASTICITY_PHENOPOWERLAW_ID)
call plastic_phenopowerlaw_results(phase_plasticityInstance(p),group)
call plastic_phenopowerlaw_results(phase_plasticityInstance(p),group)
case(PLASTICITY_KINEHARDENING_ID)
call plastic_kinehardening_results(phase_plasticityInstance(p),group)
call plastic_kinehardening_results(phase_plasticityInstance(p),group)
case(PLASTICITY_DISLOTWIN_ID)
call plastic_dislotwin_results(phase_plasticityInstance(p),group)
call plastic_dislotwin_results(phase_plasticityInstance(p),group)
case(PLASTICITY_DISLOUCLA_ID)
call plastic_disloUCLA_results(phase_plasticityInstance(p),group)
call plastic_disloUCLA_results(phase_plasticityInstance(p),group)
case(PLASTICITY_NONLOCAL_ID)
call plastic_nonlocal_results(phase_plasticityInstance(p),group)
call plastic_nonlocal_results(phase_plasticityInstance(p),group)
end select
enddo
enddo
end subroutine constitutive_results

25
src/constitutive_plastic_nonlocal.f90
View File

@ -1324,7 +1324,7 @@ end subroutine plastic_nonlocal_deltaState
!---------------------------------------------------------------------------------------------------
!> @brief calculates the rate of change of microstructure
!---------------------------------------------------------------------------------------------------
module subroutine plastic_nonlocal_dotState(Mp, Fe, Fp, Temperature, &
module subroutine plastic_nonlocal_dotState(Mp, F, Fp, Temperature, &
timestep,ip,el)
integer, intent(in) :: &
@ -1336,7 +1336,7 @@ module subroutine plastic_nonlocal_dotState(Mp, Fe, Fp, Temperature, &
real(pReal), dimension(3,3), intent(in) ::&
Mp !< MandelStress
real(pReal), dimension(3,3,homogenization_maxNgrains,discretization_nIP,discretization_nElem), intent(in) :: &
Fe, & !< elastic deformation gradient
F, & !< elastic deformation gradient
Fp !< plastic deformation gradient
integer :: &
@ -1370,7 +1370,7 @@ module subroutine plastic_nonlocal_dotState(Mp, Fe, Fp, Temperature, &
rhoDotThermalAnnihilation !< density evolution by thermal annihilation
real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phaseAt(1,el))),8) :: &
rhoSgl, & !< current single dislocation densities (positive/negative screw and edge without dipoles)
neighbor_rhoSgl, & !< 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_rhoSgl !< single dislocation densities of central ip (positive/negative screw and edge without dipoles)
real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phaseAt(1,el))),4) :: &
v, & !< current dislocation glide velocity
@ -1529,8 +1529,8 @@ module subroutine plastic_nonlocal_dotState(Mp, Fe, Fp, Temperature, &
m(1:3,1:ns,3) = -prm%slip_transverse
m(1:3,1:ns,4) = prm%slip_transverse
my_Fe = Fe(1:3,1:3,1,ip,el)
my_F = matmul(my_Fe, Fp(1:3,1:3,1,ip,el))
my_F = F(1:3,1:3,1,ip,el)
my_Fe = matmul(my_F, math_inv33(Fp(1:3,1:3,1,ip,el)))
neighbors: do n = 1,nIPneighbors
@ -1547,8 +1547,8 @@ module subroutine plastic_nonlocal_dotState(Mp, Fe, Fp, Temperature, &
if (neighbor_n > 0) then ! if neighbor exists, average deformation gradient
neighbor_instance = phase_plasticityInstance(material_phaseAt(1,neighbor_el))
neighbor_Fe = Fe(1:3,1:3,1,neighbor_ip,neighbor_el)
neighbor_F = matmul(neighbor_Fe, Fp(1:3,1:3,1,neighbor_ip,neighbor_el))
neighbor_F = F(1:3,1:3,1,neighbor_ip,neighbor_el)
neighbor_Fe = matmul(neighbor_F, math_inv33(Fp(1:3,1:3,1,neighbor_ip,neighbor_el)))
Favg = 0.5_pReal * (my_F + neighbor_F)
else ! if no neighbor, take my value as average
Favg = my_F
@ -1566,7 +1566,6 @@ module subroutine plastic_nonlocal_dotState(Mp, Fe, Fp, Temperature, &
considerEnteringFlux = .false.
neighbor_v0 = 0.0_pReal ! needed for check of sign change in flux density below
neighbor_rhoSgl = 0.0_pReal
if (neighbor_n > 0) then
if (phase_plasticity(material_phaseAt(1,neighbor_el)) == PLASTICITY_NONLOCAL_ID &
.and. any(compatibility(:,:,:,n,ip,el) > 0.0_pReal)) &
@ -1576,13 +1575,13 @@ module subroutine plastic_nonlocal_dotState(Mp, Fe, Fp, Temperature, &
enteringFlux: if (considerEnteringFlux) then
forall (s = 1:ns, t = 1:4)
neighbor_v0(s,t) = plasticState(np)%state0(iV (s,t,neighbor_instance),no)
neighbor_rhoSgl(s,t) = max(plasticState(np)%state(iRhoU(s,t,neighbor_instance),no), &
neighbor_rhoSgl0(s,t) = max(plasticState(np)%state0(iRhoU(s,t,neighbor_instance),no), &
0.0_pReal)
endforall
where (neighbor_rhoSgl * IPvolume(neighbor_ip,neighbor_el) ** 0.667_pReal < prm%significantN &
.or. neighbor_rhoSgl < prm%significantRho) &
neighbor_rhoSgl = 0.0_pReal
where (neighbor_rhoSgl0 * IPvolume(neighbor_ip,neighbor_el) ** 0.667_pReal < prm%significantN &
.or. neighbor_rhoSgl0 < prm%significantRho) &
neighbor_rhoSgl0 = 0.0_pReal
normal_neighbor2me_defConf = math_det33(Favg) * matmul(math_inv33(transpose(Favg)), &
IPareaNormal(1:3,neighbor_n,neighbor_ip,neighbor_el)) ! calculate the normal of the interface in (average) deformed configuration (now pointing from my neighbor to me!!!)
normal_neighbor2me = matmul(transpose(neighbor_Fe), normal_neighbor2me_defConf) &
@ -1595,7 +1594,7 @@ module subroutine plastic_nonlocal_dotState(Mp, Fe, Fp, Temperature, &
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 me == entering flux for me
.and. v0(s,t) * neighbor_v0(s,t) >= 0.0_pReal ) then ! ... only if no sign change in flux density
lineLength = neighbor_rhoSgl(s,t) * neighbor_v0(s,t) &
lineLength = neighbor_rhoSgl0(s,t) * neighbor_v0(s,t) &
* math_inner(m(1:3,s,t), normal_neighbor2me) * area ! positive line length that wants to enter through this interface
where (compatibility(c,1:ns,s,n,ip,el) > 0.0_pReal) & ! positive compatibility...
rhoDotFlux(1:ns,t) = rhoDotFlux(1:ns,t) &

234
src/crystallite.f90
View File

@ -22,10 +22,10 @@ module crystallite
use discretization
use lattice
use results
implicit none
implicit none
private
real(pReal), dimension(:,:,:), allocatable, public :: &
crystallite_dt !< requested time increment of each grain
real(pReal), dimension(:,:,:), allocatable :: &
@ -33,7 +33,7 @@ module crystallite
crystallite_subFrac, & !< already calculated fraction of increment
crystallite_subStep !< size of next integration step
type(rotation), dimension(:,:,:), allocatable :: &
crystallite_orientation !< current orientation
crystallite_orientation !< current orientation
real(pReal), dimension(:,:,:,:,:), allocatable, public, protected :: &
crystallite_Fe, & !< current "elastic" def grad (end of converged time step)
crystallite_P !< 1st Piola-Kirchhoff stress per grain
@ -74,33 +74,33 @@ module crystallite
crystallite_converged, & !< convergence flag
crystallite_todo, & !< flag to indicate need for further computation
crystallite_localPlasticity !< indicates this grain to have purely local constitutive law
type :: tOutput !< new requested output (per phase)
character(len=pStringLen), allocatable, dimension(:) :: &
label
end type tOutput
type(tOutput), allocatable, dimension(:) :: output_constituent
type :: tNumerics
integer :: &
iJacoLpresiduum, & !< frequency of Jacobian update of residuum in Lp
nState, & !< state loop limit
nStress !< stress loop limit
nStress !< stress loop limit
real(pReal) :: &
subStepMinCryst, & !< minimum (relative) size of sub-step allowed during cutback
subStepSizeCryst, & !< size of first substep when cutback
subStepSizeLp, & !< size of first substep when cutback in Lp calculation
subStepSizeLi, & !< size of first substep when cutback in Li calculation
stepIncreaseCryst, & !< increase of next substep size when previous substep converged
rTol_crystalliteState, & !< relative tolerance in state loop
rTol_crystalliteState, & !< relative tolerance in state loop
rTol_crystalliteStress, & !< relative tolerance in stress loop
aTol_crystalliteStress !< absolute tolerance in stress loop
end type tNumerics
type(tNumerics) :: num ! numerics parameters. Better name?
procedure(), pointer :: integrateState
public :: &
crystallite_init, &
crystallite_stress, &
@ -116,7 +116,7 @@ contains
!> @brief allocates and initialize per grain variables
!--------------------------------------------------------------------------------------------------
subroutine crystallite_init
logical, dimension(discretization_nIP,discretization_nElem) :: devNull
integer :: &
c, & !< counter in integration point component loop
@ -126,13 +126,13 @@ subroutine crystallite_init
iMax, & !< maximum number of integration points
eMax, & !< maximum number of elements
myNcomponents !< number of components at current IP
write(6,'(/,a)') ' <<<+- crystallite init -+>>>'
cMax = homogenization_maxNgrains
iMax = discretization_nIP
eMax = discretization_nElem
allocate(crystallite_S0(3,3,cMax,iMax,eMax), source=0.0_pReal)
allocate(crystallite_partionedS0(3,3,cMax,iMax,eMax), source=0.0_pReal)
allocate(crystallite_S(3,3,cMax,iMax,eMax), source=0.0_pReal)
@ -172,23 +172,23 @@ subroutine crystallite_init
allocate(crystallite_requested(cMax,iMax,eMax), source=.false.)
allocate(crystallite_todo(cMax,iMax,eMax), source=.false.)
allocate(crystallite_converged(cMax,iMax,eMax), source=.true.)
num%subStepMinCryst = config_numerics%getFloat('substepmincryst', defaultVal=1.0e-3_pReal)
num%subStepSizeCryst = config_numerics%getFloat('substepsizecryst', defaultVal=0.25_pReal)
num%stepIncreaseCryst = config_numerics%getFloat('stepincreasecryst', defaultVal=1.5_pReal)
num%subStepSizeLp = config_numerics%getFloat('substepsizelp', defaultVal=0.5_pReal)
num%subStepSizeLi = config_numerics%getFloat('substepsizeli', defaultVal=0.5_pReal)
num%rTol_crystalliteState = config_numerics%getFloat('rtol_crystallitestate', defaultVal=1.0e-6_pReal)
num%rTol_crystalliteStress = config_numerics%getFloat('rtol_crystallitestress',defaultVal=1.0e-6_pReal)
num%aTol_crystalliteStress = config_numerics%getFloat('atol_crystallitestress',defaultVal=1.0e-8_pReal)
num%iJacoLpresiduum = config_numerics%getInt ('ijacolpresiduum', defaultVal=1)
num%nState = config_numerics%getInt ('nstate', defaultVal=20)
num%nStress = config_numerics%getInt ('nstress', defaultVal=40)
if(num%subStepMinCryst <= 0.0_pReal) call IO_error(301,ext_msg='subStepMinCryst')
if(num%subStepSizeCryst <= 0.0_pReal) call IO_error(301,ext_msg='subStepSizeCryst')
if(num%stepIncreaseCryst <= 0.0_pReal) call IO_error(301,ext_msg='stepIncreaseCryst')
@ -199,12 +199,12 @@ subroutine crystallite_init
if(num%rTol_crystalliteState <= 0.0_pReal) call IO_error(301,ext_msg='rTol_crystalliteState')
if(num%rTol_crystalliteStress <= 0.0_pReal) call IO_error(301,ext_msg='rTol_crystalliteStress')
if(num%aTol_crystalliteStress <= 0.0_pReal) call IO_error(301,ext_msg='aTol_crystalliteStress')
if(num%iJacoLpresiduum < 1) call IO_error(301,ext_msg='iJacoLpresiduum')
if(num%nState < 1) call IO_error(301,ext_msg='nState')
if(num%nStress< 1) call IO_error(301,ext_msg='nStress')
select case(numerics_integrator)
case(1)
integrateState => integrateStateFPI
@ -217,11 +217,11 @@ subroutine crystallite_init
case(5)
integrateState => integrateStateRKCK45
end select
allocate(output_constituent(size(config_phase)))
do c = 1, size(config_phase)
#if defined(__GFORTRAN__)
allocate(output_constituent(c)%label(1))
allocate(output_constituent(c)%label(1))
output_constituent(c)%label(1)= 'GfortranBug86277'
output_constituent(c)%label = config_phase(c)%getStrings('(output)',defaultVal=output_constituent(c)%label )
if (output_constituent(c)%label (1) == 'GfortranBug86277') output_constituent(c)%label = [character(len=pStringLen)::]
@ -250,16 +250,16 @@ subroutine crystallite_init
enddo; enddo
enddo
!$OMP END PARALLEL DO
if(any(.not. crystallite_localPlasticity) .and. .not. usePingPong) call IO_error(601) ! exit if nonlocal but no ping-pong ToDo: Why not check earlier? or in nonlocal?
crystallite_partionedFp0 = crystallite_Fp0
crystallite_partionedFi0 = crystallite_Fi0
crystallite_partionedF0 = crystallite_F0
crystallite_partionedF = crystallite_F0
call crystallite_orientations()
!$OMP PARALLEL DO
do e = FEsolving_execElem(1),FEsolving_execElem(2)
do i = FEsolving_execIP(1),FEsolving_execIP(2)
@ -271,7 +271,7 @@ subroutine crystallite_init
enddo
enddo
!$OMP END PARALLEL DO
devNull = crystallite_stress()
call crystallite_stressTangent
@ -283,7 +283,7 @@ subroutine crystallite_init
write(6,'(a42,1x,i10)') ' # of nonlocal constituents: ',count(.not. crystallite_localPlasticity)
flush(6)
endif
call debug_info
call debug_reset
#endif
@ -295,7 +295,7 @@ end subroutine crystallite_init
!> @brief calculate stress (P)
!--------------------------------------------------------------------------------------------------
function crystallite_stress(dummyArgumentToPreventInternalCompilerErrorWithGCC)
logical, dimension(discretization_nIP,discretization_nElem) :: crystallite_stress
real(pReal), intent(in), optional :: &
dummyArgumentToPreventInternalCompilerErrorWithGCC
@ -308,7 +308,7 @@ function crystallite_stress(dummyArgumentToPreventInternalCompilerErrorWithGCC)
e, & !< counter in element loop
startIP, endIP, &
s
#ifdef DEBUG
if (iand(debug_level(debug_crystallite),debug_levelSelective) /= 0 &
.and. FEsolving_execElem(1) <= debug_e &
@ -494,7 +494,7 @@ function crystallite_stress(dummyArgumentToPreventInternalCompilerErrorWithGCC)
crystallite_stress = .false.
elementLooping5: do e = FEsolving_execElem(1),FEsolving_execElem(2)
do i = FEsolving_execIP(1),FEsolving_execIP(2)
crystallite_stress(i,e) = all(crystallite_converged(:,i,e))
crystallite_stress(i,e) = all(crystallite_converged(:,i,e))
enddo
enddo elementLooping5
@ -675,12 +675,12 @@ end subroutine crystallite_stressTangent
!> @brief calculates orientations
!--------------------------------------------------------------------------------------------------
subroutine crystallite_orientations
integer &
c, & !< counter in integration point component loop
i, & !< counter in integration point loop
e !< counter in element loop
!$OMP PARALLEL DO
do e = FEsolving_execElem(1),FEsolving_execElem(2)
do i = FEsolving_execIP(1),FEsolving_execIP(2)
@ -688,7 +688,7 @@ subroutine crystallite_orientations
call crystallite_orientation(c,i,e)%fromMatrix(transpose(math_rotationalPart33(crystallite_Fe(1:3,1:3,c,i,e))))
enddo; enddo; enddo
!$OMP END PARALLEL DO
nonlocalPresent: if (any(plasticState%nonLocal)) then
!$OMP PARALLEL DO
do e = FEsolving_execElem(1),FEsolving_execElem(2)
@ -706,7 +706,7 @@ end subroutine crystallite_orientations
!> @brief Map 2nd order tensor to reference config
!--------------------------------------------------------------------------------------------------
function crystallite_push33ToRef(ipc,ip,el, tensor33)
real(pReal), dimension(3,3) :: crystallite_push33ToRef
real(pReal), dimension(3,3), intent(in) :: tensor33
real(pReal), dimension(3,3) :: T
@ -714,7 +714,7 @@ function crystallite_push33ToRef(ipc,ip,el, tensor33)
el, &
ip, &
ipc
T = matmul(material_orientation0(ipc,ip,el)%asMatrix(), & ! ToDo: initial orientation correct?
transpose(math_inv33(crystallite_subF(1:3,1:3,ipc,ip,el))))
crystallite_push33ToRef = matmul(transpose(T),matmul(tensor33,T))
@ -731,11 +731,11 @@ subroutine crystallite_results
real(pReal), allocatable, dimension(:,:,:) :: selected_tensors
type(rotation), allocatable, dimension(:) :: selected_rotations
character(len=256) :: group,lattice_label
do p=1,size(config_name_phase)
group = trim('current/constituent')//'/'//trim(config_name_phase(p))//'/generic'
call results_closeGroup(results_addGroup(group))
call results_closeGroup(results_addGroup(group))
do o = 1, size(output_constituent(p)%label)
select case (output_constituent(p)%label(o))
@ -792,19 +792,19 @@ subroutine crystallite_results
end select
enddo
enddo
contains
!------------------------------------------------------------------------------------------------
!> @brief select tensors for output
!------------------------------------------------------------------------------------------------
function select_tensors(dataset,instance)
integer, intent(in) :: instance
real(pReal), dimension(:,:,:,:,:), intent(in) :: dataset
real(pReal), allocatable, dimension(:,:,:) :: select_tensors
integer :: e,i,c,j
allocate(select_tensors(3,3,count(material_phaseAt==instance)*homogenization_maxNgrains*discretization_nIP))
j=0
@ -818,20 +818,20 @@ subroutine crystallite_results
enddo
enddo
enddo
end function select_tensors
!--------------------------------------------------------------------------------------------------
!> @brief select rotations for output
!--------------------------------------------------------------------------------------------------
!--------------------------------------------------------------------------------------------------
function select_rotations(dataset,instance)
integer, intent(in) :: instance
type(rotation), dimension(:,:,:), intent(in) :: dataset
type(rotation), allocatable, dimension(:) :: select_rotations
integer :: e,i,c,j
allocate(select_rotations(count(material_phaseAt==instance)*homogenization_maxNgrains*discretization_nIP))
j=0
@ -845,7 +845,7 @@ subroutine crystallite_results
enddo
enddo
enddo
end function select_rotations
end subroutine crystallite_results
@ -856,12 +856,12 @@ end subroutine crystallite_results
!> intermediate acceleration of the Newton-Raphson correction
!--------------------------------------------------------------------------------------------------
logical function integrateStress(ipc,ip,el,timeFraction)
integer, intent(in):: el, & ! element index
ip, & ! integration point index
ipc ! grain index
real(pReal), optional, intent(in) :: timeFraction ! fraction of timestep
real(pReal), dimension(3,3):: Fg_new, & ! deformation gradient at end of timestep
Fp_new, & ! plastic deformation gradient at end of timestep
Fe_new, & ! elastic deformation gradient at end of timestep
@ -916,7 +916,7 @@ logical function integrateStress(ipc,ip,el,timeFraction)
jacoCounterLi ! counters to check for Jacobian update
external :: &
dgesv
!* be pessimistic
integrateStress = .false.
#ifdef DEBUG
@ -938,7 +938,7 @@ logical function integrateStress(ipc,ip,el,timeFraction)
Lpguess = crystallite_Lp(1:3,1:3,ipc,ip,el) ! take as first guess
Liguess = crystallite_Li(1:3,1:3,ipc,ip,el) ! take as first guess
Liguess_old = Liguess
invFp_current = math_inv33(crystallite_subFp0(1:3,1:3,ipc,ip,el))
failedInversionFp: if (all(dEq0(invFp_current))) then
#ifdef DEBUG
@ -964,13 +964,13 @@ logical function integrateStress(ipc,ip,el,timeFraction)
#endif
return
endif failedInversionFi
!* start Li loop with normal step length
NiterationStressLi = 0
jacoCounterLi = 0
steplengthLi = 1.0_pReal
residuumLi_old = 0.0_pReal
LiLoop: do
NiterationStressLi = NiterationStressLi + 1
LiLoopLimit: if (NiterationStressLi > num%nStress) then
@ -981,18 +981,18 @@ logical function integrateStress(ipc,ip,el,timeFraction)
#endif
return
endif LiLoopLimit
invFi_new = matmul(invFi_current,math_I3 - dt*Liguess)
Fi_new = math_inv33(invFi_new)
detInvFi = math_det33(invFi_new)
!* start Lp loop with normal step length
NiterationStressLp = 0
jacoCounterLp = 0
steplengthLp = 1.0_pReal
residuumLp_old = 0.0_pReal
Lpguess_old = Lpguess
LpLoop: do
NiterationStressLp = NiterationStressLp + 1
LpLoopLimit: if (NiterationStressLp > num%nStress) then
@ -1003,18 +1003,18 @@ logical function integrateStress(ipc,ip,el,timeFraction)
#endif
return
endif LpLoopLimit
!* calculate (elastic) 2nd Piola--Kirchhoff stress tensor and its tangent from constitutive law
B = math_I3 - dt*Lpguess
Fe = matmul(matmul(A,B), invFi_new)
call constitutive_SandItsTangents(S, dS_dFe, dS_dFi, &
Fe, Fi_new, ipc, ip, el) ! call constitutive law to calculate 2nd Piola-Kirchhoff stress and its derivative in unloaded configuration
!* calculate plastic velocity gradient and its tangent from constitutive law
call constitutive_LpAndItsTangents(Lp_constitutive, dLp_dS, dLp_dFi, &
S, Fi_new, ipc, ip, el)
#ifdef DEBUG
if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0 &
.and. ((el == debug_e .and. ip == debug_i .and. ipc == debug_g) &
@ -1032,7 +1032,7 @@ logical function integrateStress(ipc,ip,el,timeFraction)
aTolLp = max(num%rTol_crystalliteStress * max(norm2(Lpguess),norm2(Lp_constitutive)), & ! absolute tolerance from largest acceptable relative error
num%aTol_crystalliteStress) ! minimum lower cutoff
residuumLp = Lpguess - Lp_constitutive
if (any(IEEE_is_NaN(residuumLp))) then
#ifdef DEBUG
if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0) &
@ -1160,7 +1160,7 @@ logical function integrateStress(ipc,ip,el,timeFraction)
#endif
cycle LiLoop
endif
!* calculate Jacobian for correction term
if (mod(jacoCounterLi, num%iJacoLpresiduum) == 0) then
temp_33 = matmul(matmul(A,B),invFi_current)
@ -1197,11 +1197,11 @@ logical function integrateStress(ipc,ip,el,timeFraction)
#endif
return
endif
deltaLi = - math_9to33(work)
endif
jacoCounterLi = jacoCounterLi + 1
Liguess = Liguess + steplengthLi * deltaLi
#ifdef DEBUG
if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0 &
@ -1211,7 +1211,7 @@ logical function integrateStress(ipc,ip,el,timeFraction)
endif
#endif
enddo LiLoop
!* calculate new plastic and elastic deformation gradient
invFp_new = matmul(invFp_current,B)
invFp_new = invFp_new / math_det33(invFp_new)**(1.0_pReal/3.0_pReal) ! regularize
@ -1302,7 +1302,7 @@ subroutine integrateStateFPI
#endif
! store previousDotState and previousDotState2
!$OMP PARALLEL DO PRIVATE(p,c)
do e = FEsolving_execElem(1),FEsolving_execElem(2)
do i = FEsolving_execIP(1),FEsolving_execIP(2)
@ -1329,7 +1329,7 @@ subroutine integrateStateFPI
call update_dependentState
call update_stress(1.0_pReal)
call update_dotState(1.0_pReal)
!$OMP PARALLEL
!$OMP DO PRIVATE(sizeDotState,residuum_plastic,residuum_source,zeta,p,c)
do e = FEsolving_execElem(1),FEsolving_execElem(2)
@ -1342,7 +1342,7 @@ subroutine integrateStateFPI
zeta = damper(plasticState(p)%dotState (:,c), &
plasticState(p)%previousDotState (:,c), &
plasticState(p)%previousDotState2(:,c))
residuum_plastic(1:SizeDotState) = plasticState(p)%state (1:sizeDotState,c) &
- plasticState(p)%subState0(1:sizeDotState,c) &
- ( plasticState(p)%dotState (:,c) * zeta &
@ -1350,18 +1350,18 @@ subroutine integrateStateFPI
) * crystallite_subdt(g,i,e)
plasticState(p)%state(1:sizeDotState,c) = plasticState(p)%state(1:sizeDotState,c) &
- residuum_plastic(1:sizeDotState)
- residuum_plastic(1:sizeDotState)
plasticState(p)%dotState(:,c) = plasticState(p)%dotState(:,c) * zeta &
+ plasticState(p)%previousDotState(:,c) * (1.0_pReal - zeta)
crystallite_converged(g,i,e) = converged(residuum_plastic(1:sizeDotState), &
plasticState(p)%state(1:sizeDotState,c), &
plasticState(p)%aTolState(1:sizeDotState))
do s = 1, phase_Nsources(p)
sizeDotState = sourceState(p)%p(s)%sizeDotState
zeta = damper(sourceState(p)%p(s)%dotState (:,c), &
sourceState(p)%p(s)%previousDotState (:,c), &
sourceState(p)%p(s)%previousDotState2(:,c))
@ -1432,12 +1432,12 @@ subroutine integrateStateFPI
!> @brief calculate the damping for correction of state and dot state
!--------------------------------------------------------------------------------------------------
real(pReal) pure function damper(current,previous,previous2)
real(pReal), dimension(:), intent(in) ::&
current, previous, previous2
real(pReal) :: dot_prod12, dot_prod22
dot_prod12 = dot_product(current - previous, previous - previous2)
dot_prod22 = dot_product(previous - previous2, previous - previous2)
if ((dot_product(current,previous) < 0.0_pReal .or. dot_prod12 < 0.0_pReal) .and. dot_prod22 > 0.0_pReal) then
@ -1445,7 +1445,7 @@ subroutine integrateStateFPI
else
damper = 1.0_pReal
endif
end function damper
end subroutine integrateStateFPI
@ -1480,7 +1480,7 @@ subroutine integrateStateAdaptiveEuler
c, &
s, &
sizeDotState
! ToDo: MD: once all constitutives use allocate state, attach residuum arrays to the state in case of adaptive Euler
real(pReal), dimension(constitutive_plasticity_maxSizeDotState, &
homogenization_maxNgrains,discretization_nIP,discretization_nElem) :: &
@ -1501,14 +1501,14 @@ subroutine integrateStateAdaptiveEuler
if (crystallite_todo(g,i,e)) then
p = material_phaseAt(g,e); c = material_phaseMemberAt(g,i,e)
sizeDotState = plasticState(p)%sizeDotState
residuum_plastic(1:sizeDotState,g,i,e) = plasticState(p)%dotstate(1:sizeDotState,c) &
* (- 0.5_pReal * crystallite_subdt(g,i,e))
plasticState(p)%state(1:sizeDotState,c) = &
plasticState(p)%state(1:sizeDotState,c) + plasticState(p)%dotstate(1:sizeDotState,c) * crystallite_subdt(g,i,e) !ToDo: state, partitioned state?
do s = 1, phase_Nsources(p)
sizeDotState = sourceState(p)%p(s)%sizeDotState
residuum_source(1:sizeDotState,s,g,i,e) = sourceState(p)%p(s)%dotstate(1:sizeDotState,c) &
* (- 0.5_pReal * crystallite_subdt(g,i,e))
sourceState(p)%p(s)%state(1:sizeDotState,c) = &
@ -1530,17 +1530,17 @@ subroutine integrateStateAdaptiveEuler
if (crystallite_todo(g,i,e)) then
p = material_phaseAt(g,e); c = material_phaseMemberAt(g,i,e)
sizeDotState = plasticState(p)%sizeDotState
residuum_plastic(1:sizeDotState,g,i,e) = residuum_plastic(1:sizeDotState,g,i,e) &
+ 0.5_pReal * plasticState(p)%dotState(:,c) * crystallite_subdt(g,i,e)
crystallite_converged(g,i,e) = converged(residuum_plastic(1:sizeDotState,g,i,e), &
plasticState(p)%state(1:sizeDotState,c), &
plasticState(p)%aTolState(1:sizeDotState))
do s = 1, phase_Nsources(p)
sizeDotState = sourceState(p)%p(s)%sizeDotState
residuum_source(1:sizeDotState,s,g,i,e) = &
residuum_source(1:sizeDotState,s,g,i,e) + 0.5_pReal * sourceState(p)%p(s)%dotState(:,c) * crystallite_subdt(g,i,e)
@ -1549,13 +1549,13 @@ subroutine integrateStateAdaptiveEuler
sourceState(p)%p(s)%state(1:sizeDotState,c), &
sourceState(p)%p(s)%aTolState(1:sizeDotState))
enddo
endif
enddo; enddo; enddo
!$OMP END PARALLEL DO
if (any(plasticState(:)%nonlocal)) call nonlocalConvergenceCheck
end subroutine integrateStateAdaptiveEuler
@ -1720,23 +1720,23 @@ subroutine integrateStateRKCK45
do g = 1,homogenization_Ngrains(material_homogenizationAt(e))
if (crystallite_todo(g,i,e)) then
p = material_phaseAt(g,e); cc = material_phaseMemberAt(g,i,e)
sizeDotState = plasticState(p)%sizeDotState
plasticState(p)%RKCK45dotState(6,:,cc) = plasticState (p)%dotState(:,cc)
residuum_plastic(1:sizeDotState,g,i,e) = &
matmul(transpose(plasticState(p)%RKCK45dotState(1:6,1:sizeDotState,cc)),DB) & ! why transpose? Better to transpose constant DB
* crystallite_subdt(g,i,e)
plasticState(p)%dotState(:,cc) = &
matmul(transpose(plasticState(p)%RKCK45dotState(1:6,1:sizeDotState,cc)), B) ! why transpose? Better to transpose constant B
do s = 1, phase_Nsources(p)
sizeDotState = sourceState(p)%p(s)%sizeDotState
sourceState(p)%p(s)%RKCK45dotState(6,:,cc) = sourceState(p)%p(s)%dotState(:,cc)
residuum_source(1:sizeDotState,s,g,i,e) = &
matmul(transpose(sourceState(p)%p(s)%RKCK45dotState(1:6,1:sizeDotState,cc)),DB) &
* crystallite_subdt(g,i,e)
@ -1744,13 +1744,13 @@ subroutine integrateStateRKCK45
sourceState(p)%p(s)%dotState(:,cc) = &
matmul(transpose(sourceState(p)%p(s)%RKCK45dotState(1:6,1:sizeDotState,cc)),B)
enddo
endif
enddo; enddo; enddo
!$OMP END PARALLEL DO
call update_state(1.0_pReal)
! --- relative residui and state convergence ---
!$OMP PARALLEL DO PRIVATE(sizeDotState,p,cc)
@ -1759,16 +1759,16 @@ subroutine integrateStateRKCK45
do g = 1,homogenization_Ngrains(material_homogenizationAt(e))
if (crystallite_todo(g,i,e)) then
p = material_phaseAt(g,e); cc = material_phaseMemberAt(g,i,e)
sizeDotState = plasticState(p)%sizeDotState
crystallite_todo(g,i,e) = converged(residuum_plastic(1:sizeDotState,g,i,e), &
plasticState(p)%state(1:sizeDotState,cc), &
plasticState(p)%aTolState(1:sizeDotState))
do s = 1, phase_Nsources(p)
sizeDotState = sourceState(p)%p(s)%sizeDotState
crystallite_todo(g,i,e) = &
crystallite_todo(g,i,e) .and. converged(residuum_source(1:sizeDotState,s,g,i,e), &
sourceState(p)%p(s)%state(1:sizeDotState,cc), &
@ -1783,7 +1783,7 @@ subroutine integrateStateRKCK45
call update_stress(1.0_pReal)
call setConvergenceFlag
if (any(plasticState(:)%nonlocal)) call nonlocalConvergenceCheck
end subroutine integrateStateRKCK45
@ -1792,7 +1792,7 @@ end subroutine integrateStateRKCK45
!> @detail one non-converged nonlocal sets all other nonlocals to non-converged to trigger cut back
!--------------------------------------------------------------------------------------------------
subroutine nonlocalConvergenceCheck
if (any(.not. crystallite_converged .and. .not. crystallite_localPlasticity)) & ! any non-local not yet converged (or broken)...
where( .not. crystallite_localPlasticity) crystallite_converged = .false.
@ -1810,7 +1810,7 @@ subroutine setConvergenceFlag
e, & !< element index in element loop
i, & !< integration point index in ip loop
g !< grain index in grain loop
!OMP DO PARALLEL PRIVATE
do e = FEsolving_execElem(1),FEsolving_execElem(2)
do i = FEsolving_execIP(1),FEsolving_execIP(2)
@ -1826,7 +1826,7 @@ end subroutine setConvergenceFlag
!> @brief determines whether a point is converged
!--------------------------------------------------------------------------------------------------
logical pure function converged(residuum,state,aTol)
real(pReal), intent(in), dimension(:) ::&
residuum, state, aTol
real(pReal) :: &
@ -1949,11 +1949,11 @@ subroutine update_dotState(timeFraction)
g, & !< grain index in grain loop
p, &
c, &
s
s
logical :: &
NaN, &
nonlocalStop
nonlocalStop = .false.
!$OMP PARALLEL DO PRIVATE (p,c,NaN)
@ -1963,9 +1963,9 @@ subroutine update_dotState(timeFraction)
!$OMP FLUSH(nonlocalStop)
if ((crystallite_todo(g,i,e) .and. .not. crystallite_converged(g,i,e)) .and. .not. nonlocalStop) then
call constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), &
crystallite_Fe, &
crystallite_partionedF0, &
crystallite_Fi(1:3,1:3,g,i,e), &
crystallite_Fp, &
crystallite_partionedFp0, &
crystallite_subdt(g,i,e)*timeFraction, g,i,e)
p = material_phaseAt(g,e); c = material_phaseMemberAt(g,i,e)
NaN = any(IEEE_is_NaN(plasticState(p)%dotState(:,c)))
@ -1980,7 +1980,7 @@ subroutine update_dotState(timeFraction)
enddo; enddo; enddo
!$OMP END PARALLEL DO
if (nonlocalStop) crystallite_todo = crystallite_todo .and. crystallite_localPlasticity
if (nonlocalStop) crystallite_todo = crystallite_todo .and. crystallite_localPlasticity
end subroutine update_DotState
@ -1995,11 +1995,11 @@ subroutine update_deltaState
mySize, &
myOffset, &
c, &
s
s
logical :: &
NaN, &
nonlocalStop
nonlocalStop = .false.
!$OMP PARALLEL DO PRIVATE(p,c,myOffset,mySize,NaN)
@ -2016,23 +2016,23 @@ subroutine update_deltaState
myOffset = plasticState(p)%offsetDeltaState
mySize = plasticState(p)%sizeDeltaState
NaN = any(IEEE_is_NaN(plasticState(p)%deltaState(1:mySize,c)))
if (.not. NaN) then
plasticState(p)%state(myOffset + 1: myOffset + mySize,c) = &
plasticState(p)%state(myOffset + 1: myOffset + mySize,c) + plasticState(p)%deltaState(1:mySize,c)
do s = 1, phase_Nsources(p)
myOffset = sourceState(p)%p(s)%offsetDeltaState
mySize = sourceState(p)%p(s)%sizeDeltaState
NaN = NaN .or. any(IEEE_is_NaN(sourceState(p)%p(s)%deltaState(1:mySize,c)))
if (.not. NaN) then
sourceState(p)%p(s)%state(myOffset + 1:myOffset + mySize,c) = &
sourceState(p)%p(s)%state(myOffset + 1:myOffset + mySize,c) + sourceState(p)%p(s)%deltaState(1:mySize,c)
endif
enddo
endif
crystallite_todo(g,i,e) = .not. NaN
if (.not. crystallite_todo(g,i,e)) then ! if state jump fails, then convergence is broken
crystallite_converged(g,i,e) = .false.
@ -2042,7 +2042,7 @@ subroutine update_deltaState
enddo; enddo; enddo
!$OMP END PARALLEL DO
if (nonlocalStop) crystallite_todo = crystallite_todo .and. crystallite_localPlasticity
end subroutine update_deltaState