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

only one damage mechanism per phase 

material.yaml specification is designed to allow more than one, but that
requires to have two phase fields etc.
For the moment, keep it as simple as possible.
This commit is contained in:
Martin Diehl 2021-02-13 08:32:43 +01:00
parent 6e3515982d
commit b3dde6d722
8 changed files with 97 additions and 162 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
PRIVATE

@ -1 +1 @@
Subproject commit 2a0f30ec47473f42e0da922055b72b527730378d Subproject commit 571b06ac4b42ed217eacfa5c71281e571bdab1df

43
src/phase.f90
View File

@ -70,7 +70,7 @@ module phase
type(tPlasticState), allocatable, dimension(:), public :: & type(tPlasticState), allocatable, dimension(:), public :: &
plasticState plasticState
type(tSourceState), allocatable, dimension(:), public :: & type(tState), allocatable, dimension(:), public :: &
damageState damageState
@ -350,14 +350,12 @@ subroutine phase_init
PhaseLoop2:do ph = 1,phases%length PhaseLoop2:do ph = 1,phases%length
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! partition and initialize state ! partition and initialize state
plasticState(ph)%state = plasticState(ph)%state0 plasticState(ph)%state = plasticState(ph)%state0
forall(so = 1:phase_Nsources(ph)) if(phase_Nsources(ph) > 0) &
damageState(ph)%p(so)%state = damageState(ph)%p(so)%state0 damageState(ph)%state = damageState(ph)%state0
end forall
phase_source_maxSizeDotState = max(phase_source_maxSizeDotState, &
maxval(damageState(ph)%p%sizeDotState))
enddo PhaseLoop2 enddo PhaseLoop2
phase_source_maxSizeDotState = maxval(damageState%sizeDotState)
phase_plasticity_maxSizeDotState = maxval(plasticState%sizeDotState) phase_plasticity_maxSizeDotState = maxval(plasticState%sizeDotState)
end subroutine phase_init end subroutine phase_init
@ -404,15 +402,13 @@ subroutine phase_restore(ce,includeL)
integer, intent(in) :: ce integer, intent(in) :: ce
integer :: & integer :: &
co, & !< constituent number co
so
do co = 1,homogenization_Nconstituents(material_homogenizationAt2(ce)) do co = 1,homogenization_Nconstituents(material_homogenizationAt2(ce))
do so = 1, phase_Nsources(material_phaseAt2(co,ce)) if (phase_Nsources(material_phaseAt2(co,ce)) > 0) &
damageState(material_phaseAt2(co,ce))%p(so)%state( :,material_phasememberAt2(co,ce)) = & damageState(material_phaseAt2(co,ce))%state( :,material_phasememberAt2(co,ce)) = &
damageState(material_phaseAt2(co,ce))%p(so)%state0(:,material_phasememberAt2(co,ce)) damageState(material_phaseAt2(co,ce))%state0(:,material_phasememberAt2(co,ce))
enddo
enddo enddo
call mechanical_restore(ce,includeL) call mechanical_restore(ce,includeL)
@ -426,16 +422,16 @@ end subroutine phase_restore
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine phase_forward() subroutine phase_forward()
integer :: ph, so integer :: ph
call mechanical_forward() call mechanical_forward()
call thermal_forward() call thermal_forward()
do ph = 1, size(damageState) do ph = 1, size(damageState)
do so = 1,phase_Nsources(ph) if (phase_Nsources(ph) > 0) &
damageState(ph)%p(so)%state0 = damageState(ph)%p(so)%state damageState(ph)%state0 = damageState(ph)%state
enddo; enddo enddo
end subroutine phase_forward end subroutine phase_forward
@ -531,9 +527,8 @@ subroutine crystallite_init()
phases => config_material%get('phase') phases => config_material%get('phase')
do ph = 1, phases%length do ph = 1, phases%length
do so = 1, phase_Nsources(ph) if (phase_Nsources(ph) > 0) &
allocate(damageState(ph)%p(so)%subState0,source=damageState(ph)%p(so)%state0) ! ToDo: hack allocate(damageState(ph)%subState0,source=damageState(ph)%state0) ! ToDo: hack
enddo
enddo enddo
print'(a42,1x,i10)', ' # of elements: ', eMax print'(a42,1x,i10)', ' # of elements: ', eMax
@ -579,9 +574,9 @@ subroutine phase_windForward(ip,el)
call mechanical_windForward(ph,me) call mechanical_windForward(ph,me)
do so = 1, phase_Nsources(material_phaseAt(co,el)) if (phase_Nsources(ph) > 0) &
damageState(ph)%p(so)%state0(:,me) = damageState(ph)%p(so)%state(:,me) damageState(ph)%state0(:,me) = damageState(ph)%state(:,me)
enddo
enddo enddo

