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

Merge commit 'v3.0.0-alpha2-340-g014378d49'

This commit is contained in:
Test User 2021-01-22 18:21:50 +01:00
parent c3ca61c3af 014378d49c
commit 7cf15c0aab
16 changed files with 296 additions and 207 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
PRIVATE

@ -1 +1 @@
Subproject commit 72661176e9055cf6ec106bb61c50482c5bc78de8 Subproject commit 8b9836aaf5798727d96a316a2eb27df03bbd2d07

1
src/commercialFEM_fileList.f90
View File

@ -52,4 +52,5 @@
#include "homogenization_mech_isostrain.f90" #include "homogenization_mech_isostrain.f90"
#include "homogenization_mech_RGC.f90" #include "homogenization_mech_RGC.f90"
#include "homogenization_thermal.f90" #include "homogenization_thermal.f90"
#include "homogenization_damage.f90"
#include "CPFEM.f90" #include "CPFEM.f90"

54
src/constitutive.f90
View File

@ -188,6 +188,11 @@ module constitutive
real(pReal), dimension(3,3) :: P real(pReal), dimension(3,3) :: P
end function constitutive_mech_getP end function constitutive_mech_getP
module function constitutive_damage_get_phi(co,ip,el) result(phi)
integer, intent(in) :: co, ip, el
real(pReal) :: phi
end function constitutive_damage_get_phi
module function thermal_T(ph,me) result(T) module function thermal_T(ph,me) result(T)
integer, intent(in) :: ph,me integer, intent(in) :: ph,me
real(pReal) :: T real(pReal) :: T
@ -203,6 +208,11 @@ module constitutive
real(pReal), intent(in) :: T real(pReal), intent(in) :: T
integer, intent(in) :: co, ce integer, intent(in) :: co, ce
end subroutine constitutive_thermal_setT end subroutine constitutive_thermal_setT
module subroutine constitutive_damage_set_phi(phi,co,ce)
real(pReal), intent(in) :: phi
integer, intent(in) :: co, ce
end subroutine constitutive_damage_set_phi
! == cleaned:end =================================================================================== ! == cleaned:end ===================================================================================
@ -229,39 +239,14 @@ module constitutive
logical :: converged_ logical :: converged_
end function crystallite_stress end function crystallite_stress
module function constitutive_homogenizedC(co,ip,el) result(C) module function constitutive_homogenizedC(ph,me) result(C)
integer, intent(in) :: co, ip, el integer, intent(in) :: ph, me
real(pReal), dimension(6,6) :: C real(pReal), dimension(6,6) :: C
end function constitutive_homogenizedC end function constitutive_homogenizedC
module subroutine source_damage_anisoBrittle_dotState(S, co, ip, el)
integer, intent(in) :: &
co, & !< component-ID of integration point
ip, & !< integration point
el !< element
real(pReal), intent(in), dimension(3,3) :: &
S
end subroutine source_damage_anisoBrittle_dotState
module subroutine source_damage_anisoDuctile_dotState(co, ip, el)
integer, intent(in) :: &
co, & !< component-ID of integration point
ip, & !< integration point
el !< element
end subroutine source_damage_anisoDuctile_dotState
module subroutine source_damage_isoDuctile_dotState(co, ip, el)
integer, intent(in) :: &
co, & !< component-ID of integration point
ip, & !< integration point
el !< element
end subroutine source_damage_isoDuctile_dotState
module subroutine source_thermal_externalheat_dotState(phase, of)
integer, intent(in) :: &
phase, &
of
end subroutine source_thermal_externalheat_dotState
module subroutine constitutive_damage_getRateAndItsTangents(phiDot, dPhiDot_dPhi, phi, ip, el) module subroutine constitutive_damage_getRateAndItsTangents(phiDot, dPhiDot_dPhi, phi, ip, el)
integer, intent(in) :: & integer, intent(in) :: &
@ -343,19 +328,6 @@ module constitutive
dLi_dTstar !< derivative of Li with respect to Tstar (4th-order tensor defined to be zero) dLi_dTstar !< derivative of Li with respect to Tstar (4th-order tensor defined to be zero)
end subroutine kinematics_thermal_expansion_LiAndItsTangent end subroutine kinematics_thermal_expansion_LiAndItsTangent
module subroutine source_damage_isoBrittle_deltaState(C, Fe, co, ip, el)
integer, intent(in) :: &
co, & !< component-ID of integration point
ip, & !< integration point
el !< element
real(pReal), intent(in), dimension(3,3) :: &
Fe
real(pReal), intent(in), dimension(6,6) :: &
C
end subroutine source_damage_isoBrittle_deltaState
module subroutine constitutive_plastic_dependentState(co,ip,el) module subroutine constitutive_plastic_dependentState(co,ip,el)
integer, intent(in) :: & integer, intent(in) :: &
co, & !< component-ID of integration point co, & !< component-ID of integration point
@ -391,6 +363,8 @@ module constitutive
constitutive_restartRead, & constitutive_restartRead, &
integrateDamageState, & integrateDamageState, &
constitutive_thermal_setT, & constitutive_thermal_setT, &
constitutive_damage_set_phi, &
constitutive_damage_get_phi, &
constitutive_mech_getP, & constitutive_mech_getP, &
constitutive_mech_setF, & constitutive_mech_setF, &
constitutive_mech_getF, & constitutive_mech_getF, &

103
src/constitutive_damage.f90
View File

