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

Merge remote-tracking branch 'origin/separate-damage' into development

This commit is contained in:
Franz Roters 2021-01-22 14:25:32 +01:00
parent 75024433aa 32bb0d8c6e
commit 6fe8077b81
15 changed files with 295 additions and 206 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
src/commercialFEM_fileList.f90
View File

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

54
src/constitutive.f90
View File

@ -188,6 +188,11 @@ module constitutive
real(pReal), dimension(3,3) :: P
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)
integer, intent(in) :: ph,me
real(pReal) :: T
@ -203,6 +208,11 @@ module constitutive
real(pReal), intent(in) :: T
integer, intent(in) :: co, ce
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 ===================================================================================
@ -229,39 +239,14 @@ module constitutive
logical :: converged_
end function crystallite_stress
module function constitutive_homogenizedC(co,ip,el) result(C)
integer, intent(in) :: co, ip, el
module function constitutive_homogenizedC(ph,me) result(C)
integer, intent(in) :: ph, me
real(pReal), dimension(6,6) :: C
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)
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)
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)
integer, intent(in) :: &
co, & !< component-ID of integration point
@ -391,6 +363,8 @@ module constitutive
constitutive_restartRead, &
integrateDamageState, &
constitutive_thermal_setT, &
constitutive_damage_set_phi, &
constitutive_damage_get_phi, &
constitutive_mech_getP, &
constitutive_mech_setF, &
constitutive_mech_getF, &

103
src/constitutive_damage.f90
View File

@ -10,9 +10,16 @@ submodule(constitutive) constitutive_damage
DAMAGE_ANISODUCTILE_ID
end enum
type :: tDataContainer
real(pReal), dimension(:), allocatable :: phi, d_phi_d_dot_phi
end type tDataContainer
integer(kind(DAMAGE_UNDEFINED_ID)), dimension(:,:), allocatable :: &
phase_source !< active sources mechanisms of each phase
type(tDataContainer), dimension(:), allocatable :: current
interface
module function source_damage_anisoBrittle_init(source_length) result(mySources)
@ -45,6 +52,38 @@ submodule(constitutive) constitutive_damage
logical, dimension(:,:), allocatable :: myKinematics
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)
integer, intent(in) :: &
@ -120,7 +159,8 @@ contains
module subroutine damage_init
integer :: &
ph !< counter in phase loop
ph, & !< counter in phase loop
Nconstituents
class(tNode), pointer :: &
phases, &
phase, &
@ -129,10 +169,18 @@ module subroutine damage_init
phases => config_material%get('phase')
allocate(current(phases%length))
allocate(damageState (phases%length))
allocate(phase_Nsources(phases%length),source = 0) ! same for kinematics
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)
sources => phase%get('source',defaultVal=emptyList)
phase_Nsources(ph) = sources%length
@ -295,7 +343,7 @@ module function integrateDamageState(dt,co,ip,el) result(broken)
enddo
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
endif
@ -368,14 +416,14 @@ end subroutine damage_results
!--------------------------------------------------------------------------------------------------
!> @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) :: &
co, & !< component-ID of integration point
co, & !< component-ID me integration point
ip, & !< integration point
el, & !< element
ph, &
of
me
integer :: &
so !< counter in source loop
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)
case (DAMAGE_ANISOBRITTLE_ID) sourceType
call source_damage_anisoBrittle_dotState(mech_S(material_phaseAt(co,el),material_phaseMemberAt(co,ip,el)),&
co, ip, el) ! correct stress?
call source_damage_anisoBrittle_dotState(mech_S(ph,me),co, ip, el) ! correct stress?
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
@ -411,14 +458,11 @@ end function constitutive_damage_collectDotState
!> @brief for constitutive models having an instantaneous change of state
!> 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) :: &
co, & !< component-ID of integration point
ip, & !< integration point
el, & !< element
ph, &
of
me
real(pReal), intent(in), dimension(3,3) :: &
Fe !< elastic deformation gradient
integer :: &
@ -436,14 +480,13 @@ function constitutive_damage_deltaState(Fe, co, ip, el, ph, of) result(broken)
sourceType: select case (phase_source(so,ph))
case (DAMAGE_ISOBRITTLE_ID) sourceType
call source_damage_isoBrittle_deltaState (constitutive_homogenizedC(co,ip,el), Fe, &
co, ip, el)
broken = any(IEEE_is_NaN(damageState(ph)%p(so)%deltaState(:,of)))
call source_damage_isoBrittle_deltaState(constitutive_homogenizedC(ph,me), Fe, ph,me)
broken = any(IEEE_is_NaN(damageState(ph)%p(so)%deltaState(:,me)))
if(.not. broken) then
myOffset = damageState(ph)%p(so)%offsetDeltaState
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,of) + damageState(ph)%p(so)%deltaState(1:mySize,of)
damageState(ph)%p(so)%state(myOffset + 1: myOffset + mySize,me) = &
damageState(ph)%p(so)%state(myOffset + 1: myOffset + mySize,me) + damageState(ph)%p(so)%deltaState(1:mySize,me)
endif
end select sourceType
@ -484,4 +527,28 @@ function source_active(source_label,src_length) result(active_source)
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

