Merge branch 'no-partitioned-state' into 'development'

No partitioned state

See merge request damask/DAMASK!327
This commit is contained in:
Franz Roters 2021-01-26 19:37:10 +01:00
commit a0a7009143
12 changed files with 260 additions and 329 deletions

View File

@ -115,12 +115,6 @@ module constitutive
integer, intent(in) :: ph
end subroutine damage_results
module subroutine mech_initializeRestorationPoints(ph,me)
integer, intent(in) :: ph, me
end subroutine mech_initializeRestorationPoints
module subroutine mech_windForward(ph,me)
integer, intent(in) :: ph, me
end subroutine mech_windForward
@ -359,7 +353,6 @@ module constitutive
constitutive_mech_getP, &
constitutive_mech_setF, &
constitutive_mech_getF, &
constitutive_initializeRestorationPoints, &
constitutive_windForward, &
KINEMATICS_UNDEFINED_ID ,&
KINEMATICS_CLEAVAGE_OPENING_ID, &
@ -404,11 +397,9 @@ subroutine constitutive_init
PhaseLoop2:do ph = 1,phases%length
!--------------------------------------------------------------------------------------------------
! partition and initialize state
plasticState(ph)%partitionedState0 = plasticState(ph)%state0
plasticState(ph)%state = plasticState(ph)%partitionedState0
plasticState(ph)%state = plasticState(ph)%state0
forall(so = 1:phase_Nsources(ph))
damageState(ph)%p(so)%partitionedState0 = damageState(ph)%p(so)%state0
damageState(ph)%p(so)%state = damageState(ph)%p(so)%partitionedState0
damageState(ph)%p(so)%state = damageState(ph)%p(so)%state0
end forall
constitutive_source_maxSizeDotState = max(constitutive_source_maxSizeDotState, &
@ -473,7 +464,6 @@ subroutine constitutive_allocateState(state, &
allocate(state%atol (sizeState), source=0.0_pReal)
allocate(state%state0 (sizeState,Nconstituents), source=0.0_pReal)
allocate(state%partitionedState0(sizeState,Nconstituents), source=0.0_pReal)
allocate(state%state (sizeState,Nconstituents), source=0.0_pReal)
allocate(state%dotState (sizeDotState,Nconstituents), source=0.0_pReal)
@ -499,8 +489,8 @@ subroutine constitutive_restore(ce,includeL)
do co = 1,homogenization_Nconstituents(material_homogenizationAt2(ce))
do so = 1, phase_Nsources(material_phaseAt2(co,ce))
damageState(material_phaseAt2(co,ce))%p(so)%state( :,material_phasememberAt2(co,ce)) = &
damageState(material_phaseAt2(co,ce))%p(so)%partitionedState0(:,material_phasememberAt2(co,ce))
damageState(material_phaseAt2(co,ce))%p(so)%state( :,material_phasememberAt2(co,ce)) = &
damageState(material_phaseAt2(co,ce))%p(so)%state0(:,material_phasememberAt2(co,ce))
enddo
enddo
@ -569,7 +559,6 @@ subroutine crystallite_init()
iMax, & !< maximum number of integration points
eMax !< maximum number of elements
class(tNode), pointer :: &
num_crystallite, &
debug_crystallite, & ! pointer to debug options for crystallite
@ -651,33 +640,6 @@ subroutine crystallite_init()
end subroutine crystallite_init
!--------------------------------------------------------------------------------------------------
!> @brief Backup data for homog cutback.
!--------------------------------------------------------------------------------------------------
subroutine constitutive_initializeRestorationPoints(ip,el)
integer, intent(in) :: &
ip, & !< integration point number
el !< element number
integer :: &
co, & !< constituent number
so,ph, me
do co = 1,homogenization_Nconstituents(material_homogenizationAt(el))
ph = material_phaseAt(co,el)
me = material_phaseMemberAt(co,ip,el)
call mech_initializeRestorationPoints(ph,me)
do so = 1, size(damageState(ph)%p)
damageState(ph)%p(so)%partitionedState0(:,me) = damageState(ph)%p(so)%state0(:,me)
enddo
enddo
end subroutine constitutive_initializeRestorationPoints
!--------------------------------------------------------------------------------------------------
!> @brief Wind homog inc forward.
!--------------------------------------------------------------------------------------------------
@ -699,7 +661,7 @@ subroutine constitutive_windForward(ip,el)
call mech_windForward(ph,me)
do so = 1, phase_Nsources(material_phaseAt(co,el))
damageState(ph)%p(so)%partitionedState0(:,me) = damageState(ph)%p(so)%state(:,me)
damageState(ph)%p(so)%state0(:,me) = damageState(ph)%p(so)%state(:,me)
enddo
enddo

View File

@ -39,16 +39,7 @@ submodule(constitutive) constitutive_mech
constitutive_mech_F0, &
constitutive_mech_Li0, &
constitutive_mech_Lp0, &
constitutive_mech_S0, &
! converged value at end of last homogenization increment (RGC only)
constitutive_mech_partitionedFi0, &
constitutive_mech_partitionedFp0, &
constitutive_mech_partitionedF0, &
constitutive_mech_partitionedLi0, &
constitutive_mech_partitionedLp0, &
constitutive_mech_partitionedS0
constitutive_mech_S0
integer(kind(PLASTICITY_undefined_ID)), dimension(:), allocatable :: &
@ -361,23 +352,17 @@ module subroutine mech_init(phases)
allocate(constitutive_mech_Fe(phases%length))
allocate(constitutive_mech_Fi(phases%length))
allocate(constitutive_mech_Fi0(phases%length))
allocate(constitutive_mech_partitionedFi0(phases%length))
allocate(constitutive_mech_Fp(phases%length))
allocate(constitutive_mech_Fp0(phases%length))
allocate(constitutive_mech_partitionedFp0(phases%length))
allocate(constitutive_mech_F(phases%length))
allocate(constitutive_mech_F0(phases%length))
allocate(constitutive_mech_partitionedF0(phases%length))
allocate(constitutive_mech_Li(phases%length))
allocate(constitutive_mech_Li0(phases%length))
allocate(constitutive_mech_partitionedLi0(phases%length))
allocate(constitutive_mech_partitionedLp0(phases%length))
allocate(constitutive_mech_Lp0(phases%length))
allocate(constitutive_mech_Lp(phases%length))
allocate(constitutive_mech_S(phases%length))
allocate(constitutive_mech_P(phases%length))
allocate(constitutive_mech_S0(phases%length))
allocate(constitutive_mech_partitionedS0(phases%length))
do ph = 1, phases%length
Nconstituents = count(material_phaseAt == ph) * discretization_nIPs
@ -385,23 +370,17 @@ module subroutine mech_init(phases)
allocate(constitutive_mech_Fi(ph)%data(3,3,Nconstituents))
allocate(constitutive_mech_Fe(ph)%data(3,3,Nconstituents))
allocate(constitutive_mech_Fi0(ph)%data(3,3,Nconstituents))
allocate(constitutive_mech_partitionedFi0(ph)%data(3,3,Nconstituents))
allocate(constitutive_mech_Fp(ph)%data(3,3,Nconstituents))
allocate(constitutive_mech_Fp0(ph)%data(3,3,Nconstituents))
allocate(constitutive_mech_partitionedFp0(ph)%data(3,3,Nconstituents))
allocate(constitutive_mech_Li(ph)%data(3,3,Nconstituents))
allocate(constitutive_mech_Li0(ph)%data(3,3,Nconstituents))
allocate(constitutive_mech_partitionedLi0(ph)%data(3,3,Nconstituents))
allocate(constitutive_mech_partitionedLp0(ph)%data(3,3,Nconstituents))
allocate(constitutive_mech_Lp0(ph)%data(3,3,Nconstituents))
allocate(constitutive_mech_Lp(ph)%data(3,3,Nconstituents))
allocate(constitutive_mech_S(ph)%data(3,3,Nconstituents),source=0.0_pReal)
allocate(constitutive_mech_P(ph)%data(3,3,Nconstituents),source=0.0_pReal)
allocate(constitutive_mech_S0(ph)%data(3,3,Nconstituents),source=0.0_pReal)
allocate(constitutive_mech_partitionedS0(ph)%data(3,3,Nconstituents),source=0.0_pReal)
allocate(constitutive_mech_F(ph)%data(3,3,Nconstituents))
allocate(constitutive_mech_F0(ph)%data(3,3,Nconstituents))
allocate(constitutive_mech_partitionedF0(ph)%data(3,3,Nconstituents))
phase => phases%get(ph)
mech => phase%get('mechanics')
@ -454,10 +433,6 @@ module subroutine mech_init(phases)
constitutive_mech_Fi(ph)%data(1:3,1:3,me) = constitutive_mech_Fi0(ph)%data(1:3,1:3,me)
constitutive_mech_F(ph)%data(1:3,1:3,me) = constitutive_mech_F0(ph)%data(1:3,1:3,me)
constitutive_mech_partitionedFi0(ph)%data(1:3,1:3,me) = constitutive_mech_Fi0(ph)%data(1:3,1:3,me)
constitutive_mech_partitionedFp0(ph)%data(1:3,1:3,me) = constitutive_mech_Fp0(ph)%data(1:3,1:3,me)
constitutive_mech_partitionedF0(ph)%data(1:3,1:3,me) = constitutive_mech_F0(ph)%data(1:3,1:3,me)
enddo
enddo; enddo
!$OMP END PARALLEL DO
@ -1464,26 +1439,6 @@ subroutine crystallite_results(group,ph)
end subroutine crystallite_results
!--------------------------------------------------------------------------------------------------
!> @brief Backup data for homog cutback.
!--------------------------------------------------------------------------------------------------
module subroutine mech_initializeRestorationPoints(ph,me)
integer, intent(in) :: ph, me
constitutive_mech_partitionedFi0(ph)%data(1:3,1:3,me) = constitutive_mech_Fi0(ph)%data(1:3,1:3,me)
constitutive_mech_partitionedFp0(ph)%data(1:3,1:3,me) = constitutive_mech_Fp0(ph)%data(1:3,1:3,me)
constitutive_mech_partitionedF0(ph)%data(1:3,1:3,me) = constitutive_mech_F0(ph)%data(1:3,1:3,me)
constitutive_mech_partitionedLi0(ph)%data(1:3,1:3,me) = constitutive_mech_Li0(ph)%data(1:3,1:3,me)
constitutive_mech_partitionedLp0(ph)%data(1:3,1:3,me) = constitutive_mech_Lp0(ph)%data(1:3,1:3,me)
constitutive_mech_partitionedS0(ph)%data(1:3,1:3,me) = constitutive_mech_S0(ph)%data(1:3,1:3,me)
plasticState(ph)%partitionedState0(:,me) = plasticState(ph)%state0(:,me)
end subroutine mech_initializeRestorationPoints
!--------------------------------------------------------------------------------------------------
!> @brief Wind homog inc forward.
!--------------------------------------------------------------------------------------------------
@ -1492,14 +1447,14 @@ module subroutine mech_windForward(ph,me)
integer, intent(in) :: ph, me
constitutive_mech_partitionedFp0(ph)%data(1:3,1:3,me) = constitutive_mech_Fp(ph)%data(1:3,1:3,me)
constitutive_mech_partitionedFi0(ph)%data(1:3,1:3,me) = constitutive_mech_Fi(ph)%data(1:3,1:3,me)
constitutive_mech_partitionedF0(ph)%data(1:3,1:3,me) = constitutive_mech_F(ph)%data(1:3,1:3,me)
constitutive_mech_partitionedLi0(ph)%data(1:3,1:3,me) = constitutive_mech_Li(ph)%data(1:3,1:3,me)
constitutive_mech_partitionedLp0(ph)%data(1:3,1:3,me) = constitutive_mech_Lp(ph)%data(1:3,1:3,me)
constitutive_mech_partitionedS0(ph)%data(1:3,1:3,me) = constitutive_mech_S(ph)%data(1:3,1:3,me)
constitutive_mech_Fp0(ph)%data(1:3,1:3,me) = constitutive_mech_Fp(ph)%data(1:3,1:3,me)
constitutive_mech_Fi0(ph)%data(1:3,1:3,me) = constitutive_mech_Fi(ph)%data(1:3,1:3,me)
constitutive_mech_F0(ph)%data(1:3,1:3,me) = constitutive_mech_F(ph)%data(1:3,1:3,me)
constitutive_mech_Li0(ph)%data(1:3,1:3,me) = constitutive_mech_Li(ph)%data(1:3,1:3,me)
constitutive_mech_Lp0(ph)%data(1:3,1:3,me) = constitutive_mech_Lp(ph)%data(1:3,1:3,me)
constitutive_mech_S0(ph)%data(1:3,1:3,me) = constitutive_mech_S(ph)%data(1:3,1:3,me)
plasticState(ph)%partitionedState0(:,me) = plasticState(ph)%state(:,me)
plasticState(ph)%State0(:,me) = plasticState(ph)%state(:,me)
end subroutine mech_windForward
@ -1578,17 +1533,17 @@ module function crystallite_stress(dt,co,ip,el) result(converged_)
me = material_phaseMemberAt(co,ip,el)
sizeDotState = plasticState(ph)%sizeDotState
subLi0 = constitutive_mech_partitionedLi0(ph)%data(1:3,1:3,me)
subLp0 = constitutive_mech_partitionedLp0(ph)%data(1:3,1:3,me)
subState0 = plasticState(ph)%partitionedState0(:,me)
subLi0 = constitutive_mech_Li0(ph)%data(1:3,1:3,me)
subLp0 = constitutive_mech_Lp0(ph)%data(1:3,1:3,me)
subState0 = plasticState(ph)%State0(:,me)
do so = 1, phase_Nsources(ph)
damageState(ph)%p(so)%subState0(:,me) = damageState(ph)%p(so)%partitionedState0(:,me)
damageState(ph)%p(so)%subState0(:,me) = damageState(ph)%p(so)%state0(:,me)
enddo
subFp0 = constitutive_mech_partitionedFp0(ph)%data(1:3,1:3,me)
subFi0 = constitutive_mech_partitionedFi0(ph)%data(1:3,1:3,me)
subF0 = constitutive_mech_partitionedF0(ph)%data(1:3,1:3,me)
subFp0 = constitutive_mech_Fp0(ph)%data(1:3,1:3,me)
subFi0 = constitutive_mech_Fi0(ph)%data(1:3,1:3,me)
subF0 = constitutive_mech_F0(ph)%data(1:3,1:3,me)
subFrac = 0.0_pReal
subStep = 1.0_pReal/num%subStepSizeCryst
todo = .true.
@ -1638,7 +1593,7 @@ module function crystallite_stress(dt,co,ip,el) result(converged_)
! prepare for integration
if (todo) then
subF = subF0 &
+ subStep * (constitutive_mech_F(ph)%data(1:3,1:3,me) - constitutive_mech_partitionedF0(ph)%data(1:3,1:3,me))
+ subStep * (constitutive_mech_F(ph)%data(1:3,1:3,me) - constitutive_mech_F0(ph)%data(1:3,1:3,me))
constitutive_mech_Fe(ph)%data(1:3,1:3,me) = matmul(subF,math_inv33(matmul(constitutive_mech_Fi(ph)%data(1:3,1:3,me), &
constitutive_mech_Fp(ph)%data(1:3,1:3,me))))
converged_ = .not. integrateState(subF0,subF,subFp0,subFi0,subState0(1:sizeDotState),subStep * dt,co,ip,el)
@ -1667,15 +1622,15 @@ module subroutine mech_restore(ce,includeL)
ph = material_phaseAt2(co,ce)
me = material_phaseMemberAt2(co,ce)
if (includeL) then
constitutive_mech_Lp(ph)%data(1:3,1:3,me) = constitutive_mech_partitionedLp0(ph)%data(1:3,1:3,me)
constitutive_mech_Li(ph)%data(1:3,1:3,me) = constitutive_mech_partitionedLi0(ph)%data(1:3,1:3,me)
constitutive_mech_Lp(ph)%data(1:3,1:3,me) = constitutive_mech_Lp0(ph)%data(1:3,1:3,me)
constitutive_mech_Li(ph)%data(1:3,1:3,me) = constitutive_mech_Li0(ph)%data(1:3,1:3,me)
endif ! maybe protecting everything from overwriting makes more sense
constitutive_mech_Fp(ph)%data(1:3,1:3,me) = constitutive_mech_partitionedFp0(ph)%data(1:3,1:3,me)
constitutive_mech_Fi(ph)%data(1:3,1:3,me) = constitutive_mech_partitionedFi0(ph)%data(1:3,1:3,me)
constitutive_mech_S(ph)%data(1:3,1:3,me) = constitutive_mech_partitionedS0(ph)%data(1:3,1:3,me)
constitutive_mech_Fp(ph)%data(1:3,1:3,me) = constitutive_mech_Fp0(ph)%data(1:3,1:3,me)
constitutive_mech_Fi(ph)%data(1:3,1:3,me) = constitutive_mech_Fi0(ph)%data(1:3,1:3,me)
constitutive_mech_S(ph)%data(1:3,1:3,me) = constitutive_mech_S0(ph)%data(1:3,1:3,me)
plasticState(ph)%state(:,me) = plasticState(ph)%partitionedState0(:,me)
plasticState(ph)%state(:,me) = plasticState(ph)%State0(:,me)
enddo
end subroutine mech_restore
@ -1727,8 +1682,8 @@ module function constitutive_mech_dPdF(dt,co,ip,el) result(dPdF)
invFp = math_inv33(constitutive_mech_Fp(ph)%data(1:3,1:3,me))
invFi = math_inv33(constitutive_mech_Fi(ph)%data(1:3,1:3,me))
invSubFp0 = math_inv33(constitutive_mech_partitionedFp0(ph)%data(1:3,1:3,me))
invSubFi0 = math_inv33(constitutive_mech_partitionedFi0(ph)%data(1:3,1:3,me))
invSubFp0 = math_inv33(constitutive_mech_Fp0(ph)%data(1:3,1:3,me))
invSubFi0 = math_inv33(constitutive_mech_Fi0(ph)%data(1:3,1:3,me))
if (sum(abs(dLidS)) < tol_math_check) then
dFidS = 0.0_pReal

View File

@ -116,12 +116,11 @@ module subroutine thermal_init(phases)
PhaseLoop2:do ph = 1,phases%length
do so = 1,thermal_Nsources(ph)
deallocate(thermalState(ph)%p(so)%partitionedState0)
thermalState(ph)%p(so)%state = thermalState(ph)%p(so)%state0
thermalState(ph)%p(so)%state = thermalState(ph)%p(so)%state0
enddo
thermal_source_maxSizeDotState = max(thermal_source_maxSizeDotState, &
maxval(thermalState(ph)%p%sizeDotState))
thermal_source_maxSizeDotState = max(thermal_source_maxSizeDotState, &
maxval(thermalState(ph)%p%sizeDotState))
enddo PhaseLoop2
!--------------------------------------------------------------------------------------------------

View File

@ -27,7 +27,6 @@ subroutine damage_none_init
Nmaterialpoints = count(material_homogenizationAt == h)
damageState_h(h)%sizeState = 0
allocate(damageState_h(h)%state0 (0,Nmaterialpoints))
allocate(damageState_h(h)%subState0(0,Nmaterialpoints))
allocate(damageState_h(h)%state (0,Nmaterialpoints))
allocate (damage(h)%p(Nmaterialpoints), source=1.0_pReal)

View File

@ -14,28 +14,17 @@ module damage_nonlocal
implicit none
private
type :: tParameters
character(len=pStringLen), allocatable, dimension(:) :: &
output
end type tParameters
type, private :: tNumerics
real(pReal) :: &
charLength !< characteristic length scale for gradient problems
end type tNumerics
type(tparameters), dimension(:), allocatable :: &
param
type(tNumerics), private :: &
num
public :: &
damage_nonlocal_init, &
damage_nonlocal_getSourceAndItsTangent, &
damage_nonlocal_getDiffusion, &
damage_nonlocal_getMobility, &
damage_nonlocal_putNonLocalDamage, &
damage_nonlocal_results
damage_nonlocal_getDiffusion
contains
@ -48,9 +37,7 @@ subroutine damage_nonlocal_init
integer :: Ninstances,Nmaterialpoints,h
class(tNode), pointer :: &
num_generic, &
material_homogenization, &
homog, &
homogDamage
material_homogenization
print'(/,a)', ' <<<+- damage_nonlocal init -+>>>'; flush(6)
@ -60,58 +47,23 @@ subroutine damage_nonlocal_init
num%charLength = num_generic%get_asFloat('charLength',defaultVal=1.0_pReal)
Ninstances = count(damage_type == DAMAGE_nonlocal_ID)
allocate(param(Ninstances))
material_homogenization => config_material%get('homogenization')
do h = 1, material_homogenization%length
if (damage_type(h) /= DAMAGE_NONLOCAL_ID) cycle
homog => material_homogenization%get(h)
homogDamage => homog%get('damage')
associate(prm => param(damage_typeInstance(h)))
#if defined (__GFORTRAN__)
prm%output = output_asStrings(homogDamage)
#else
prm%output = homogDamage%get_asStrings('output',defaultVal=emptyStringArray)
#endif
Nmaterialpoints = count(material_homogenizationAt == h)
damageState_h(h)%sizeState = 1
allocate(damageState_h(h)%state0 (1,Nmaterialpoints), source=1.0_pReal)
allocate(damageState_h(h)%subState0(1,Nmaterialpoints), source=1.0_pReal)
allocate(damageState_h(h)%state (1,Nmaterialpoints), source=1.0_pReal)
damage(h)%p => damageState_h(h)%state(1,:)
end associate
enddo
end subroutine damage_nonlocal_init
!--------------------------------------------------------------------------------------------------
!> @brief calculates homogenized damage driving forces
!--------------------------------------------------------------------------------------------------
subroutine damage_nonlocal_getSourceAndItsTangent(phiDot, dPhiDot_dPhi, phi, ip, el)
integer, intent(in) :: &
ip, & !< integration point number
el !< element number
real(pReal), intent(in) :: &
phi
real(pReal) :: &
phiDot, dPhiDot_dPhi
phiDot = 0.0_pReal
dPhiDot_dPhi = 0.0_pReal
call constitutive_damage_getRateAndItsTangents(phiDot, dPhiDot_dPhi, phi, ip, el)
phiDot = phiDot/real(homogenization_Nconstituents(material_homogenizationAt(el)),pReal)
dPhiDot_dPhi = dPhiDot_dPhi/real(homogenization_Nconstituents(material_homogenizationAt(el)),pReal)
end subroutine damage_nonlocal_getSourceAndItsTangent
!--------------------------------------------------------------------------------------------------
!> @brief returns homogenized non local damage diffusion tensor in reference configuration
!--------------------------------------------------------------------------------------------------
@ -139,70 +91,4 @@ function damage_nonlocal_getDiffusion(ip,el)
end function damage_nonlocal_getDiffusion
!--------------------------------------------------------------------------------------------------
!> @brief Returns homogenized nonlocal damage mobility
!--------------------------------------------------------------------------------------------------
real(pReal) function damage_nonlocal_getMobility(ip,el)
integer, intent(in) :: &
ip, & !< integration point number
el !< element number
integer :: &
co
damage_nonlocal_getMobility = 0.0_pReal
do co = 1, homogenization_Nconstituents(material_homogenizationAt(el))
damage_nonlocal_getMobility = damage_nonlocal_getMobility + lattice_M(material_phaseAt(co,el))
enddo
damage_nonlocal_getMobility = damage_nonlocal_getMobility/&
real(homogenization_Nconstituents(material_homogenizationAt(el)),pReal)
end function damage_nonlocal_getMobility
!--------------------------------------------------------------------------------------------------
!> @brief updated nonlocal damage field with solution from damage phase field PDE
!--------------------------------------------------------------------------------------------------
subroutine damage_nonlocal_putNonLocalDamage(phi,ip,el)
integer, intent(in) :: &
ip, & !< integration point number
el !< element number
real(pReal), intent(in) :: &
phi
integer :: &
homog, &
offset
homog = material_homogenizationAt(el)
offset = material_homogenizationMemberAt(ip,el)
damage(homog)%p(offset) = phi
end subroutine damage_nonlocal_putNonLocalDamage
!--------------------------------------------------------------------------------------------------
!> @brief writes results to HDF5 output file
!--------------------------------------------------------------------------------------------------
subroutine damage_nonlocal_results(homog,group)
integer, intent(in) :: homog
character(len=*), intent(in) :: group
integer :: o
associate(prm => param(damage_typeInstance(homog)))
outputsLoop: do o = 1,size(prm%output)
select case(prm%output(o))
case ('phi')
call results_writeDataset(group,damage(homog)%p,prm%output(o),&
'damage indicator','-')
end select
enddo outputsLoop
end associate
end subroutine damage_nonlocal_results
end module damage_nonlocal

View File

@ -198,7 +198,6 @@ function grid_damage_spectral_solution(timeinc) result(solution)
do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1)
cell = cell + 1
call damage_nonlocal_putNonLocalDamage(phi_current(i,j,k),1,cell)
homogenization_phi(cell) = phi_current(i,j,k)
enddo; enddo; enddo
call VecMin(solution_vec,devNull,phi_min,ierr); CHKERRQ(ierr)
@ -236,7 +235,6 @@ subroutine grid_damage_spectral_forward(cutBack)
do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1)
cell = cell + 1
call damage_nonlocal_putNonLocalDamage(phi_current(i,j,k),1,cell)
homogenization_phi(cell) = phi_current(i,j,k)
enddo; enddo; enddo
else
phi_lastInc = phi_current

