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

Merge remote-tracking branch 'remotes/origin/explicitNonlocal' into development

This commit is contained in:
Franz Roters 2020-02-11 10:32:53 +01:00
parent ee11c7ea4b 61810e6a08
commit 81ae66860a
4 changed files with 587 additions and 556 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
PRIVATE

@ -1 +1 @@
Subproject commit 64432754ce3c590c882cf4987695539cee524ee8 Subproject commit 9dc7065beedf2a097aa60a656bbfa52e55d7147c

204
src/constitutive.f90
View File

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

699
src/constitutive_plastic_nonlocal.f90
View File

File diff suppressed because it is too large Load Diff

238
src/crystallite.f90
View File

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