71
src/phase_damage.f90
View File

@ -176,8 +176,8 @@ module subroutine damage_init
phase => phases%get(ph) phase => phases%get(ph)
sources => phase%get('damage',defaultVal=emptyList) sources => phase%get('damage',defaultVal=emptyList)
phase_Nsources(ph) = sources%length if (sources%length > 1) error stop
allocate(damageState(ph)%p(phase_Nsources(ph)))
enddo enddo
allocate(phase_source(maxval(phase_Nsources),phases%length), source = DAMAGE_UNDEFINED_ID) allocate(phase_source(maxval(phase_Nsources),phases%length), source = DAMAGE_UNDEFINED_ID)
@ -267,57 +267,48 @@ module function integrateDamageState(dt,co,ip,el) result(broken)
NiterationState, & !< number of iterations in state loop NiterationState, & !< number of iterations in state loop
ph, & ph, &
me, & me, &
so
integer, dimension(maxval(phase_Nsources)) :: &
size_so size_so
real(pReal) :: & real(pReal) :: &
zeta zeta
real(pReal), dimension(phase_source_maxSizeDotState) :: & real(pReal), dimension(phase_source_maxSizeDotState) :: &
r ! state residuum r ! state residuum
real(pReal), dimension(phase_source_maxSizeDotState,2,maxval(phase_Nsources)) :: source_dotState real(pReal), dimension(phase_source_maxSizeDotState,2) :: source_dotState
logical :: & logical :: &
converged_ converged_
ph = material_phaseAt(co,el) ph = material_phaseAt(co,el)
if (phase_Nsources(ph) == 0) return
me = material_phaseMemberAt(co,ip,el) me = material_phaseMemberAt(co,ip,el)
converged_ = .true. converged_ = .true.
broken = phase_damage_collectDotState(co,ip,el,ph,me) broken = phase_damage_collectDotState(co,ip,el,ph,me)
if(broken) return if(broken) return
do so = 1, phase_Nsources(ph) size_so = damageState(ph)%sizeDotState
size_so(so) = damageState(ph)%p(so)%sizeDotState damageState(ph)%state(1:size_so,me) = damageState(ph)%subState0(1:size_so,me) &
damageState(ph)%p(so)%state(1:size_so(so),me) = damageState(ph)%p(so)%subState0(1:size_so(so),me) & + damageState(ph)%dotState (1:size_so,me) * dt
+ damageState(ph)%p(so)%dotState (1:size_so(so),me) * dt source_dotState(1:size_so,2) = 0.0_pReal
source_dotState(1:size_so(so),2,so) = 0.0_pReal
enddo
iteration: do NiterationState = 1, num%nState iteration: do NiterationState = 1, num%nState
do so = 1, phase_Nsources(ph) if(nIterationState > 1) source_dotState(1:size_so,2) = source_dotState(1:size_so,1)
if(nIterationState > 1) source_dotState(1:size_so(so),2,so) = source_dotState(1:size_so(so),1,so) source_dotState(1:size_so,1) = damageState(ph)%dotState(:,me)
source_dotState(1:size_so(so),1,so) = damageState(ph)%p(so)%dotState(:,me)
enddo
broken = phase_damage_collectDotState(co,ip,el,ph,me) broken = phase_damage_collectDotState(co,ip,el,ph,me)
if(broken) exit iteration if(broken) exit iteration
do so = 1, phase_Nsources(ph)
zeta = damper(damageState(ph)%p(so)%dotState(:,me), & zeta = damper(damageState(ph)%dotState(:,me),source_dotState(1:size_so,1),source_dotState(1:size_so,2))
source_dotState(1:size_so(so),1,so),& damageState(ph)%dotState(:,me) = damageState(ph)%dotState(:,me) * zeta &
source_dotState(1:size_so(so),2,so)) + source_dotState(1:size_so,1)* (1.0_pReal - zeta)
damageState(ph)%p(so)%dotState(:,me) = damageState(ph)%p(so)%dotState(:,me) * zeta & r(1:size_so) = damageState(ph)%state (1:size_so,me) &
+ source_dotState(1:size_so(so),1,so)* (1.0_pReal - zeta) - damageState(ph)%subState0(1:size_so,me) &
r(1:size_so(so)) = damageState(ph)%p(so)%state (1:size_so(so),me) & - damageState(ph)%dotState (1:size_so,me) * dt
- damageState(ph)%p(so)%subState0(1:size_so(so),me) & damageState(ph)%state(1:size_so,me) = damageState(ph)%state(1:size_so,me) - r(1:size_so)
- damageState(ph)%p(so)%dotState (1:size_so(so),me) * dt converged_ = converged_ .and. converged(r(1:size_so), &
damageState(ph)%p(so)%state(1:size_so(so),me) = damageState(ph)%p(so)%state(1:size_so(so),me) & damageState(ph)%state(1:size_so,me), &
- r(1:size_so(so)) damageState(ph)%atol(1:size_so))
converged_ = converged_ .and. converged(r(1:size_so(so)), &
damageState(ph)%p(so)%state(1:size_so(so),me), &
damageState(ph)%p(so)%atol(1:size_so(so)))
enddo
if(converged_) then if(converged_) then
broken = phase_damage_deltaState(mechanical_F_e(ph,me),ph,me) broken = phase_damage_deltaState(mechanical_F_e(ph,me),ph,me)
@ -423,7 +414,7 @@ function phase_damage_collectDotState(co,ip,el,ph,me) result(broken)
end select sourceType end select sourceType
broken = broken .or. any(IEEE_is_NaN(damageState(ph)%p(so)%dotState(:,me))) broken = broken .or. any(IEEE_is_NaN(damageState(ph)%dotState(:,me)))
enddo SourceLoop enddo SourceLoop
@ -443,7 +434,6 @@ function phase_damage_deltaState(Fe, ph, me) result(broken)
real(pReal), intent(in), dimension(3,3) :: & real(pReal), intent(in), dimension(3,3) :: &
Fe !< elastic deformation gradient Fe !< elastic deformation gradient
integer :: & integer :: &
so, &
myOffset, & myOffset, &
mySize mySize
logical :: & logical :: &
@ -452,23 +442,22 @@ function phase_damage_deltaState(Fe, ph, me) result(broken)
broken = .false. broken = .false.
sourceLoop: do so = 1, phase_Nsources(ph) if (phase_Nsources(ph) == 0) return
sourceType: select case (phase_source(so,ph)) sourceType: select case (phase_source(1,ph))
case (DAMAGE_ISOBRITTLE_ID) sourceType case (DAMAGE_ISOBRITTLE_ID) sourceType
call source_damage_isoBrittle_deltaState(phase_homogenizedC(ph,me), Fe, ph,me) call source_damage_isoBrittle_deltaState(phase_homogenizedC(ph,me), Fe, ph,me)
broken = any(IEEE_is_NaN(damageState(ph)%p(so)%deltaState(:,me))) broken = any(IEEE_is_NaN(damageState(ph)%deltaState(:,me)))
if(.not. broken) then if(.not. broken) then
myOffset = damageState(ph)%p(so)%offsetDeltaState myOffset = damageState(ph)%offsetDeltaState
mySize = damageState(ph)%p(so)%sizeDeltaState mySize = damageState(ph)%sizeDeltaState
damageState(ph)%p(so)%state(myOffset + 1: myOffset + mySize,me) = & damageState(ph)%state(myOffset + 1: myOffset + mySize,me) = &
damageState(ph)%p(so)%state(myOffset + 1: myOffset + mySize,me) + damageState(ph)%p(so)%deltaState(1:mySize,me) damageState(ph)%state(myOffset + 1: myOffset + mySize,me) + damageState(ph)%deltaState(1:mySize,me)
endif endif
end select sourceType end select sourceType
enddo SourceLoop
end function phase_damage_deltaState end function phase_damage_deltaState