31
src/constitutive_mech.f90
View File

@ -307,13 +307,9 @@ submodule(constitutive) constitutive_mech
character(len=*), intent(in) :: group
end subroutine plastic_nonlocal_results
module function plastic_dislotwin_homogenizedC(co,ip,el) result(homogenizedC)
real(pReal), dimension(6,6) :: &
homogenizedC
integer, intent(in) :: &
co, & !< component-ID of integration point
ip, & !< integration point
el !< element
module function plastic_dislotwin_homogenizedC(ph,me) result(homogenizedC)
real(pReal), dimension(6,6) :: homogenizedC
integer, intent(in) :: ph,me
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) :: &
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
real(pReal), dimension(3,3) :: E
real(pReal), dimension(3,3,3,3) :: C
integer :: &
ho, & !< homogenization
d !< counter in degradation loop
integer :: &
i, j
d, & !< counter in degradation loop
i, j, ph, me
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))
degradationType: select case(phase_stiffnessDegradation(d,material_phaseAt(co,el)))
@ -1535,19 +1531,16 @@ end subroutine mech_forward
!> @brief returns the homogenize elasticity matrix
!> 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
integer, intent(in) :: &
co, & !< component-ID of integration point
ip, & !< integration point
el !< element
integer, intent(in) :: ph, me
plasticityType: select case (phase_plasticity(material_phaseAt(co,el)))
plasticityType: select case (phase_plasticity(ph))
case (PLASTICITY_DISLOTWIN_ID) plasticityType
C = plastic_dislotwin_homogenizedC(co,ip,el)
C = plastic_dislotwin_homogenizedC(ph,me)
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 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)
print'(a,i2)', ' # instances: ',Ninstances; flush(IO_STDOUT)
if(Ninstances == 0) return
print*, 'Ma and Roters, Acta Materialia 52(12):3603–3612, 2004'
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.
!--------------------------------------------------------------------------------------------------
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) :: &
homogenizedC
integer, intent(in) :: &
co, & !< component-ID of integration point
ip, & !< integration point
el !< element
integer :: i, &
of
integer :: i
real(pReal) :: f_unrotated
of = material_phasememberAt(co,ip,el)
associate(prm => param(phase_plasticityInstance(material_phaseAt(co,el))),&
stt => state(phase_plasticityInstance(material_phaseAT(co,el))))
associate(prm => param(phase_plasticityInstance(ph)),&
stt => state(phase_plasticityInstance(ph)))
f_unrotated = 1.0_pReal &
- sum(stt%f_tw(1:prm%sum_N_tw,of)) &
- sum(stt%f_tr(1:prm%sum_N_tr,of))
- sum(stt%f_tw(1:prm%sum_N_tw,me)) &
- sum(stt%f_tr(1:prm%sum_N_tr,me))
homogenizedC = f_unrotated * prm%C66
do i=1,prm%sum_N_tw
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
do i=1,prm%sum_N_tr
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
end associate

13
src/constitutive_thermal.f90
View File