View File

@ -25,9 +25,6 @@ module homogenization
!--------------------------------------------------------------------------------------------------
! General variables for the homogenization at a material point
real(pReal), dimension(:), allocatable, public :: &
homogenization_phi, &
homogenization_dot_phi
real(pReal), dimension(:,:,:), allocatable, public :: &
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
@ -41,10 +38,6 @@ module homogenization
type :: tNumerics
integer :: &
nMPstate !< materialpoint state loop limit
real(pReal) :: &
subStepMinHomog, & !< minimum (relative) size of sub-step allowed during cutback in homogenization
subStepSizeHomog, & !< size of first substep when cutback in homogenization
stepIncreaseHomog !< increase of next substep size when previous substep converged in homogenization
end type tNumerics
type(tNumerics) :: num
@ -75,8 +68,7 @@ module homogenization
integer, intent(in) :: ce
end subroutine thermal_partition
module subroutine damage_partition(phi,ce)
real(pReal), intent(in) :: phi
module subroutine damage_partition(ce)
integer, intent(in) :: ce
end subroutine damage_partition
@ -149,6 +141,43 @@ module homogenization
real(pReal), intent(out) :: Tdot
end subroutine thermal_conduction_getSource
module function damage_nonlocal_getMobility(ip,el) result(M)
integer, intent(in) :: &
ip, & !< integration point number
el !< element number
integer :: &
co
real(pReal) :: M
end function damage_nonlocal_getMobility
module subroutine damage_nonlocal_getSourceAndItsTangent(phiDot, dPhiDot_dPhi, phi, ip, el)
integer, intent(in) :: &
ip, & !< integration point number
el !< element number
real(pReal), intent(in) :: &
phi
real(pReal) :: &
phiDot, dPhiDot_dPhi
end subroutine damage_nonlocal_getSourceAndItsTangent
module subroutine damage_nonlocal_putNonLocalDamage(phi,ip,el)
integer, intent(in) :: &
ip, & !< integration point number
el !< element number
real(pReal), intent(in) :: &
phi
end subroutine damage_nonlocal_putNonLocalDamage
module subroutine damage_nonlocal_results(homog,group)
integer, intent(in) :: homog
character(len=*), intent(in) :: group
end subroutine damage_nonlocal_results
end interface
public :: &
@ -158,6 +187,9 @@ module homogenization
thermal_conduction_getConductivity, &
thermal_conduction_getMassDensity, &
thermal_conduction_getSource, &
damage_nonlocal_getMobility, &
damage_nonlocal_getSourceAndItsTangent, &
damage_nonlocal_putNonLocalDamage, &
homogenization_thermal_setfield, &
homogenization_thermal_T, &
homogenization_forward, &
@ -183,14 +215,7 @@ subroutine homogenization_init()
num_homogGeneric => num_homog%get('generic',defaultVal=emptyDict)
num%nMPstate = num_homogGeneric%get_asInt ('nMPstate', defaultVal=10)
num%subStepMinHomog = num_homogGeneric%get_asFloat('subStepMin', defaultVal=1.0e-3_pReal)
num%subStepSizeHomog = num_homogGeneric%get_asFloat('subStepSize', defaultVal=0.25_pReal)
num%stepIncreaseHomog = num_homogGeneric%get_asFloat('stepIncrease', defaultVal=1.5_pReal)
if (num%nMPstate < 1) call IO_error(301,ext_msg='nMPstate')
if (num%subStepMinHomog <= 0.0_pReal) call IO_error(301,ext_msg='subStepMinHomog')
if (num%subStepSizeHomog <= 0.0_pReal) call IO_error(301,ext_msg='subStepSizeHomog')
if (num%stepIncreaseHomog <= 0.0_pReal) call IO_error(301,ext_msg='stepIncreaseHomog')
call mech_init(num_homog)
@ -216,16 +241,13 @@ subroutine materialpoint_stressAndItsTangent(dt,FEsolving_execIP,FEsolving_execE
ip, & !< integration point number
el, & !< element number
myNgrains, co, ce, ho, me, ph
real(pReal) :: &
subFrac, &
subStep
logical :: &
converged
logical, dimension(2) :: &
doneAndHappy
!$OMP PARALLEL
!$OMP DO PRIVATE(ce,me,ho,myNgrains,NiterationMPstate,subFrac,converged,subStep,doneAndHappy)
!$OMP DO PRIVATE(ce,me,ho,myNgrains,NiterationMPstate,converged,doneAndHappy)
do el = FEsolving_execElem(1),FEsolving_execElem(2)
ho = material_homogenizationAt(el)
myNgrains = homogenization_Nconstituents(ho)
@ -233,78 +255,39 @@ subroutine materialpoint_stressAndItsTangent(dt,FEsolving_execIP,FEsolving_execE
ce = (el-1)*discretization_nIPs + ip
me = material_homogenizationMemberAt2(ce)
call constitutive_initializeRestorationPoints(ip,el)
call constitutive_restore(ce,.false.) ! wrong name (is more a forward function)
subFrac = 0.0_pReal
converged = .false. ! pretend failed step ...
subStep = 1.0_pReal/num%subStepSizeHomog ! ... larger then the requested calculation
if(homogState(ho)%sizeState > 0) homogState(ho)%State(:,me) = homogState(ho)%State0(:,me)
if(damageState_h(ho)%sizeState > 0) damageState_h(ho)%State(:,me) = damageState_h(ho)%State0(:,me)
if (homogState(ho)%sizeState > 0) homogState(ho)%subState0(:,me) = homogState(ho)%State0(:,me)
if (damageState_h(ho)%sizeState > 0) damageState_h(ho)%subState0(:,me) = damageState_h(ho)%State0(:,me)
doneAndHappy = [.false.,.true.]
cutBackLooping: do while (.not. terminallyIll .and. subStep > num%subStepMinHomog)
NiterationMPstate = 0
convergenceLooping: do while (.not. (terminallyIll .or. doneAndHappy(1)) &
.and. NiterationMPstate < num%nMPstate)
NiterationMPstate = NiterationMPstate + 1
if (converged) then
subFrac = subFrac + subStep
subStep = min(1.0_pReal-subFrac,num%stepIncreaseHomog*subStep) ! introduce flexibility for step increase/acceleration
steppingNeeded: if (subStep > num%subStepMinHomog) then
if (.not. doneAndHappy(1)) then
call mech_partition(homogenization_F(1:3,1:3,ce),ip,el)
converged = .true.
do co = 1, myNgrains
converged = converged .and. crystallite_stress(dt,co,ip,el)
enddo
! wind forward grain starting point
call constitutive_windForward(ip,el)
if(homogState(ho)%sizeState > 0) homogState(ho)%subState0(:,me) = homogState(ho)%State(:,me)
if(damageState_h(ho)%sizeState > 0) damageState_h(ho)%subState0(:,me) = damageState_h(ho)%State(:,me)
endif steppingNeeded
elseif ( (myNgrains == 1 .and. subStep <= 1.0 ) .or. & ! single grain already tried internal subStepping in crystallite
num%subStepSizeHomog * subStep <= num%subStepMinHomog ) then ! would require too small subStep
! cutback makes no sense
if (.not. terminallyIll) & ! so first signals terminally ill...
print*, ' Integration point ', ip,' at element ', el, ' terminally ill'
terminallyIll = .true. ! ...and kills all others
else ! cutback makes sense
subStep = num%subStepSizeHomog * subStep ! crystallite had severe trouble, so do a significant cutback
call constitutive_restore(ce,subStep < 1.0_pReal)
if(homogState(ho)%sizeState > 0) homogState(ho)%State(:,me) = homogState(ho)%subState0(:,me)
if(damageState_h(ho)%sizeState > 0) damageState_h(ho)%State(:,me) = damageState_h(ho)%subState0(:,me)
if (.not. converged) then
doneAndHappy = [.true.,.false.]
else
doneAndHappy = mech_updateState(dt,homogenization_F(1:3,1:3,ce),ip,el)
converged = all(doneAndHappy)
endif
endif
if (subStep > num%subStepMinHomog) doneAndHappy = [.false.,.true.]
NiterationMPstate = 0
convergenceLooping: do while (.not. (terminallyIll .or. doneAndHappy(1)) &
.and. NiterationMPstate < num%nMPstate)
NiterationMPstate = NiterationMPstate + 1
!--------------------------------------------------------------------------------------------------
! deformation partitioning
if (.not. doneAndHappy(1)) then
call mech_partition( homogenization_F0(1:3,1:3,ce) &
+ (homogenization_F(1:3,1:3,ce)-homogenization_F0(1:3,1:3,ce))*(subStep+subFrac), &
ip,el)
converged = .true.
do co = 1, myNgrains
converged = converged .and. crystallite_stress(dt*subStep,co,ip,el)
enddo
if (.not. converged) then
doneAndHappy = [.true.,.false.]
else
doneAndHappy = mech_updateState(dt*subStep, &
homogenization_F0(1:3,1:3,ce) &
+ (homogenization_F(1:3,1:3,ce)-homogenization_F0(1:3,1:3,ce)) &
*(subStep+subFrac), &
ip,el)
converged = all(doneAndHappy)
endif
endif
enddo convergenceLooping
enddo cutBackLooping
enddo convergenceLooping
if (.not. converged) then
if (.not. terminallyIll) print*, ' Integration point ', ip,' at element ', el, ' terminally ill'
terminallyIll = .true.
endif
enddo
enddo
!$OMP END DO
@ -330,6 +313,26 @@ subroutine materialpoint_stressAndItsTangent(dt,FEsolving_execIP,FEsolving_execE
enddo
!$OMP END DO
! !$OMP DO PRIVATE(ho,ph,ce)
! do el = FEsolving_execElem(1),FEsolving_execElem(2)
! if (terminallyIll) continue
! ho = material_homogenizationAt(el)
! do ip = FEsolving_execIP(1),FEsolving_execIP(2)
! ce = (el-1)*discretization_nIPs + ip
! call damage_partition(ce)
! do co = 1, homogenization_Nconstituents(ho)
! ph = material_phaseAt(co,el)
! if (.not. thermal_stress(dt,ph,material_phaseMemberAt(co,ip,el))) then
! if (.not. terminallyIll) & ! so first signals terminally ill...
! print*, ' Integration point ', ip,' at element ', el, ' terminally ill'
! terminallyIll = .true. ! ...and kills all others
! endif
! call thermal_homogenize(ip,el)
! enddo
! enddo
! enddo
! !$OMP END DO
!$OMP DO PRIVATE(ho)
elementLooping3: do el = FEsolving_execElem(1),FEsolving_execElem(2)
ho = material_homogenizationAt(el)

View File

@ -3,19 +3,60 @@
!--------------------------------------------------------------------------------------------------
submodule(homogenization) homogenization_damage
use lattice
type :: tDataContainer
real(pReal), dimension(:), allocatable :: phi
end type tDataContainer
type(tDataContainer), dimension(:), allocatable :: current
type :: tParameters
character(len=pStringLen), allocatable, dimension(:) :: &
output
end type tParameters
type(tparameters), dimension(:), allocatable :: &
param
contains
!--------------------------------------------------------------------------------------------------
!> @brief Allocate variables and set parameters.
!--------------------------------------------------------------------------------------------------
module subroutine damage_init()
class(tNode), pointer :: &
configHomogenizations, &
configHomogenization, &
configHomogenizationDamage
integer :: ho
print'(/,a)', ' <<<+- homogenization_damage init -+>>>'
allocate(homogenization_phi(discretization_nIPs*discretization_Nelems))
allocate(homogenization_dot_phi(discretization_nIPs*discretization_Nelems))
configHomogenizations => config_material%get('homogenization')
allocate(param(configHomogenizations%length))
allocate(current(configHomogenizations%length))
do ho = 1, configHomogenizations%length
allocate(current(ho)%phi(count(material_homogenizationAt2==ho)), source=1.0_pReal)
configHomogenization => configHomogenizations%get(ho)
associate(prm => param(ho))
if (configHomogenization%contains('damage')) then
configHomogenizationDamage => configHomogenization%get('damage')
#if defined (__GFORTRAN__)
prm%output = output_asStrings(configHomogenizationDamage)
#else
prm%output = configHomogenizationDamage%get_asStrings('output',defaultVal=emptyStringArray)
#endif
else
prm%output = emptyStringArray
endif
end associate
enddo
end subroutine damage_init
@ -23,13 +64,15 @@ end subroutine damage_init
!--------------------------------------------------------------------------------------------------
!> @brief Partition temperature onto the individual constituents.
!--------------------------------------------------------------------------------------------------
module subroutine damage_partition(phi,ce)
module subroutine damage_partition(ce)
real(pReal), intent(in) :: phi
real(pReal) :: phi
integer, intent(in) :: ce
integer :: co
phi = current(material_homogenizationAt2(ce))%phi(material_homogenizationMemberAt2(ce))
do co = 1, homogenization_Nconstituents(material_homogenizationAt2(ce))
call constitutive_damage_set_phi(phi,co,ce)
enddo
@ -37,4 +80,94 @@ module subroutine damage_partition(phi,ce)
end subroutine damage_partition
!--------------------------------------------------------------------------------------------------
!> @brief Returns homogenized nonlocal damage mobility
!--------------------------------------------------------------------------------------------------
module function damage_nonlocal_getMobility(ip,el) result(M)
integer, intent(in) :: &
ip, & !< integration point number
el !< element number
integer :: &
co
real(pReal) :: M
M = 0.0_pReal
do co = 1, homogenization_Nconstituents(material_homogenizationAt(el))
M = M + lattice_M(material_phaseAt(co,el))
enddo
M = M/real(homogenization_Nconstituents(material_homogenizationAt(el)),pReal)
end function damage_nonlocal_getMobility
!--------------------------------------------------------------------------------------------------
!> @brief calculates homogenized damage driving forces
!--------------------------------------------------------------------------------------------------
module subroutine damage_nonlocal_getSourceAndItsTangent(phiDot, dPhiDot_dPhi, phi, ip, el)
integer, intent(in) :: &
ip, & !< integration point number
el !< element number
real(pReal), intent(in) :: &
phi
real(pReal) :: &
phiDot, dPhiDot_dPhi
phiDot = 0.0_pReal
dPhiDot_dPhi = 0.0_pReal
call constitutive_damage_getRateAndItsTangents(phiDot, dPhiDot_dPhi, phi, ip, el)
phiDot = phiDot/real(homogenization_Nconstituents(material_homogenizationAt(el)),pReal)
dPhiDot_dPhi = dPhiDot_dPhi/real(homogenization_Nconstituents(material_homogenizationAt(el)),pReal)
end subroutine damage_nonlocal_getSourceAndItsTangent
!--------------------------------------------------------------------------------------------------
!> @brief updated nonlocal damage field with solution from damage phase field PDE
!--------------------------------------------------------------------------------------------------
module subroutine damage_nonlocal_putNonLocalDamage(phi,ip,el)
integer, intent(in) :: &
ip, & !< integration point number
el !< element number
real(pReal), intent(in) :: &
phi
integer :: &
homog, &
offset
homog = material_homogenizationAt(el)
offset = material_homogenizationMemberAt(ip,el)
damage(homog)%p(offset) = phi
end subroutine damage_nonlocal_putNonLocalDamage
!--------------------------------------------------------------------------------------------------
!> @brief writes results to HDF5 output file
!--------------------------------------------------------------------------------------------------
module subroutine damage_nonlocal_results(homog,group)
integer, intent(in) :: homog
character(len=*), intent(in) :: group
integer :: o
associate(prm => param(damage_typeInstance(homog)))
outputsLoop: do o = 1,size(prm%output)
select case(prm%output(o))
case ('phi')
call results_writeDataset(group,damage(homog)%p,prm%output(o),&
'damage indicator','-')
end select
enddo outputsLoop
end associate
end subroutine damage_nonlocal_results
end submodule homogenization_damage

View File

@ -171,7 +171,6 @@ module subroutine mech_RGC_init(num_homogMech)
homogState(h)%sizeState = sizeState
allocate(homogState(h)%state0 (sizeState,Nmaterialpoints), source=0.0_pReal)
allocate(homogState(h)%subState0(sizeState,Nmaterialpoints), source=0.0_pReal)
allocate(homogState(h)%state (sizeState,Nmaterialpoints), source=0.0_pReal)
stt%relaxationVector => homogState(h)%state(1:nIntFaceTot,:)

View File

@ -64,7 +64,6 @@ module subroutine mech_isostrain_init
Nmaterialpoints = count(material_homogenizationAt == h)
homogState(h)%sizeState = 0
allocate(homogState(h)%state0 (0,Nmaterialpoints))
allocate(homogState(h)%subState0(0,Nmaterialpoints))
allocate(homogState(h)%state (0,Nmaterialpoints))
end associate

View File

@ -32,7 +32,6 @@ module subroutine mech_none_init
Nmaterialpoints = count(material_homogenizationAt == h)
homogState(h)%sizeState = 0
allocate(homogState(h)%state0 (0,Nmaterialpoints))
allocate(homogState(h)%subState0(0,Nmaterialpoints))
allocate(homogState(h)%state (0,Nmaterialpoints))
enddo

View File

@ -32,13 +32,13 @@ module prec
real(pReal), dimension(:), pointer :: p
end type group_float
! http://stackoverflow.com/questions/3948210/can-i-have-a-pointer-to-an-item-in-an-allocatable-array
type :: tState
integer :: &
sizeState = 0, & !< size of state
sizeDotState = 0, & !< size of dot state, i.e. state(1:sizeDot) follows time evolution by dotState rates
offsetDeltaState = 0, & !< index offset of delta state
sizeDeltaState = 0 !< size of delta state, i.e. state(offset+1:offset+sizeDelta) follows time evolution by deltaState increments
! http://stackoverflow.com/questions/3948210
real(pReal), pointer, dimension(:), contiguous :: &
atol
real(pReal), pointer, dimension(:,:), contiguous :: & ! a pointer is needed here because we might point to state/doState. However, they will never point to something, but are rather allocated and, hence, contiguous
@ -47,7 +47,6 @@ module prec
dotState, & !< rate of state change
deltaState !< increment of state change
real(pReal), allocatable, dimension(:,:) :: &
partitionedState0, &
subState0
end type