30
src/phase_damage_anisobrittle.f90
View File

@ -6,8 +6,7 @@
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
submodule (phase:damagee) anisobrittle submodule (phase:damagee) anisobrittle
integer, dimension(:), allocatable :: & integer, dimension(:), allocatable :: &
source_damage_anisoBrittle_offset, & !< which source is my current source mechanism?
source_damage_anisoBrittle_instance !< instance of source mechanism source_damage_anisoBrittle_instance !< instance of source mechanism
type :: tParameters !< container type for internal constitutive parameters type :: tParameters !< container type for internal constitutive parameters
@ -58,7 +57,6 @@ module function anisobrittle_init(source_length) result(mySources)
phases => config_material%get('phase') phases => config_material%get('phase')
allocate(param(Ninstances)) allocate(param(Ninstances))
allocate(source_damage_anisoBrittle_offset (phases%length), source=0)
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
@ -68,7 +66,6 @@ module function anisobrittle_init(source_length) result(mySources)
sources => phase%get('damage') sources => phase%get('damage')
do sourceOffset = 1, sources%length do sourceOffset = 1, sources%length
if(mySources(sourceOffset,p)) then if(mySources(sourceOffset,p)) then
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)
@ -101,9 +98,9 @@ module function anisobrittle_init(source_length) result(mySources)
if (any(prm%s_crit < 0.0_pReal)) extmsg = trim(extmsg)//' s_crit' if (any(prm%s_crit < 0.0_pReal)) extmsg = trim(extmsg)//' s_crit'
Nconstituents = count(material_phaseAt==p) * discretization_nIPs Nconstituents = count(material_phaseAt==p) * discretization_nIPs
call phase_allocateState(damageState(p)%p(sourceOffset),Nconstituents,1,1,0) call phase_allocateState(damageState(p),Nconstituents,1,1,0)
damageState(p)%p(sourceOffset)%atol = src%get_asFloat('anisobrittle_atol',defaultVal=1.0e-3_pReal) damageState(p)%atol = src%get_asFloat('anisobrittle_atol',defaultVal=1.0e-3_pReal)
if(any(damageState(p)%p(sourceOffset)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' anisobrittle_atol' if(any(damageState(p)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' anisobrittle_atol'
end associate end associate
@ -139,14 +136,13 @@ module subroutine 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
ph = material_phaseAt(co,el) ph = material_phaseAt(co,el)
me = material_phasememberAt(co,ip,el) me = material_phasememberAt(co,ip,el)
sourceOffset = source_damage_anisoBrittle_offset(ph)
homog = material_homogenizationAt(el)
damageOffset = material_homogenizationMemberAt(ip,el)
associate(prm => param(source_damage_anisoBrittle_instance(ph))) associate(prm => param(source_damage_anisoBrittle_instance(ph)))
damageState(ph)%p(sourceOffset)%dotState(1,me) = 0.0_pReal damageState(ph)%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 +150,8 @@ module subroutine anisobrittle_dotState(S, co, ip, el)
traction_crit = prm%g_crit(i)*phase_damage_get_phi(co,ip,el)**2.0_pReal traction_crit = prm%g_crit(i)*phase_damage_get_phi(co,ip,el)**2.0_pReal
damageState(ph)%p(sourceOffset)%dotState(1,me) & damageState(ph)%dotState(1,me) &
= damageState(ph)%p(sourceOffset)%dotState(1,me) & = damageState(ph)%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 + &
@ -180,12 +176,8 @@ module subroutine source_damage_anisobrittle_getRateAndItsTangent(localphiDot, d
localphiDot, & localphiDot, &
dLocalphiDot_dPhi dLocalphiDot_dPhi
integer :: &
sourceOffset
sourceOffset = source_damage_anisoBrittle_offset(phase) dLocalphiDot_dPhi = -damageState(phase)%state(1,constituent)
dLocalphiDot_dPhi = -damageState(phase)%p(sourceOffset)%state(1,constituent)
localphiDot = 1.0_pReal & localphiDot = 1.0_pReal &
+ dLocalphiDot_dPhi*phi + dLocalphiDot_dPhi*phi
@ -204,7 +196,7 @@ module subroutine anisobrittle_results(phase,group)
integer :: o integer :: o
associate(prm => param(source_damage_anisoBrittle_instance(phase)), & associate(prm => param(source_damage_anisoBrittle_instance(phase)), &
stt => damageState(phase)%p(source_damage_anisoBrittle_offset(phase))%state) stt => damageState(phase)%state)
outputsLoop: do o = 1,size(prm%output) outputsLoop: do o = 1,size(prm%output)
select case(trim(prm%output(o))) select case(trim(prm%output(o)))
case ('f_phi') case ('f_phi')

29
src/phase_damage_anisoductile.f90
View File

@ -7,7 +7,6 @@
submodule(phase:damagee) anisoductile submodule(phase:damagee) anisoductile
integer, dimension(:), allocatable :: & integer, dimension(:), allocatable :: &
source_damage_anisoDuctile_offset, & !< which source is my current damage mechanism?
source_damage_anisoDuctile_instance !< instance of damage source mechanism source_damage_anisoDuctile_instance !< instance of damage source mechanism
type :: tParameters !< container type for internal constitutive parameters type :: tParameters !< container type for internal constitutive parameters
@ -53,7 +52,6 @@ module function anisoductile_init(source_length) result(mySources)
phases => config_material%get('phase') phases => config_material%get('phase')
allocate(param(Ninstances)) allocate(param(Ninstances))
allocate(source_damage_anisoDuctile_offset (phases%length), source=0)
allocate(source_damage_anisoDuctile_instance(phases%length), source=0) allocate(source_damage_anisoDuctile_instance(phases%length), source=0)
do p = 1, phases%length do p = 1, phases%length
@ -65,7 +63,6 @@ module function anisoductile_init(source_length) result(mySources)
sources => phase%get('source') sources => phase%get('source')
do sourceOffset = 1, sources%length do sourceOffset = 1, sources%length
if(mySources(sourceOffset,p)) then if(mySources(sourceOffset,p)) then
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)
@ -87,9 +84,9 @@ module function anisoductile_init(source_length) result(mySources)
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'
Nconstituents=count(material_phaseAt==p) * discretization_nIPs Nconstituents=count(material_phaseAt==p) * discretization_nIPs
call phase_allocateState(damageState(p)%p(sourceOffset),Nconstituents,1,1,0) call phase_allocateState(damageState(p),Nconstituents,1,1,0)
damageState(p)%p(sourceOffset)%atol = src%get_asFloat('anisoDuctile_atol',defaultVal=1.0e-3_pReal) damageState(p)%atol = src%get_asFloat('anisoDuctile_atol',defaultVal=1.0e-3_pReal)
if(any(damageState(p)%p(sourceOffset)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' anisoductile_atol' if(any(damageState(p)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' anisoductile_atol'
end associate end associate
@ -116,20 +113,14 @@ module subroutine anisoductile_dotState(co, ip, el)
integer :: & integer :: &
ph, & ph, &
me, & me
sourceOffset, &
damageOffset, &
homog
ph = material_phaseAt(co,el) ph = material_phaseAt(co,el)
me = material_phasememberAt(co,ip,el) me = material_phasememberAt(co,ip,el)
sourceOffset = source_damage_anisoDuctile_offset(ph)
homog = material_homogenizationAt(el)
damageOffset = material_homogenizationMemberAt(ip,el)
associate(prm => param(source_damage_anisoDuctile_instance(ph))) associate(prm => param(source_damage_anisoDuctile_instance(ph)))
damageState(ph)%p(sourceOffset)%dotState(1,me) & damageState(ph)%dotState(1,me) = sum(plasticState(ph)%slipRate(:,me)/(phase_damage_get_phi(co,ip,el)**prm%q)/prm%gamma_crit)
= sum(plasticState(ph)%slipRate(:,me)/(phase_damage_get_phi(co,ip,el)**prm%q)/prm%gamma_crit)
end associate end associate
end subroutine anisoductile_dotState end subroutine anisoductile_dotState
@ -149,12 +140,8 @@ module subroutine source_damage_anisoDuctile_getRateAndItsTangent(localphiDot, d
localphiDot, & localphiDot, &
dLocalphiDot_dPhi dLocalphiDot_dPhi
integer :: &
sourceOffset
sourceOffset = source_damage_anisoDuctile_offset(phase) dLocalphiDot_dPhi = -damageState(phase)%state(1,constituent)
dLocalphiDot_dPhi = -damageState(phase)%p(sourceOffset)%state(1,constituent)
localphiDot = 1.0_pReal & localphiDot = 1.0_pReal &
+ dLocalphiDot_dPhi*phi + dLocalphiDot_dPhi*phi
@ -173,7 +160,7 @@ module subroutine anisoductile_results(phase,group)
integer :: o integer :: o
associate(prm => param(source_damage_anisoDuctile_instance(phase)), & associate(prm => param(source_damage_anisoDuctile_instance(phase)), &
stt => damageState(phase)%p(source_damage_anisoDuctile_offset(phase))%state) stt => damageState(phase)%state)
outputsLoop: do o = 1,size(prm%output) outputsLoop: do o = 1,size(prm%output)
select case(trim(prm%output(o))) select case(trim(prm%output(o)))
case ('f_phi') case ('f_phi')

37
src/phase_damage_isobrittle.f90
View File

@ -7,7 +7,6 @@
submodule(phase:damagee) isobrittle submodule(phase:damagee) isobrittle
integer, dimension(:), allocatable :: & integer, dimension(:), allocatable :: &
source_damage_isoBrittle_offset, &
source_damage_isoBrittle_instance source_damage_isoBrittle_instance
type :: tParameters !< container type for internal constitutive parameters type :: tParameters !< container type for internal constitutive parameters
@ -48,7 +47,6 @@ module function isobrittle_init(source_length) result(mySources)
phases => config_material%get('phase') phases => config_material%get('phase')
allocate(param(Ninstances)) allocate(param(Ninstances))
allocate(source_damage_isoBrittle_offset (phases%length), source=0)
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
@ -58,7 +56,6 @@ module function isobrittle_init(source_length) result(mySources)
sources => phase%get('damage') sources => phase%get('damage')
do sourceOffset = 1, sources%length do sourceOffset = 1, sources%length
if(mySources(sourceOffset,p)) then if(mySources(sourceOffset,p)) then
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)
@ -74,9 +71,9 @@ module function isobrittle_init(source_length) result(mySources)
if (prm%W_crit <= 0.0_pReal) extmsg = trim(extmsg)//' W_crit' if (prm%W_crit <= 0.0_pReal) extmsg = trim(extmsg)//' W_crit'
Nconstituents = count(material_phaseAt==p) * discretization_nIPs Nconstituents = count(material_phaseAt==p) * discretization_nIPs
call phase_allocateState(damageState(p)%p(sourceOffset),Nconstituents,1,1,1) call phase_allocateState(damageState(p),Nconstituents,1,1,1)
damageState(p)%p(sourceOffset)%atol = src%get_asFloat('isoBrittle_atol',defaultVal=1.0e-3_pReal) damageState(p)%atol = src%get_asFloat('isoBrittle_atol',defaultVal=1.0e-3_pReal)
if(any(damageState(p)%p(sourceOffset)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' isobrittle_atol' if(any(damageState(p)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' isobrittle_atol'
end associate end associate
@ -102,14 +99,11 @@ module subroutine source_damage_isoBrittle_deltaState(C, Fe, ph,me)
real(pReal), intent(in), dimension(6,6) :: & real(pReal), intent(in), dimension(6,6) :: &
C C
integer :: &
sourceOffset
real(pReal), dimension(6) :: & real(pReal), dimension(6) :: &
strain strain
real(pReal) :: & real(pReal) :: &
strainenergy strainenergy
sourceOffset = source_damage_isoBrittle_offset(ph)
strain = 0.5_pReal*math_sym33to6(matmul(transpose(Fe),Fe)-math_I3) strain = 0.5_pReal*math_sym33to6(matmul(transpose(Fe),Fe)-math_I3)
@ -117,14 +111,11 @@ module subroutine source_damage_isoBrittle_deltaState(C, Fe, ph,me)
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(ph)%p(sourceOffset)%subState0(1,me)) then if (strainenergy > damageState(ph)%subState0(1,me)) then
damageState(ph)%p(sourceOffset)%deltaState(1,me) = & damageState(ph)%deltaState(1,me) = strainenergy - damageState(ph)%state(1,me)
strainenergy - damageState(ph)%p(sourceOffset)%state(1,me) else
else damageState(ph)%deltaState(1,me) = damageState(ph)%subState0(1,me) - damageState(ph)%state(1,me)
damageState(ph)%p(sourceOffset)%deltaState(1,me) = & endif
damageState(ph)%p(sourceOffset)%subState0(1,me) - &
damageState(ph)%p(sourceOffset)%state(1,me)
endif
end associate end associate
end subroutine source_damage_isoBrittle_deltaState end subroutine source_damage_isoBrittle_deltaState
@ -144,15 +135,11 @@ module subroutine source_damage_isoBrittle_getRateAndItsTangent(localphiDot, dLo
localphiDot, & localphiDot, &
dLocalphiDot_dPhi dLocalphiDot_dPhi
integer :: &
sourceOffset
sourceOffset = source_damage_isoBrittle_offset(phase)
associate(prm => param(source_damage_isoBrittle_instance(phase))) associate(prm => param(source_damage_isoBrittle_instance(phase)))
localphiDot = 1.0_pReal & localphiDot = 1.0_pReal &
- phi*damageState(phase)%p(sourceOffset)%state(1,constituent) - phi*damageState(phase)%state(1,constituent)
dLocalphiDot_dPhi = - damageState(phase)%p(sourceOffset)%state(1,constituent) dLocalphiDot_dPhi = - damageState(phase)%state(1,constituent)
end associate end associate
end subroutine source_damage_isoBrittle_getRateAndItsTangent end subroutine source_damage_isoBrittle_getRateAndItsTangent
@ -169,7 +156,7 @@ module subroutine isobrittle_results(phase,group)
integer :: o integer :: o
associate(prm => param(source_damage_isoBrittle_instance(phase)), & associate(prm => param(source_damage_isoBrittle_instance(phase)), &
stt => damageState(phase)%p(source_damage_isoBrittle_offset(phase))%state) stt => damageState(phase)%state)
outputsLoop: do o = 1,size(prm%output) outputsLoop: do o = 1,size(prm%output)
select case(trim(prm%output(o))) select case(trim(prm%output(o)))
case ('f_phi') case ('f_phi')

29
src/phase_damage_isoductile.f90
View File

@ -7,7 +7,6 @@
submodule(phase:damagee) isoductile submodule(phase:damagee) isoductile
integer, dimension(:), allocatable :: & integer, dimension(:), allocatable :: &
source_damage_isoDuctile_offset, & !< which source is my current damage mechanism?
source_damage_isoDuctile_instance !< instance of damage source mechanism source_damage_isoDuctile_instance !< instance of damage source mechanism
type:: tParameters !< container type for internal constitutive parameters type:: tParameters !< container type for internal constitutive parameters
@ -50,7 +49,6 @@ module function isoductile_init(source_length) result(mySources)
phases => config_material%get('phase') phases => config_material%get('phase')
allocate(param(Ninstances)) allocate(param(Ninstances))
allocate(source_damage_isoDuctile_offset (phases%length), source=0)
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
@ -60,7 +58,6 @@ module function isoductile_init(source_length) result(mySources)
sources => phase%get('damage') sources => phase%get('damage')
do sourceOffset = 1, sources%length do sourceOffset = 1, sources%length
if(mySources(sourceOffset,p)) then if(mySources(sourceOffset,p)) then
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)
@ -78,9 +75,9 @@ module function isoductile_init(source_length) result(mySources)
if (prm%gamma_crit <= 0.0_pReal) extmsg = trim(extmsg)//' gamma_crit' if (prm%gamma_crit <= 0.0_pReal) extmsg = trim(extmsg)//' gamma_crit'
Nconstituents=count(material_phaseAt==p) * discretization_nIPs Nconstituents=count(material_phaseAt==p) * discretization_nIPs
call phase_allocateState(damageState(p)%p(sourceOffset),Nconstituents,1,1,0) call phase_allocateState(damageState(p),Nconstituents,1,1,0)
damageState(p)%p(sourceOffset)%atol = src%get_asFloat('isoDuctile_atol',defaultVal=1.0e-3_pReal) damageState(p)%atol = src%get_asFloat('isoDuctile_atol',defaultVal=1.0e-3_pReal)
if(any(damageState(p)%p(sourceOffset)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' isoductile_atol' if(any(damageState(p)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' isoductile_atol'
end associate end associate
@ -107,20 +104,14 @@ module subroutine isoductile_dotState(co, ip, el)
integer :: & integer :: &
ph, & ph, &
me, & me
sourceOffset, &
damageOffset, &
homog
ph = material_phaseAt(co,el) ph = material_phaseAt(co,el)
me = material_phasememberAt(co,ip,el) me = material_phasememberAt(co,ip,el)
sourceOffset = source_damage_isoDuctile_offset(ph)
homog = material_homogenizationAt(el)
damageOffset = material_homogenizationMemberAt(ip,el)
associate(prm => param(source_damage_isoDuctile_instance(ph))) associate(prm => param(source_damage_isoDuctile_instance(ph)))
damageState(ph)%p(sourceOffset)%dotState(1,me) = & damageState(ph)%dotState(1,me) = sum(plasticState(ph)%slipRate(:,me))/(phase_damage_get_phi(co,ip,el)**prm%q)/prm%gamma_crit
sum(plasticState(ph)%slipRate(:,me))/(phase_damage_get_phi(co,ip,el)**prm%q)/prm%gamma_crit
end associate end associate
end subroutine isoductile_dotState end subroutine isoductile_dotState
@ -140,12 +131,8 @@ module subroutine source_damage_isoDuctile_getRateAndItsTangent(localphiDot, dLo
localphiDot, & localphiDot, &
dLocalphiDot_dPhi dLocalphiDot_dPhi
integer :: &
sourceOffset
sourceOffset = source_damage_isoDuctile_offset(phase) dLocalphiDot_dPhi = -damageState(phase)%state(1,constituent)
dLocalphiDot_dPhi = -damageState(phase)%p(sourceOffset)%state(1,constituent)
localphiDot = 1.0_pReal & localphiDot = 1.0_pReal &
+ dLocalphiDot_dPhi*phi + dLocalphiDot_dPhi*phi
@ -164,7 +151,7 @@ module subroutine isoductile_results(phase,group)
integer :: o integer :: o
associate(prm => param(source_damage_isoDuctile_instance(phase)), & associate(prm => param(source_damage_isoDuctile_instance(phase)), &
stt => damageState(phase)%p(source_damage_isoDuctile_offset(phase))%state) stt => damageState(phase)%state)
outputsLoop: do o = 1,size(prm%output) outputsLoop: do o = 1,size(prm%output)
select case(trim(prm%output(o))) select case(trim(prm%output(o)))
case ('f_phi') case ('f_phi')

18
src/phase_mechanical.f90
View File

@ -1138,7 +1138,7 @@ module function crystallite_stress(dt,co,ip,el) result(converged_)
real(pReal) :: & real(pReal) :: &
formerSubStep formerSubStep
integer :: & integer :: &
so, ph, me, sizeDotState ph, me, sizeDotState
logical :: todo logical :: todo
real(pReal) :: subFrac,subStep real(pReal) :: subFrac,subStep
real(pReal), dimension(3,3) :: & real(pReal), dimension(3,3) :: &
@ -1159,10 +1159,9 @@ module function crystallite_stress(dt,co,ip,el) result(converged_)
subLp0 = phase_mechanical_Lp0(ph)%data(1:3,1:3,me) subLp0 = phase_mechanical_Lp0(ph)%data(1:3,1:3,me)
subState0 = plasticState(ph)%State0(:,me) subState0 = plasticState(ph)%State0(:,me)
if (phase_Nsources(ph) > 0) &
damageState(ph)%subState0(:,me) = damageState(ph)%state0(:,me)
do so = 1, phase_Nsources(ph)
damageState(ph)%p(so)%subState0(:,me) = damageState(ph)%p(so)%state0(:,me)
enddo
subFp0 = phase_mechanical_Fp0(ph)%data(1:3,1:3,me) subFp0 = phase_mechanical_Fp0(ph)%data(1:3,1:3,me)
subFi0 = phase_mechanical_Fi0(ph)%data(1:3,1:3,me) subFi0 = phase_mechanical_Fi0(ph)%data(1:3,1:3,me)
subF0 = phase_mechanical_F0(ph)%data(1:3,1:3,me) subF0 = phase_mechanical_F0(ph)%data(1:3,1:3,me)
@ -1188,9 +1187,9 @@ module function crystallite_stress(dt,co,ip,el) result(converged_)
subFp0 = phase_mechanical_Fp(ph)%data(1:3,1:3,me) subFp0 = phase_mechanical_Fp(ph)%data(1:3,1:3,me)
subFi0 = phase_mechanical_Fi(ph)%data(1:3,1:3,me) subFi0 = phase_mechanical_Fi(ph)%data(1:3,1:3,me)
subState0 = plasticState(ph)%state(:,me) subState0 = plasticState(ph)%state(:,me)
do so = 1, phase_Nsources(ph) if (phase_Nsources(ph) > 0) &
damageState(ph)%p(so)%subState0(:,me) = damageState(ph)%p(so)%state(:,me) damageState(ph)%subState0(:,me) = damageState(ph)%state(:,me)
enddo
endif endif
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! cut back (reduced time and restore) ! cut back (reduced time and restore)
@ -1204,9 +1203,8 @@ module function crystallite_stress(dt,co,ip,el) result(converged_)
phase_mechanical_Li(ph)%data(1:3,1:3,me) = subLi0 phase_mechanical_Li(ph)%data(1:3,1:3,me) = subLi0
endif endif
plasticState(ph)%state(:,me) = subState0 plasticState(ph)%state(:,me) = subState0
do so = 1, phase_Nsources(ph) if (phase_Nsources(ph) > 0) &
damageState(ph)%p(so)%state(:,me) = damageState(ph)%p(so)%subState0(:,me) damageState(ph)%state(:,me) = damageState(ph)%subState0(:,me)
enddo
todo = subStep > num%subStepMinCryst ! still on track or already done (beyond repair) todo = subStep > num%subStepMinCryst ! still on track or already done (beyond repair)
endif endif