@ -36,6 +36,13 @@ submodule(constitutive) constitutive_thermal
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)
integer, intent(in) :: &
phase !< phase ID of element
@ -47,10 +54,10 @@ submodule(constitutive) constitutive_thermal
TDot
end subroutine thermal_dissipation_getRate
module subroutine thermal_externalheat_getRate(TDot, phase,of)
module subroutine thermal_externalheat_getRate(TDot, ph,me)
integer, intent(in) :: &
phase, &
of
ph, &
me
real(pReal), intent(out) :: &
TDot
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
!> @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) :: &
phase, &
of
ph, &
me
integer :: &
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
@ -104,11 +104,11 @@ end subroutine source_thermal_externalheat_dotState
!--------------------------------------------------------------------------------------------------
!> @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) :: &
phase, &
of
ph, &
me
real(pReal), intent(out) :: &
TDot
@ -117,18 +117,18 @@ module subroutine thermal_externalheat_getRate(TDot, phase, of)
real(pReal) :: &
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
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
if ( (frac_time < 0.0_pReal .and. interval == 1) &
.or. (frac_time >= 1.0_pReal .and. interval == prm%nIntervals) &
.or. (frac_time >= 0.0_pReal .and. frac_time < 1.0_pReal) ) &
TDot = prm%f_T(interval ) * (1.0_pReal - frac_time) + &
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
end associate

69
src/grid/grid_damage_spectral.f90
View File