@ -10,9 +10,16 @@ submodule(constitutive) constitutive_damage
DAMAGE_ANISODUCTILE_ID DAMAGE_ANISODUCTILE_ID
end enum end enum
type :: tDataContainer
real(pReal), dimension(:), allocatable :: phi, d_phi_d_dot_phi
end type tDataContainer
integer(kind(DAMAGE_UNDEFINED_ID)), dimension(:,:), allocatable :: & integer(kind(DAMAGE_UNDEFINED_ID)), dimension(:,:), allocatable :: &
phase_source !< active sources mechanisms of each phase phase_source !< active sources mechanisms of each phase
type(tDataContainer), dimension(:), allocatable :: current
interface interface
module function source_damage_anisoBrittle_init(source_length) result(mySources) module function source_damage_anisoBrittle_init(source_length) result(mySources)
@ -45,6 +52,38 @@ submodule(constitutive) constitutive_damage
logical, dimension(:,:), allocatable :: myKinematics logical, dimension(:,:), allocatable :: myKinematics
end function kinematics_slipplane_opening_init end function kinematics_slipplane_opening_init
module subroutine source_damage_isoBrittle_deltaState(C, Fe, ph, me)
integer, intent(in) :: ph,me
real(pReal), intent(in), dimension(3,3) :: &
Fe
real(pReal), intent(in), dimension(6,6) :: &
C
end subroutine source_damage_isoBrittle_deltaState
module subroutine source_damage_anisoBrittle_dotState(S, co, ip, el)
integer, intent(in) :: &
co, & !< component-ID of integration point
ip, & !< integration point
el !< element
real(pReal), intent(in), dimension(3,3) :: &
S
end subroutine source_damage_anisoBrittle_dotState
module subroutine source_damage_anisoDuctile_dotState(co, ip, el)
integer, intent(in) :: &
co, & !< component-ID of integration point
ip, & !< integration point
el !< element
end subroutine source_damage_anisoDuctile_dotState
module subroutine source_damage_isoDuctile_dotState(co, ip, el)
integer, intent(in) :: &
co, & !< component-ID of integration point
ip, & !< integration point
el !< element
end subroutine source_damage_isoDuctile_dotState
module subroutine source_damage_anisobrittle_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent) module subroutine source_damage_anisobrittle_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent)
integer, intent(in) :: & integer, intent(in) :: &
@ -120,7 +159,8 @@ contains
module subroutine damage_init module subroutine damage_init
integer :: & integer :: &
ph !< counter in phase loop ph, & !< counter in phase loop
Nconstituents
class(tNode), pointer :: & class(tNode), pointer :: &
phases, & phases, &
phase, & phase, &
@ -129,10 +169,18 @@ module subroutine damage_init
phases => config_material%get('phase') phases => config_material%get('phase')
allocate(current(phases%length))
allocate(damageState (phases%length)) allocate(damageState (phases%length))
allocate(phase_Nsources(phases%length),source = 0) ! same for kinematics allocate(phase_Nsources(phases%length),source = 0) ! same for kinematics
do ph = 1,phases%length do ph = 1,phases%length
Nconstituents = count(material_phaseAt == ph) * discretization_nIPs
allocate(current(ph)%phi(Nconstituents),source=1.0_pReal)
allocate(current(ph)%d_phi_d_dot_phi(Nconstituents),source=0.0_pReal)
phase => phases%get(ph) phase => phases%get(ph)
sources => phase%get('source',defaultVal=emptyList) sources => phase%get('source',defaultVal=emptyList)
phase_Nsources(ph) = sources%length phase_Nsources(ph) = sources%length
@ -295,7 +343,7 @@ module function integrateDamageState(dt,co,ip,el) result(broken)
enddo enddo
if(converged_) then if(converged_) then
broken = constitutive_damage_deltaState(mech_F_e(ph,me),co,ip,el,ph,me) broken = constitutive_damage_deltaState(mech_F_e(ph,me),ph,me)
exit iteration exit iteration
endif endif
@ -368,14 +416,14 @@ end subroutine damage_results
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @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
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function constitutive_damage_collectDotState(co,ip,el,ph,of) result(broken) function constitutive_damage_collectDotState(co,ip,el,ph,me) result(broken)
integer, intent(in) :: & integer, intent(in) :: &
co, & !< component-ID of integration point co, & !< component-ID me integration point
ip, & !< integration point ip, & !< integration point
el, & !< element el, & !< element
ph, & ph, &
of me
integer :: & integer :: &
so !< counter in source loop so !< counter in source loop
logical :: broken logical :: broken
@ -394,12 +442,11 @@ function constitutive_damage_collectDotState(co,ip,el,ph,of) result(broken)
call source_damage_anisoDuctile_dotState(co, ip, el) call source_damage_anisoDuctile_dotState(co, ip, el)
case (DAMAGE_ANISOBRITTLE_ID) sourceType case (DAMAGE_ANISOBRITTLE_ID) sourceType
call source_damage_anisoBrittle_dotState(mech_S(material_phaseAt(co,el),material_phaseMemberAt(co,ip,el)),& call source_damage_anisoBrittle_dotState(mech_S(ph,me),co, ip, el) ! correct stress?
co, ip, el) ! correct stress?
end select sourceType end select sourceType
broken = broken .or. any(IEEE_is_NaN(damageState(ph)%p(so)%dotState(:,of))) broken = broken .or. any(IEEE_is_NaN(damageState(ph)%p(so)%dotState(:,me)))
enddo SourceLoop enddo SourceLoop
@ -411,14 +458,11 @@ end function constitutive_damage_collectDotState
!> @brief for constitutive models having an instantaneous change of state !> @brief for constitutive models having an instantaneous change of state
!> will return false if delta state is not needed/supported by the constitutive model !> will return false if delta state is not needed/supported by the constitutive model
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function constitutive_damage_deltaState(Fe, co, ip, el, ph, of) result(broken) function constitutive_damage_deltaState(Fe, ph, me) result(broken)
integer, intent(in) :: & integer, intent(in) :: &
co, & !< component-ID of integration point
ip, & !< integration point
el, & !< element
ph, & ph, &
of me
real(pReal), intent(in), dimension(3,3) :: & real(pReal), intent(in), dimension(3,3) :: &
Fe !< elastic deformation gradient Fe !< elastic deformation gradient
integer :: & integer :: &
@ -436,14 +480,13 @@ function constitutive_damage_deltaState(Fe, co, ip, el, ph, of) result(broken)
sourceType: select case (phase_source(so,ph)) sourceType: select case (phase_source(so,ph))
case (DAMAGE_ISOBRITTLE_ID) sourceType case (DAMAGE_ISOBRITTLE_ID) sourceType
call source_damage_isoBrittle_deltaState (constitutive_homogenizedC(co,ip,el), Fe, & call source_damage_isoBrittle_deltaState(constitutive_homogenizedC(ph,me), Fe, ph,me)
co, ip, el) broken = any(IEEE_is_NaN(damageState(ph)%p(so)%deltaState(:,me)))
broken = any(IEEE_is_NaN(damageState(ph)%p(so)%deltaState(:,of)))
if(.not. broken) then if(.not. broken) then
myOffset = damageState(ph)%p(so)%offsetDeltaState myOffset = damageState(ph)%p(so)%offsetDeltaState
mySize = damageState(ph)%p(so)%sizeDeltaState mySize = damageState(ph)%p(so)%sizeDeltaState
damageState(ph)%p(so)%state(myOffset + 1: myOffset + mySize,of) = & damageState(ph)%p(so)%state(myOffset + 1: myOffset + mySize,me) = &
damageState(ph)%p(so)%state(myOffset + 1: myOffset + mySize,of) + damageState(ph)%p(so)%deltaState(1:mySize,of) damageState(ph)%p(so)%state(myOffset + 1: myOffset + mySize,me) + damageState(ph)%p(so)%deltaState(1:mySize,me)
endif endif
end select sourceType end select sourceType
@ -484,4 +527,28 @@ function source_active(source_label,src_length) result(active_source)
end function source_active end function source_active
!----------------------------------------------------------------------------------------------
!< @brief Set damage parameter
!----------------------------------------------------------------------------------------------
module subroutine constitutive_damage_set_phi(phi,co,ce)
real(pReal), intent(in) :: phi
integer, intent(in) :: ce, co
current(material_phaseAt2(co,ce))%phi(material_phaseMemberAt2(co,ce)) = phi
end subroutine constitutive_damage_set_phi
module function constitutive_damage_get_phi(co,ip,el) result(phi)
integer, intent(in) :: co, ip, el
real(pReal) :: phi
phi = current(material_phaseAt(co,el))%phi(material_phaseMemberAt(co,ip,el))
end function constitutive_damage_get_phi
end submodule constitutive_damage end submodule constitutive_damage

31
src/constitutive_mech.f90
View File