@ -17,8 +17,9 @@ module grid_damage_spectral
use discretization_grid
use damage_nonlocal
use YAML_types
use homogenization
use config
implicit none
private
@ -43,13 +44,13 @@ module grid_damage_spectral
phi_current, & !< field of current damage
phi_lastInc, & !< field of previous 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
real(pReal), dimension(3,3) :: K_ref
real(pReal) :: mu_ref
public :: &
grid_damage_spectral_init, &
grid_damage_spectral_solution, &
@ -62,8 +63,8 @@ contains
! ToDo: Restart not implemented
!--------------------------------------------------------------------------------------------------
subroutine grid_damage_spectral_init
PetscInt, dimension(0:worldsize-1) :: localK
PetscInt, dimension(0:worldsize-1) :: localK
DM :: damage_grid
Vec :: uBound, lBound
PetscErrorCode :: ierr
@ -72,22 +73,22 @@ subroutine grid_damage_spectral_init
num_generic
character(len=pStringLen) :: &
snes_type
print'(/,a)', ' <<<+- grid_spectral_damage init -+>>>'
print*, 'Shanthraj et al., Handbook of Mechanics of Materials, 2019'
print*, 'https://doi.org/10.1007/978-981-10-6855-3_80'
!-------------------------------------------------------------------------------------------------
! read numerical parameters and do sanity checks
num_grid => config_numerics%get('grid',defaultVal=emptyDict)
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_rtol = num_grid%get_asFloat ('eps_damage_rtol',defaultVal=1.0e-6_pReal)
num_generic => config_numerics%get('generic',defaultVal=emptyDict)
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%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')
@ -104,7 +105,7 @@ subroutine grid_damage_spectral_init
!--------------------------------------------------------------------------------------------------
! initialize solver specific parts of PETSc
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(worldrank) = grid3
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 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
CHKERRQ(ierr)
CHKERRQ(ierr)
call SNESSetFromOptions(damage_snes,ierr); CHKERRQ(ierr) ! pull it all together with additional CLI arguments
call SNESGetType(damage_snes,snes_type,ierr); CHKERRQ(ierr)
if (trim(snes_type) == 'vinewtonrsls' .or. &
@ -137,12 +138,12 @@ subroutine grid_damage_spectral_init
endif
!--------------------------------------------------------------------------------------------------
! init fields
! init fields
call DMDAGetCorners(damage_grid,xstart,ystart,zstart,xend,yend,zend,ierr)
CHKERRQ(ierr)
xend = xstart + xend - 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_lastInc(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
!--------------------------------------------------------------------------------------------------
function grid_damage_spectral_solution(timeinc) result(solution)
real(pReal), intent(in) :: &
timeinc !< increment in time for current solution
integer :: i, j, k, cell
type(tSolutionState) :: solution
PetscInt :: devNull
PetscReal :: phi_min, phi_max, stagNorm, solnNorm
PetscErrorCode :: ierr
SNESConvergedReason :: reason
solution%converged =.false.
!--------------------------------------------------------------------------------------------------
! set module wide availabe data
! set module wide availabe data
params%timeinc = timeinc
call SNESSolve(damage_snes,PETSC_NULL_VEC,solution_vec,ierr); CHKERRQ(ierr)
call SNESGetConvergedReason(damage_snes,reason,ierr); CHKERRQ(ierr)
if (reason < 1) then
solution%converged = .false.
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)
!--------------------------------------------------------------------------------------------------
! updating damage state
! updating damage state
cell = 0
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)
call VecMax(solution_vec,devNull,phi_max,ierr); CHKERRQ(ierr)
if (solution%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)', ' ==========================================================================='
flush(IO_STDOUT)
flush(IO_STDOUT)
end function grid_damage_spectral_solution
@ -214,31 +216,32 @@ end function grid_damage_spectral_solution
!> @brief spectral damage forwarding routine
!--------------------------------------------------------------------------------------------------
subroutine grid_damage_spectral_forward(cutBack)
logical, intent(in) :: cutBack
integer :: i, j, k, cell
DM :: dm_local
PetscScalar, dimension(:,:,:), pointer :: x_scal
PetscErrorCode :: ierr
if (cutBack) then
if (cutBack) then
phi_current = phi_lastInc
phi_stagInc = phi_lastInc
!--------------------------------------------------------------------------------------------------
! reverting damage field state
! reverting damage field state
cell = 0
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
x_scal(xstart:xend,ystart:yend,zstart:zend) = phi_current
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)
cell = cell + 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
call updateReference
endif
endif
end subroutine grid_damage_spectral_forward
@ -247,7 +250,7 @@ end subroutine grid_damage_spectral_forward
!> @brief forms the spectral damage residual vector
!--------------------------------------------------------------------------------------------------
subroutine formResidual(in,x_scal,f_scal,dummy,ierr)
DMDALocalInfo, dimension(DMDA_LOCAL_INFO_SIZE) :: &
in
PetscScalar, dimension( &
@ -261,11 +264,11 @@ subroutine formResidual(in,x_scal,f_scal,dummy,ierr)
integer :: i, j, k, cell
real(pReal) :: phiDot, dPhiDot_dPhi, mobility
phi_current = x_scal
phi_current = x_scal
!--------------------------------------------------------------------------------------------------
! evaluate polarization field
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_fourierScalarGradient !< calculate gradient of damage field
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
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
!--------------------------------------------------------------------------------------------------
! constructing residual
f_scal = scalarField_real(1:grid(1),1:grid(2),1:grid3) - phi_current
@ -311,7 +314,7 @@ end subroutine formResidual
subroutine updateReference
integer :: i,j,k,cell,ierr
cell = 0
K_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
real(pReal), dimension(:), allocatable, public :: &
homogenization_T, &
homogenization_dot_T
homogenization_dot_T, &
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
@ -62,6 +64,9 @@ module homogenization
module subroutine thermal_init
end subroutine thermal_init
module subroutine damage_init
end subroutine damage_init
module subroutine mech_partition(subF,ip,el)
real(pReal), intent(in), dimension(3,3) :: &
subF
@ -75,6 +80,11 @@ module homogenization
integer, intent(in) :: ce
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)
integer, intent(in) :: ip,el
end subroutine thermal_homogenize
@ -118,7 +128,7 @@ contains
!--------------------------------------------------------------------------------------------------
!> @brief module initialization
!--------------------------------------------------------------------------------------------------
subroutine homogenization_init
subroutine homogenization_init()
class (tNode) , pointer :: &
num_homog, &
@ -142,6 +152,7 @@ subroutine homogenization_init
call mech_init(num_homog)
call thermal_init()
call damage_init()
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)

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
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

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
q !< damage rate sensitivity
real(pReal), dimension(:), allocatable :: &
s_crit, & !< critical displacement
s_crit, & !< critical displacement
g_crit !< critical load
real(pReal), dimension(:,:,:,:), allocatable :: &
cleavage_systems
@ -37,7 +37,7 @@ contains
!--------------------------------------------------------------------------------------------------
module function source_damage_anisoBrittle_init(source_length) result(mySources)
integer, intent(in) :: source_length
integer, intent(in) :: source_length
logical, dimension(:,:), allocatable :: mySources
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)
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(count(mySources(:,p)) == 0) cycle
sources => phase%get('source')
@ -70,20 +70,20 @@ module function source_damage_anisoBrittle_init(source_length) result(mySources)
if(mySources(sourceOffset,p)) then
source_damage_anisoBrittle_offset(p) = sourceOffset
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)
prm%sum_N_cl = sum(abs(N_cl))
prm%q = src%get_asFloat('q')
prm%dot_o = src%get_asFloat('dot_o')
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%cleavage_systems = lattice_SchmidMatrix_cleavage(N_cl,phase%get_asString('lattice'),&
phase%get_asFloat('c/a',defaultVal=0.0_pReal))
! expand: family => system
prm%s_crit = math_expand(prm%s_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
prm%output = src%get_asStrings('output',defaultVal=emptyStringArray)
#endif
! sanity checks
if (prm%q <= 0.0_pReal) extmsg = trim(extmsg)//' q'
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
integer :: &
phase, &
constituent, &
ph, &
me, &
sourceOffset, &
damageOffset, &
homog, &
@ -139,14 +139,14 @@ module subroutine source_damage_anisoBrittle_dotState(S, co, ip, el)
real(pReal) :: &
traction_d, traction_t, traction_n, traction_crit
phase = material_phaseAt(co,el)
constituent = material_phasememberAt(co,ip,el)
sourceOffset = source_damage_anisoBrittle_offset(phase)
ph = material_phaseAt(co,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(phase)))
damageState(phase)%p(sourceOffset)%dotState(1,constituent) = 0.0_pReal
associate(prm => param(source_damage_anisoBrittle_instance(ph)))
damageState(ph)%p(sourceOffset)%dotState(1,me) = 0.0_pReal
do i = 1, prm%sum_N_cl
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))
@ -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
damageState(phase)%p(sourceOffset)%dotState(1,constituent) &
= damageState(phase)%p(sourceOffset)%dotState(1,constituent) &
damageState(ph)%p(sourceOffset)%dotState(1,me) &
= damageState(ph)%p(sourceOffset)%dotState(1,me) &
+ 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_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)
integer, intent(in) :: source_length
integer, intent(in) :: source_length
logical, dimension(:,:), allocatable :: mySources
class(tNode), pointer :: &
@ -67,7 +67,7 @@ module function source_damage_anisoDuctile_init(source_length) result(mySources)
if(mySources(sourceOffset,p)) then
source_damage_anisoDuctile_offset(p) = sourceOffset
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)
prm%q = src%get_asFloat('q')
@ -81,7 +81,7 @@ module function source_damage_anisoDuctile_init(source_length) result(mySources)
#else
prm%output = src%get_asStrings('output',defaultVal=emptyStringArray)
#endif
! sanity checks
if (prm%q <= 0.0_pReal) extmsg = trim(extmsg)//' q'
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
integer :: &
phase, &
constituent, &
ph, &
me, &
sourceOffset, &
damageOffset, &
homog
phase = material_phaseAt(co,el)
constituent = material_phasememberAt(co,ip,el)
sourceOffset = source_damage_anisoDuctile_offset(phase)
ph = material_phaseAt(co,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(phase)))
damageState(phase)%p(sourceOffset)%dotState(1,constituent) &
= sum(plasticState(phase)%slipRate(:,constituent)/(damage(homog)%p(damageOffset)**prm%q)/prm%gamma_crit)
associate(prm => param(source_damage_anisoDuctile_instance(ph)))
damageState(ph)%p(sourceOffset)%dotState(1,me) &
= sum(plasticState(ph)%slipRate(:,me)/(damage(homog)%p(damageOffset)**prm%q)/prm%gamma_crit)
end associate
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)
integer, intent(in) :: source_length
integer, intent(in) :: source_length
logical, dimension(:,:), allocatable :: mySources
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)
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(count(mySources(:,p)) == 0) cycle
sources => phase%get('source')
@ -60,7 +60,7 @@ module function source_damage_isoBrittle_init(source_length) result(mySources)
if(mySources(sourceOffset,p)) then
source_damage_isoBrittle_offset(p) = sourceOffset
associate(prm => param(source_damage_isoBrittle_instance(p)))
src => sources%get(sourceOffset)
src => sources%get(sourceOffset)
prm%W_crit = src%get_asFloat('W_crit')
@ -69,7 +69,7 @@ module function source_damage_isoBrittle_init(source_length) result(mySources)
#else
prm%output = src%get_asStrings('output',defaultVal=emptyStringArray)
#endif
! sanity checks
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
!--------------------------------------------------------------------------------------------------
module subroutine source_damage_isoBrittle_deltaState(C, Fe, co, ip, el)
module subroutine source_damage_isoBrittle_deltaState(C, Fe, ph,me)
integer, intent(in) :: &
co, & !< component-ID of integration point
ip, & !< integration point
el !< element
integer, intent(in) :: ph,me
real(pReal), intent(in), dimension(3,3) :: &
Fe
real(pReal), intent(in), dimension(6,6) :: &
C
integer :: &
phase, &
constituent, &
sourceOffset
real(pReal), dimension(6) :: &
strain
real(pReal) :: &
strainenergy
phase = material_phaseAt(co,el) !< phase ID at co,ip,el
constituent = material_phasememberAt(co,ip,el) !< state array offset for phase ID at co,ip,el
sourceOffset = source_damage_isoBrittle_offset(phase)
sourceOffset = source_damage_isoBrittle_offset(ph)
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
! 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
damageState(phase)%p(sourceOffset)%deltaState(1,constituent) = &
strainenergy - damageState(phase)%p(sourceOffset)%state(1,constituent)
if (strainenergy > damageState(ph)%p(sourceOffset)%subState0(1,me)) then
damageState(ph)%p(sourceOffset)%deltaState(1,me) = &
strainenergy - damageState(ph)%p(sourceOffset)%state(1,me)
else
damageState(phase)%p(sourceOffset)%deltaState(1,constituent) = &
damageState(phase)%p(sourceOffset)%subState0(1,constituent) - &
damageState(phase)%p(sourceOffset)%state(1,constituent)
damageState(ph)%p(sourceOffset)%deltaState(1,me) = &
damageState(ph)%p(sourceOffset)%subState0(1,me) - &
damageState(ph)%p(sourceOffset)%state(1,me)
endif
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)
integer, intent(in) :: source_length
integer, intent(in) :: source_length
logical, dimension(:,:), allocatable :: mySources
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)
do p = 1, phases%length
phase => phases%get(p)
phase => phases%get(p)
if(count(mySources(:,p)) == 0) cycle
if(any(mySources(:,p))) source_damage_isoDuctile_instance(p) = count(mySources(:,1:p))
sources => phase%get('source')
@ -62,7 +62,7 @@ module function source_damage_isoDuctile_init(source_length) result(mySources)
if(mySources(sourceOffset,p)) then
source_damage_isoDuctile_offset(p) = sourceOffset
associate(prm => param(source_damage_isoDuctile_instance(p)))
src => sources%get(sourceOffset)
src => sources%get(sourceOffset)
prm%q = src%get_asFloat('q')
prm%gamma_crit = src%get_asFloat('gamma_crit')
@ -72,7 +72,7 @@ module function source_damage_isoDuctile_init(source_length) result(mySources)
#else
prm%output = src%get_asStrings('output',defaultVal=emptyStringArray)
#endif
! sanity checks
if (prm%q <= 0.0_pReal) extmsg = trim(extmsg)//' q'
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
integer :: &
phase, &
constituent, &
ph, &
me, &
sourceOffset, &
damageOffset, &
homog
phase = material_phaseAt(co,el)
constituent = material_phasememberAt(co,ip,el)
sourceOffset = source_damage_isoDuctile_offset(phase)
ph = material_phaseAt(co,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(phase)))
damageState(phase)%p(sourceOffset)%dotState(1,constituent) = &
sum(plasticState(phase)%slipRate(:,constituent))/(damage(homog)%p(damageOffset)**prm%q)/prm%gamma_crit
associate(prm => param(source_damage_isoDuctile_instance(ph)))
damageState(ph)%p(sourceOffset)%dotState(1,me) = &
sum(plasticState(ph)%slipRate(:,me))/(damage(homog)%p(damageOffset)**prm%q)/prm%gamma_crit
end associate
end subroutine source_damage_isoDuctile_dotState