@ -307,13 +307,9 @@ submodule(constitutive) constitutive_mech
character(len=*), intent(in) :: group character(len=*), intent(in) :: group
end subroutine plastic_nonlocal_results end subroutine plastic_nonlocal_results
module function plastic_dislotwin_homogenizedC(co,ip,el) result(homogenizedC) module function plastic_dislotwin_homogenizedC(ph,me) result(homogenizedC)
real(pReal), dimension(6,6) :: & real(pReal), dimension(6,6) :: homogenizedC
homogenizedC integer, intent(in) :: ph,me
integer, intent(in) :: &
co, & !< component-ID of integration point
ip, & !< integration point
el !< element
end function plastic_dislotwin_homogenizedC end function plastic_dislotwin_homogenizedC
@ -559,16 +555,16 @@ subroutine constitutive_hooke_SandItsTangents(S, dS_dFe, dS_dFi, &
real(pReal), intent(out), dimension(3,3,3,3) :: & real(pReal), intent(out), dimension(3,3,3,3) :: &
dS_dFe, & !< derivative of 2nd P-K stress with respect to elastic deformation gradient dS_dFe, & !< derivative of 2nd P-K stress with respect to elastic deformation gradient
dS_dFi !< derivative of 2nd P-K stress with respect to intermediate deformation gradient dS_dFi !< derivative of 2nd P-K stress with respect to intermediate deformation gradient
real(pReal), dimension(3,3) :: E real(pReal), dimension(3,3) :: E
real(pReal), dimension(3,3,3,3) :: C real(pReal), dimension(3,3,3,3) :: C
integer :: & integer :: &
ho, & !< homogenization ho, & !< homogenization
d !< counter in degradation loop d, & !< counter in degradation loop
integer :: & i, j, ph, me
i, j
ho = material_homogenizationAt(el) ho = material_homogenizationAt(el)
C = math_66toSym3333(constitutive_homogenizedC(co,ip,el)) C = math_66toSym3333(constitutive_homogenizedC(material_phaseAt(co,el),material_phaseMemberAt(co,ip,el)))
DegradationLoop: do d = 1, phase_NstiffnessDegradations(material_phaseAt(co,el)) DegradationLoop: do d = 1, phase_NstiffnessDegradations(material_phaseAt(co,el))
degradationType: select case(phase_stiffnessDegradation(d,material_phaseAt(co,el))) degradationType: select case(phase_stiffnessDegradation(d,material_phaseAt(co,el)))
@ -1535,19 +1531,16 @@ end subroutine mech_forward
!> @brief returns the homogenize elasticity matrix !> @brief returns the homogenize elasticity matrix
!> ToDo: homogenizedC66 would be more consistent !> ToDo: homogenizedC66 would be more consistent
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module function constitutive_homogenizedC(co,ip,el) result(C) module function constitutive_homogenizedC(ph,me) result(C)
real(pReal), dimension(6,6) :: C real(pReal), dimension(6,6) :: C
integer, intent(in) :: & integer, intent(in) :: ph, me
co, & !< component-ID of integration point
ip, & !< integration point
el !< element
plasticityType: select case (phase_plasticity(material_phaseAt(co,el))) plasticityType: select case (phase_plasticity(ph))
case (PLASTICITY_DISLOTWIN_ID) plasticityType case (PLASTICITY_DISLOTWIN_ID) plasticityType
C = plastic_dislotwin_homogenizedC(co,ip,el) C = plastic_dislotwin_homogenizedC(ph,me)
case default plasticityType case default plasticityType
C = lattice_C66(1:6,1:6,material_phaseAt(co,el)) C = lattice_C66(1:6,1:6,ph)
end select plasticityType end select plasticityType
end function constitutive_homogenizedC end function constitutive_homogenizedC

27
src/constitutive_plastic_dislotwin.f90
View File

@ -150,7 +150,7 @@ module function plastic_dislotwin_init() result(myPlasticity)
Ninstances = count(myPlasticity) Ninstances = count(myPlasticity)
print'(a,i2)', ' # instances: ',Ninstances; flush(IO_STDOUT) print'(a,i2)', ' # instances: ',Ninstances; flush(IO_STDOUT)
if(Ninstances == 0) return if(Ninstances == 0) return
print*, 'Ma and Roters, Acta Materialia 52(12):3603–3612, 2004' print*, 'Ma and Roters, Acta Materialia 52(12):3603–3612, 2004'
print*, 'https://doi.org/10.1016/j.actamat.2004.04.012'//IO_EOL print*, 'https://doi.org/10.1016/j.actamat.2004.04.012'//IO_EOL
@ -485,35 +485,32 @@ end function plastic_dislotwin_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Return the homogenized elasticity matrix. !> @brief Return the homogenized elasticity matrix.
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module function plastic_dislotwin_homogenizedC(co,ip,el) result(homogenizedC) module function plastic_dislotwin_homogenizedC(ph,me) result(homogenizedC)
integer, intent(in) :: &
ph, me
real(pReal), dimension(6,6) :: & real(pReal), dimension(6,6) :: &
homogenizedC homogenizedC
integer, intent(in) :: &
co, & !< component-ID of integration point
ip, & !< integration point
el !< element
integer :: i, & integer :: i
of
real(pReal) :: f_unrotated real(pReal) :: f_unrotated
of = material_phasememberAt(co,ip,el)
associate(prm => param(phase_plasticityInstance(material_phaseAt(co,el))),& associate(prm => param(phase_plasticityInstance(ph)),&
stt => state(phase_plasticityInstance(material_phaseAT(co,el)))) stt => state(phase_plasticityInstance(ph)))
f_unrotated = 1.0_pReal & f_unrotated = 1.0_pReal &
- sum(stt%f_tw(1:prm%sum_N_tw,of)) & - sum(stt%f_tw(1:prm%sum_N_tw,me)) &
- sum(stt%f_tr(1:prm%sum_N_tr,of)) - sum(stt%f_tr(1:prm%sum_N_tr,me))
homogenizedC = f_unrotated * prm%C66 homogenizedC = f_unrotated * prm%C66
do i=1,prm%sum_N_tw do i=1,prm%sum_N_tw
homogenizedC = homogenizedC & homogenizedC = homogenizedC &
+ stt%f_tw(i,of)*prm%C66_tw(1:6,1:6,i) + stt%f_tw(i,me)*prm%C66_tw(1:6,1:6,i)
enddo enddo
do i=1,prm%sum_N_tr do i=1,prm%sum_N_tr
homogenizedC = homogenizedC & homogenizedC = homogenizedC &
+ stt%f_tr(i,of)*prm%C66_tr(1:6,1:6,i) + stt%f_tr(i,me)*prm%C66_tr(1:6,1:6,i)
enddo enddo
end associate end associate

13
src/constitutive_thermal.f90
View File

@ -36,6 +36,13 @@ submodule(constitutive) constitutive_thermal
end function kinematics_thermal_expansion_init end function kinematics_thermal_expansion_init
module subroutine source_thermal_externalheat_dotState(ph, me)
integer, intent(in) :: &
ph, &
me
end subroutine source_thermal_externalheat_dotState
module subroutine thermal_dissipation_getRate(TDot, Tstar,Lp,phase) module subroutine thermal_dissipation_getRate(TDot, Tstar,Lp,phase)
integer, intent(in) :: & integer, intent(in) :: &
phase !< phase ID of element phase !< phase ID of element
@ -47,10 +54,10 @@ submodule(constitutive) constitutive_thermal
TDot TDot
end subroutine thermal_dissipation_getRate end subroutine thermal_dissipation_getRate
module subroutine thermal_externalheat_getRate(TDot, phase,of) module subroutine thermal_externalheat_getRate(TDot, ph,me)
integer, intent(in) :: & integer, intent(in) :: &
phase, & ph, &
of me
real(pReal), intent(out) :: & real(pReal), intent(out) :: &
TDot TDot
end subroutine thermal_externalheat_getRate end subroutine thermal_externalheat_getRate

24
src/constitutive_thermal_externalheat.f90
View File

@ -85,18 +85,18 @@ end function source_thermal_externalheat_init
!> @brief rate of change of state !> @brief rate of change of state
!> @details state only contains current time to linearly interpolate given heat powers !> @details state only contains current time to linearly interpolate given heat powers
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module subroutine source_thermal_externalheat_dotState(phase, of) module subroutine source_thermal_externalheat_dotState(ph, me)
integer, intent(in) :: & integer, intent(in) :: &
phase, & ph, &
of me
integer :: & integer :: &
sourceOffset sourceOffset
sourceOffset = source_thermal_externalheat_offset(phase) sourceOffset = source_thermal_externalheat_offset(ph)
thermalState(phase)%p(sourceOffset)%dotState(1,of) = 1.0_pReal ! state is current time thermalState(ph)%p(sourceOffset)%dotState(1,me) = 1.0_pReal ! state is current time
end subroutine source_thermal_externalheat_dotState end subroutine source_thermal_externalheat_dotState
@ -104,11 +104,11 @@ end subroutine source_thermal_externalheat_dotState
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief returns local heat generation rate !> @brief returns local heat generation rate
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module subroutine thermal_externalheat_getRate(TDot, phase, of) module subroutine thermal_externalheat_getRate(TDot, ph, me)
integer, intent(in) :: & integer, intent(in) :: &
phase, & ph, &
of me
real(pReal), intent(out) :: & real(pReal), intent(out) :: &
TDot TDot
@ -117,18 +117,18 @@ module subroutine thermal_externalheat_getRate(TDot, phase, of)
real(pReal) :: & real(pReal) :: &
frac_time frac_time
sourceOffset = source_thermal_externalheat_offset(phase) sourceOffset = source_thermal_externalheat_offset(ph)
associate(prm => param(source_thermal_externalheat_instance(phase))) associate(prm => param(source_thermal_externalheat_instance(ph)))
do interval = 1, prm%nIntervals ! scan through all rate segments do interval = 1, prm%nIntervals ! scan through all rate segments
frac_time = (thermalState(phase)%p(sourceOffset)%state(1,of) - prm%t_n(interval)) & frac_time = (thermalState(ph)%p(sourceOffset)%state(1,me) - prm%t_n(interval)) &
/ (prm%t_n(interval+1) - prm%t_n(interval)) ! fractional time within segment / (prm%t_n(interval+1) - prm%t_n(interval)) ! fractional time within segment
if ( (frac_time < 0.0_pReal .and. interval == 1) & if ( (frac_time < 0.0_pReal .and. interval == 1) &
.or. (frac_time >= 1.0_pReal .and. interval == prm%nIntervals) & .or. (frac_time >= 1.0_pReal .and. interval == prm%nIntervals) &
.or. (frac_time >= 0.0_pReal .and. frac_time < 1.0_pReal) ) & .or. (frac_time >= 0.0_pReal .and. frac_time < 1.0_pReal) ) &
TDot = prm%f_T(interval ) * (1.0_pReal - frac_time) + & TDot = prm%f_T(interval ) * (1.0_pReal - frac_time) + &
prm%f_T(interval+1) * frac_time ! interpolate heat rate between segment boundaries... prm%f_T(interval+1) * frac_time ! interpolate heat rate between segment boundaries...
! ...or extrapolate if outside of bounds ! ...or extrapolate if outside me bounds
enddo enddo
end associate end associate

69
src/grid/grid_damage_spectral.f90
View File

@ -17,8 +17,9 @@ module grid_damage_spectral
use discretization_grid use discretization_grid
use damage_nonlocal use damage_nonlocal
use YAML_types use YAML_types
use homogenization
use config use config
implicit none implicit none
private private
@ -43,13 +44,13 @@ module grid_damage_spectral
phi_current, & !< field of current damage phi_current, & !< field of current damage
phi_lastInc, & !< field of previous damage phi_lastInc, & !< field of previous damage
phi_stagInc !< field of staggered damage phi_stagInc !< field of staggered damage
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! reference diffusion tensor, mobility etc. ! reference diffusion tensor, mobility etc.
integer :: totalIter = 0 !< total iteration in current increment integer :: totalIter = 0 !< total iteration in current increment
real(pReal), dimension(3,3) :: K_ref real(pReal), dimension(3,3) :: K_ref
real(pReal) :: mu_ref real(pReal) :: mu_ref
public :: & public :: &
grid_damage_spectral_init, & grid_damage_spectral_init, &
grid_damage_spectral_solution, & grid_damage_spectral_solution, &
@ -62,8 +63,8 @@ contains
! ToDo: Restart not implemented ! ToDo: Restart not implemented
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine grid_damage_spectral_init subroutine grid_damage_spectral_init
PetscInt, dimension(0:worldsize-1) :: localK PetscInt, dimension(0:worldsize-1) :: localK
DM :: damage_grid DM :: damage_grid
Vec :: uBound, lBound Vec :: uBound, lBound
PetscErrorCode :: ierr PetscErrorCode :: ierr
@ -72,22 +73,22 @@ subroutine grid_damage_spectral_init
num_generic num_generic
character(len=pStringLen) :: & character(len=pStringLen) :: &
snes_type snes_type
print'(/,a)', ' <<<+- grid_spectral_damage init -+>>>' print'(/,a)', ' <<<+- grid_spectral_damage init -+>>>'
print*, 'Shanthraj et al., Handbook of Mechanics of Materials, 2019' print*, 'Shanthraj et al., Handbook of Mechanics of Materials, 2019'
print*, 'https://doi.org/10.1007/978-981-10-6855-3_80' print*, 'https://doi.org/10.1007/978-981-10-6855-3_80'
!------------------------------------------------------------------------------------------------- !-------------------------------------------------------------------------------------------------
! read numerical parameters and do sanity checks ! read numerical parameters and do sanity checks
num_grid => config_numerics%get('grid',defaultVal=emptyDict) num_grid => config_numerics%get('grid',defaultVal=emptyDict)
num%itmax = num_grid%get_asInt ('itmax',defaultVal=250) num%itmax = num_grid%get_asInt ('itmax',defaultVal=250)
num%eps_damage_atol = num_grid%get_asFloat ('eps_damage_atol',defaultVal=1.0e-2_pReal) num%eps_damage_atol = num_grid%get_asFloat ('eps_damage_atol',defaultVal=1.0e-2_pReal)
num%eps_damage_rtol = num_grid%get_asFloat ('eps_damage_rtol',defaultVal=1.0e-6_pReal) num%eps_damage_rtol = num_grid%get_asFloat ('eps_damage_rtol',defaultVal=1.0e-6_pReal)
num_generic => config_numerics%get('generic',defaultVal=emptyDict) num_generic => config_numerics%get('generic',defaultVal=emptyDict)
num%residualStiffness = num_generic%get_asFloat('residualStiffness', defaultVal=1.0e-6_pReal) num%residualStiffness = num_generic%get_asFloat('residualStiffness', defaultVal=1.0e-6_pReal)
if (num%residualStiffness < 0.0_pReal) call IO_error(301,ext_msg='residualStiffness') if (num%residualStiffness < 0.0_pReal) call IO_error(301,ext_msg='residualStiffness')
if (num%itmax <= 1) call IO_error(301,ext_msg='itmax') if (num%itmax <= 1) call IO_error(301,ext_msg='itmax')
if (num%eps_damage_atol <= 0.0_pReal) call IO_error(301,ext_msg='eps_damage_atol') if (num%eps_damage_atol <= 0.0_pReal) call IO_error(301,ext_msg='eps_damage_atol')
@ -104,7 +105,7 @@ subroutine grid_damage_spectral_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! initialize solver specific parts of PETSc ! initialize solver specific parts of PETSc
call SNESCreate(PETSC_COMM_WORLD,damage_snes,ierr); CHKERRQ(ierr) call SNESCreate(PETSC_COMM_WORLD,damage_snes,ierr); CHKERRQ(ierr)
call SNESSetOptionsPrefix(damage_snes,'damage_',ierr);CHKERRQ(ierr) call SNESSetOptionsPrefix(damage_snes,'damage_',ierr);CHKERRQ(ierr)
localK = 0 localK = 0
localK(worldrank) = grid3 localK(worldrank) = grid3
call MPI_Allreduce(MPI_IN_PLACE,localK,worldsize,MPI_INTEGER,MPI_SUM,PETSC_COMM_WORLD,ierr) call MPI_Allreduce(MPI_IN_PLACE,localK,worldsize,MPI_INTEGER,MPI_SUM,PETSC_COMM_WORLD,ierr)
@ -122,7 +123,7 @@ subroutine grid_damage_spectral_init
call DMsetUp(damage_grid,ierr); CHKERRQ(ierr) call DMsetUp(damage_grid,ierr); CHKERRQ(ierr)
call DMCreateGlobalVector(damage_grid,solution_vec,ierr); CHKERRQ(ierr) ! global solution vector (grid x 1, i.e. every def grad tensor) call DMCreateGlobalVector(damage_grid,solution_vec,ierr); CHKERRQ(ierr) ! global solution vector (grid x 1, i.e. every def grad tensor)
call DMDASNESSetFunctionLocal(damage_grid,INSERT_VALUES,formResidual,PETSC_NULL_SNES,ierr) ! residual vector of same shape as solution vector call DMDASNESSetFunctionLocal(damage_grid,INSERT_VALUES,formResidual,PETSC_NULL_SNES,ierr) ! residual vector of same shape as solution vector
CHKERRQ(ierr) CHKERRQ(ierr)
call SNESSetFromOptions(damage_snes,ierr); CHKERRQ(ierr) ! pull it all together with additional CLI arguments call SNESSetFromOptions(damage_snes,ierr); CHKERRQ(ierr) ! pull it all together with additional CLI arguments
call SNESGetType(damage_snes,snes_type,ierr); CHKERRQ(ierr) call SNESGetType(damage_snes,snes_type,ierr); CHKERRQ(ierr)
if (trim(snes_type) == 'vinewtonrsls' .or. & if (trim(snes_type) == 'vinewtonrsls' .or. &
@ -137,12 +138,12 @@ subroutine grid_damage_spectral_init
endif endif
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! init fields ! init fields
call DMDAGetCorners(damage_grid,xstart,ystart,zstart,xend,yend,zend,ierr) call DMDAGetCorners(damage_grid,xstart,ystart,zstart,xend,yend,zend,ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
xend = xstart + xend - 1 xend = xstart + xend - 1
yend = ystart + yend - 1 yend = ystart + yend - 1
zend = zstart + zend - 1 zend = zstart + zend - 1
allocate(phi_current(grid(1),grid(2),grid3), source=1.0_pReal) allocate(phi_current(grid(1),grid(2),grid3), source=1.0_pReal)
allocate(phi_lastInc(grid(1),grid(2),grid3), source=1.0_pReal) allocate(phi_lastInc(grid(1),grid(2),grid3), source=1.0_pReal)
allocate(phi_stagInc(grid(1),grid(2),grid3), source=1.0_pReal) allocate(phi_stagInc(grid(1),grid(2),grid3), source=1.0_pReal)
@ -157,26 +158,26 @@ end subroutine grid_damage_spectral_init
!> @brief solution for the spectral damage scheme with internal iterations !> @brief solution for the spectral damage scheme with internal iterations
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function grid_damage_spectral_solution(timeinc) result(solution) function grid_damage_spectral_solution(timeinc) result(solution)
real(pReal), intent(in) :: & real(pReal), intent(in) :: &
timeinc !< increment in time for current solution timeinc !< increment in time for current solution
integer :: i, j, k, cell integer :: i, j, k, cell
type(tSolutionState) :: solution type(tSolutionState) :: solution
PetscInt :: devNull PetscInt :: devNull
PetscReal :: phi_min, phi_max, stagNorm, solnNorm PetscReal :: phi_min, phi_max, stagNorm, solnNorm
PetscErrorCode :: ierr PetscErrorCode :: ierr
SNESConvergedReason :: reason SNESConvergedReason :: reason
solution%converged =.false. solution%converged =.false.
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! set module wide availabe data ! set module wide availabe data
params%timeinc = timeinc params%timeinc = timeinc
call SNESSolve(damage_snes,PETSC_NULL_VEC,solution_vec,ierr); CHKERRQ(ierr) call SNESSolve(damage_snes,PETSC_NULL_VEC,solution_vec,ierr); CHKERRQ(ierr)
call SNESGetConvergedReason(damage_snes,reason,ierr); CHKERRQ(ierr) call SNESGetConvergedReason(damage_snes,reason,ierr); CHKERRQ(ierr)
if (reason < 1) then if (reason < 1) then
solution%converged = .false. solution%converged = .false.
solution%iterationsNeeded = num%itmax solution%iterationsNeeded = num%itmax
@ -192,20 +193,21 @@ function grid_damage_spectral_solution(timeinc) result(solution)
solution%stagConverged = stagNorm < max(num%eps_damage_atol, num%eps_damage_rtol*solnNorm) solution%stagConverged = stagNorm < max(num%eps_damage_atol, num%eps_damage_rtol*solnNorm)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! updating damage state ! updating damage state
cell = 0 cell = 0
do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1)
cell = cell + 1 cell = cell + 1
call damage_nonlocal_putNonLocalDamage(phi_current(i,j,k),1,cell) call damage_nonlocal_putNonLocalDamage(phi_current(i,j,k),1,cell)
homogenization_phi(cell) = phi_current(i,j,k)
enddo; enddo; enddo enddo; enddo; enddo
call VecMin(solution_vec,devNull,phi_min,ierr); CHKERRQ(ierr) call VecMin(solution_vec,devNull,phi_min,ierr); CHKERRQ(ierr)
call VecMax(solution_vec,devNull,phi_max,ierr); CHKERRQ(ierr) call VecMax(solution_vec,devNull,phi_max,ierr); CHKERRQ(ierr)
if (solution%converged) & if (solution%converged) &
print'(/,a)', ' ... nonlocal damage converged .....................................' print'(/,a)', ' ... nonlocal damage converged .....................................'
print'(/,a,f8.6,2x,f8.6,2x,e11.4)', ' Minimum|Maximum|Delta Damage = ', phi_min, phi_max, stagNorm print'(/,a,f8.6,2x,f8.6,2x,e11.4)', ' Minimum|Maximum|Delta Damage = ', phi_min, phi_max, stagNorm
print'(/,a)', ' ===========================================================================' print'(/,a)', ' ==========================================================================='
flush(IO_STDOUT) flush(IO_STDOUT)
end function grid_damage_spectral_solution end function grid_damage_spectral_solution
@ -214,31 +216,32 @@ end function grid_damage_spectral_solution
!> @brief spectral damage forwarding routine !> @brief spectral damage forwarding routine
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine grid_damage_spectral_forward(cutBack) subroutine grid_damage_spectral_forward(cutBack)
logical, intent(in) :: cutBack logical, intent(in) :: cutBack
integer :: i, j, k, cell integer :: i, j, k, cell
DM :: dm_local DM :: dm_local
PetscScalar, dimension(:,:,:), pointer :: x_scal PetscScalar, dimension(:,:,:), pointer :: x_scal
PetscErrorCode :: ierr PetscErrorCode :: ierr
if (cutBack) then if (cutBack) then
phi_current = phi_lastInc phi_current = phi_lastInc
phi_stagInc = phi_lastInc phi_stagInc = phi_lastInc
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! reverting damage field state ! reverting damage field state
cell = 0 cell = 0
call SNESGetDM(damage_snes,dm_local,ierr); CHKERRQ(ierr) call SNESGetDM(damage_snes,dm_local,ierr); CHKERRQ(ierr)
call DMDAVecGetArrayF90(dm_local,solution_vec,x_scal,ierr); CHKERRQ(ierr) !< get the data out of PETSc to work with call DMDAVecGetArrayF90(dm_local,solution_vec,x_scal,ierr); CHKERRQ(ierr) !< get the data out of PETSc to work with
x_scal(xstart:xend,ystart:yend,zstart:zend) = phi_current x_scal(xstart:xend,ystart:yend,zstart:zend) = phi_current
call DMDAVecRestoreArrayF90(dm_local,solution_vec,x_scal,ierr); CHKERRQ(ierr) call DMDAVecRestoreArrayF90(dm_local,solution_vec,x_scal,ierr); CHKERRQ(ierr)
do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1)
cell = cell + 1 cell = cell + 1
call damage_nonlocal_putNonLocalDamage(phi_current(i,j,k),1,cell) call damage_nonlocal_putNonLocalDamage(phi_current(i,j,k),1,cell)
homogenization_phi(cell) = phi_current(i,j,k)
enddo; enddo; enddo enddo; enddo; enddo
else else
phi_lastInc = phi_current phi_lastInc = phi_current
call updateReference call updateReference
endif endif
end subroutine grid_damage_spectral_forward end subroutine grid_damage_spectral_forward
@ -247,7 +250,7 @@ end subroutine grid_damage_spectral_forward
!> @brief forms the spectral damage residual vector !> @brief forms the spectral damage residual vector
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine formResidual(in,x_scal,f_scal,dummy,ierr) subroutine formResidual(in,x_scal,f_scal,dummy,ierr)
DMDALocalInfo, dimension(DMDA_LOCAL_INFO_SIZE) :: & DMDALocalInfo, dimension(DMDA_LOCAL_INFO_SIZE) :: &
in in
PetscScalar, dimension( & PetscScalar, dimension( &
@ -261,11 +264,11 @@ subroutine formResidual(in,x_scal,f_scal,dummy,ierr)
integer :: i, j, k, cell integer :: i, j, k, cell
real(pReal) :: phiDot, dPhiDot_dPhi, mobility real(pReal) :: phiDot, dPhiDot_dPhi, mobility
phi_current = x_scal phi_current = x_scal
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! evaluate polarization field ! evaluate polarization field
scalarField_real = 0.0_pReal scalarField_real = 0.0_pReal
scalarField_real(1:grid(1),1:grid(2),1:grid3) = phi_current scalarField_real(1:grid(1),1:grid(2),1:grid3) = phi_current
call utilities_FFTscalarForward call utilities_FFTscalarForward
call utilities_fourierScalarGradient !< calculate gradient of damage field call utilities_fourierScalarGradient !< calculate gradient of damage field
call utilities_FFTvectorBackward call utilities_FFTvectorBackward
@ -297,7 +300,7 @@ subroutine formResidual(in,x_scal,f_scal,dummy,ierr)
scalarField_real(1:grid(1),1:grid(2),1:grid3) = phi_lastInc scalarField_real(1:grid(1),1:grid(2),1:grid3) = phi_lastInc
where(scalarField_real(1:grid(1),1:grid(2),1:grid3) < num%residualStiffness) & where(scalarField_real(1:grid(1),1:grid(2),1:grid3) < num%residualStiffness) &
scalarField_real(1:grid(1),1:grid(2),1:grid3) = num%residualStiffness scalarField_real(1:grid(1),1:grid(2),1:grid3) = num%residualStiffness
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! constructing residual ! constructing residual
f_scal = scalarField_real(1:grid(1),1:grid(2),1:grid3) - phi_current f_scal = scalarField_real(1:grid(1),1:grid(2),1:grid3) - phi_current
@ -311,7 +314,7 @@ end subroutine formResidual
subroutine updateReference subroutine updateReference
integer :: i,j,k,cell,ierr integer :: i,j,k,cell,ierr
cell = 0 cell = 0
K_ref = 0.0_pReal K_ref = 0.0_pReal
mu_ref = 0.0_pReal mu_ref = 0.0_pReal

15
src/homogenization.f90
View File

@ -29,7 +29,9 @@ module homogenization
! General variables for the homogenization at a material point ! General variables for the homogenization at a material point
real(pReal), dimension(:), allocatable, public :: & real(pReal), dimension(:), allocatable, public :: &
homogenization_T, & homogenization_T, &
homogenization_dot_T homogenization_dot_T, &
homogenization_phi, &
homogenization_dot_phi
real(pReal), dimension(:,:,:), allocatable, public :: & real(pReal), dimension(:,:,:), allocatable, public :: &
homogenization_F0, & !< def grad of IP at start of FE increment homogenization_F0, & !< def grad of IP at start of FE increment
homogenization_F !< def grad of IP to be reached at end of FE increment homogenization_F !< def grad of IP to be reached at end of FE increment
@ -62,6 +64,9 @@ module homogenization
module subroutine thermal_init module subroutine thermal_init
end subroutine thermal_init end subroutine thermal_init
module subroutine damage_init
end subroutine damage_init
module subroutine mech_partition(subF,ip,el) module subroutine mech_partition(subF,ip,el)
real(pReal), intent(in), dimension(3,3) :: & real(pReal), intent(in), dimension(3,3) :: &
subF subF
@ -75,6 +80,11 @@ module homogenization
integer, intent(in) :: ce integer, intent(in) :: ce
end subroutine thermal_partition end subroutine thermal_partition
module subroutine damage_partition(phi,ce)
real(pReal), intent(in) :: phi
integer, intent(in) :: ce
end subroutine damage_partition
module subroutine thermal_homogenize(ip,el) module subroutine thermal_homogenize(ip,el)
integer, intent(in) :: ip,el integer, intent(in) :: ip,el
end subroutine thermal_homogenize end subroutine thermal_homogenize
@ -118,7 +128,7 @@ contains
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief module initialization !> @brief module initialization
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine homogenization_init subroutine homogenization_init()
class (tNode) , pointer :: & class (tNode) , pointer :: &
num_homog, & num_homog, &
@ -142,6 +152,7 @@ subroutine homogenization_init
call mech_init(num_homog) call mech_init(num_homog)
call thermal_init() call thermal_init()
call damage_init()
if (any(thermal_type == THERMAL_isothermal_ID)) call thermal_isothermal_init(homogenization_T) if (any(thermal_type == THERMAL_isothermal_ID)) call thermal_isothermal_init(homogenization_T)
if (any(thermal_type == THERMAL_conduction_ID)) call thermal_conduction_init(homogenization_T) if (any(thermal_type == THERMAL_conduction_ID)) call thermal_conduction_init(homogenization_T)

40
src/homogenization_damage.f90 Normal file
View File

@ -0,0 +1,40 @@
!--------------------------------------------------------------------------------------------------
!> @author Martin Diehl, KU Leuven
!--------------------------------------------------------------------------------------------------
submodule(homogenization) homogenization_damage
contains
!--------------------------------------------------------------------------------------------------
!> @brief Allocate variables and set parameters.
!--------------------------------------------------------------------------------------------------
module subroutine damage_init()
print'(/,a)', ' <<<+- homogenization_damage init -+>>>'
allocate(homogenization_phi(discretization_nIPs*discretization_Nelems))
allocate(homogenization_dot_phi(discretization_nIPs*discretization_Nelems))
end subroutine damage_init
!--------------------------------------------------------------------------------------------------
!> @brief Partition temperature onto the individual constituents.
!--------------------------------------------------------------------------------------------------
module subroutine damage_partition(phi,ce)
real(pReal), intent(in) :: phi
integer, intent(in) :: ce
integer :: co
do co = 1, homogenization_Nconstituents(material_homogenizationAt2(ce))
call constitutive_damage_set_phi(phi,co,ce)
enddo
end subroutine damage_partition
end submodule homogenization_damage

5
src/homogenization_mech_RGC.f90
View File

@ -656,8 +656,11 @@ module function mech_RGC_updateState(P,F,avgF,dt,dPdF,ip,el) result(doneAndHappy
ip, & !< integration point number ip, & !< integration point number
el !< element number el !< element number
real(pReal), dimension(6,6) :: C
equivalentMu = lattice_equivalent_mu(constitutive_homogenizedC(grainID,ip,el),'voigt')
C = constitutive_homogenizedC(material_phaseAt(grainID,el),material_phaseMemberAt(grainID,ip,el))
equivalentMu = lattice_equivalent_mu(C,'voigt')
end function equivalentMu end function equivalentMu

38
src/source_damage_anisoBrittle.f90
View File

@ -15,7 +15,7 @@ submodule (constitutive:constitutive_damage) source_damage_anisoBrittle
dot_o, & !< opening rate of cleavage planes dot_o, & !< opening rate of cleavage planes
q !< damage rate sensitivity q !< damage rate sensitivity
real(pReal), dimension(:), allocatable :: & real(pReal), dimension(:), allocatable :: &
s_crit, & !< critical displacement s_crit, & !< critical displacement
g_crit !< critical load g_crit !< critical load
real(pReal), dimension(:,:,:,:), allocatable :: & real(pReal), dimension(:,:,:,:), allocatable :: &
cleavage_systems cleavage_systems
@ -37,7 +37,7 @@ contains
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module function source_damage_anisoBrittle_init(source_length) result(mySources) module function source_damage_anisoBrittle_init(source_length) result(mySources)
integer, intent(in) :: source_length integer, intent(in) :: source_length
logical, dimension(:,:), allocatable :: mySources logical, dimension(:,:), allocatable :: mySources
class(tNode), pointer :: & class(tNode), pointer :: &
@ -62,7 +62,7 @@ module function source_damage_anisoBrittle_init(source_length) result(mySources)
allocate(source_damage_anisoBrittle_instance(phases%length), source=0) allocate(source_damage_anisoBrittle_instance(phases%length), source=0)
do p = 1, phases%length do p = 1, phases%length
phase => phases%get(p) phase => phases%get(p)
if(any(mySources(:,p))) source_damage_anisoBrittle_instance(p) = count(mySources(:,1:p)) if(any(mySources(:,p))) source_damage_anisoBrittle_instance(p) = count(mySources(:,1:p))
if(count(mySources(:,p)) == 0) cycle if(count(mySources(:,p)) == 0) cycle
sources => phase%get('source') sources => phase%get('source')
@ -70,20 +70,20 @@ module function source_damage_anisoBrittle_init(source_length) result(mySources)
if(mySources(sourceOffset,p)) then if(mySources(sourceOffset,p)) then
source_damage_anisoBrittle_offset(p) = sourceOffset source_damage_anisoBrittle_offset(p) = sourceOffset
associate(prm => param(source_damage_anisoBrittle_instance(p))) associate(prm => param(source_damage_anisoBrittle_instance(p)))
src => sources%get(sourceOffset) src => sources%get(sourceOffset)
N_cl = src%get_asInts('N_cl',defaultVal=emptyIntArray) N_cl = src%get_asInts('N_cl',defaultVal=emptyIntArray)
prm%sum_N_cl = sum(abs(N_cl)) prm%sum_N_cl = sum(abs(N_cl))
prm%q = src%get_asFloat('q') prm%q = src%get_asFloat('q')
prm%dot_o = src%get_asFloat('dot_o') prm%dot_o = src%get_asFloat('dot_o')
prm%s_crit = src%get_asFloats('s_crit', requiredSize=size(N_cl)) prm%s_crit = src%get_asFloats('s_crit', requiredSize=size(N_cl))
prm%g_crit = src%get_asFloats('g_crit', requiredSize=size(N_cl)) prm%g_crit = src%get_asFloats('g_crit', requiredSize=size(N_cl))
prm%cleavage_systems = lattice_SchmidMatrix_cleavage(N_cl,phase%get_asString('lattice'),& prm%cleavage_systems = lattice_SchmidMatrix_cleavage(N_cl,phase%get_asString('lattice'),&
phase%get_asFloat('c/a',defaultVal=0.0_pReal)) phase%get_asFloat('c/a',defaultVal=0.0_pReal))
! expand: family => system ! expand: family => system
prm%s_crit = math_expand(prm%s_crit,N_cl) prm%s_crit = math_expand(prm%s_crit,N_cl)
prm%g_crit = math_expand(prm%g_crit,N_cl) prm%g_crit = math_expand(prm%g_crit,N_cl)
@ -93,7 +93,7 @@ module function source_damage_anisoBrittle_init(source_length) result(mySources)
#else #else
prm%output = src%get_asStrings('output',defaultVal=emptyStringArray) prm%output = src%get_asStrings('output',defaultVal=emptyStringArray)
#endif #endif
! sanity checks ! sanity checks
if (prm%q <= 0.0_pReal) extmsg = trim(extmsg)//' q' if (prm%q <= 0.0_pReal) extmsg = trim(extmsg)//' q'
if (prm%dot_o <= 0.0_pReal) extmsg = trim(extmsg)//' dot_o' if (prm%dot_o <= 0.0_pReal) extmsg = trim(extmsg)//' dot_o'
@ -130,8 +130,8 @@ module subroutine source_damage_anisoBrittle_dotState(S, co, ip, el)
S S
integer :: & integer :: &
phase, & ph, &
constituent, & me, &
sourceOffset, & sourceOffset, &
damageOffset, & damageOffset, &
homog, & homog, &
@ -139,14 +139,14 @@ module subroutine source_damage_anisoBrittle_dotState(S, co, ip, el)
real(pReal) :: & real(pReal) :: &
traction_d, traction_t, traction_n, traction_crit traction_d, traction_t, traction_n, traction_crit
phase = material_phaseAt(co,el) ph = material_phaseAt(co,el)
constituent = material_phasememberAt(co,ip,el) me = material_phasememberAt(co,ip,el)
sourceOffset = source_damage_anisoBrittle_offset(phase) sourceOffset = source_damage_anisoBrittle_offset(ph)
homog = material_homogenizationAt(el) homog = material_homogenizationAt(el)
damageOffset = material_homogenizationMemberAt(ip,el) damageOffset = material_homogenizationMemberAt(ip,el)
associate(prm => param(source_damage_anisoBrittle_instance(phase))) associate(prm => param(source_damage_anisoBrittle_instance(ph)))
damageState(phase)%p(sourceOffset)%dotState(1,constituent) = 0.0_pReal damageState(ph)%p(sourceOffset)%dotState(1,me) = 0.0_pReal
do i = 1, prm%sum_N_cl do i = 1, prm%sum_N_cl
traction_d = math_tensordot(S,prm%cleavage_systems(1:3,1:3,1,i)) traction_d = math_tensordot(S,prm%cleavage_systems(1:3,1:3,1,i))
traction_t = math_tensordot(S,prm%cleavage_systems(1:3,1:3,2,i)) traction_t = math_tensordot(S,prm%cleavage_systems(1:3,1:3,2,i))
@ -154,8 +154,8 @@ module subroutine source_damage_anisoBrittle_dotState(S, co, ip, el)
traction_crit = prm%g_crit(i)*damage(homog)%p(damageOffset)**2.0_pReal traction_crit = prm%g_crit(i)*damage(homog)%p(damageOffset)**2.0_pReal
damageState(phase)%p(sourceOffset)%dotState(1,constituent) & damageState(ph)%p(sourceOffset)%dotState(1,me) &
= damageState(phase)%p(sourceOffset)%dotState(1,constituent) & = damageState(ph)%p(sourceOffset)%dotState(1,me) &
+ prm%dot_o / prm%s_crit(i) & + prm%dot_o / prm%s_crit(i) &
* ((max(0.0_pReal, abs(traction_d) - traction_crit)/traction_crit)**prm%q + & * ((max(0.0_pReal, abs(traction_d) - traction_crit)/traction_crit)**prm%q + &
(max(0.0_pReal, abs(traction_t) - traction_crit)/traction_crit)**prm%q + & (max(0.0_pReal, abs(traction_t) - traction_crit)/traction_crit)**prm%q + &

22
src/source_damage_anisoDuctile.f90
View File

@ -30,7 +30,7 @@ contains
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module function source_damage_anisoDuctile_init(source_length) result(mySources) module function source_damage_anisoDuctile_init(source_length) result(mySources)
integer, intent(in) :: source_length integer, intent(in) :: source_length
logical, dimension(:,:), allocatable :: mySources logical, dimension(:,:), allocatable :: mySources
class(tNode), pointer :: & class(tNode), pointer :: &
@ -67,7 +67,7 @@ module function source_damage_anisoDuctile_init(source_length) result(mySources)
if(mySources(sourceOffset,p)) then if(mySources(sourceOffset,p)) then
source_damage_anisoDuctile_offset(p) = sourceOffset source_damage_anisoDuctile_offset(p) = sourceOffset
associate(prm => param(source_damage_anisoDuctile_instance(p))) associate(prm => param(source_damage_anisoDuctile_instance(p)))
src => sources%get(sourceOffset) src => sources%get(sourceOffset)
N_sl = pl%get_asInts('N_sl',defaultVal=emptyIntArray) N_sl = pl%get_asInts('N_sl',defaultVal=emptyIntArray)
prm%q = src%get_asFloat('q') prm%q = src%get_asFloat('q')
@ -81,7 +81,7 @@ module function source_damage_anisoDuctile_init(source_length) result(mySources)
#else #else
prm%output = src%get_asStrings('output',defaultVal=emptyStringArray) prm%output = src%get_asStrings('output',defaultVal=emptyStringArray)
#endif #endif
! sanity checks ! sanity checks
if (prm%q <= 0.0_pReal) extmsg = trim(extmsg)//' q' if (prm%q <= 0.0_pReal) extmsg = trim(extmsg)//' q'
if (any(prm%gamma_crit < 0.0_pReal)) extmsg = trim(extmsg)//' gamma_crit' if (any(prm%gamma_crit < 0.0_pReal)) extmsg = trim(extmsg)//' gamma_crit'
@ -115,21 +115,21 @@ module subroutine source_damage_anisoDuctile_dotState(co, ip, el)
el !< element el !< element
integer :: & integer :: &
phase, & ph, &
constituent, & me, &
sourceOffset, & sourceOffset, &
damageOffset, & damageOffset, &
homog homog
phase = material_phaseAt(co,el) ph = material_phaseAt(co,el)
constituent = material_phasememberAt(co,ip,el) me = material_phasememberAt(co,ip,el)
sourceOffset = source_damage_anisoDuctile_offset(phase) sourceOffset = source_damage_anisoDuctile_offset(ph)
homog = material_homogenizationAt(el) homog = material_homogenizationAt(el)
damageOffset = material_homogenizationMemberAt(ip,el) damageOffset = material_homogenizationMemberAt(ip,el)
associate(prm => param(source_damage_anisoDuctile_instance(phase))) associate(prm => param(source_damage_anisoDuctile_instance(ph)))
damageState(phase)%p(sourceOffset)%dotState(1,constituent) & damageState(ph)%p(sourceOffset)%dotState(1,me) &
= sum(plasticState(phase)%slipRate(:,constituent)/(damage(homog)%p(damageOffset)**prm%q)/prm%gamma_crit) = sum(plasticState(ph)%slipRate(:,me)/(damage(homog)%p(damageOffset)**prm%q)/prm%gamma_crit)
end associate end associate
end subroutine source_damage_anisoDuctile_dotState end subroutine source_damage_anisoDuctile_dotState

35
src/source_damage_isoBrittle.f90
View File

@ -28,7 +28,7 @@ contains
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module function source_damage_isoBrittle_init(source_length) result(mySources) module function source_damage_isoBrittle_init(source_length) result(mySources)
integer, intent(in) :: source_length integer, intent(in) :: source_length
logical, dimension(:,:), allocatable :: mySources logical, dimension(:,:), allocatable :: mySources
class(tNode), pointer :: & class(tNode), pointer :: &
@ -52,7 +52,7 @@ module function source_damage_isoBrittle_init(source_length) result(mySources)
allocate(source_damage_isoBrittle_instance(phases%length), source=0) allocate(source_damage_isoBrittle_instance(phases%length), source=0)
do p = 1, phases%length do p = 1, phases%length
phase => phases%get(p) phase => phases%get(p)
if(any(mySources(:,p))) source_damage_isoBrittle_instance(p) = count(mySources(:,1:p)) if(any(mySources(:,p))) source_damage_isoBrittle_instance(p) = count(mySources(:,1:p))
if(count(mySources(:,p)) == 0) cycle if(count(mySources(:,p)) == 0) cycle
sources => phase%get('source') sources => phase%get('source')
@ -60,7 +60,7 @@ module function source_damage_isoBrittle_init(source_length) result(mySources)
if(mySources(sourceOffset,p)) then if(mySources(sourceOffset,p)) then
source_damage_isoBrittle_offset(p) = sourceOffset source_damage_isoBrittle_offset(p) = sourceOffset
associate(prm => param(source_damage_isoBrittle_instance(p))) associate(prm => param(source_damage_isoBrittle_instance(p)))
src => sources%get(sourceOffset) src => sources%get(sourceOffset)
prm%W_crit = src%get_asFloat('W_crit') prm%W_crit = src%get_asFloat('W_crit')
@ -69,7 +69,7 @@ module function source_damage_isoBrittle_init(source_length) result(mySources)
#else #else
prm%output = src%get_asStrings('output',defaultVal=emptyStringArray) prm%output = src%get_asStrings('output',defaultVal=emptyStringArray)
#endif #endif
! sanity checks ! sanity checks
if (prm%W_crit <= 0.0_pReal) extmsg = trim(extmsg)//' W_crit' if (prm%W_crit <= 0.0_pReal) extmsg = trim(extmsg)//' W_crit'
@ -94,43 +94,36 @@ end function source_damage_isoBrittle_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief calculates derived quantities from state !> @brief calculates derived quantities from state
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module subroutine source_damage_isoBrittle_deltaState(C, Fe, co, ip, el) module subroutine source_damage_isoBrittle_deltaState(C, Fe, ph,me)
integer, intent(in) :: & integer, intent(in) :: ph,me
co, & !< component-ID of integration point
ip, & !< integration point
el !< element
real(pReal), intent(in), dimension(3,3) :: & real(pReal), intent(in), dimension(3,3) :: &
Fe Fe
real(pReal), intent(in), dimension(6,6) :: & real(pReal), intent(in), dimension(6,6) :: &
C C
integer :: & integer :: &
phase, &
constituent, &
sourceOffset sourceOffset
real(pReal), dimension(6) :: & real(pReal), dimension(6) :: &
strain strain
real(pReal) :: & real(pReal) :: &
strainenergy strainenergy
phase = material_phaseAt(co,el) !< phase ID at co,ip,el sourceOffset = source_damage_isoBrittle_offset(ph)
constituent = material_phasememberAt(co,ip,el) !< state array offset for phase ID at co,ip,el
sourceOffset = source_damage_isoBrittle_offset(phase)
strain = 0.5_pReal*math_sym33to6(matmul(transpose(Fe),Fe)-math_I3) strain = 0.5_pReal*math_sym33to6(matmul(transpose(Fe),Fe)-math_I3)
associate(prm => param(source_damage_isoBrittle_instance(phase))) associate(prm => param(source_damage_isoBrittle_instance(ph)))
strainenergy = 2.0_pReal*sum(strain*matmul(C,strain))/prm%W_crit strainenergy = 2.0_pReal*sum(strain*matmul(C,strain))/prm%W_crit
! ToDo: check strainenergy = 2.0_pReal*dot_product(strain,matmul(C,strain))/prm%W_crit ! ToDo: check strainenergy = 2.0_pReal*dot_product(strain,matmul(C,strain))/prm%W_crit
if (strainenergy > damageState(phase)%p(sourceOffset)%subState0(1,constituent)) then if (strainenergy > damageState(ph)%p(sourceOffset)%subState0(1,me)) then
damageState(phase)%p(sourceOffset)%deltaState(1,constituent) = & damageState(ph)%p(sourceOffset)%deltaState(1,me) = &
strainenergy - damageState(phase)%p(sourceOffset)%state(1,constituent) strainenergy - damageState(ph)%p(sourceOffset)%state(1,me)
else else
damageState(phase)%p(sourceOffset)%deltaState(1,constituent) = & damageState(ph)%p(sourceOffset)%deltaState(1,me) = &
damageState(phase)%p(sourceOffset)%subState0(1,constituent) - & damageState(ph)%p(sourceOffset)%subState0(1,me) - &
damageState(phase)%p(sourceOffset)%state(1,constituent) damageState(ph)%p(sourceOffset)%state(1,me)
endif endif
end associate end associate

24
src/source_damage_isoDuctile.f90
View File

@ -30,7 +30,7 @@ contains
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module function source_damage_isoDuctile_init(source_length) result(mySources) module function source_damage_isoDuctile_init(source_length) result(mySources)
integer, intent(in) :: source_length integer, intent(in) :: source_length
logical, dimension(:,:), allocatable :: mySources logical, dimension(:,:), allocatable :: mySources
class(tNode), pointer :: & class(tNode), pointer :: &
@ -54,7 +54,7 @@ module function source_damage_isoDuctile_init(source_length) result(mySources)
allocate(source_damage_isoDuctile_instance(phases%length), source=0) allocate(source_damage_isoDuctile_instance(phases%length), source=0)
do p = 1, phases%length do p = 1, phases%length
phase => phases%get(p) phase => phases%get(p)
if(count(mySources(:,p)) == 0) cycle if(count(mySources(:,p)) == 0) cycle
if(any(mySources(:,p))) source_damage_isoDuctile_instance(p) = count(mySources(:,1:p)) if(any(mySources(:,p))) source_damage_isoDuctile_instance(p) = count(mySources(:,1:p))
sources => phase%get('source') sources => phase%get('source')
@ -62,7 +62,7 @@ module function source_damage_isoDuctile_init(source_length) result(mySources)
if(mySources(sourceOffset,p)) then if(mySources(sourceOffset,p)) then
source_damage_isoDuctile_offset(p) = sourceOffset source_damage_isoDuctile_offset(p) = sourceOffset
associate(prm => param(source_damage_isoDuctile_instance(p))) associate(prm => param(source_damage_isoDuctile_instance(p)))
src => sources%get(sourceOffset) src => sources%get(sourceOffset)
prm%q = src%get_asFloat('q') prm%q = src%get_asFloat('q')
prm%gamma_crit = src%get_asFloat('gamma_crit') prm%gamma_crit = src%get_asFloat('gamma_crit')
@ -72,7 +72,7 @@ module function source_damage_isoDuctile_init(source_length) result(mySources)
#else #else
prm%output = src%get_asStrings('output',defaultVal=emptyStringArray) prm%output = src%get_asStrings('output',defaultVal=emptyStringArray)
#endif #endif
! sanity checks ! sanity checks
if (prm%q <= 0.0_pReal) extmsg = trim(extmsg)//' q' if (prm%q <= 0.0_pReal) extmsg = trim(extmsg)//' q'
if (prm%gamma_crit <= 0.0_pReal) extmsg = trim(extmsg)//' gamma_crit' if (prm%gamma_crit <= 0.0_pReal) extmsg = trim(extmsg)//' gamma_crit'
@ -106,21 +106,21 @@ module subroutine source_damage_isoDuctile_dotState(co, ip, el)
el !< element el !< element
integer :: & integer :: &
phase, & ph, &
constituent, & me, &
sourceOffset, & sourceOffset, &
damageOffset, & damageOffset, &
homog homog
phase = material_phaseAt(co,el) ph = material_phaseAt(co,el)
constituent = material_phasememberAt(co,ip,el) me = material_phasememberAt(co,ip,el)
sourceOffset = source_damage_isoDuctile_offset(phase) sourceOffset = source_damage_isoDuctile_offset(ph)
homog = material_homogenizationAt(el) homog = material_homogenizationAt(el)
damageOffset = material_homogenizationMemberAt(ip,el) damageOffset = material_homogenizationMemberAt(ip,el)
associate(prm => param(source_damage_isoDuctile_instance(phase))) associate(prm => param(source_damage_isoDuctile_instance(ph)))
damageState(phase)%p(sourceOffset)%dotState(1,constituent) = & damageState(ph)%p(sourceOffset)%dotState(1,me) = &
sum(plasticState(phase)%slipRate(:,constituent))/(damage(homog)%p(damageOffset)**prm%q)/prm%gamma_crit sum(plasticState(ph)%slipRate(:,me))/(damage(homog)%p(damageOffset)**prm%q)/prm%gamma_crit
end associate end associate
end subroutine source_damage_isoDuctile_dotState end subroutine source_damage_isoDuctile_dotState