From 4515920b695e01aa6b0a343eb3ffa41887e7fc4d Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 11 Feb 2021 23:54:43 +0100 Subject: [PATCH 01/25] not needed --- src/homogenization.f90 | 20 -------------------- 1 file changed, 20 deletions(-) diff --git a/src/homogenization.f90 b/src/homogenization.f90 index 62ac52f06..14dd26726 100644 --- a/src/homogenization.f90 +++ b/src/homogenization.f90 @@ -311,26 +311,6 @@ 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) From f6be3fe0b74be7b4ee2e3d728fee9b75d0327f42 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 12 Feb 2021 00:12:57 +0100 Subject: [PATCH 02/25] no need for pointer --- src/damage_none.f90 | 4 ---- src/damage_nonlocal.f90 | 2 -- src/homogenization_damage.f90 | 4 ++-- src/material.f90 | 6 ------ src/phase_damage_anisobrittle.f90 | 2 +- src/phase_damage_anisoductile.f90 | 2 +- src/phase_damage_isoductile.f90 | 2 +- src/phase_mechanical.f90 | 2 +- src/phase_mechanical_eigen_cleavageopening.f90 | 2 +- src/phase_mechanical_eigen_slipplaneopening.f90 | 2 +- 10 files changed, 8 insertions(+), 20 deletions(-) diff --git a/src/damage_none.f90 b/src/damage_none.f90 index 680110d55..c4b1c1589 100644 --- a/src/damage_none.f90 +++ b/src/damage_none.f90 @@ -3,8 +3,6 @@ !> @brief material subroutine for constant damage field !-------------------------------------------------------------------------------------------------- module damage_none - use prec - use config use material implicit none @@ -29,8 +27,6 @@ subroutine damage_none_init allocate(damageState_h(h)%state0 (0,Nmaterialpoints)) allocate(damageState_h(h)%state (0,Nmaterialpoints)) - allocate (damage(h)%p(Nmaterialpoints), source=1.0_pReal) - enddo end subroutine damage_none_init diff --git a/src/damage_nonlocal.f90 b/src/damage_nonlocal.f90 index 807231889..dd50cf4cf 100644 --- a/src/damage_nonlocal.f90 +++ b/src/damage_nonlocal.f90 @@ -57,8 +57,6 @@ subroutine damage_nonlocal_init allocate(damageState_h(h)%state0 (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,:) - enddo end subroutine damage_nonlocal_init diff --git a/src/homogenization_damage.f90 b/src/homogenization_damage.f90 index 8b73ed54a..4f23af3be 100644 --- a/src/homogenization_damage.f90 +++ b/src/homogenization_damage.f90 @@ -143,7 +143,7 @@ module subroutine damage_nonlocal_putNonLocalDamage(phi,ip,el) homog = material_homogenizationAt(el) offset = material_homogenizationMemberAt(ip,el) - damage(homog)%p(offset) = phi + damagestate_h(homog)%state(1,offset) = phi end subroutine damage_nonlocal_putNonLocalDamage @@ -162,7 +162,7 @@ module subroutine damage_nonlocal_results(homog,group) outputsLoop: do o = 1,size(prm%output) select case(prm%output(o)) case ('phi') - call results_writeDataset(group,damage(homog)%p,prm%output(o),& + call results_writeDataset(group,damagestate_h(homog)%state(1,:),prm%output(o),& 'damage indicator','-') end select enddo outputsLoop diff --git a/src/material.f90 b/src/material.f90 index ccc6898be..3689cd0b0 100644 --- a/src/material.f90 +++ b/src/material.f90 @@ -70,9 +70,6 @@ module material type(Rotation), dimension(:,:,:), allocatable, public, protected :: & material_orientation0 !< initial orientation of each grain,IP,element - type(group_float), allocatable, dimension(:), public :: & - damage !< damage field - public :: & material_init, & THERMAL_ISOTHERMAL_ID, & @@ -105,9 +102,6 @@ subroutine material_init(restart) allocate(homogState (size(material_name_homogenization))) allocate(damageState_h (size(material_name_homogenization))) - allocate(damage (size(material_name_homogenization))) - - if (.not. restart) then call results_openJobFile call results_mapping_phase(material_phaseAt,material_phaseMemberAt,material_name_phase) diff --git a/src/phase_damage_anisobrittle.f90 b/src/phase_damage_anisobrittle.f90 index 51ca7786a..b328dfed7 100644 --- a/src/phase_damage_anisobrittle.f90 +++ b/src/phase_damage_anisobrittle.f90 @@ -152,7 +152,7 @@ module subroutine anisobrittle_dotState(S, co, ip, el) traction_t = math_tensordot(S,prm%cleavage_systems(1:3,1:3,2,i)) traction_n = math_tensordot(S,prm%cleavage_systems(1:3,1:3,3,i)) - traction_crit = prm%g_crit(i)*damage(homog)%p(damageOffset)**2.0_pReal + traction_crit = prm%g_crit(i)*damagestate_h(homog)%state(1,damageOffset)**2.0_pReal damageState(ph)%p(sourceOffset)%dotState(1,me) & = damageState(ph)%p(sourceOffset)%dotState(1,me) & diff --git a/src/phase_damage_anisoductile.f90 b/src/phase_damage_anisoductile.f90 index e54eff201..18c29fdfc 100644 --- a/src/phase_damage_anisoductile.f90 +++ b/src/phase_damage_anisoductile.f90 @@ -129,7 +129,7 @@ module subroutine anisoductile_dotState(co, ip, el) 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) + = sum(plasticState(ph)%slipRate(:,me)/(damagestate_h(homog)%state(1,damageOffset)**prm%q)/prm%gamma_crit) end associate end subroutine anisoductile_dotState diff --git a/src/phase_damage_isoductile.f90 b/src/phase_damage_isoductile.f90 index 71f99078b..9dac17cd8 100644 --- a/src/phase_damage_isoductile.f90 +++ b/src/phase_damage_isoductile.f90 @@ -120,7 +120,7 @@ module subroutine isoductile_dotState(co, ip, el) 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 + sum(plasticState(ph)%slipRate(:,me))/(damagestate_h(homog)%state(1,damageOffset)**prm%q)/prm%gamma_crit end associate end subroutine isoductile_dotState diff --git a/src/phase_mechanical.f90 b/src/phase_mechanical.f90 index efd4d51d2..f2f3333d3 100644 --- a/src/phase_mechanical.f90 +++ b/src/phase_mechanical.f90 @@ -381,7 +381,7 @@ subroutine phase_hooke_SandItsTangents(S, dS_dFe, dS_dFi, & DegradationLoop: do d = 1, phase_NstiffnessDegradations(material_phaseAt(co,el)) degradationType: select case(phase_stiffnessDegradation(d,material_phaseAt(co,el))) case (STIFFNESS_DEGRADATION_damage_ID) degradationType - C = C * damage(ho)%p(material_homogenizationMemberAt(ip,el))**2 + C = C * damagestate_h(ho)%state(1,material_homogenizationMemberAt(ip,el))**2 end select degradationType enddo DegradationLoop diff --git a/src/phase_mechanical_eigen_cleavageopening.f90 b/src/phase_mechanical_eigen_cleavageopening.f90 index 59be7837a..b41cb1e9c 100644 --- a/src/phase_mechanical_eigen_cleavageopening.f90 +++ b/src/phase_mechanical_eigen_cleavageopening.f90 @@ -126,7 +126,7 @@ module subroutine kinematics_cleavage_opening_LiAndItsTangent(Ld, dLd_dTstar, S, dLd_dTstar = 0.0_pReal associate(prm => param(kinematics_cleavage_opening_instance(material_phaseAt(co,el)))) do i = 1,prm%sum_N_cl - traction_crit = prm%g_crit(i)* damage(homog)%p(damageOffset)**2.0_pReal + traction_crit = prm%g_crit(i)*damagestate_h(homog)%state(1,damageOffset)**2.0_pReal traction_d = math_tensordot(S,prm%cleavage_systems(1:3,1:3,1,i)) if (abs(traction_d) > traction_crit + tol_math_check) then diff --git a/src/phase_mechanical_eigen_slipplaneopening.f90 b/src/phase_mechanical_eigen_slipplaneopening.f90 index a144d39a6..682f0b825 100644 --- a/src/phase_mechanical_eigen_slipplaneopening.f90 +++ b/src/phase_mechanical_eigen_slipplaneopening.f90 @@ -152,7 +152,7 @@ module subroutine kinematics_slipplane_opening_LiAndItsTangent(Ld, dLd_dTstar, S traction_t = math_tensordot(S,prm%P_t(1:3,1:3,i)) traction_n = math_tensordot(S,prm%P_n(1:3,1:3,i)) - traction_crit = prm%g_crit(i)* damage(homog)%p(damageOffset) ! degrading critical load carrying capacity by damage + traction_crit = prm%g_crit(i)* damagestate_h(homog)%state(1,damageOffset) ! degrading critical load carrying capacity by damage udotd = sign(1.0_pReal,traction_d)* prm%dot_o* ( abs(traction_d)/traction_crit & - abs(traction_d)/prm%g_crit(i))**prm%q From 462ca1a30b2893fb0897b26eb1eeb5ac95fbf07d Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 12 Feb 2021 08:10:25 +0100 Subject: [PATCH 03/25] not needed --- src/commercialFEM_fileList.f90 | 1 - src/damage_none.f90 | 34 ---------------------------------- src/homogenization.f90 | 2 -- 3 files changed, 37 deletions(-) delete mode 100644 src/damage_none.f90 diff --git a/src/commercialFEM_fileList.f90 b/src/commercialFEM_fileList.f90 index a191298a3..cf48e19ce 100644 --- a/src/commercialFEM_fileList.f90 +++ b/src/commercialFEM_fileList.f90 @@ -41,7 +41,6 @@ #include "phase_damage_isoductile.f90" #include "phase_damage_anisobrittle.f90" #include "phase_damage_anisoductile.f90" -#include "damage_none.f90" #include "damage_nonlocal.f90" #include "homogenization.f90" #include "homogenization_mechanical.f90" diff --git a/src/damage_none.f90 b/src/damage_none.f90 deleted file mode 100644 index c4b1c1589..000000000 --- a/src/damage_none.f90 +++ /dev/null @@ -1,34 +0,0 @@ -!-------------------------------------------------------------------------------------------------- -!> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH -!> @brief material subroutine for constant damage field -!-------------------------------------------------------------------------------------------------- -module damage_none - use material - - implicit none - public - -contains - -!-------------------------------------------------------------------------------------------------- -!> @brief allocates all neccessary fields, reads information from material configuration file -!-------------------------------------------------------------------------------------------------- -subroutine damage_none_init - - integer :: h,Nmaterialpoints - - print'(/,a)', ' <<<+- damage_none init -+>>>'; flush(6) - - do h = 1, size(material_name_homogenization) - if (damage_type(h) /= DAMAGE_NONE_ID) cycle - - Nmaterialpoints = count(material_homogenizationAt == h) - damageState_h(h)%sizeState = 0 - allocate(damageState_h(h)%state0 (0,Nmaterialpoints)) - allocate(damageState_h(h)%state (0,Nmaterialpoints)) - - enddo - -end subroutine damage_none_init - -end module damage_none diff --git a/src/homogenization.f90 b/src/homogenization.f90 index 14dd26726..87b15c43e 100644 --- a/src/homogenization.f90 +++ b/src/homogenization.f90 @@ -12,7 +12,6 @@ module homogenization use material use phase use discretization - use damage_none use damage_nonlocal use HDF5_utilities use results @@ -220,7 +219,6 @@ subroutine homogenization_init() call thermal_init() call damage_init() - if (any(damage_type == DAMAGE_none_ID)) call damage_none_init if (any(damage_type == DAMAGE_nonlocal_ID)) call damage_nonlocal_init From 7bec3e03630d8daa0c957b7cb5714fecf1274d5f Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 12 Feb 2021 09:31:02 +0100 Subject: [PATCH 04/25] use partitioned damage --- src/homogenization.f90 | 6 ++++-- src/homogenization_damage.f90 | 3 ++- src/material.f90 | 2 +- src/phase_damage_anisobrittle.f90 | 2 +- src/phase_damage_anisoductile.f90 | 2 +- src/phase_damage_isoductile.f90 | 2 +- src/phase_mechanical.f90 | 2 +- src/phase_mechanical_eigen_cleavageopening.f90 | 2 +- src/phase_mechanical_eigen_slipplaneopening.f90 | 2 +- 9 files changed, 13 insertions(+), 10 deletions(-) diff --git a/src/homogenization.f90 b/src/homogenization.f90 index 87b15c43e..3d0a4bb7a 100644 --- a/src/homogenization.f90 +++ b/src/homogenization.f90 @@ -219,7 +219,7 @@ subroutine homogenization_init() call thermal_init() call damage_init() - if (any(damage_type == DAMAGE_nonlocal_ID)) call damage_nonlocal_init + call damage_nonlocal_init end subroutine homogenization_init @@ -255,6 +255,7 @@ subroutine materialpoint_stressAndItsTangent(dt,FEsolving_execIP,FEsolving_execE 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) + call damage_partition(ce) doneAndHappy = [.false.,.true.] @@ -375,7 +376,8 @@ subroutine homogenization_forward do ho = 1, size(material_name_homogenization) homogState (ho)%state0 = homogState (ho)%state - damageState_h(ho)%state0 = damageState_h(ho)%state + if(damageState_h(ho)%sizeState > 0) & + damageState_h(ho)%state0 = damageState_h(ho)%state enddo end subroutine homogenization_forward diff --git a/src/homogenization_damage.f90 b/src/homogenization_damage.f90 index 4f23af3be..502174309 100644 --- a/src/homogenization_damage.f90 +++ b/src/homogenization_damage.f90 @@ -72,7 +72,8 @@ module subroutine damage_partition(ce) integer :: co - phi = current(material_homogenizationAt2(ce))%phi(material_homogenizationMemberAt2(ce)) + if(damageState_h(material_homogenizationAt2(ce))%sizeState < 1) return + phi = damagestate_h(material_homogenizationAt2(ce))%state(1,material_homogenizationMemberAt2(ce)) do co = 1, homogenization_Nconstituents(material_homogenizationAt2(ce)) call phase_damage_set_phi(phi,co,ce) enddo diff --git a/src/material.f90 b/src/material.f90 index 3689cd0b0..b76b2b296 100644 --- a/src/material.f90 +++ b/src/material.f90 @@ -100,7 +100,7 @@ subroutine material_init(restart) allocate(homogState (size(material_name_homogenization))) - allocate(damageState_h (size(material_name_homogenization))) + allocate(damageState_h (size(material_name_homogenization))) if (.not. restart) then call results_openJobFile diff --git a/src/phase_damage_anisobrittle.f90 b/src/phase_damage_anisobrittle.f90 index b328dfed7..699e9508b 100644 --- a/src/phase_damage_anisobrittle.f90 +++ b/src/phase_damage_anisobrittle.f90 @@ -152,7 +152,7 @@ module subroutine anisobrittle_dotState(S, co, ip, el) traction_t = math_tensordot(S,prm%cleavage_systems(1:3,1:3,2,i)) traction_n = math_tensordot(S,prm%cleavage_systems(1:3,1:3,3,i)) - traction_crit = prm%g_crit(i)*damagestate_h(homog)%state(1,damageOffset)**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)%p(sourceOffset)%dotState(1,me) & diff --git a/src/phase_damage_anisoductile.f90 b/src/phase_damage_anisoductile.f90 index 18c29fdfc..b22f5b3d6 100644 --- a/src/phase_damage_anisoductile.f90 +++ b/src/phase_damage_anisoductile.f90 @@ -129,7 +129,7 @@ module subroutine anisoductile_dotState(co, ip, el) associate(prm => param(source_damage_anisoDuctile_instance(ph))) damageState(ph)%p(sourceOffset)%dotState(1,me) & - = sum(plasticState(ph)%slipRate(:,me)/(damagestate_h(homog)%state(1,damageOffset)**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 subroutine anisoductile_dotState diff --git a/src/phase_damage_isoductile.f90 b/src/phase_damage_isoductile.f90 index 9dac17cd8..e9aec1944 100644 --- a/src/phase_damage_isoductile.f90 +++ b/src/phase_damage_isoductile.f90 @@ -120,7 +120,7 @@ module subroutine isoductile_dotState(co, ip, el) associate(prm => param(source_damage_isoDuctile_instance(ph))) damageState(ph)%p(sourceOffset)%dotState(1,me) = & - sum(plasticState(ph)%slipRate(:,me))/(damagestate_h(homog)%state(1,damageOffset)**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 subroutine isoductile_dotState diff --git a/src/phase_mechanical.f90 b/src/phase_mechanical.f90 index f2f3333d3..61151d33b 100644 --- a/src/phase_mechanical.f90 +++ b/src/phase_mechanical.f90 @@ -381,7 +381,7 @@ subroutine phase_hooke_SandItsTangents(S, dS_dFe, dS_dFi, & DegradationLoop: do d = 1, phase_NstiffnessDegradations(material_phaseAt(co,el)) degradationType: select case(phase_stiffnessDegradation(d,material_phaseAt(co,el))) case (STIFFNESS_DEGRADATION_damage_ID) degradationType - C = C * damagestate_h(ho)%state(1,material_homogenizationMemberAt(ip,el))**2 + C = C * phase_damage_get_phi(co,ip,el)**2 end select degradationType enddo DegradationLoop diff --git a/src/phase_mechanical_eigen_cleavageopening.f90 b/src/phase_mechanical_eigen_cleavageopening.f90 index b41cb1e9c..55102e196 100644 --- a/src/phase_mechanical_eigen_cleavageopening.f90 +++ b/src/phase_mechanical_eigen_cleavageopening.f90 @@ -126,7 +126,7 @@ module subroutine kinematics_cleavage_opening_LiAndItsTangent(Ld, dLd_dTstar, S, dLd_dTstar = 0.0_pReal associate(prm => param(kinematics_cleavage_opening_instance(material_phaseAt(co,el)))) do i = 1,prm%sum_N_cl - traction_crit = prm%g_crit(i)*damagestate_h(homog)%state(1,damageOffset)**2.0_pReal + traction_crit = prm%g_crit(i)*phase_damage_get_phi(co,ip,el)**2.0_pReal traction_d = math_tensordot(S,prm%cleavage_systems(1:3,1:3,1,i)) if (abs(traction_d) > traction_crit + tol_math_check) then diff --git a/src/phase_mechanical_eigen_slipplaneopening.f90 b/src/phase_mechanical_eigen_slipplaneopening.f90 index 682f0b825..cf5e918b3 100644 --- a/src/phase_mechanical_eigen_slipplaneopening.f90 +++ b/src/phase_mechanical_eigen_slipplaneopening.f90 @@ -152,7 +152,7 @@ module subroutine kinematics_slipplane_opening_LiAndItsTangent(Ld, dLd_dTstar, S traction_t = math_tensordot(S,prm%P_t(1:3,1:3,i)) traction_n = math_tensordot(S,prm%P_n(1:3,1:3,i)) - traction_crit = prm%g_crit(i)* damagestate_h(homog)%state(1,damageOffset) ! degrading critical load carrying capacity by damage + traction_crit = prm%g_crit(i)* phase_damage_get_phi(co,ip,el) ! degrading critical load carrying capacity by damage udotd = sign(1.0_pReal,traction_d)* prm%dot_o* ( abs(traction_d)/traction_crit & - abs(traction_d)/prm%g_crit(i))**prm%q From a09989fe0b53946e7607a2ddb309db6af31c6bb4 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 12 Feb 2021 12:35:50 +0100 Subject: [PATCH 05/25] homogenized damage only needed in homogenization --- src/damage_nonlocal.f90 | 69 ----------------------------------- src/homogenization.f90 | 79 +++++++++++++++++++++++++++++++++++++++++ src/material.f90 | 7 ++-- 3 files changed, 81 insertions(+), 74 deletions(-) diff --git a/src/damage_nonlocal.f90 b/src/damage_nonlocal.f90 index dd50cf4cf..fab08d500 100644 --- a/src/damage_nonlocal.f90 +++ b/src/damage_nonlocal.f90 @@ -14,79 +14,10 @@ module damage_nonlocal implicit none private - type, private :: tNumerics - real(pReal) :: & - charLength !< characteristic length scale for gradient problems - end type tNumerics - type(tNumerics), private :: & - num - - public :: & - damage_nonlocal_init, & - damage_nonlocal_getDiffusion contains -!-------------------------------------------------------------------------------------------------- -!> @brief module initialization -!> @details reads in material parameters, allocates arrays, and does sanity checks -!-------------------------------------------------------------------------------------------------- -subroutine damage_nonlocal_init - - integer :: Ninstances,Nmaterialpoints,h - class(tNode), pointer :: & - num_generic, & - material_homogenization - - print'(/,a)', ' <<<+- damage_nonlocal init -+>>>'; flush(6) - -!------------------------------------------------------------------------------------ -! read numerics parameter - num_generic => config_numerics%get('generic',defaultVal= emptyDict) - num%charLength = num_generic%get_asFloat('charLength',defaultVal=1.0_pReal) - - Ninstances = count(damage_type == DAMAGE_nonlocal_ID) - - material_homogenization => config_material%get('homogenization') - do h = 1, material_homogenization%length - if (damage_type(h) /= DAMAGE_NONLOCAL_ID) cycle - - 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)%state (1,Nmaterialpoints), source=1.0_pReal) - - enddo - -end subroutine damage_nonlocal_init - - -!-------------------------------------------------------------------------------------------------- -!> @brief returns homogenized non local damage diffusion tensor in reference configuration -!-------------------------------------------------------------------------------------------------- -function damage_nonlocal_getDiffusion(ip,el) - - integer, intent(in) :: & - ip, & !< integration point number - el !< element number - real(pReal), dimension(3,3) :: & - damage_nonlocal_getDiffusion - integer :: & - homog, & - grain - - homog = material_homogenizationAt(el) - damage_nonlocal_getDiffusion = 0.0_pReal - do grain = 1, homogenization_Nconstituents(homog) - damage_nonlocal_getDiffusion = damage_nonlocal_getDiffusion + & - crystallite_push33ToRef(grain,ip,el,lattice_D(1:3,1:3,material_phaseAt(grain,el))) - enddo - - damage_nonlocal_getDiffusion = & - num%charLength**2*damage_nonlocal_getDiffusion/real(homogenization_Nconstituents(homog),pReal) - -end function damage_nonlocal_getDiffusion end module damage_nonlocal diff --git a/src/homogenization.f90 b/src/homogenization.f90 index 3d0a4bb7a..36c03a314 100644 --- a/src/homogenization.f90 +++ b/src/homogenization.f90 @@ -15,10 +15,20 @@ module homogenization use damage_nonlocal use HDF5_utilities use results + use lattice implicit none private + type, private :: tNumerics_damage + real(pReal) :: & + charLength !< characteristic length scale for gradient problems + end type tNumerics_damage + + type(tNumerics_damage), private :: & + num_damage + + logical, public :: & terminallyIll = .false. !< at least one material point is terminally ill @@ -194,6 +204,10 @@ module homogenization homogenization_restartRead, & homogenization_restartWrite + public :: & + damage_nonlocal_init, & + damage_nonlocal_getDiffusion + contains @@ -208,6 +222,9 @@ subroutine homogenization_init() print'(/,a)', ' <<<+- homogenization init -+>>>'; flush(IO_STDOUT) + call material_parseHomogenization + print*, 'Homogenization parsed' + num_homog => config_numerics%get('homogenization',defaultVal=emptyDict) num_homogGeneric => num_homog%get('generic',defaultVal=emptyDict) @@ -437,4 +454,66 @@ subroutine homogenization_restartRead(fileHandle) end subroutine homogenization_restartRead + +!-------------------------------------------------------------------------------------------------- +!> @brief module initialization +!> @details reads in material parameters, allocates arrays, and does sanity checks +!-------------------------------------------------------------------------------------------------- +subroutine damage_nonlocal_init + + integer :: Ninstances,Nmaterialpoints,h + class(tNode), pointer :: & + num_generic, & + material_homogenization + + print'(/,a)', ' <<<+- damage_nonlocal init -+>>>'; flush(6) + +!------------------------------------------------------------------------------------ +! read numerics parameter + num_generic => config_numerics%get('generic',defaultVal= emptyDict) + num_damage%charLength = num_generic%get_asFloat('charLength',defaultVal=1.0_pReal) + + Ninstances = count(damage_type == DAMAGE_nonlocal_ID) + + material_homogenization => config_material%get('homogenization') + do h = 1, material_homogenization%length + if (damage_type(h) /= DAMAGE_NONLOCAL_ID) cycle + + 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)%state (1,Nmaterialpoints), source=1.0_pReal) + + enddo + +end subroutine damage_nonlocal_init + + +!-------------------------------------------------------------------------------------------------- +!> @brief returns homogenized non local damage diffusion tensor in reference configuration +!-------------------------------------------------------------------------------------------------- +function damage_nonlocal_getDiffusion(ip,el) + + integer, intent(in) :: & + ip, & !< integration point number + el !< element number + real(pReal), dimension(3,3) :: & + damage_nonlocal_getDiffusion + integer :: & + homog, & + grain + + homog = material_homogenizationAt(el) + damage_nonlocal_getDiffusion = 0.0_pReal + do grain = 1, homogenization_Nconstituents(homog) + damage_nonlocal_getDiffusion = damage_nonlocal_getDiffusion + & + crystallite_push33ToRef(grain,ip,el,lattice_D(1:3,1:3,material_phaseAt(grain,el))) + enddo + + damage_nonlocal_getDiffusion = & + num_damage%charLength**2*damage_nonlocal_getDiffusion/real(homogenization_Nconstituents(homog),pReal) + +end function damage_nonlocal_getDiffusion + + end module homogenization diff --git a/src/material.f90 b/src/material.f90 index b76b2b296..1c8292e5f 100644 --- a/src/material.f90 +++ b/src/material.f90 @@ -78,7 +78,8 @@ module material DAMAGE_NONLOCAL_ID, & HOMOGENIZATION_NONE_ID, & HOMOGENIZATION_ISOSTRAIN_ID, & - HOMOGENIZATION_RGC_ID + HOMOGENIZATION_RGC_ID, & + material_parseHomogenization contains @@ -95,10 +96,6 @@ subroutine material_init(restart) call material_parseMaterial print*, 'Material parsed' - call material_parseHomogenization - print*, 'Homogenization parsed' - - allocate(homogState (size(material_name_homogenization))) allocate(damageState_h (size(material_name_homogenization))) From 4eb2a981cab505e3a53be2a907d9e7d4aa3bf204 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 12 Feb 2021 15:31:43 +0100 Subject: [PATCH 06/25] keeping variables local --- src/commercialFEM_fileList.f90 | 1 - src/damage_nonlocal.f90 | 23 ------ src/grid/grid_damage_spectral.f90 | 1 - src/homogenization.f90 | 116 +++++++++++++++++++++++++++- src/homogenization_mechanical.f90 | 4 +- src/material.f90 | 123 ++---------------------------- 6 files changed, 120 insertions(+), 148 deletions(-) delete mode 100644 src/damage_nonlocal.f90 diff --git a/src/commercialFEM_fileList.f90 b/src/commercialFEM_fileList.f90 index cf48e19ce..97e7520f9 100644 --- a/src/commercialFEM_fileList.f90 +++ b/src/commercialFEM_fileList.f90 @@ -41,7 +41,6 @@ #include "phase_damage_isoductile.f90" #include "phase_damage_anisobrittle.f90" #include "phase_damage_anisoductile.f90" -#include "damage_nonlocal.f90" #include "homogenization.f90" #include "homogenization_mechanical.f90" #include "homogenization_mechanical_pass.f90" diff --git a/src/damage_nonlocal.f90 b/src/damage_nonlocal.f90 deleted file mode 100644 index fab08d500..000000000 --- a/src/damage_nonlocal.f90 +++ /dev/null @@ -1,23 +0,0 @@ -!-------------------------------------------------------------------------------------------------- -!> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH -!> @brief material subroutine for non-locally evolving damage field -!-------------------------------------------------------------------------------------------------- -module damage_nonlocal - use prec - use material - use config - use YAML_types - use lattice - use phase - use results - - implicit none - private - - - -contains - - - -end module damage_nonlocal diff --git a/src/grid/grid_damage_spectral.f90 b/src/grid/grid_damage_spectral.f90 index dc78d9e44..afbae027f 100644 --- a/src/grid/grid_damage_spectral.f90 +++ b/src/grid/grid_damage_spectral.f90 @@ -15,7 +15,6 @@ module grid_damage_spectral use IO use spectral_utilities use discretization_grid - use damage_nonlocal use YAML_types use homogenization use config diff --git a/src/homogenization.f90 b/src/homogenization.f90 index 36c03a314..073482c9a 100644 --- a/src/homogenization.f90 +++ b/src/homogenization.f90 @@ -12,7 +12,6 @@ module homogenization use material use phase use discretization - use damage_nonlocal use HDF5_utilities use results use lattice @@ -20,6 +19,37 @@ module homogenization implicit none private + + enum, bind(c); enumerator :: & + THERMAL_ISOTHERMAL_ID, & + THERMAL_CONDUCTION_ID, & + DAMAGE_NONE_ID, & + DAMAGE_NONLOCAL_ID, & + HOMOGENIZATION_UNDEFINED_ID, & + HOMOGENIZATION_NONE_ID, & + HOMOGENIZATION_ISOSTRAIN_ID, & + HOMOGENIZATION_RGC_ID + end enum + + type(tState), allocatable, dimension(:), public :: & + homogState, & + damageState_h + + integer, dimension(:), allocatable, public, protected :: & + homogenization_typeInstance, & !< instance of particular type of each homogenization + thermal_typeInstance, & !< instance of particular type of each thermal transport + damage_typeInstance !< instance of particular type of each nonlocal damage + + real(pReal), dimension(:), allocatable, public, protected :: & + thermal_initialT + + integer(kind(THERMAL_isothermal_ID)), dimension(:), allocatable, public, protected :: & + thermal_type !< thermal transport model + integer(kind(DAMAGE_none_ID)), dimension(:), allocatable, public, protected :: & + damage_type !< nonlocal damage model + integer(kind(HOMOGENIZATION_undefined_ID)), dimension(:), allocatable, public, protected :: & + homogenization_type !< type of each homogenization + type, private :: tNumerics_damage real(pReal) :: & charLength !< characteristic length scale for gradient problems @@ -202,7 +232,9 @@ module homogenization homogenization_forward, & homogenization_results, & homogenization_restartRead, & - homogenization_restartWrite + homogenization_restartWrite, & + THERMAL_CONDUCTION_ID, & + DAMAGE_NONLOCAL_ID public :: & damage_nonlocal_init, & @@ -222,6 +254,9 @@ subroutine homogenization_init() print'(/,a)', ' <<<+- homogenization init -+>>>'; flush(IO_STDOUT) + + allocate(homogState (size(material_name_homogenization))) + allocate(damageState_h (size(material_name_homogenization))) call material_parseHomogenization print*, 'Homogenization parsed' @@ -516,4 +551,81 @@ function damage_nonlocal_getDiffusion(ip,el) end function damage_nonlocal_getDiffusion +!-------------------------------------------------------------------------------------------------- +!> @brief parses the homogenization part from the material configuration +! ToDo: This should be done in homogenization +!-------------------------------------------------------------------------------------------------- +subroutine material_parseHomogenization + + class(tNode), pointer :: & + material_homogenization, & + homog, & + homogMech, & + homogThermal, & + homogDamage + + integer :: h + + material_homogenization => config_material%get('homogenization') + + allocate(homogenization_type(size(material_name_homogenization)), source=HOMOGENIZATION_undefined_ID) + allocate(thermal_type(size(material_name_homogenization)), source=THERMAL_isothermal_ID) + allocate(damage_type (size(material_name_homogenization)), source=DAMAGE_none_ID) + allocate(homogenization_typeInstance(size(material_name_homogenization)), source=0) + allocate(thermal_typeInstance(size(material_name_homogenization)), source=0) + allocate(damage_typeInstance(size(material_name_homogenization)), source=0) + allocate(thermal_initialT(size(material_name_homogenization)), source=300.0_pReal) + + do h=1, size(material_name_homogenization) + homog => material_homogenization%get(h) + homogMech => homog%get('mechanics') + select case (homogMech%get_asString('type')) + case('pass') + homogenization_type(h) = HOMOGENIZATION_NONE_ID + case('isostrain') + homogenization_type(h) = HOMOGENIZATION_ISOSTRAIN_ID + case('RGC') + homogenization_type(h) = HOMOGENIZATION_RGC_ID + case default + call IO_error(500,ext_msg=homogMech%get_asString('type')) + end select + + homogenization_typeInstance(h) = count(homogenization_type==homogenization_type(h)) + + if(homog%contains('thermal')) then + homogThermal => homog%get('thermal') + thermal_initialT(h) = homogThermal%get_asFloat('T_0',defaultVal=300.0_pReal) + + select case (homogThermal%get_asString('type')) + case('isothermal') + thermal_type(h) = THERMAL_isothermal_ID + case('conduction') + thermal_type(h) = THERMAL_conduction_ID + case default + call IO_error(500,ext_msg=homogThermal%get_asString('type')) + end select + endif + + if(homog%contains('damage')) then + homogDamage => homog%get('damage') + select case (homogDamage%get_asString('type')) + case('none') + damage_type(h) = DAMAGE_none_ID + case('nonlocal') + damage_type(h) = DAMAGE_nonlocal_ID + case default + call IO_error(500,ext_msg=homogDamage%get_asString('type')) + end select + endif + enddo + + do h=1, size(material_name_homogenization) + homogenization_typeInstance(h) = count(homogenization_type(1:h) == homogenization_type(h)) + thermal_typeInstance(h) = count(thermal_type (1:h) == thermal_type (h)) + damage_typeInstance(h) = count(damage_type (1:h) == damage_type (h)) + enddo + +end subroutine material_parseHomogenization + + end module homogenization diff --git a/src/homogenization_mechanical.f90 b/src/homogenization_mechanical.f90 index 91e1d3d84..c05b576da 100644 --- a/src/homogenization_mechanical.f90 +++ b/src/homogenization_mechanical.f90 @@ -225,8 +225,6 @@ end function mechanical_updateState !> @brief Write results to file. !-------------------------------------------------------------------------------------------------- module subroutine mechanical_results(group_base,h) - use material, only: & - material_homogenization_type => homogenization_type character(len=*), intent(in) :: group_base integer, intent(in) :: h @@ -236,7 +234,7 @@ module subroutine mechanical_results(group_base,h) group = trim(group_base)//'/mech' call results_closeGroup(results_addGroup(group)) - select case(material_homogenization_type(h)) + select case(homogenization_type(h)) case(HOMOGENIZATION_rgc_ID) call mechanical_RGC_results(homogenization_typeInstance(h),group) diff --git a/src/material.f90 b/src/material.f90 index 1c8292e5f..df6dc6dd1 100644 --- a/src/material.f90 +++ b/src/material.f90 @@ -6,50 +6,26 @@ !-------------------------------------------------------------------------------------------------- module material use prec - use math use config use results use IO use rotations use discretization + use YAML_types implicit none private - enum, bind(c); enumerator :: & - THERMAL_ISOTHERMAL_ID, & - THERMAL_CONDUCTION_ID, & - DAMAGE_NONE_ID, & - DAMAGE_NONLOCAL_ID, & - HOMOGENIZATION_UNDEFINED_ID, & - HOMOGENIZATION_NONE_ID, & - HOMOGENIZATION_ISOSTRAIN_ID, & - HOMOGENIZATION_RGC_ID - end enum + integer, dimension(:), allocatable, public, protected :: & + homogenization_Nconstituents !< number of grains in each homogenization character(len=:), public, protected, allocatable, dimension(:) :: & material_name_phase, & !< name of each phase material_name_homogenization !< name of each homogenization - integer(kind(THERMAL_isothermal_ID)), dimension(:), allocatable, public, protected :: & - thermal_type !< thermal transport model - integer(kind(DAMAGE_none_ID)), dimension(:), allocatable, public, protected :: & - damage_type !< nonlocal damage model - integer(kind(HOMOGENIZATION_undefined_ID)), dimension(:), allocatable, public, protected :: & - homogenization_type !< type of each homogenization - integer, public, protected :: & homogenization_maxNconstituents !< max number of grains in any USED homogenization - integer, dimension(:), allocatable, public, protected :: & - homogenization_Nconstituents, & !< number of grains in each homogenization - homogenization_typeInstance, & !< instance of particular type of each homogenization - thermal_typeInstance, & !< instance of particular type of each thermal transport - damage_typeInstance !< instance of particular type of each nonlocal damage - - real(pReal), dimension(:), allocatable, public, protected :: & - thermal_initialT !< initial temperature per each homogenization - integer, dimension(:), allocatable, public, protected :: & ! (elem) material_homogenizationAt, & !< homogenization ID of each element material_homogenizationAt2, & !< per cell @@ -63,23 +39,11 @@ module material integer, dimension(:,:,:), allocatable, public, protected :: & ! (constituent,IP,elem) material_phaseMemberAt !< position of the element within its phase instance - type(tState), allocatable, dimension(:), public :: & - homogState, & - damageState_h - type(Rotation), dimension(:,:,:), allocatable, public, protected :: & material_orientation0 !< initial orientation of each grain,IP,element public :: & - material_init, & - THERMAL_ISOTHERMAL_ID, & - THERMAL_CONDUCTION_ID, & - DAMAGE_NONE_ID, & - DAMAGE_NONLOCAL_ID, & - HOMOGENIZATION_NONE_ID, & - HOMOGENIZATION_ISOSTRAIN_ID, & - HOMOGENIZATION_RGC_ID, & - material_parseHomogenization + material_init contains @@ -90,14 +54,13 @@ subroutine material_init(restart) logical, intent(in) :: restart + print'(/,a)', ' <<<+- material init -+>>>'; flush(IO_STDOUT) call material_parseMaterial print*, 'Material parsed' - allocate(homogState (size(material_name_homogenization))) - allocate(damageState_h (size(material_name_homogenization))) if (.not. restart) then call results_openJobFile @@ -109,82 +72,6 @@ subroutine material_init(restart) end subroutine material_init -!-------------------------------------------------------------------------------------------------- -!> @brief parses the homogenization part from the material configuration -! ToDo: This should be done in homogenization -!-------------------------------------------------------------------------------------------------- -subroutine material_parseHomogenization - - class(tNode), pointer :: & - material_homogenization, & - homog, & - homogMech, & - homogThermal, & - homogDamage - - integer :: h - - material_homogenization => config_material%get('homogenization') - - allocate(homogenization_type(size(material_name_homogenization)), source=HOMOGENIZATION_undefined_ID) - allocate(thermal_type(size(material_name_homogenization)), source=THERMAL_isothermal_ID) - allocate(damage_type (size(material_name_homogenization)), source=DAMAGE_none_ID) - allocate(homogenization_typeInstance(size(material_name_homogenization)), source=0) - allocate(thermal_typeInstance(size(material_name_homogenization)), source=0) - allocate(damage_typeInstance(size(material_name_homogenization)), source=0) - allocate(thermal_initialT(size(material_name_homogenization)), source=300.0_pReal) - - do h=1, size(material_name_homogenization) - homog => material_homogenization%get(h) - homogMech => homog%get('mechanics') - select case (homogMech%get_asString('type')) - case('pass') - homogenization_type(h) = HOMOGENIZATION_NONE_ID - case('isostrain') - homogenization_type(h) = HOMOGENIZATION_ISOSTRAIN_ID - case('RGC') - homogenization_type(h) = HOMOGENIZATION_RGC_ID - case default - call IO_error(500,ext_msg=homogMech%get_asString('type')) - end select - - homogenization_typeInstance(h) = count(homogenization_type==homogenization_type(h)) - - if(homog%contains('thermal')) then - homogThermal => homog%get('thermal') - thermal_initialT(h) = homogThermal%get_asFloat('T_0',defaultVal=300.0_pReal) - - select case (homogThermal%get_asString('type')) - case('isothermal') - thermal_type(h) = THERMAL_isothermal_ID - case('conduction') - thermal_type(h) = THERMAL_conduction_ID - case default - call IO_error(500,ext_msg=homogThermal%get_asString('type')) - end select - endif - - if(homog%contains('damage')) then - homogDamage => homog%get('damage') - select case (homogDamage%get_asString('type')) - case('none') - damage_type(h) = DAMAGE_none_ID - case('nonlocal') - damage_type(h) = DAMAGE_nonlocal_ID - case default - call IO_error(500,ext_msg=homogDamage%get_asString('type')) - end select - endif - enddo - - do h=1, size(material_name_homogenization) - homogenization_typeInstance(h) = count(homogenization_type(1:h) == homogenization_type(h)) - thermal_typeInstance(h) = count(thermal_type (1:h) == thermal_type (h)) - damage_typeInstance(h) = count(damage_type (1:h) == damage_type (h)) - enddo - -end subroutine material_parseHomogenization - !-------------------------------------------------------------------------------------------------- !> @brief parses the material part in the material configuration file From 830d00fa675b960c6d5225638c48585915e2c04b Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 12 Feb 2021 22:41:35 +0100 Subject: [PATCH 07/25] simplified structure choice of damage model triggers eigendeformation, no repeated variables This implementation is the most ugly hack I could imagine. I just serves the purpose of having a stable material.yaml --- PRIVATE | 2 +- src/phase_damage.f90 | 4 +- src/phase_damage_anisobrittle.f90 | 4 +- src/phase_damage_isobrittle.f90 | 4 +- src/phase_damage_isoductile.f90 | 4 +- src/phase_mechanical_eigen.f90 | 44 ++++++++++++++++++- ...phase_mechanical_eigen_cleavageopening.f90 | 4 +- ...hase_mechanical_eigen_slipplaneopening.f90 | 4 +- 8 files changed, 55 insertions(+), 15 deletions(-) diff --git a/PRIVATE b/PRIVATE index e74cf0062..07d125fe4 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit e74cf00628285a587ced1e551cc9837c1011ca1c +Subproject commit 07d125fe476621224d9f2ef73b0cf88f3b07bc60 diff --git a/src/phase_damage.f90 b/src/phase_damage.f90 index 57145c550..6063f3583 100644 --- a/src/phase_damage.f90 +++ b/src/phase_damage.f90 @@ -175,7 +175,7 @@ module subroutine damage_init allocate(current(ph)%d_phi_d_dot_phi(Nconstituents),source=0.0_pReal) phase => phases%get(ph) - sources => phase%get('source',defaultVal=emptyList) + sources => phase%get('damage',defaultVal=emptyList) phase_Nsources(ph) = sources%length allocate(damageState(ph)%p(phase_Nsources(ph))) enddo @@ -493,7 +493,7 @@ function source_active(source_label,src_length) result(active_source) allocate(active_source(src_length,phases%length), source = .false. ) do p = 1, phases%length phase => phases%get(p) - sources => phase%get('source',defaultVal=emptyList) + sources => phase%get('damage',defaultVal=emptyList) do s = 1, sources%length src => sources%get(s) if(src%get_asString('type') == source_label) active_source(s,p) = .true. diff --git a/src/phase_damage_anisobrittle.f90 b/src/phase_damage_anisobrittle.f90 index 699e9508b..ad7f80604 100644 --- a/src/phase_damage_anisobrittle.f90 +++ b/src/phase_damage_anisobrittle.f90 @@ -51,7 +51,7 @@ module function anisobrittle_init(source_length) result(mySources) print'(/,a)', ' <<<+- phase:damage:anisobrittle init -+>>>' - mySources = source_active('damage_anisoBrittle',source_length) + mySources = source_active('anisobrittle',source_length) Ninstances = count(mySources) print'(a,i2)', ' # instances: ',Ninstances; flush(IO_STDOUT) if(Ninstances == 0) return @@ -65,7 +65,7 @@ module function anisobrittle_init(source_length) result(mySources) 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') + sources => phase%get('damage') do sourceOffset = 1, sources%length if(mySources(sourceOffset,p)) then source_damage_anisoBrittle_offset(p) = sourceOffset diff --git a/src/phase_damage_isobrittle.f90 b/src/phase_damage_isobrittle.f90 index 529ecd404..a2ae7dfb7 100644 --- a/src/phase_damage_isobrittle.f90 +++ b/src/phase_damage_isobrittle.f90 @@ -41,7 +41,7 @@ module function isobrittle_init(source_length) result(mySources) print'(/,a)', ' <<<+- phase:damage:isobrittle init -+>>>' - mySources = source_active('damage_isoBrittle',source_length) + mySources = source_active('isobrittle',source_length) Ninstances = count(mySources) print'(a,i2)', ' # instances: ',Ninstances; flush(IO_STDOUT) if(Ninstances == 0) return @@ -55,7 +55,7 @@ module function isobrittle_init(source_length) result(mySources) 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') + sources => phase%get('damage') do sourceOffset = 1, sources%length if(mySources(sourceOffset,p)) then source_damage_isoBrittle_offset(p) = sourceOffset diff --git a/src/phase_damage_isoductile.f90 b/src/phase_damage_isoductile.f90 index e9aec1944..f378b7ad2 100644 --- a/src/phase_damage_isoductile.f90 +++ b/src/phase_damage_isoductile.f90 @@ -43,7 +43,7 @@ module function isoductile_init(source_length) result(mySources) print'(/,a)', ' <<<+- phase:damage:isoductile init -+>>>' - mySources = source_active('damage_isoDuctile',source_length) + mySources = source_active('isoductile',source_length) Ninstances = count(mySources) print'(a,i2)', ' # instances: ',Ninstances; flush(IO_STDOUT) if(Ninstances == 0) return @@ -57,7 +57,7 @@ module function isoductile_init(source_length) result(mySources) 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') + sources => phase%get('damage') do sourceOffset = 1, sources%length if(mySources(sourceOffset,p)) then source_damage_isoDuctile_offset(p) = sourceOffset diff --git a/src/phase_mechanical_eigen.f90 b/src/phase_mechanical_eigen.f90 index 6df993565..3fb2d1bf2 100644 --- a/src/phase_mechanical_eigen.f90 +++ b/src/phase_mechanical_eigen.f90 @@ -66,7 +66,8 @@ module subroutine eigendeformation_init(phases) ph class(tNode), pointer :: & phase, & - kinematics + kinematics, & + damage print'(/,a)', ' <<<+- phase:mechanics:eigendeformation init -+>>>' @@ -77,6 +78,12 @@ module subroutine eigendeformation_init(phases) phase => phases%get(ph) kinematics => phase%get('kinematics',defaultVal=emptyList) phase_Nkinematics(ph) = kinematics%length + kinematics => phase%get('damage',defaultVal=emptyList) + if(kinematics%length >0) then + damage => kinematics%get(1) + if(damage%get_asString('type') == 'anisobrittle') phase_Nkinematics(ph) = phase_Nkinematics(ph) +1 + if(damage%get_asString('type') == 'isoductile') phase_Nkinematics(ph) = phase_Nkinematics(ph) +1 + endif enddo allocate(phase_kinematics(maxval(phase_Nkinematics),phases%length), source = KINEMATICS_undefined_ID) @@ -90,7 +97,6 @@ module subroutine eigendeformation_init(phases) end subroutine eigendeformation_init - !-------------------------------------------------------------------------------------------------- !> @brief checks if a kinematic mechanism is active or not !-------------------------------------------------------------------------------------------------- @@ -122,6 +128,38 @@ function kinematics_active(kinematics_label,kinematics_length) result(active_ki end function kinematics_active + +!-------------------------------------------------------------------------------------------------- +!> @brief checks if a kinematic mechanism is active or not +!-------------------------------------------------------------------------------------------------- +function kinematics_active2(kinematics_label,kinematics_length) result(active_kinematics) + + character(len=*), intent(in) :: kinematics_label !< name of kinematic mechanism + integer, intent(in) :: kinematics_length !< max. number of kinematics in system + logical, dimension(:,:), allocatable :: active_kinematics + + class(tNode), pointer :: & + phases, & + phase, & + kinematics, & + kinematics_type + integer :: p,k + + phases => config_material%get('phase') + allocate(active_kinematics(kinematics_length,phases%length), source = .false. ) + do p = 1, phases%length + phase => phases%get(p) + kinematics => phase%get('damage',defaultVal=emptyList) + do k = 1, kinematics%length + kinematics_type => kinematics%get(k) + if(kinematics_type%get_asString('type') == kinematics_label) active_kinematics(k,p) = .true. + enddo + enddo + + +end function kinematics_active2 + + !-------------------------------------------------------------------------------------------------- !> @brief contains the constitutive equation for calculating the velocity gradient ! ToDo: MD: S is Mi? @@ -175,8 +213,10 @@ module subroutine phase_LiAndItsTangents(Li, dLi_dS, dLi_dFi, & KinematicsLoop: do k = 1, phase_Nkinematics(material_phaseAt(co,el)) kinematicsType: select case (phase_kinematics(k,material_phaseAt(co,el))) case (KINEMATICS_cleavage_opening_ID) kinematicsType + print*, 'clea' call kinematics_cleavage_opening_LiAndItsTangent(my_Li, my_dLi_dS, S, co, ip, el) case (KINEMATICS_slipplane_opening_ID) kinematicsType + print*, 'slipp' call kinematics_slipplane_opening_LiAndItsTangent(my_Li, my_dLi_dS, S, co, ip, el) case (KINEMATICS_thermal_expansion_ID) kinematicsType me = material_phaseMemberAt(co,ip,el) diff --git a/src/phase_mechanical_eigen_cleavageopening.f90 b/src/phase_mechanical_eigen_cleavageopening.f90 index 55102e196..f1774c78e 100644 --- a/src/phase_mechanical_eigen_cleavageopening.f90 +++ b/src/phase_mechanical_eigen_cleavageopening.f90 @@ -46,7 +46,7 @@ module function kinematics_cleavage_opening_init(kinematics_length) result(myKin print'(/,a)', ' <<<+- phase:mechanics:eigendeformation:cleavageopening init -+>>>' - myKinematics = kinematics_active('cleavage_opening',kinematics_length) + myKinematics = kinematics_active2('anisobrittle',kinematics_length) Ninstances = count(myKinematics) print'(a,i2)', ' # instances: ',Ninstances; flush(IO_STDOUT) if(Ninstances == 0) return @@ -59,7 +59,7 @@ module function kinematics_cleavage_opening_init(kinematics_length) result(myKin if(any(myKinematics(:,p))) kinematics_cleavage_opening_instance(p) = count(myKinematics(:,1:p)) phase => phases%get(p) if(count(myKinematics(:,p)) == 0) cycle - kinematics => phase%get('kinematics') + kinematics => phase%get('damage') do k = 1, kinematics%length if(myKinematics(k,p)) then associate(prm => param(kinematics_cleavage_opening_instance(p))) diff --git a/src/phase_mechanical_eigen_slipplaneopening.f90 b/src/phase_mechanical_eigen_slipplaneopening.f90 index cf5e918b3..dd5f172d9 100644 --- a/src/phase_mechanical_eigen_slipplaneopening.f90 +++ b/src/phase_mechanical_eigen_slipplaneopening.f90 @@ -51,7 +51,7 @@ module function kinematics_slipplane_opening_init(kinematics_length) result(myKi print'(/,a)', ' <<<+- phase:mechanics:eigendeformation:slipplaneopening init -+>>>' - myKinematics = kinematics_active('slipplane_opening',kinematics_length) + myKinematics = kinematics_active2('isoductile',kinematics_length) Ninstances = count(myKinematics) print'(a,i2)', ' # instances: ',Ninstances; flush(IO_STDOUT) if(Ninstances == 0) return @@ -66,7 +66,7 @@ module function kinematics_slipplane_opening_init(kinematics_length) result(myKi mech => phase%get('mechanics') pl => mech%get('plasticity') if(count(myKinematics(:,p)) == 0) cycle - kinematics => phase%get('kinematics') + kinematics => phase%get('damage') do k = 1, kinematics%length if(myKinematics(k,p)) then associate(prm => param(kinematics_slipplane_opening_instance(p))) From 2b0b1aeffe0c18e0350ebcf2e9039cae98151724 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 12 Feb 2021 23:35:06 +0100 Subject: [PATCH 08/25] jobname.yaml not supported anymore --- PRIVATE | 2 +- src/config.f90 | 23 +++++++---------------- 2 files changed, 8 insertions(+), 17 deletions(-) diff --git a/PRIVATE b/PRIVATE index 07d125fe4..2a0f30ec4 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit 07d125fe476621224d9f2ef73b0cf88f3b07bc60 +Subproject commit 2a0f30ec47473f42e0da922055b72b527730378d diff --git a/src/config.f90 b/src/config.f90 index b10edf013..02b16f2a2 100644 --- a/src/config.f90 +++ b/src/config.f90 @@ -5,16 +5,10 @@ !! precedence over material.yaml. !-------------------------------------------------------------------------------------------------- module config - use prec - use DAMASK_interface use IO use YAML_parse use YAML_types -#ifdef PETSc -#include - use petscsys -#endif implicit none private @@ -50,17 +44,12 @@ end subroutine config_init subroutine parse_material logical :: fileExists - character(len=:), allocatable :: fname - fname = getSolverJobName()//'.yaml' - inquire(file=fname,exist=fileExists) - if(.not. fileExists) then - fname = 'material.yaml' - inquire(file=fname,exist=fileExists) - if(.not. fileExists) call IO_error(100,ext_msg=fname) - endif - print*, 'reading '//fname; flush(IO_STDOUT) - config_material => YAML_parse_file(fname) + + inquire(file='material.yaml',exist=fileExists) + if(.not. fileExists) call IO_error(100,ext_msg='material.yaml') + print*, 'reading material.yaml'; flush(IO_STDOUT) + config_material => YAML_parse_file('material.yaml') end subroutine parse_material @@ -72,6 +61,7 @@ subroutine parse_numerics logical :: fexist + config_numerics => emptyDict inquire(file='numerics.yaml', exist=fexist) if (fexist) then @@ -89,6 +79,7 @@ subroutine parse_debug logical :: fexist + config_debug => emptyDict inquire(file='debug.yaml', exist=fexist) fileExists: if (fexist) then From c790c82a427de2a024e935af18488ee9dcc6c85d Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 13 Feb 2021 06:59:18 +0100 Subject: [PATCH 09/25] separating by instance just complicates things --- src/material.f90 | 2 +- src/phase_mechanical_eigen_cleavageopening.f90 | 15 ++++----------- src/phase_mechanical_eigen_slipplaneopening.f90 | 14 +++----------- 3 files changed, 8 insertions(+), 23 deletions(-) diff --git a/src/material.f90 b/src/material.f90 index df6dc6dd1..3656f0f0a 100644 --- a/src/material.f90 +++ b/src/material.f90 @@ -48,7 +48,7 @@ module material contains !-------------------------------------------------------------------------------------------------- -!> @brief parses material configuration file +!> @brief Parse material configuration file (material.yaml). !-------------------------------------------------------------------------------------------------- subroutine material_init(restart) diff --git a/src/phase_mechanical_eigen_cleavageopening.f90 b/src/phase_mechanical_eigen_cleavageopening.f90 index f1774c78e..8878c10bd 100644 --- a/src/phase_mechanical_eigen_cleavageopening.f90 +++ b/src/phase_mechanical_eigen_cleavageopening.f90 @@ -6,8 +6,6 @@ !-------------------------------------------------------------------------------------------------- submodule(phase:eigendeformation) cleavageopening - integer, dimension(:), allocatable :: kinematics_cleavage_opening_instance - type :: tParameters !< container type for internal constitutive parameters integer :: & sum_N_cl !< total number of cleavage planes @@ -20,7 +18,7 @@ submodule(phase:eigendeformation) cleavageopening cleavage_systems end type tParameters - type(tParameters), dimension(:), allocatable :: param !< containers of constitutive parameters (len Ninstances) + type(tParameters), dimension(:), allocatable :: param !< containers of constitutive parameters contains @@ -52,17 +50,15 @@ module function kinematics_cleavage_opening_init(kinematics_length) result(myKin if(Ninstances == 0) return phases => config_material%get('phase') - allocate(param(Ninstances)) - allocate(kinematics_cleavage_opening_instance(phases%length), source=0) + allocate(param(phases%length)) do p = 1, phases%length - if(any(myKinematics(:,p))) kinematics_cleavage_opening_instance(p) = count(myKinematics(:,1:p)) phase => phases%get(p) if(count(myKinematics(:,p)) == 0) cycle kinematics => phase%get('damage') do k = 1, kinematics%length if(myKinematics(k,p)) then - associate(prm => param(kinematics_cleavage_opening_instance(p))) + associate(prm => param(p)) kinematic_type => kinematics%get(k) N_cl = kinematic_type%get_asInts('N_cl') @@ -113,18 +109,15 @@ module subroutine kinematics_cleavage_opening_LiAndItsTangent(Ld, dLd_dTstar, S, dLd_dTstar !< derivative of Ld with respect to Tstar (4th-order tensor) integer :: & - homog, damageOffset, & i, k, l, m, n real(pReal) :: & traction_d, traction_t, traction_n, traction_crit, & udotd, dudotd_dt, udott, dudott_dt, udotn, dudotn_dt - homog = material_homogenizationAt(el) - damageOffset = material_homogenizationMemberAt(ip,el) Ld = 0.0_pReal dLd_dTstar = 0.0_pReal - associate(prm => param(kinematics_cleavage_opening_instance(material_phaseAt(co,el)))) + associate(prm => param(material_phaseAt(co,el))) do i = 1,prm%sum_N_cl traction_crit = prm%g_crit(i)*phase_damage_get_phi(co,ip,el)**2.0_pReal diff --git a/src/phase_mechanical_eigen_slipplaneopening.f90 b/src/phase_mechanical_eigen_slipplaneopening.f90 index dd5f172d9..e8410a935 100644 --- a/src/phase_mechanical_eigen_slipplaneopening.f90 +++ b/src/phase_mechanical_eigen_slipplaneopening.f90 @@ -57,11 +57,9 @@ module function kinematics_slipplane_opening_init(kinematics_length) result(myKi if(Ninstances == 0) return phases => config_material%get('phase') - allocate(kinematics_slipplane_opening_instance(phases%length), source=0) - allocate(param(Ninstances)) + allocate(param(phases%length)) do p = 1, phases%length - if(any(myKinematics(:,p))) kinematics_slipplane_opening_instance(p) = count(myKinematics(:,1:p)) phase => phases%get(p) mech => phase%get('mechanics') pl => mech%get('plasticity') @@ -69,7 +67,7 @@ module function kinematics_slipplane_opening_init(kinematics_length) result(myKi kinematics => phase%get('damage') do k = 1, kinematics%length if(myKinematics(k,p)) then - associate(prm => param(kinematics_slipplane_opening_instance(p))) + associate(prm => param(p)) kinematic_type => kinematics%get(k) prm%dot_o = kinematic_type%get_asFloat('dot_o') @@ -131,19 +129,13 @@ module subroutine kinematics_slipplane_opening_LiAndItsTangent(Ld, dLd_dTstar, S dLd_dTstar !< derivative of Ld with respect to Tstar (4th-order tensor) integer :: & - instance, phase, & - homog, damageOffset, & i, k, l, m, n real(pReal) :: & traction_d, traction_t, traction_n, traction_crit, & udotd, dudotd_dt, udott, dudott_dt, udotn, dudotn_dt - phase = material_phaseAt(co,el) - instance = kinematics_slipplane_opening_instance(phase) - homog = material_homogenizationAt(el) - damageOffset = material_homogenizationMemberAt(ip,el) - associate(prm => param(instance)) + associate(prm => param(material_phaseAt(co,el))) Ld = 0.0_pReal dLd_dTstar = 0.0_pReal do i = 1, prm%sum_N_sl From f95e3bc08d58f7434a23c730c02c1585bc30035a Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 13 Feb 2021 07:15:57 +0100 Subject: [PATCH 10/25] simplified access pattern --- src/phase.f90 | 5 +++ src/phase_damage.f90 | 11 +++++ src/phase_mechanical.f90 | 11 ++--- src/phase_mechanical_eigen.f90 | 43 ++++++------------- ...phase_mechanical_eigen_cleavageopening.f90 | 10 ++--- ...hase_mechanical_eigen_slipplaneopening.f90 | 10 ++--- 6 files changed, 42 insertions(+), 48 deletions(-) diff --git a/src/phase.f90 b/src/phase.f90 index 628bfc443..82b82d71d 100644 --- a/src/phase.f90 +++ b/src/phase.f90 @@ -181,6 +181,11 @@ module phase real(pReal) :: dot_T end function thermal_dot_T + module function damage_phi(ph,me) result(phi) + integer, intent(in) :: ph,me + real(pReal) :: phi + end function damage_phi + module subroutine phase_mechanical_setF(F,co,ip,el) real(pReal), dimension(3,3), intent(in) :: F diff --git a/src/phase_damage.f90 b/src/phase_damage.f90 index 6063f3583..8defcf3a1 100644 --- a/src/phase_damage.f90 +++ b/src/phase_damage.f90 @@ -528,4 +528,15 @@ module function phase_damage_get_phi(co,ip,el) result(phi) end function phase_damage_get_phi +module function damage_phi(ph,me) result(phi) + + integer, intent(in) :: ph, me + real(pReal) :: phi + + + phi = current(ph)%phi(me) + +end function damage_phi + + end submodule damagee diff --git a/src/phase_mechanical.f90 b/src/phase_mechanical.f90 index 61151d33b..f7631a38d 100644 --- a/src/phase_mechanical.f90 +++ b/src/phase_mechanical.f90 @@ -98,12 +98,9 @@ submodule(phase) mechanics end function plastic_deltaState module subroutine phase_LiAndItsTangents(Li, dLi_dS, dLi_dFi, & - S, Fi, co, ip, el) - + S, Fi, ph,me) integer, intent(in) :: & - co, & !< component-ID of integration point - ip, & !< integration point - el !< element + ph,me real(pReal), intent(in), dimension(3,3) :: & S !< 2nd Piola-Kirchhoff stress real(pReal), intent(in), dimension(3,3) :: & @@ -583,7 +580,7 @@ function integrateStress(F,subFp0,subFi0,Delta_t,co,ip,el) result(broken) enddo LpLoop call phase_LiAndItsTangents(Li_constitutive, dLi_dS, dLi_dFi, & - S, Fi_new, co, ip, el) + S, Fi_new, ph,me) !* update current residuum and check for convergence of loop atol_Li = max(num%rtol_crystalliteStress * max(norm2(Liguess),norm2(Li_constitutive)), & ! absolute tolerance from largest acceptable relative error @@ -1303,7 +1300,7 @@ module function phase_mechanical_dPdF(dt,co,ip,el) result(dPdF) call phase_LiAndItsTangents(devNull,dLidS,dLidFi, & phase_mechanical_S(ph)%data(1:3,1:3,me), & phase_mechanical_Fi(ph)%data(1:3,1:3,me), & - co,ip,el) + ph,me) invFp = math_inv33(phase_mechanical_Fp(ph)%data(1:3,1:3,me)) invFi = math_inv33(phase_mechanical_Fi(ph)%data(1:3,1:3,me)) diff --git a/src/phase_mechanical_eigen.f90 b/src/phase_mechanical_eigen.f90 index 3fb2d1bf2..9eb3c5767 100644 --- a/src/phase_mechanical_eigen.f90 +++ b/src/phase_mechanical_eigen.f90 @@ -16,11 +16,8 @@ submodule(phase:mechanics) eigendeformation logical, dimension(:,:), allocatable :: myKinematics end function kinematics_thermal_expansion_init - module subroutine kinematics_cleavage_opening_LiAndItsTangent(Ld, dLd_dTstar, S, co, ip, el) - integer, intent(in) :: & - co, & !< grain number - ip, & !< integration point number - el !< element number + module subroutine kinematics_cleavage_opening_LiAndItsTangent(Ld, dLd_dTstar, S, ph,me) + integer, intent(in) :: ph, me real(pReal), intent(in), dimension(3,3) :: & S real(pReal), intent(out), dimension(3,3) :: & @@ -29,11 +26,8 @@ submodule(phase:mechanics) eigendeformation dLd_dTstar !< derivative of Ld with respect to Tstar (4th-order tensor) end subroutine kinematics_cleavage_opening_LiAndItsTangent - module subroutine kinematics_slipplane_opening_LiAndItsTangent(Ld, dLd_dTstar, S, co, ip, el) - integer, intent(in) :: & - co, & !< grain number - ip, & !< integration point number - el !< element number + module subroutine kinematics_slipplane_opening_LiAndItsTangent(Ld, dLd_dTstar, S, ph,me) + integer, intent(in) :: ph, me real(pReal), intent(in), dimension(3,3) :: & S real(pReal), intent(out), dimension(3,3) :: & @@ -44,7 +38,6 @@ submodule(phase:mechanics) eigendeformation module subroutine thermalexpansion_LiAndItsTangent(Li, dLi_dTstar, ph,me) integer, intent(in) :: ph, me - !< element number real(pReal), intent(out), dimension(3,3) :: & Li !< thermal velocity gradient real(pReal), intent(out), dimension(3,3,3,3) :: & @@ -165,12 +158,10 @@ end function kinematics_active2 ! ToDo: MD: S is Mi? !-------------------------------------------------------------------------------------------------- module subroutine phase_LiAndItsTangents(Li, dLi_dS, dLi_dFi, & - S, Fi, co, ip, el) + S, Fi, ph,me) integer, intent(in) :: & - co, & !< component-ID of integration point - ip, & !< integration point - el !< element + ph,me real(pReal), intent(in), dimension(3,3) :: & S !< 2nd Piola-Kirchhoff stress real(pReal), intent(in), dimension(3,3) :: & @@ -190,18 +181,16 @@ module subroutine phase_LiAndItsTangents(Li, dLi_dS, dLi_dFi, & real(pReal) :: & detFi integer :: & - k, i, j, & - instance, of, me, ph + k, i, j Li = 0.0_pReal dLi_dS = 0.0_pReal dLi_dFi = 0.0_pReal - plasticType: select case (phase_plasticity(material_phaseAt(co,el))) + + plasticType: select case (phase_plasticity(ph)) case (PLASTICITY_isotropic_ID) plasticType - of = material_phasememberAt(co,ip,el) - instance = phase_plasticInstance(material_phaseAt(co,el)) - call plastic_isotropic_LiAndItsTangent(my_Li, my_dLi_dS, S ,instance,of) + call plastic_isotropic_LiAndItsTangent(my_Li, my_dLi_dS, S ,phase_plasticInstance(ph),me) case default plasticType my_Li = 0.0_pReal my_dLi_dS = 0.0_pReal @@ -210,17 +199,13 @@ module subroutine phase_LiAndItsTangents(Li, dLi_dS, dLi_dFi, & Li = Li + my_Li dLi_dS = dLi_dS + my_dLi_dS - KinematicsLoop: do k = 1, phase_Nkinematics(material_phaseAt(co,el)) - kinematicsType: select case (phase_kinematics(k,material_phaseAt(co,el))) + KinematicsLoop: do k = 1, phase_Nkinematics(ph) + kinematicsType: select case (phase_kinematics(k,ph)) case (KINEMATICS_cleavage_opening_ID) kinematicsType - print*, 'clea' - call kinematics_cleavage_opening_LiAndItsTangent(my_Li, my_dLi_dS, S, co, ip, el) + call kinematics_cleavage_opening_LiAndItsTangent(my_Li, my_dLi_dS, S, ph, me) case (KINEMATICS_slipplane_opening_ID) kinematicsType - print*, 'slipp' - call kinematics_slipplane_opening_LiAndItsTangent(my_Li, my_dLi_dS, S, co, ip, el) + call kinematics_slipplane_opening_LiAndItsTangent(my_Li, my_dLi_dS, S, ph, me) case (KINEMATICS_thermal_expansion_ID) kinematicsType - me = material_phaseMemberAt(co,ip,el) - ph = material_phaseAt(co,el) call thermalexpansion_LiAndItsTangent(my_Li, my_dLi_dS, ph,me) case default kinematicsType my_Li = 0.0_pReal diff --git a/src/phase_mechanical_eigen_cleavageopening.f90 b/src/phase_mechanical_eigen_cleavageopening.f90 index 8878c10bd..5bb9847f4 100644 --- a/src/phase_mechanical_eigen_cleavageopening.f90 +++ b/src/phase_mechanical_eigen_cleavageopening.f90 @@ -95,12 +95,10 @@ end function kinematics_cleavage_opening_init !-------------------------------------------------------------------------------------------------- !> @brief contains the constitutive equation for calculating the velocity gradient !-------------------------------------------------------------------------------------------------- -module subroutine kinematics_cleavage_opening_LiAndItsTangent(Ld, dLd_dTstar, S, co, ip, el) +module subroutine kinematics_cleavage_opening_LiAndItsTangent(Ld, dLd_dTstar, S, ph,me) integer, intent(in) :: & - co, & !< grain number - ip, & !< integration point number - el !< element number + ph,me real(pReal), intent(in), dimension(3,3) :: & S real(pReal), intent(out), dimension(3,3) :: & @@ -117,9 +115,9 @@ module subroutine kinematics_cleavage_opening_LiAndItsTangent(Ld, dLd_dTstar, S, Ld = 0.0_pReal dLd_dTstar = 0.0_pReal - associate(prm => param(material_phaseAt(co,el))) + associate(prm => param(ph)) do i = 1,prm%sum_N_cl - traction_crit = prm%g_crit(i)*phase_damage_get_phi(co,ip,el)**2.0_pReal + traction_crit = prm%g_crit(i)*damage_phi(ph,me)**2.0_pReal traction_d = math_tensordot(S,prm%cleavage_systems(1:3,1:3,1,i)) if (abs(traction_d) > traction_crit + tol_math_check) then diff --git a/src/phase_mechanical_eigen_slipplaneopening.f90 b/src/phase_mechanical_eigen_slipplaneopening.f90 index e8410a935..64988d87f 100644 --- a/src/phase_mechanical_eigen_slipplaneopening.f90 +++ b/src/phase_mechanical_eigen_slipplaneopening.f90 @@ -115,12 +115,10 @@ end function kinematics_slipplane_opening_init !-------------------------------------------------------------------------------------------------- !> @brief contains the constitutive equation for calculating the velocity gradient !-------------------------------------------------------------------------------------------------- -module subroutine kinematics_slipplane_opening_LiAndItsTangent(Ld, dLd_dTstar, S, co, ip, el) +module subroutine kinematics_slipplane_opening_LiAndItsTangent(Ld, dLd_dTstar, S, ph,me) integer, intent(in) :: & - co, & !< grain number - ip, & !< integration point number - el !< element number + ph, me real(pReal), intent(in), dimension(3,3) :: & S real(pReal), intent(out), dimension(3,3) :: & @@ -135,7 +133,7 @@ module subroutine kinematics_slipplane_opening_LiAndItsTangent(Ld, dLd_dTstar, S udotd, dudotd_dt, udott, dudott_dt, udotn, dudotn_dt - associate(prm => param(material_phaseAt(co,el))) + associate(prm => param(ph)) Ld = 0.0_pReal dLd_dTstar = 0.0_pReal do i = 1, prm%sum_N_sl @@ -144,7 +142,7 @@ module subroutine kinematics_slipplane_opening_LiAndItsTangent(Ld, dLd_dTstar, S traction_t = math_tensordot(S,prm%P_t(1:3,1:3,i)) traction_n = math_tensordot(S,prm%P_n(1:3,1:3,i)) - traction_crit = prm%g_crit(i)* phase_damage_get_phi(co,ip,el) ! degrading critical load carrying capacity by damage + traction_crit = prm%g_crit(i)* damage_phi(ph,me) udotd = sign(1.0_pReal,traction_d)* prm%dot_o* ( abs(traction_d)/traction_crit & - abs(traction_d)/prm%g_crit(i))**prm%q From 6e3515982d6fe4c727731297577d4aafe1930665 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 13 Feb 2021 07:55:32 +0100 Subject: [PATCH 11/25] not needed outside of thermal --- src/phase.f90 | 2 +- src/phase_thermal.f90 | 3 +++ 2 files changed, 4 insertions(+), 1 deletion(-) diff --git a/src/phase.f90 b/src/phase.f90 index 82b82d71d..1d08f680d 100644 --- a/src/phase.f90 +++ b/src/phase.f90 @@ -71,7 +71,7 @@ module phase type(tPlasticState), allocatable, dimension(:), public :: & plasticState type(tSourceState), allocatable, dimension(:), public :: & - damageState, thermalState + damageState integer, public, protected :: & diff --git a/src/phase_thermal.f90 b/src/phase_thermal.f90 index 85ca2f7a5..5d3d89ccb 100644 --- a/src/phase_thermal.f90 +++ b/src/phase_thermal.f90 @@ -3,6 +3,9 @@ !---------------------------------------------------------------------------------------------------- submodule(phase) thermal + type(tSourceState), allocatable, dimension(:) :: & + thermalState + enum, bind(c); enumerator :: & THERMAL_UNDEFINED_ID ,& THERMAL_DISSIPATION_ID, & From b3dde6d7227b60a830c2486d068ca8b6b5b0343c Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 13 Feb 2021 08:32:43 +0100 Subject: [PATCH 12/25] 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. --- PRIVATE | 2 +- src/phase.f90 | 43 +++++++++---------- src/phase_damage.f90 | 71 +++++++++++++------------------ src/phase_damage_anisobrittle.f90 | 30 +++++-------- src/phase_damage_anisoductile.f90 | 29 ++++--------- src/phase_damage_isobrittle.f90 | 37 ++++++---------- src/phase_damage_isoductile.f90 | 29 ++++--------- src/phase_mechanical.f90 | 18 ++++---- 8 files changed, 97 insertions(+), 162 deletions(-) diff --git a/PRIVATE b/PRIVATE index 2a0f30ec4..571b06ac4 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit 2a0f30ec47473f42e0da922055b72b527730378d +Subproject commit 571b06ac4b42ed217eacfa5c71281e571bdab1df diff --git a/src/phase.f90 b/src/phase.f90 index 1d08f680d..8ea1fa5f7 100644 --- a/src/phase.f90 +++ b/src/phase.f90 @@ -70,7 +70,7 @@ module phase type(tPlasticState), allocatable, dimension(:), public :: & plasticState - type(tSourceState), allocatable, dimension(:), public :: & + type(tState), allocatable, dimension(:), public :: & damageState @@ -350,14 +350,12 @@ subroutine phase_init PhaseLoop2:do ph = 1,phases%length !-------------------------------------------------------------------------------------------------- ! partition and initialize state - plasticState(ph)%state = plasticState(ph)%state0 - forall(so = 1:phase_Nsources(ph)) - damageState(ph)%p(so)%state = damageState(ph)%p(so)%state0 - end forall - - phase_source_maxSizeDotState = max(phase_source_maxSizeDotState, & - maxval(damageState(ph)%p%sizeDotState)) + plasticState(ph)%state = plasticState(ph)%state0 + if(phase_Nsources(ph) > 0) & + damageState(ph)%state = damageState(ph)%state0 enddo PhaseLoop2 + + phase_source_maxSizeDotState = maxval(damageState%sizeDotState) phase_plasticity_maxSizeDotState = maxval(plasticState%sizeDotState) end subroutine phase_init @@ -404,15 +402,13 @@ subroutine phase_restore(ce,includeL) integer, intent(in) :: ce integer :: & - co, & !< constituent number - so + co 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)%state0(:,material_phasememberAt2(co,ce)) - enddo + if (phase_Nsources(material_phaseAt2(co,ce)) > 0) & + damageState(material_phaseAt2(co,ce))%state( :,material_phasememberAt2(co,ce)) = & + damageState(material_phaseAt2(co,ce))%state0(:,material_phasememberAt2(co,ce)) enddo call mechanical_restore(ce,includeL) @@ -426,16 +422,16 @@ end subroutine phase_restore !-------------------------------------------------------------------------------------------------- subroutine phase_forward() - integer :: ph, so + integer :: ph call mechanical_forward() call thermal_forward() do ph = 1, size(damageState) - do so = 1,phase_Nsources(ph) - damageState(ph)%p(so)%state0 = damageState(ph)%p(so)%state - enddo; enddo + if (phase_Nsources(ph) > 0) & + damageState(ph)%state0 = damageState(ph)%state + enddo end subroutine phase_forward @@ -531,9 +527,8 @@ subroutine crystallite_init() phases => config_material%get('phase') do ph = 1, phases%length - do so = 1, phase_Nsources(ph) - allocate(damageState(ph)%p(so)%subState0,source=damageState(ph)%p(so)%state0) ! ToDo: hack - enddo + if (phase_Nsources(ph) > 0) & + allocate(damageState(ph)%subState0,source=damageState(ph)%state0) ! ToDo: hack enddo print'(a42,1x,i10)', ' # of elements: ', eMax @@ -579,9 +574,9 @@ subroutine phase_windForward(ip,el) call mechanical_windForward(ph,me) - do so = 1, phase_Nsources(material_phaseAt(co,el)) - damageState(ph)%p(so)%state0(:,me) = damageState(ph)%p(so)%state(:,me) - enddo + if (phase_Nsources(ph) > 0) & + damageState(ph)%state0(:,me) = damageState(ph)%state(:,me) + enddo diff --git a/src/phase_damage.f90 b/src/phase_damage.f90 index 8defcf3a1..5f323c3a2 100644 --- a/src/phase_damage.f90 +++ b/src/phase_damage.f90 @@ -176,8 +176,8 @@ module subroutine damage_init phase => phases%get(ph) sources => phase%get('damage',defaultVal=emptyList) - phase_Nsources(ph) = sources%length - allocate(damageState(ph)%p(phase_Nsources(ph))) + if (sources%length > 1) error stop + enddo 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 ph, & me, & - so - integer, dimension(maxval(phase_Nsources)) :: & size_so real(pReal) :: & zeta real(pReal), dimension(phase_source_maxSizeDotState) :: & 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 :: & converged_ - ph = material_phaseAt(co,el) + if (phase_Nsources(ph) == 0) return me = material_phaseMemberAt(co,ip,el) converged_ = .true. broken = phase_damage_collectDotState(co,ip,el,ph,me) if(broken) return - do so = 1, phase_Nsources(ph) - size_so(so) = damageState(ph)%p(so)%sizeDotState - damageState(ph)%p(so)%state(1:size_so(so),me) = damageState(ph)%p(so)%subState0(1:size_so(so),me) & - + damageState(ph)%p(so)%dotState (1:size_so(so),me) * dt - source_dotState(1:size_so(so),2,so) = 0.0_pReal - enddo + size_so = damageState(ph)%sizeDotState + damageState(ph)%state(1:size_so,me) = damageState(ph)%subState0(1:size_so,me) & + + damageState(ph)%dotState (1:size_so,me) * dt + source_dotState(1:size_so,2) = 0.0_pReal iteration: do NiterationState = 1, num%nState - do so = 1, phase_Nsources(ph) - if(nIterationState > 1) source_dotState(1:size_so(so),2,so) = source_dotState(1:size_so(so),1,so) - source_dotState(1:size_so(so),1,so) = damageState(ph)%p(so)%dotState(:,me) - enddo + if(nIterationState > 1) source_dotState(1:size_so,2) = source_dotState(1:size_so,1) + source_dotState(1:size_so,1) = damageState(ph)%dotState(:,me) broken = phase_damage_collectDotState(co,ip,el,ph,me) if(broken) exit iteration - do so = 1, phase_Nsources(ph) - zeta = damper(damageState(ph)%p(so)%dotState(:,me), & - source_dotState(1:size_so(so),1,so),& - source_dotState(1:size_so(so),2,so)) - damageState(ph)%p(so)%dotState(:,me) = damageState(ph)%p(so)%dotState(:,me) * zeta & - + source_dotState(1:size_so(so),1,so)* (1.0_pReal - zeta) - r(1:size_so(so)) = damageState(ph)%p(so)%state (1:size_so(so),me) & - - damageState(ph)%p(so)%subState0(1:size_so(so),me) & - - damageState(ph)%p(so)%dotState (1:size_so(so),me) * dt - damageState(ph)%p(so)%state(1:size_so(so),me) = damageState(ph)%p(so)%state(1:size_so(so),me) & - - r(1:size_so(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 + + zeta = damper(damageState(ph)%dotState(:,me),source_dotState(1:size_so,1),source_dotState(1:size_so,2)) + damageState(ph)%dotState(:,me) = damageState(ph)%dotState(:,me) * zeta & + + source_dotState(1:size_so,1)* (1.0_pReal - zeta) + r(1:size_so) = damageState(ph)%state (1:size_so,me) & + - damageState(ph)%subState0(1:size_so,me) & + - damageState(ph)%dotState (1:size_so,me) * dt + damageState(ph)%state(1:size_so,me) = damageState(ph)%state(1:size_so,me) - r(1:size_so) + converged_ = converged_ .and. converged(r(1:size_so), & + damageState(ph)%state(1:size_so,me), & + damageState(ph)%atol(1:size_so)) + if(converged_) then 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 - 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 @@ -443,7 +434,6 @@ function phase_damage_deltaState(Fe, ph, me) result(broken) real(pReal), intent(in), dimension(3,3) :: & Fe !< elastic deformation gradient integer :: & - so, & myOffset, & mySize logical :: & @@ -452,23 +442,22 @@ function phase_damage_deltaState(Fe, ph, me) result(broken) 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 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 - myOffset = damageState(ph)%p(so)%offsetDeltaState - mySize = damageState(ph)%p(so)%sizeDeltaState - 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) + myOffset = damageState(ph)%offsetDeltaState + mySize = damageState(ph)%sizeDeltaState + damageState(ph)%state(myOffset + 1: myOffset + mySize,me) = & + damageState(ph)%state(myOffset + 1: myOffset + mySize,me) + damageState(ph)%deltaState(1:mySize,me) endif end select sourceType - enddo SourceLoop end function phase_damage_deltaState diff --git a/src/phase_damage_anisobrittle.f90 b/src/phase_damage_anisobrittle.f90 index ad7f80604..422b2a91a 100644 --- a/src/phase_damage_anisobrittle.f90 +++ b/src/phase_damage_anisobrittle.f90 @@ -6,8 +6,7 @@ !-------------------------------------------------------------------------------------------------- submodule (phase:damagee) anisobrittle - integer, dimension(:), allocatable :: & - source_damage_anisoBrittle_offset, & !< which source is my current source mechanism? + integer, dimension(:), allocatable :: & source_damage_anisoBrittle_instance !< instance of source mechanism 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') allocate(param(Ninstances)) - allocate(source_damage_anisoBrittle_offset (phases%length), source=0) allocate(source_damage_anisoBrittle_instance(phases%length), source=0) do p = 1, phases%length @@ -68,7 +66,6 @@ module function anisobrittle_init(source_length) result(mySources) sources => phase%get('damage') do sourceOffset = 1, sources%length if(mySources(sourceOffset,p)) then - source_damage_anisoBrittle_offset(p) = sourceOffset associate(prm => param(source_damage_anisoBrittle_instance(p))) 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' Nconstituents = count(material_phaseAt==p) * discretization_nIPs - call phase_allocateState(damageState(p)%p(sourceOffset),Nconstituents,1,1,0) - damageState(p)%p(sourceOffset)%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' + call phase_allocateState(damageState(p),Nconstituents,1,1,0) + damageState(p)%atol = src%get_asFloat('anisobrittle_atol',defaultVal=1.0e-3_pReal) + if(any(damageState(p)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' anisobrittle_atol' end associate @@ -139,14 +136,13 @@ module subroutine anisobrittle_dotState(S, co, ip, el) real(pReal) :: & traction_d, traction_t, traction_n, traction_crit + 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(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 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 +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 - damageState(ph)%p(sourceOffset)%dotState(1,me) & - = damageState(ph)%p(sourceOffset)%dotState(1,me) & + damageState(ph)%dotState(1,me) & + = damageState(ph)%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 + & @@ -180,12 +176,8 @@ module subroutine source_damage_anisobrittle_getRateAndItsTangent(localphiDot, d localphiDot, & dLocalphiDot_dPhi - integer :: & - sourceOffset - sourceOffset = source_damage_anisoBrittle_offset(phase) - - dLocalphiDot_dPhi = -damageState(phase)%p(sourceOffset)%state(1,constituent) + dLocalphiDot_dPhi = -damageState(phase)%state(1,constituent) localphiDot = 1.0_pReal & + dLocalphiDot_dPhi*phi @@ -204,7 +196,7 @@ module subroutine anisobrittle_results(phase,group) integer :: o 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) select case(trim(prm%output(o))) case ('f_phi') diff --git a/src/phase_damage_anisoductile.f90 b/src/phase_damage_anisoductile.f90 index b22f5b3d6..bafcf80a2 100644 --- a/src/phase_damage_anisoductile.f90 +++ b/src/phase_damage_anisoductile.f90 @@ -7,7 +7,6 @@ submodule(phase:damagee) anisoductile integer, dimension(:), allocatable :: & - source_damage_anisoDuctile_offset, & !< which source is my current damage mechanism? source_damage_anisoDuctile_instance !< instance of damage source mechanism 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') allocate(param(Ninstances)) - allocate(source_damage_anisoDuctile_offset (phases%length), source=0) allocate(source_damage_anisoDuctile_instance(phases%length), source=0) do p = 1, phases%length @@ -65,7 +63,6 @@ module function anisoductile_init(source_length) result(mySources) sources => phase%get('source') do sourceOffset = 1, sources%length if(mySources(sourceOffset,p)) then - source_damage_anisoDuctile_offset(p) = sourceOffset associate(prm => param(source_damage_anisoDuctile_instance(p))) 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' Nconstituents=count(material_phaseAt==p) * discretization_nIPs - call phase_allocateState(damageState(p)%p(sourceOffset),Nconstituents,1,1,0) - damageState(p)%p(sourceOffset)%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' + call phase_allocateState(damageState(p),Nconstituents,1,1,0) + damageState(p)%atol = src%get_asFloat('anisoDuctile_atol',defaultVal=1.0e-3_pReal) + if(any(damageState(p)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' anisoductile_atol' end associate @@ -116,20 +113,14 @@ module subroutine anisoductile_dotState(co, ip, el) integer :: & ph, & - me, & - sourceOffset, & - damageOffset, & - homog + me 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(ph))) - damageState(ph)%p(sourceOffset)%dotState(1,me) & - = sum(plasticState(ph)%slipRate(:,me)/(phase_damage_get_phi(co,ip,el)**prm%q)/prm%gamma_crit) + damageState(ph)%dotState(1,me) = sum(plasticState(ph)%slipRate(:,me)/(phase_damage_get_phi(co,ip,el)**prm%q)/prm%gamma_crit) end associate end subroutine anisoductile_dotState @@ -149,12 +140,8 @@ module subroutine source_damage_anisoDuctile_getRateAndItsTangent(localphiDot, d localphiDot, & dLocalphiDot_dPhi - integer :: & - sourceOffset - sourceOffset = source_damage_anisoDuctile_offset(phase) - - dLocalphiDot_dPhi = -damageState(phase)%p(sourceOffset)%state(1,constituent) + dLocalphiDot_dPhi = -damageState(phase)%state(1,constituent) localphiDot = 1.0_pReal & + dLocalphiDot_dPhi*phi @@ -173,7 +160,7 @@ module subroutine anisoductile_results(phase,group) integer :: o 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) select case(trim(prm%output(o))) case ('f_phi') diff --git a/src/phase_damage_isobrittle.f90 b/src/phase_damage_isobrittle.f90 index a2ae7dfb7..0a87493a6 100644 --- a/src/phase_damage_isobrittle.f90 +++ b/src/phase_damage_isobrittle.f90 @@ -7,7 +7,6 @@ submodule(phase:damagee) isobrittle integer, dimension(:), allocatable :: & - source_damage_isoBrittle_offset, & source_damage_isoBrittle_instance 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') allocate(param(Ninstances)) - allocate(source_damage_isoBrittle_offset (phases%length), source=0) allocate(source_damage_isoBrittle_instance(phases%length), source=0) do p = 1, phases%length @@ -58,7 +56,6 @@ module function isobrittle_init(source_length) result(mySources) sources => phase%get('damage') do sourceOffset = 1, sources%length if(mySources(sourceOffset,p)) then - source_damage_isoBrittle_offset(p) = sourceOffset associate(prm => param(source_damage_isoBrittle_instance(p))) 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' Nconstituents = count(material_phaseAt==p) * discretization_nIPs - call phase_allocateState(damageState(p)%p(sourceOffset),Nconstituents,1,1,1) - damageState(p)%p(sourceOffset)%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' + call phase_allocateState(damageState(p),Nconstituents,1,1,1) + damageState(p)%atol = src%get_asFloat('isoBrittle_atol',defaultVal=1.0e-3_pReal) + if(any(damageState(p)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' isobrittle_atol' end associate @@ -102,14 +99,11 @@ module subroutine source_damage_isoBrittle_deltaState(C, Fe, ph,me) real(pReal), intent(in), dimension(6,6) :: & C - integer :: & - sourceOffset real(pReal), dimension(6) :: & strain real(pReal) :: & strainenergy - sourceOffset = source_damage_isoBrittle_offset(ph) 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 ! 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 - damageState(ph)%p(sourceOffset)%deltaState(1,me) = & - strainenergy - damageState(ph)%p(sourceOffset)%state(1,me) - else - damageState(ph)%p(sourceOffset)%deltaState(1,me) = & - damageState(ph)%p(sourceOffset)%subState0(1,me) - & - damageState(ph)%p(sourceOffset)%state(1,me) - endif + if (strainenergy > damageState(ph)%subState0(1,me)) then + damageState(ph)%deltaState(1,me) = strainenergy - damageState(ph)%state(1,me) + else + damageState(ph)%deltaState(1,me) = damageState(ph)%subState0(1,me) - damageState(ph)%state(1,me) + endif end associate end subroutine source_damage_isoBrittle_deltaState @@ -144,15 +135,11 @@ module subroutine source_damage_isoBrittle_getRateAndItsTangent(localphiDot, dLo localphiDot, & dLocalphiDot_dPhi - integer :: & - sourceOffset - - sourceOffset = source_damage_isoBrittle_offset(phase) associate(prm => param(source_damage_isoBrittle_instance(phase))) - localphiDot = 1.0_pReal & - - phi*damageState(phase)%p(sourceOffset)%state(1,constituent) - dLocalphiDot_dPhi = - damageState(phase)%p(sourceOffset)%state(1,constituent) + localphiDot = 1.0_pReal & + - phi*damageState(phase)%state(1,constituent) + dLocalphiDot_dPhi = - damageState(phase)%state(1,constituent) end associate end subroutine source_damage_isoBrittle_getRateAndItsTangent @@ -169,7 +156,7 @@ module subroutine isobrittle_results(phase,group) integer :: o 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) select case(trim(prm%output(o))) case ('f_phi') diff --git a/src/phase_damage_isoductile.f90 b/src/phase_damage_isoductile.f90 index f378b7ad2..657f9da39 100644 --- a/src/phase_damage_isoductile.f90 +++ b/src/phase_damage_isoductile.f90 @@ -7,7 +7,6 @@ submodule(phase:damagee) isoductile integer, dimension(:), allocatable :: & - source_damage_isoDuctile_offset, & !< which source is my current damage mechanism? source_damage_isoDuctile_instance !< instance of damage source mechanism 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') allocate(param(Ninstances)) - allocate(source_damage_isoDuctile_offset (phases%length), source=0) allocate(source_damage_isoDuctile_instance(phases%length), source=0) do p = 1, phases%length @@ -60,7 +58,6 @@ module function isoductile_init(source_length) result(mySources) sources => phase%get('damage') do sourceOffset = 1, sources%length if(mySources(sourceOffset,p)) then - source_damage_isoDuctile_offset(p) = sourceOffset associate(prm => param(source_damage_isoDuctile_instance(p))) 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' Nconstituents=count(material_phaseAt==p) * discretization_nIPs - call phase_allocateState(damageState(p)%p(sourceOffset),Nconstituents,1,1,0) - damageState(p)%p(sourceOffset)%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' + call phase_allocateState(damageState(p),Nconstituents,1,1,0) + damageState(p)%atol = src%get_asFloat('isoDuctile_atol',defaultVal=1.0e-3_pReal) + if(any(damageState(p)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' isoductile_atol' end associate @@ -107,20 +104,14 @@ module subroutine isoductile_dotState(co, ip, el) integer :: & ph, & - me, & - sourceOffset, & - damageOffset, & - homog + me 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(ph))) - damageState(ph)%p(sourceOffset)%dotState(1,me) = & - sum(plasticState(ph)%slipRate(:,me))/(phase_damage_get_phi(co,ip,el)**prm%q)/prm%gamma_crit + damageState(ph)%dotState(1,me) = sum(plasticState(ph)%slipRate(:,me))/(phase_damage_get_phi(co,ip,el)**prm%q)/prm%gamma_crit end associate end subroutine isoductile_dotState @@ -140,12 +131,8 @@ module subroutine source_damage_isoDuctile_getRateAndItsTangent(localphiDot, dLo localphiDot, & dLocalphiDot_dPhi - integer :: & - sourceOffset - sourceOffset = source_damage_isoDuctile_offset(phase) - - dLocalphiDot_dPhi = -damageState(phase)%p(sourceOffset)%state(1,constituent) + dLocalphiDot_dPhi = -damageState(phase)%state(1,constituent) localphiDot = 1.0_pReal & + dLocalphiDot_dPhi*phi @@ -164,7 +151,7 @@ module subroutine isoductile_results(phase,group) integer :: o 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) select case(trim(prm%output(o))) case ('f_phi') diff --git a/src/phase_mechanical.f90 b/src/phase_mechanical.f90 index f7631a38d..bc4b2cd47 100644 --- a/src/phase_mechanical.f90 +++ b/src/phase_mechanical.f90 @@ -1138,7 +1138,7 @@ module function crystallite_stress(dt,co,ip,el) result(converged_) real(pReal) :: & formerSubStep integer :: & - so, ph, me, sizeDotState + ph, me, sizeDotState logical :: todo real(pReal) :: subFrac,subStep 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) 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) subFi0 = phase_mechanical_Fi0(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) subFi0 = phase_mechanical_Fi(ph)%data(1:3,1:3,me) subState0 = plasticState(ph)%state(:,me) - do so = 1, phase_Nsources(ph) - damageState(ph)%p(so)%subState0(:,me) = damageState(ph)%p(so)%state(:,me) - enddo + if (phase_Nsources(ph) > 0) & + damageState(ph)%subState0(:,me) = damageState(ph)%state(:,me) + endif !-------------------------------------------------------------------------------------------------- ! 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 endif plasticState(ph)%state(:,me) = subState0 - do so = 1, phase_Nsources(ph) - damageState(ph)%p(so)%state(:,me) = damageState(ph)%p(so)%subState0(:,me) - enddo + if (phase_Nsources(ph) > 0) & + damageState(ph)%state(:,me) = damageState(ph)%subState0(:,me) todo = subStep > num%subStepMinCryst ! still on track or already done (beyond repair) endif From 775a51faa11ae9d75d85cdf30e583e8f5dcc4966 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 13 Feb 2021 10:11:39 +0100 Subject: [PATCH 13/25] explicit instance mapping not needed --- src/phase_damage_anisobrittle.f90 | 52 ++++++++++++++----------------- src/phase_damage_anisoductile.f90 | 29 +++++++---------- src/phase_damage_isobrittle.f90 | 32 +++++++++---------- src/phase_damage_isoductile.f90 | 28 +++++++---------- 4 files changed, 60 insertions(+), 81 deletions(-) diff --git a/src/phase_damage_anisobrittle.f90 b/src/phase_damage_anisobrittle.f90 index 422b2a91a..81072e16d 100644 --- a/src/phase_damage_anisobrittle.f90 +++ b/src/phase_damage_anisobrittle.f90 @@ -6,9 +6,6 @@ !-------------------------------------------------------------------------------------------------- submodule (phase:damagee) anisobrittle - integer, dimension(:), allocatable :: & - source_damage_anisoBrittle_instance !< instance of source mechanism - type :: tParameters !< container type for internal constitutive parameters real(pReal) :: & dot_o, & !< opening rate of cleavage planes @@ -56,17 +53,16 @@ module function anisobrittle_init(source_length) result(mySources) if(Ninstances == 0) return phases => config_material%get('phase') - allocate(param(Ninstances)) - allocate(source_damage_anisoBrittle_instance(phases%length), source=0) + allocate(param(phases%length)) + do p = 1, phases%length 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('damage') do sourceOffset = 1, sources%length if(mySources(sourceOffset,p)) then - associate(prm => param(source_damage_anisoBrittle_instance(p))) + associate(prm => param(p)) src => sources%get(sourceOffset) N_cl = src%get_asInts('N_cl',defaultVal=emptyIntArray) @@ -141,22 +137,21 @@ module subroutine anisobrittle_dotState(S, co, ip, el) me = material_phasememberAt(co,ip,el) - associate(prm => param(source_damage_anisoBrittle_instance(ph))) - damageState(ph)%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)) - traction_n = math_tensordot(S,prm%cleavage_systems(1:3,1:3,3,i)) + associate(prm => param(ph)) + damageState(ph)%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)) + traction_n = math_tensordot(S,prm%cleavage_systems(1:3,1:3,3,i)) - 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)%dotState(1,me) & - = damageState(ph)%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 + & - (max(0.0_pReal, abs(traction_n) - traction_crit)/traction_crit)**prm%q) - enddo + damageState(ph)%dotState(1,me) = damageState(ph)%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 + & + (max(0.0_pReal, abs(traction_n) - traction_crit)/traction_crit)**prm%q) + enddo end associate end subroutine anisobrittle_dotState @@ -195,14 +190,13 @@ module subroutine anisobrittle_results(phase,group) integer :: o - associate(prm => param(source_damage_anisoBrittle_instance(phase)), & - stt => damageState(phase)%state) - outputsLoop: do o = 1,size(prm%output) - select case(trim(prm%output(o))) - case ('f_phi') - call results_writeDataset(group,stt,trim(prm%output(o)),'driving force','J/m³') - end select - enddo outputsLoop + associate(prm => param(phase), stt => damageState(phase)%state) + outputsLoop: do o = 1,size(prm%output) + select case(trim(prm%output(o))) + case ('f_phi') + call results_writeDataset(group,stt,trim(prm%output(o)),'driving force','J/m³') + end select + enddo outputsLoop end associate end subroutine anisobrittle_results diff --git a/src/phase_damage_anisoductile.f90 b/src/phase_damage_anisoductile.f90 index bafcf80a2..4d2392589 100644 --- a/src/phase_damage_anisoductile.f90 +++ b/src/phase_damage_anisoductile.f90 @@ -6,9 +6,6 @@ !-------------------------------------------------------------------------------------------------- submodule(phase:damagee) anisoductile - integer, dimension(:), allocatable :: & - source_damage_anisoDuctile_instance !< instance of damage source mechanism - type :: tParameters !< container type for internal constitutive parameters real(pReal) :: & q !< damage rate sensitivity @@ -18,7 +15,7 @@ submodule(phase:damagee) anisoductile output end type tParameters - type(tParameters), dimension(:), allocatable :: param !< containers of constitutive parameters (len Ninstances) + type(tParameters), dimension(:), allocatable :: param !< containers of constitutive parameters contains @@ -51,19 +48,17 @@ module function anisoductile_init(source_length) result(mySources) if(Ninstances == 0) return phases => config_material%get('phase') - allocate(param(Ninstances)) - allocate(source_damage_anisoDuctile_instance(phases%length), source=0) + allocate(param(phases%length)) do p = 1, phases%length phase => phases%get(p) - if(any(mySources(:,p))) source_damage_anisoDuctile_instance(p) = count(mySources(:,1:p)) if(count(mySources(:,p)) == 0) cycle mech => phase%get('mechanics') pl => mech%get('plasticity') sources => phase%get('source') do sourceOffset = 1, sources%length if(mySources(sourceOffset,p)) then - associate(prm => param(source_damage_anisoDuctile_instance(p))) + associate(prm => param(p)) src => sources%get(sourceOffset) N_sl = pl%get_asInts('N_sl',defaultVal=emptyIntArray) @@ -119,8 +114,8 @@ module subroutine anisoductile_dotState(co, ip, el) me = material_phasememberAt(co,ip,el) - associate(prm => param(source_damage_anisoDuctile_instance(ph))) - damageState(ph)%dotState(1,me) = sum(plasticState(ph)%slipRate(:,me)/(phase_damage_get_phi(co,ip,el)**prm%q)/prm%gamma_crit) + associate(prm => param(ph)) + damageState(ph)%dotState(1,me) = sum(plasticState(ph)%slipRate(:,me)/(damage_phi(ph,me)**prm%q)/prm%gamma_crit) end associate end subroutine anisoductile_dotState @@ -159,14 +154,14 @@ module subroutine anisoductile_results(phase,group) integer :: o - associate(prm => param(source_damage_anisoDuctile_instance(phase)), & + associate(prm => param(phase), & stt => damageState(phase)%state) - outputsLoop: do o = 1,size(prm%output) - select case(trim(prm%output(o))) - case ('f_phi') - call results_writeDataset(group,stt,trim(prm%output(o)),'driving force','J/m³') - end select - enddo outputsLoop + outputsLoop: do o = 1,size(prm%output) + select case(trim(prm%output(o))) + case ('f_phi') + call results_writeDataset(group,stt,trim(prm%output(o)),'driving force','J/m³') + end select + enddo outputsLoop end associate end subroutine anisoductile_results diff --git a/src/phase_damage_isobrittle.f90 b/src/phase_damage_isobrittle.f90 index 0a87493a6..4e37f3ee9 100644 --- a/src/phase_damage_isobrittle.f90 +++ b/src/phase_damage_isobrittle.f90 @@ -6,9 +6,6 @@ !-------------------------------------------------------------------------------------------------- submodule(phase:damagee) isobrittle - integer, dimension(:), allocatable :: & - source_damage_isoBrittle_instance - type :: tParameters !< container type for internal constitutive parameters real(pReal) :: & W_crit !< critical elastic strain energy @@ -46,17 +43,15 @@ module function isobrittle_init(source_length) result(mySources) if(Ninstances == 0) return phases => config_material%get('phase') - allocate(param(Ninstances)) - allocate(source_damage_isoBrittle_instance(phases%length), source=0) + allocate(param(phases%length)) do p = 1, phases%length 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('damage') do sourceOffset = 1, sources%length if(mySources(sourceOffset,p)) then - associate(prm => param(source_damage_isoBrittle_instance(p))) + associate(prm => param(p)) src => sources%get(sourceOffset) prm%W_crit = src%get_asFloat('W_crit') @@ -107,9 +102,9 @@ module subroutine source_damage_isoBrittle_deltaState(C, Fe, ph,me) strain = 0.5_pReal*math_sym33to6(matmul(transpose(Fe),Fe)-math_I3) - 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 + associate(prm => param(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(ph)%subState0(1,me)) then damageState(ph)%deltaState(1,me) = strainenergy - damageState(ph)%state(1,me) @@ -136,7 +131,7 @@ module subroutine source_damage_isoBrittle_getRateAndItsTangent(localphiDot, dLo dLocalphiDot_dPhi - associate(prm => param(source_damage_isoBrittle_instance(phase))) + associate(prm => param(phase)) localphiDot = 1.0_pReal & - phi*damageState(phase)%state(1,constituent) dLocalphiDot_dPhi = - damageState(phase)%state(1,constituent) @@ -155,14 +150,15 @@ module subroutine isobrittle_results(phase,group) integer :: o - associate(prm => param(source_damage_isoBrittle_instance(phase)), & + + associate(prm => param(phase), & stt => damageState(phase)%state) - outputsLoop: do o = 1,size(prm%output) - select case(trim(prm%output(o))) - case ('f_phi') - call results_writeDataset(group,stt,trim(prm%output(o)),'driving force','J/m³') - end select - enddo outputsLoop + outputsLoop: do o = 1,size(prm%output) + select case(trim(prm%output(o))) + case ('f_phi') + call results_writeDataset(group,stt,trim(prm%output(o)),'driving force','J/m³') + end select + enddo outputsLoop end associate end subroutine isobrittle_results diff --git a/src/phase_damage_isoductile.f90 b/src/phase_damage_isoductile.f90 index 657f9da39..89211b746 100644 --- a/src/phase_damage_isoductile.f90 +++ b/src/phase_damage_isoductile.f90 @@ -6,9 +6,6 @@ !-------------------------------------------------------------------------------------------------- submodule(phase:damagee) isoductile - integer, dimension(:), allocatable :: & - source_damage_isoDuctile_instance !< instance of damage source mechanism - type:: tParameters !< container type for internal constitutive parameters real(pReal) :: & gamma_crit, & !< critical plastic strain @@ -48,17 +45,15 @@ module function isoductile_init(source_length) result(mySources) if(Ninstances == 0) return phases => config_material%get('phase') - allocate(param(Ninstances)) - allocate(source_damage_isoDuctile_instance(phases%length), source=0) + allocate(param(phases%length)) do p = 1, phases%length 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('damage') do sourceOffset = 1, sources%length if(mySources(sourceOffset,p)) then - associate(prm => param(source_damage_isoDuctile_instance(p))) + associate(prm => param(p)) src => sources%get(sourceOffset) prm%q = src%get_asFloat('q') @@ -110,8 +105,8 @@ module subroutine isoductile_dotState(co, ip, el) me = material_phasememberAt(co,ip,el) - associate(prm => param(source_damage_isoDuctile_instance(ph))) - damageState(ph)%dotState(1,me) = sum(plasticState(ph)%slipRate(:,me))/(phase_damage_get_phi(co,ip,el)**prm%q)/prm%gamma_crit + associate(prm => param(ph)) + damageState(ph)%dotState(1,me) = sum(plasticState(ph)%slipRate(:,me))/(damage_phi(ph,me)**prm%q)/prm%gamma_crit end associate end subroutine isoductile_dotState @@ -150,14 +145,13 @@ module subroutine isoductile_results(phase,group) integer :: o - associate(prm => param(source_damage_isoDuctile_instance(phase)), & - stt => damageState(phase)%state) - outputsLoop: do o = 1,size(prm%output) - select case(trim(prm%output(o))) - case ('f_phi') - call results_writeDataset(group,stt,trim(prm%output(o)),'driving force','J/m³') - end select - enddo outputsLoop + associate(prm => param(phase), stt => damageState(phase)%state) + outputsLoop: do o = 1,size(prm%output) + select case(trim(prm%output(o))) + case ('f_phi') + call results_writeDataset(group,stt,trim(prm%output(o)),'driving force','J/m³') + end select + enddo outputsLoop end associate end subroutine isoductile_results From d9699b0f2e2c5126871f62b4ad8ab3183fd14b1d Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 13 Feb 2021 11:01:08 +0100 Subject: [PATCH 14/25] simplified access pattern --- src/phase_damage.f90 | 92 ++++++++++++------------------- src/phase_damage_anisobrittle.f90 | 25 +++------ src/phase_damage_anisoductile.f90 | 24 +++----- src/phase_damage_isobrittle.f90 | 25 ++++----- src/phase_damage_isoductile.f90 | 23 +++----- 5 files changed, 71 insertions(+), 118 deletions(-) diff --git a/src/phase_damage.f90 b/src/phase_damage.f90 index 5f323c3a2..803a9251d 100644 --- a/src/phase_damage.f90 +++ b/src/phase_damage.f90 @@ -43,82 +43,65 @@ submodule(phase) damagee end function isoductile_init - module subroutine source_damage_isoBrittle_deltaState(C, Fe, ph, me) + module subroutine 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 + end subroutine isobrittle_deltaState - module subroutine anisobrittle_dotState(S, co, ip, el) - integer, intent(in) :: & - co, & !< component-ID of integration point - ip, & !< integration point - el !< element + module subroutine anisobrittle_dotState(S, ph, me) + integer, intent(in) :: ph,me real(pReal), intent(in), dimension(3,3) :: & S end subroutine anisobrittle_dotState - module subroutine anisoductile_dotState(co, ip, el) - integer, intent(in) :: & - co, & !< component-ID of integration point - ip, & !< integration point - el !< element + module subroutine anisoductile_dotState(ph,me) + integer, intent(in) :: ph,me end subroutine anisoductile_dotState - module subroutine isoductile_dotState(co, ip, el) - integer, intent(in) :: & - co, & !< component-ID of integration point - ip, & !< integration point - el !< element + module subroutine isoductile_dotState(ph,me) + integer, intent(in) :: ph,me end subroutine isoductile_dotState - module subroutine source_damage_anisobrittle_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent) - integer, intent(in) :: & - phase, & !< phase ID of element - constituent !< position of element within its phase instance + module subroutine anisobrittle_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, ph, me) + integer, intent(in) :: ph,me real(pReal), intent(in) :: & phi !< damage parameter real(pReal), intent(out) :: & localphiDot, & dLocalphiDot_dPhi - end subroutine source_damage_anisoBrittle_getRateAndItsTangent + end subroutine anisobrittle_getRateAndItsTangent - module subroutine source_damage_anisoDuctile_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent) - integer, intent(in) :: & - phase, & !< phase ID of element - constituent !< position of element within its phase instance + module subroutine anisoductile_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, ph,me) + integer, intent(in) :: ph,me real(pReal), intent(in) :: & phi !< damage parameter real(pReal), intent(out) :: & localphiDot, & dLocalphiDot_dPhi - end subroutine source_damage_anisoDuctile_getRateAndItsTangent + end subroutine anisoductile_getRateAndItsTangent - module subroutine source_damage_isoBrittle_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent) - integer, intent(in) :: & - phase, & !< phase ID of element - constituent !< position of element within its phase instance + module subroutine isobrittle_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, ph,me) + integer, intent(in) :: ph,me real(pReal), intent(in) :: & phi !< damage parameter real(pReal), intent(out) :: & localphiDot, & dLocalphiDot_dPhi - end subroutine source_damage_isoBrittle_getRateAndItsTangent + end subroutine isobrittle_getRateAndItsTangent - module subroutine source_damage_isoDuctile_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent) - integer, intent(in) :: & - phase, & !< phase ID of element - constituent !< position of element within its phase instance + module subroutine isoductile_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, ph,me) + integer, intent(in) :: ph,me real(pReal), intent(in) :: & phi !< damage parameter real(pReal), intent(out) :: & localphiDot, & dLocalphiDot_dPhi - end subroutine source_damage_isoDuctile_getRateAndItsTangent + end subroutine isoductile_getRateAndItsTangent module subroutine anisobrittle_results(phase,group) integer, intent(in) :: phase @@ -225,16 +208,16 @@ module subroutine phase_damage_getRateAndItsTangents(phiDot, dPhiDot_dPhi, phi, do so = 1, phase_Nsources(ph) select case(phase_source(so,ph)) case (DAMAGE_ISOBRITTLE_ID) - call source_damage_isobrittle_getRateAndItsTangent (localphiDot, dLocalphiDot_dPhi, phi, ph, me) + call isobrittle_getRateAndItsTangent (localphiDot, dLocalphiDot_dPhi, phi, ph, me) case (DAMAGE_ISODUCTILE_ID) - call source_damage_isoductile_getRateAndItsTangent (localphiDot, dLocalphiDot_dPhi, phi, ph, me) + call isoductile_getRateAndItsTangent (localphiDot, dLocalphiDot_dPhi, phi, ph, me) case (DAMAGE_ANISOBRITTLE_ID) - call source_damage_anisobrittle_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, ph, me) + call anisobrittle_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, ph, me) case (DAMAGE_ANISODUCTILE_ID) - call source_damage_anisoductile_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, ph, me) + call anisoductile_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, ph, me) case default localphiDot = 0.0_pReal @@ -281,7 +264,7 @@ module function integrateDamageState(dt,co,ip,el) result(broken) me = material_phaseMemberAt(co,ip,el) converged_ = .true. - broken = phase_damage_collectDotState(co,ip,el,ph,me) + broken = phase_damage_collectDotState(ph,me) if(broken) return size_so = damageState(ph)%sizeDotState @@ -294,7 +277,7 @@ module function integrateDamageState(dt,co,ip,el) result(broken) if(nIterationState > 1) source_dotState(1:size_so,2) = source_dotState(1:size_so,1) source_dotState(1:size_so,1) = damageState(ph)%dotState(:,me) - broken = phase_damage_collectDotState(co,ip,el,ph,me) + broken = phase_damage_collectDotState(ph,me) if(broken) exit iteration @@ -384,39 +367,34 @@ end subroutine damage_results !-------------------------------------------------------------------------------------------------- !> @brief contains the constitutive equation for calculating the rate of change of microstructure !-------------------------------------------------------------------------------------------------- -function phase_damage_collectDotState(co,ip,el,ph,me) result(broken) +function phase_damage_collectDotState(ph,me) result(broken) integer, intent(in) :: & - co, & !< component-ID me integration point - ip, & !< integration point - el, & !< element ph, & - me - integer :: & - so !< counter in source loop + me !< counter in source loop logical :: broken broken = .false. - SourceLoop: do so = 1, phase_Nsources(ph) + if (phase_Nsources(ph)==1) then - sourceType: select case (phase_source(so,ph)) + sourceType: select case (phase_source(1,ph)) case (DAMAGE_ISODUCTILE_ID) sourceType - call isoductile_dotState(co, ip, el) + call isoductile_dotState(ph,me) case (DAMAGE_ANISODUCTILE_ID) sourceType - call anisoductile_dotState(co, ip, el) + call anisoductile_dotState(ph,me) case (DAMAGE_ANISOBRITTLE_ID) sourceType - call anisobrittle_dotState(mechanical_S(ph,me),co, ip, el) ! correct stress? + call anisobrittle_dotState(mechanical_S(ph,me), ph,me) ! correct stress? end select sourceType broken = broken .or. any(IEEE_is_NaN(damageState(ph)%dotState(:,me))) - enddo SourceLoop + endif end function phase_damage_collectDotState @@ -447,7 +425,7 @@ function phase_damage_deltaState(Fe, ph, me) result(broken) sourceType: select case (phase_source(1,ph)) case (DAMAGE_ISOBRITTLE_ID) sourceType - call source_damage_isoBrittle_deltaState(phase_homogenizedC(ph,me), Fe, ph,me) + call isobrittle_deltaState(phase_homogenizedC(ph,me), Fe, ph,me) broken = any(IEEE_is_NaN(damageState(ph)%deltaState(:,me))) if(.not. broken) then myOffset = damageState(ph)%offsetDeltaState diff --git a/src/phase_damage_anisobrittle.f90 b/src/phase_damage_anisobrittle.f90 index 81072e16d..a121ee932 100644 --- a/src/phase_damage_anisobrittle.f90 +++ b/src/phase_damage_anisobrittle.f90 @@ -113,18 +113,14 @@ end function anisobrittle_init !-------------------------------------------------------------------------------------------------- !> @brief calculates derived quantities from state !-------------------------------------------------------------------------------------------------- -module subroutine anisobrittle_dotState(S, co, ip, el) +module subroutine anisobrittle_dotState(S, ph,me) integer, intent(in) :: & - co, & !< component-ID of integration point - ip, & !< integration point - el !< element + ph,me real(pReal), intent(in), dimension(3,3) :: & S integer :: & - ph, & - me, & sourceOffset, & damageOffset, & homog, & @@ -133,10 +129,6 @@ module subroutine anisobrittle_dotState(S, co, ip, el) traction_d, traction_t, traction_n, traction_crit - ph = material_phaseAt(co,el) - me = material_phasememberAt(co,ip,el) - - associate(prm => param(ph)) damageState(ph)%dotState(1,me) = 0.0_pReal do i = 1, prm%sum_N_cl @@ -144,7 +136,7 @@ module subroutine anisobrittle_dotState(S, co, ip, el) traction_t = math_tensordot(S,prm%cleavage_systems(1:3,1:3,2,i)) traction_n = math_tensordot(S,prm%cleavage_systems(1:3,1:3,3,i)) - traction_crit = prm%g_crit(i)*phase_damage_get_phi(co,ip,el)**2.0_pReal + traction_crit = prm%g_crit(i)*damage_phi(ph,me)**2.0_pReal damageState(ph)%dotState(1,me) = damageState(ph)%dotState(1,me) & + prm%dot_o / prm%s_crit(i) & @@ -160,11 +152,11 @@ end subroutine anisobrittle_dotState !-------------------------------------------------------------------------------------------------- !> @brief returns local part of nonlocal damage driving force !-------------------------------------------------------------------------------------------------- -module subroutine source_damage_anisobrittle_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent) +module subroutine anisobrittle_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, ph, me) integer, intent(in) :: & - phase, & - constituent + ph, & + me real(pReal), intent(in) :: & phi real(pReal), intent(out) :: & @@ -172,12 +164,12 @@ module subroutine source_damage_anisobrittle_getRateAndItsTangent(localphiDot, d dLocalphiDot_dPhi - dLocalphiDot_dPhi = -damageState(phase)%state(1,constituent) + dLocalphiDot_dPhi = -damageState(ph)%state(1,me) localphiDot = 1.0_pReal & + dLocalphiDot_dPhi*phi -end subroutine source_damage_anisoBrittle_getRateAndItsTangent +end subroutine anisobrittle_getRateAndItsTangent !-------------------------------------------------------------------------------------------------- @@ -190,6 +182,7 @@ module subroutine anisobrittle_results(phase,group) integer :: o + associate(prm => param(phase), stt => damageState(phase)%state) outputsLoop: do o = 1,size(prm%output) select case(trim(prm%output(o))) diff --git a/src/phase_damage_anisoductile.f90 b/src/phase_damage_anisoductile.f90 index 4d2392589..cf55e398b 100644 --- a/src/phase_damage_anisoductile.f90 +++ b/src/phase_damage_anisoductile.f90 @@ -99,20 +99,12 @@ end function anisoductile_init !-------------------------------------------------------------------------------------------------- !> @brief calculates derived quantities from state !-------------------------------------------------------------------------------------------------- -module subroutine anisoductile_dotState(co, ip, el) +module subroutine anisoductile_dotState(ph,me) integer, intent(in) :: & - co, & !< component-ID of integration point - ip, & !< integration point - el !< element - - integer :: & ph, & me - ph = material_phaseAt(co,el) - me = material_phasememberAt(co,ip,el) - associate(prm => param(ph)) damageState(ph)%dotState(1,me) = sum(plasticState(ph)%slipRate(:,me)/(damage_phi(ph,me)**prm%q)/prm%gamma_crit) @@ -124,11 +116,11 @@ end subroutine anisoductile_dotState !-------------------------------------------------------------------------------------------------- !> @brief returns local part of nonlocal damage driving force !-------------------------------------------------------------------------------------------------- -module subroutine source_damage_anisoDuctile_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent) +module subroutine anisoductile_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, ph,me) integer, intent(in) :: & - phase, & - constituent + ph, & + me real(pReal), intent(in) :: & phi real(pReal), intent(out) :: & @@ -136,12 +128,12 @@ module subroutine source_damage_anisoDuctile_getRateAndItsTangent(localphiDot, d dLocalphiDot_dPhi - dLocalphiDot_dPhi = -damageState(phase)%state(1,constituent) + dLocalphiDot_dPhi = -damageState(ph)%state(1,me) localphiDot = 1.0_pReal & + dLocalphiDot_dPhi*phi -end subroutine source_damage_anisoDuctile_getRateAndItsTangent +end subroutine anisoductile_getRateAndItsTangent !-------------------------------------------------------------------------------------------------- @@ -154,8 +146,8 @@ module subroutine anisoductile_results(phase,group) integer :: o - associate(prm => param(phase), & - stt => damageState(phase)%state) + + associate(prm => param(phase), stt => damageState(phase)%state) outputsLoop: do o = 1,size(prm%output) select case(trim(prm%output(o))) case ('f_phi') diff --git a/src/phase_damage_isobrittle.f90 b/src/phase_damage_isobrittle.f90 index 4e37f3ee9..e267ffcc9 100644 --- a/src/phase_damage_isobrittle.f90 +++ b/src/phase_damage_isobrittle.f90 @@ -86,7 +86,7 @@ end function isobrittle_init !-------------------------------------------------------------------------------------------------- !> @brief calculates derived quantities from state !-------------------------------------------------------------------------------------------------- -module subroutine source_damage_isoBrittle_deltaState(C, Fe, ph,me) +module subroutine isobrittle_deltaState(C, Fe, ph,me) integer, intent(in) :: ph,me real(pReal), intent(in), dimension(3,3) :: & @@ -106,24 +106,21 @@ module subroutine source_damage_isoBrittle_deltaState(C, Fe, ph,me) 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(ph)%subState0(1,me)) then - damageState(ph)%deltaState(1,me) = strainenergy - damageState(ph)%state(1,me) - else - damageState(ph)%deltaState(1,me) = damageState(ph)%subState0(1,me) - damageState(ph)%state(1,me) - endif + damageState(ph)%deltaState(1,me) = merge(strainenergy - damageState(ph)%state(1,me), & + damageState(ph)%subState0(1,me) - damageState(ph)%state(1,me), & + strainenergy > damageState(ph)%subState0(1,me)) end associate -end subroutine source_damage_isoBrittle_deltaState +end subroutine isobrittle_deltaState !-------------------------------------------------------------------------------------------------- !> @brief returns local part of nonlocal damage driving force !-------------------------------------------------------------------------------------------------- -module subroutine source_damage_isoBrittle_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent) +module subroutine isobrittle_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, ph, me) integer, intent(in) :: & - phase, & - constituent + ph, me real(pReal), intent(in) :: & phi real(pReal), intent(out) :: & @@ -131,13 +128,13 @@ module subroutine source_damage_isoBrittle_getRateAndItsTangent(localphiDot, dLo dLocalphiDot_dPhi - associate(prm => param(phase)) + associate(prm => param(ph)) localphiDot = 1.0_pReal & - - phi*damageState(phase)%state(1,constituent) - dLocalphiDot_dPhi = - damageState(phase)%state(1,constituent) + - phi*damageState(ph)%state(1,me) + dLocalphiDot_dPhi = - damageState(ph)%state(1,me) end associate -end subroutine source_damage_isoBrittle_getRateAndItsTangent +end subroutine isobrittle_getRateAndItsTangent !-------------------------------------------------------------------------------------------------- diff --git a/src/phase_damage_isoductile.f90 b/src/phase_damage_isoductile.f90 index 89211b746..f59952222 100644 --- a/src/phase_damage_isoductile.f90 +++ b/src/phase_damage_isoductile.f90 @@ -90,23 +90,16 @@ end function isoductile_init !-------------------------------------------------------------------------------------------------- !> @brief calculates derived quantities from state !-------------------------------------------------------------------------------------------------- -module subroutine isoductile_dotState(co, ip, el) +module subroutine isoductile_dotState(ph, me) integer, intent(in) :: & - co, & !< component-ID of integration point - ip, & !< integration point - el !< element - - integer :: & ph, & me - ph = material_phaseAt(co,el) - me = material_phasememberAt(co,ip,el) - associate(prm => param(ph)) - damageState(ph)%dotState(1,me) = sum(plasticState(ph)%slipRate(:,me))/(damage_phi(ph,me)**prm%q)/prm%gamma_crit + damageState(ph)%dotState(1,me) = sum(plasticState(ph)%slipRate(:,me)) & + / (prm%gamma_crit*damage_phi(ph,me)**prm%q) end associate end subroutine isoductile_dotState @@ -115,11 +108,11 @@ end subroutine isoductile_dotState !-------------------------------------------------------------------------------------------------- !> @brief returns local part of nonlocal damage driving force !-------------------------------------------------------------------------------------------------- -module subroutine source_damage_isoDuctile_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent) +module subroutine isoductile_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, ph, me) integer, intent(in) :: & - phase, & - constituent + ph, & + me real(pReal), intent(in) :: & phi real(pReal), intent(out) :: & @@ -127,12 +120,12 @@ module subroutine source_damage_isoDuctile_getRateAndItsTangent(localphiDot, dLo dLocalphiDot_dPhi - dLocalphiDot_dPhi = -damageState(phase)%state(1,constituent) + dLocalphiDot_dPhi = -damageState(ph)%state(1,me) localphiDot = 1.0_pReal & + dLocalphiDot_dPhi*phi -end subroutine source_damage_isoDuctile_getRateAndItsTangent +end subroutine isoductile_getRateAndItsTangent !-------------------------------------------------------------------------------------------------- From 570086c814c22062bfd7b9582b64d805c6ae970b Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 13 Feb 2021 12:44:47 +0100 Subject: [PATCH 15/25] hard code at max 1 damage mechanism --- src/homogenization.f90 | 4 +- src/phase_damage.f90 | 61 ++++++++++++++----------------- src/phase_damage_anisobrittle.f90 | 18 ++++----- src/phase_damage_anisoductile.f90 | 29 +++++++-------- src/phase_damage_isobrittle.f90 | 20 +++++----- src/phase_damage_isoductile.f90 | 23 +++++------- 6 files changed, 70 insertions(+), 85 deletions(-) diff --git a/src/homogenization.f90 b/src/homogenization.f90 index 073482c9a..3e1b939b5 100644 --- a/src/homogenization.f90 +++ b/src/homogenization.f90 @@ -257,8 +257,8 @@ subroutine homogenization_init() allocate(homogState (size(material_name_homogenization))) allocate(damageState_h (size(material_name_homogenization))) - call material_parseHomogenization - print*, 'Homogenization parsed' + call material_parseHomogenization + num_homog => config_numerics%get('homogenization',defaultVal=emptyDict) num_homogGeneric => num_homog%get('generic',defaultVal=emptyDict) diff --git a/src/phase_damage.f90 b/src/phase_damage.f90 index 803a9251d..ca08d8aea 100644 --- a/src/phase_damage.f90 +++ b/src/phase_damage.f90 @@ -15,31 +15,27 @@ submodule(phase) damagee real(pReal), dimension(:), allocatable :: phi, d_phi_d_dot_phi end type tDataContainer - integer(kind(DAMAGE_UNDEFINED_ID)), dimension(:,:), allocatable :: & + integer(kind(DAMAGE_UNDEFINED_ID)), dimension(:), allocatable :: & phase_source !< active sources mechanisms of each phase type(tDataContainer), dimension(:), allocatable :: current interface - module function anisobrittle_init(source_length) result(mySources) - integer, intent(in) :: source_length - logical, dimension(:,:), allocatable :: mySources + module function anisobrittle_init() result(mySources) + logical, dimension(:), allocatable :: mySources end function anisobrittle_init - module function anisoductile_init(source_length) result(mySources) - integer, intent(in) :: source_length - logical, dimension(:,:), allocatable :: mySources + module function anisoductile_init() result(mySources) + logical, dimension(:), allocatable :: mySources end function anisoductile_init - module function isobrittle_init(source_length) result(mySources) - integer, intent(in) :: source_length - logical, dimension(:,:), allocatable :: mySources + module function isobrittle_init() result(mySources) + logical, dimension(:), allocatable :: mySources end function isobrittle_init - module function isoductile_init(source_length) result(mySources) - integer, intent(in) :: source_length - logical, dimension(:,:), allocatable :: mySources + module function isoductile_init() result(mySources) + logical, dimension(:), allocatable :: mySources end function isoductile_init @@ -163,14 +159,14 @@ module subroutine damage_init enddo - allocate(phase_source(maxval(phase_Nsources),phases%length), source = DAMAGE_UNDEFINED_ID) + allocate(phase_source(phases%length), source = DAMAGE_UNDEFINED_ID) ! initialize source mechanisms if(maxval(phase_Nsources) /= 0) then - where(isobrittle_init (maxval(phase_Nsources))) phase_source = DAMAGE_ISOBRITTLE_ID - where(isoductile_init (maxval(phase_Nsources))) phase_source = DAMAGE_ISODUCTILE_ID - where(anisobrittle_init (maxval(phase_Nsources))) phase_source = DAMAGE_ANISOBRITTLE_ID - where(anisoductile_init (maxval(phase_Nsources))) phase_source = DAMAGE_ANISODUCTILE_ID + where(isobrittle_init() ) phase_source = DAMAGE_ISOBRITTLE_ID + where(isoductile_init() ) phase_source = DAMAGE_ISODUCTILE_ID + where(anisobrittle_init()) phase_source = DAMAGE_ANISOBRITTLE_ID + where(anisoductile_init()) phase_source = DAMAGE_ANISODUCTILE_ID endif end subroutine damage_init @@ -206,7 +202,7 @@ module subroutine phase_damage_getRateAndItsTangents(phiDot, dPhiDot_dPhi, phi, ph = material_phaseAt(co,el) me = material_phasememberAt(co,ip,el) do so = 1, phase_Nsources(ph) - select case(phase_source(so,ph)) + select case(phase_source(ph)) case (DAMAGE_ISOBRITTLE_ID) call isobrittle_getRateAndItsTangent (localphiDot, dLocalphiDot_dPhi, phi, ph, me) @@ -340,10 +336,10 @@ module subroutine damage_results(group,ph) sourceLoop: do so = 1, phase_Nsources(ph) - if (phase_source(so,ph) /= DAMAGE_UNDEFINED_ID) & + if (phase_source(ph) /= DAMAGE_UNDEFINED_ID) & call results_closeGroup(results_addGroup(group//'sources/')) ! should be 'damage' - sourceType: select case (phase_source(so,ph)) + sourceType: select case (phase_source(ph)) case (DAMAGE_ISOBRITTLE_ID) sourceType call isobrittle_results(ph,group//'sources/') @@ -379,7 +375,7 @@ function phase_damage_collectDotState(ph,me) result(broken) if (phase_Nsources(ph)==1) then - sourceType: select case (phase_source(1,ph)) + sourceType: select case (phase_source(ph)) case (DAMAGE_ISODUCTILE_ID) sourceType call isoductile_dotState(ph,me) @@ -422,7 +418,7 @@ function phase_damage_deltaState(Fe, ph, me) result(broken) if (phase_Nsources(ph) == 0) return - sourceType: select case (phase_source(1,ph)) + sourceType: select case (phase_source(ph)) case (DAMAGE_ISOBRITTLE_ID) sourceType call isobrittle_deltaState(phase_homogenizedC(ph,me), Fe, ph,me) @@ -443,28 +439,25 @@ end function phase_damage_deltaState !-------------------------------------------------------------------------------------------------- !> @brief checks if a source mechanism is active or not !-------------------------------------------------------------------------------------------------- -function source_active(source_label,src_length) result(active_source) +function source_active(source_label) result(active_source) character(len=*), intent(in) :: source_label !< name of source mechanism - integer, intent(in) :: src_length !< max. number of sources in system - logical, dimension(:,:), allocatable :: active_source + logical, dimension(:), allocatable :: active_source class(tNode), pointer :: & phases, & phase, & sources, & src - integer :: p,s + integer :: ph phases => config_material%get('phase') - allocate(active_source(src_length,phases%length), source = .false. ) - do p = 1, phases%length - phase => phases%get(p) + allocate(active_source(phases%length)) + do ph = 1, phases%length + phase => phases%get(ph) sources => phase%get('damage',defaultVal=emptyList) - do s = 1, sources%length - src => sources%get(s) - if(src%get_asString('type') == source_label) active_source(s,p) = .true. - enddo + src => sources%get(1) + active_source(ph) = src%get_asString('type') == source_label enddo diff --git a/src/phase_damage_anisobrittle.f90 b/src/phase_damage_anisobrittle.f90 index a121ee932..07dd42945 100644 --- a/src/phase_damage_anisobrittle.f90 +++ b/src/phase_damage_anisobrittle.f90 @@ -31,23 +31,22 @@ contains !> @brief module initialization !> @details reads in material parameters, allocates arrays, and does sanity checks !-------------------------------------------------------------------------------------------------- -module function anisobrittle_init(source_length) result(mySources) +module function anisobrittle_init() result(mySources) - integer, intent(in) :: source_length - logical, dimension(:,:), allocatable :: mySources + logical, dimension(:), allocatable :: mySources class(tNode), pointer :: & phases, & phase, & sources, & src - integer :: Ninstances,sourceOffset,Nconstituents,p + integer :: Ninstances,Nconstituents,p integer, dimension(:), allocatable :: N_cl character(len=pStringLen) :: extmsg = '' print'(/,a)', ' <<<+- phase:damage:anisobrittle init -+>>>' - mySources = source_active('anisobrittle',source_length) + mySources = source_active('anisobrittle') Ninstances = count(mySources) print'(a,i2)', ' # instances: ',Ninstances; flush(IO_STDOUT) if(Ninstances == 0) return @@ -57,13 +56,12 @@ module function anisobrittle_init(source_length) result(mySources) do p = 1, phases%length + if(mySources(p)) then phase => phases%get(p) - if(count(mySources(:,p)) == 0) cycle sources => phase%get('damage') - do sourceOffset = 1, sources%length - if(mySources(sourceOffset,p)) then + associate(prm => param(p)) - src => sources%get(sourceOffset) + src => sources%get(1) N_cl = src%get_asInts('N_cl',defaultVal=emptyIntArray) prm%sum_N_cl = sum(abs(N_cl)) @@ -104,7 +102,7 @@ module function anisobrittle_init(source_length) result(mySources) ! exit if any parameter is out of range if (extmsg /= '') call IO_error(211,ext_msg=trim(extmsg)//'(damage_anisoBrittle)') endif - enddo + enddo end function anisobrittle_init diff --git a/src/phase_damage_anisoductile.f90 b/src/phase_damage_anisoductile.f90 index cf55e398b..aad361238 100644 --- a/src/phase_damage_anisoductile.f90 +++ b/src/phase_damage_anisoductile.f90 @@ -24,10 +24,9 @@ contains !> @brief module initialization !> @details reads in material parameters, allocates arrays, and does sanity checks !-------------------------------------------------------------------------------------------------- -module function anisoductile_init(source_length) result(mySources) +module function anisoductile_init() result(mySources) - integer, intent(in) :: source_length - logical, dimension(:,:), allocatable :: mySources + logical, dimension(:), allocatable :: mySources class(tNode), pointer :: & phases, & @@ -36,13 +35,13 @@ module function anisoductile_init(source_length) result(mySources) pl, & sources, & src - integer :: Ninstances,sourceOffset,Nconstituents,p + integer :: Ninstances,Nconstituents,p integer, dimension(:), allocatable :: N_sl character(len=pStringLen) :: extmsg = '' print'(/,a)', ' <<<+- phase:damage:anisoductile init -+>>>' - mySources = source_active('damage_anisoDuctile',source_length) + mySources = source_active('damage_anisoDuctile') Ninstances = count(mySources) print'(a,i2)', ' # instances: ',Ninstances; flush(IO_STDOUT) if(Ninstances == 0) return @@ -51,15 +50,15 @@ module function anisoductile_init(source_length) result(mySources) allocate(param(phases%length)) do p = 1, phases%length - phase => phases%get(p) - if(count(mySources(:,p)) == 0) cycle - mech => phase%get('mechanics') - pl => mech%get('plasticity') - sources => phase%get('source') - do sourceOffset = 1, sources%length - if(mySources(sourceOffset,p)) then + if(mySources(p)) then + phase => phases%get(p) + mech => phase%get('mechanics') + pl => mech%get('plasticity') + sources => phase%get('damage') + + associate(prm => param(p)) - src => sources%get(sourceOffset) + src => sources%get(1) N_sl = pl%get_asInts('N_sl',defaultVal=emptyIntArray) prm%q = src%get_asFloat('q') @@ -78,7 +77,7 @@ module function anisoductile_init(source_length) result(mySources) if (prm%q <= 0.0_pReal) extmsg = trim(extmsg)//' q' if (any(prm%gamma_crit < 0.0_pReal)) extmsg = trim(extmsg)//' gamma_crit' - Nconstituents=count(material_phaseAt==p) * discretization_nIPs + Nconstituents=count(material_phaseAt2==p) call phase_allocateState(damageState(p),Nconstituents,1,1,0) damageState(p)%atol = src%get_asFloat('anisoDuctile_atol',defaultVal=1.0e-3_pReal) if(any(damageState(p)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' anisoductile_atol' @@ -89,7 +88,7 @@ module function anisoductile_init(source_length) result(mySources) ! exit if any parameter is out of range if (extmsg /= '') call IO_error(211,ext_msg=trim(extmsg)//'(damage_anisoDuctile)') endif - enddo + enddo diff --git a/src/phase_damage_isobrittle.f90 b/src/phase_damage_isobrittle.f90 index e267ffcc9..973817ec6 100644 --- a/src/phase_damage_isobrittle.f90 +++ b/src/phase_damage_isobrittle.f90 @@ -22,22 +22,21 @@ contains !> @brief module initialization !> @details reads in material parameters, allocates arrays, and does sanity checks !-------------------------------------------------------------------------------------------------- -module function isobrittle_init(source_length) result(mySources) +module function isobrittle_init() result(mySources) - integer, intent(in) :: source_length - logical, dimension(:,:), allocatable :: mySources + logical, dimension(:), allocatable :: mySources class(tNode), pointer :: & phases, & phase, & sources, & src - integer :: Ninstances,sourceOffset,Nconstituents,p + integer :: Ninstances,Nconstituents,p character(len=pStringLen) :: extmsg = '' print'(/,a)', ' <<<+- phase:damage:isobrittle init -+>>>' - mySources = source_active('isobrittle',source_length) + mySources = source_active('isobrittle') Ninstances = count(mySources) print'(a,i2)', ' # instances: ',Ninstances; flush(IO_STDOUT) if(Ninstances == 0) return @@ -46,13 +45,12 @@ module function isobrittle_init(source_length) result(mySources) allocate(param(phases%length)) do p = 1, phases%length + if(mySources(p)) then phase => phases%get(p) - if(count(mySources(:,p)) == 0) cycle sources => phase%get('damage') - do sourceOffset = 1, sources%length - if(mySources(sourceOffset,p)) then + associate(prm => param(p)) - src => sources%get(sourceOffset) + src => sources%get(1) prm%W_crit = src%get_asFloat('W_crit') @@ -65,7 +63,7 @@ module function isobrittle_init(source_length) result(mySources) ! sanity checks if (prm%W_crit <= 0.0_pReal) extmsg = trim(extmsg)//' W_crit' - Nconstituents = count(material_phaseAt==p) * discretization_nIPs + Nconstituents = count(material_phaseAt2==p) call phase_allocateState(damageState(p),Nconstituents,1,1,1) damageState(p)%atol = src%get_asFloat('isoBrittle_atol',defaultVal=1.0e-3_pReal) if(any(damageState(p)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' isobrittle_atol' @@ -76,7 +74,7 @@ module function isobrittle_init(source_length) result(mySources) ! exit if any parameter is out of range if (extmsg /= '') call IO_error(211,ext_msg=trim(extmsg)//'(damage_isoBrittle)') endif - enddo + enddo diff --git a/src/phase_damage_isoductile.f90 b/src/phase_damage_isoductile.f90 index f59952222..abcef1236 100644 --- a/src/phase_damage_isoductile.f90 +++ b/src/phase_damage_isoductile.f90 @@ -24,22 +24,21 @@ contains !> @brief module initialization !> @details reads in material parameters, allocates arrays, and does sanity checks !-------------------------------------------------------------------------------------------------- -module function isoductile_init(source_length) result(mySources) +module function isoductile_init() result(mySources) - integer, intent(in) :: source_length - logical, dimension(:,:), allocatable :: mySources + logical, dimension(:), allocatable :: mySources class(tNode), pointer :: & phases, & phase, & sources, & src - integer :: Ninstances,sourceOffset,Nconstituents,p + integer :: Ninstances,Nconstituents,p character(len=pStringLen) :: extmsg = '' print'(/,a)', ' <<<+- phase:damage:isoductile init -+>>>' - mySources = source_active('isoductile',source_length) + mySources = source_active('isoductile') Ninstances = count(mySources) print'(a,i2)', ' # instances: ',Ninstances; flush(IO_STDOUT) if(Ninstances == 0) return @@ -48,13 +47,12 @@ module function isoductile_init(source_length) result(mySources) allocate(param(phases%length)) do p = 1, phases%length - phase => phases%get(p) - if(count(mySources(:,p)) == 0) cycle - sources => phase%get('damage') - do sourceOffset = 1, sources%length - if(mySources(sourceOffset,p)) then + if(mySources(p)) then + phase => phases%get(p) + sources => phase%get('damage') + associate(prm => param(p)) - src => sources%get(sourceOffset) + src => sources%get(1) prm%q = src%get_asFloat('q') prm%gamma_crit = src%get_asFloat('gamma_crit') @@ -69,7 +67,7 @@ module function isoductile_init(source_length) result(mySources) if (prm%q <= 0.0_pReal) extmsg = trim(extmsg)//' q' if (prm%gamma_crit <= 0.0_pReal) extmsg = trim(extmsg)//' gamma_crit' - Nconstituents=count(material_phaseAt==p) * discretization_nIPs + Nconstituents=count(material_phaseAt2==p) call phase_allocateState(damageState(p),Nconstituents,1,1,0) damageState(p)%atol = src%get_asFloat('isoDuctile_atol',defaultVal=1.0e-3_pReal) if(any(damageState(p)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' isoductile_atol' @@ -80,7 +78,6 @@ module function isoductile_init(source_length) result(mySources) ! exit if any parameter is out of range if (extmsg /= '') call IO_error(211,ext_msg=trim(extmsg)//'(damage_isoDuctile)') endif - enddo enddo From ab202b8e7314ca2a1049a2e4195744d9eae8de98 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 13 Feb 2021 13:07:35 +0100 Subject: [PATCH 16/25] less verbose reporting --- src/phase_damage_anisobrittle.f90 | 11 +++++---- src/phase_damage_anisoductile.f90 | 11 +++++---- src/phase_damage_isobrittle.f90 | 11 +++++---- src/phase_damage_isoductile.f90 | 9 +++---- ...phase_mechanical_plastic_dislotungsten.f90 | 7 +++--- src/phase_mechanical_plastic_dislotwin.f90 | 6 ++--- src/phase_mechanical_plastic_isotropic.f90 | 7 +++--- ...phase_mechanical_plastic_kinehardening.f90 | 7 +++--- src/phase_mechanical_plastic_none.f90 | 24 +++++-------------- src/phase_mechanical_plastic_nonlocal.f90 | 6 ++--- ...phase_mechanical_plastic_phenopowerlaw.f90 | 6 +++-- 11 files changed, 51 insertions(+), 54 deletions(-) diff --git a/src/phase_damage_anisobrittle.f90 b/src/phase_damage_anisobrittle.f90 index 07dd42945..bfd069031 100644 --- a/src/phase_damage_anisobrittle.f90 +++ b/src/phase_damage_anisobrittle.f90 @@ -40,16 +40,17 @@ module function anisobrittle_init() result(mySources) phase, & sources, & src - integer :: Ninstances,Nconstituents,p + integer :: Nconstituents,p integer, dimension(:), allocatable :: N_cl character(len=pStringLen) :: extmsg = '' - print'(/,a)', ' <<<+- phase:damage:anisobrittle init -+>>>' mySources = source_active('anisobrittle') - Ninstances = count(mySources) - print'(a,i2)', ' # instances: ',Ninstances; flush(IO_STDOUT) - if(Ninstances == 0) return + if(count(mySources) == 0) return + + print'(/,a)', ' <<<+- phase:damage:anisobrittle init -+>>>' + print'(a,i0)', ' # phases: ',count(mySources); flush(IO_STDOUT) + phases => config_material%get('phase') allocate(param(phases%length)) diff --git a/src/phase_damage_anisoductile.f90 b/src/phase_damage_anisoductile.f90 index aad361238..f50c05f07 100644 --- a/src/phase_damage_anisoductile.f90 +++ b/src/phase_damage_anisoductile.f90 @@ -39,12 +39,13 @@ module function anisoductile_init() result(mySources) integer, dimension(:), allocatable :: N_sl character(len=pStringLen) :: extmsg = '' - print'(/,a)', ' <<<+- phase:damage:anisoductile init -+>>>' - mySources = source_active('damage_anisoDuctile') - Ninstances = count(mySources) - print'(a,i2)', ' # instances: ',Ninstances; flush(IO_STDOUT) - if(Ninstances == 0) return + mySources = source_active('anisoductile') + if(count(mySources) == 0) return + + print'(/,a)', ' <<<+- phase:damage:anisoductile init -+>>>' + print'(a,i0)', ' # phases: ',count(mySources); flush(IO_STDOUT) + phases => config_material%get('phase') allocate(param(phases%length)) diff --git a/src/phase_damage_isobrittle.f90 b/src/phase_damage_isobrittle.f90 index 973817ec6..1cc0d7900 100644 --- a/src/phase_damage_isobrittle.f90 +++ b/src/phase_damage_isobrittle.f90 @@ -31,15 +31,16 @@ module function isobrittle_init() result(mySources) phase, & sources, & src - integer :: Ninstances,Nconstituents,p + integer :: Nconstituents,p character(len=pStringLen) :: extmsg = '' - print'(/,a)', ' <<<+- phase:damage:isobrittle init -+>>>' mySources = source_active('isobrittle') - Ninstances = count(mySources) - print'(a,i2)', ' # instances: ',Ninstances; flush(IO_STDOUT) - if(Ninstances == 0) return + if(count(mySources) == 0) return + + print'(/,a)', ' <<<+- phase:damage:isobrittle init -+>>>' + print'(a,i0)', ' # phases: ',count(mySources); flush(IO_STDOUT) + phases => config_material%get('phase') allocate(param(phases%length)) diff --git a/src/phase_damage_isoductile.f90 b/src/phase_damage_isoductile.f90 index abcef1236..247cd715a 100644 --- a/src/phase_damage_isoductile.f90 +++ b/src/phase_damage_isoductile.f90 @@ -36,12 +36,13 @@ module function isoductile_init() result(mySources) integer :: Ninstances,Nconstituents,p character(len=pStringLen) :: extmsg = '' - print'(/,a)', ' <<<+- phase:damage:isoductile init -+>>>' mySources = source_active('isoductile') - Ninstances = count(mySources) - print'(a,i2)', ' # instances: ',Ninstances; flush(IO_STDOUT) - if(Ninstances == 0) return + if(count(mySources) == 0) return + + print'(/,a)', ' <<<+- phase:damage:isoductile init -+>>>' + print'(a,i0)', ' # phases: ',count(mySources); flush(IO_STDOUT) + phases => config_material%get('phase') allocate(param(phases%length)) diff --git a/src/phase_mechanical_plastic_dislotungsten.f90 b/src/phase_mechanical_plastic_dislotungsten.f90 index a9946b6e7..4db546b67 100644 --- a/src/phase_mechanical_plastic_dislotungsten.f90 +++ b/src/phase_mechanical_plastic_dislotungsten.f90 @@ -97,13 +97,14 @@ module function plastic_dislotungsten_init() result(myPlasticity) mech, & pl - print'(/,a)', ' <<<+- phase:mechanics:plastic:dislotungsten init -+>>>' - myPlasticity = plastic_active('dislotungsten') Ninstances = count(myPlasticity) - print'(a,i2)', ' # instances: ',Ninstances; flush(IO_STDOUT) if(Ninstances == 0) return + print'(/,a)', ' <<<+- phase:mechanics:plastic:dislotungsten init -+>>>' + print'(a,i0)', ' # phases: ',Ninstances; flush(IO_STDOUT) + + print*, 'Cereceda et al., International Journal of Plasticity 78:242–256, 2016' print*, 'https://dx.doi.org/10.1016/j.ijplas.2015.09.002' diff --git a/src/phase_mechanical_plastic_dislotwin.f90 b/src/phase_mechanical_plastic_dislotwin.f90 index e1259fbd2..2a011dfda 100644 --- a/src/phase_mechanical_plastic_dislotwin.f90 +++ b/src/phase_mechanical_plastic_dislotwin.f90 @@ -144,13 +144,13 @@ module function plastic_dislotwin_init() result(myPlasticity) mech, & pl - print'(/,a)', ' <<<+- phase:mechanics:plastic:dislotwin init -+>>>' - myPlasticity = plastic_active('dislotwin') Ninstances = count(myPlasticity) - print'(a,i2)', ' # instances: ',Ninstances; flush(IO_STDOUT) if(Ninstances == 0) return + print'(/,a)', ' <<<+- phase:mechanics:plastic:dislotwin init -+>>>' + print'(a,i0)', ' # phases: ',Ninstances; flush(IO_STDOUT) + print*, 'Ma and Roters, Acta Materialia 52(12):3603–3612, 2004' print*, 'https://doi.org/10.1016/j.actamat.2004.04.012'//IO_EOL diff --git a/src/phase_mechanical_plastic_isotropic.f90 b/src/phase_mechanical_plastic_isotropic.f90 index 1a4e1449a..4cb1bf28c 100644 --- a/src/phase_mechanical_plastic_isotropic.f90 +++ b/src/phase_mechanical_plastic_isotropic.f90 @@ -68,13 +68,14 @@ module function plastic_isotropic_init() result(myPlasticity) mech, & pl - print'(/,a)', ' <<<+- phase:mechanics:plastic:isotropic init -+>>>' - myPlasticity = plastic_active('isotropic') Ninstances = count(myPlasticity) - print'(a,i2)', ' # instances: ',Ninstances; flush(IO_STDOUT) if(Ninstances == 0) return + print'(/,a)', ' <<<+- phase:mechanics:plastic:isotropic init -+>>>' + print'(a,i0)', ' # phses: ',Ninstances; flush(IO_STDOUT) + + print*, 'Maiti and Eisenlohr, Scripta Materialia 145:37–40, 2018' print*, 'https://doi.org/10.1016/j.scriptamat.2017.09.047' diff --git a/src/phase_mechanical_plastic_kinehardening.f90 b/src/phase_mechanical_plastic_kinehardening.f90 index f00e4171b..0ef762ffe 100644 --- a/src/phase_mechanical_plastic_kinehardening.f90 +++ b/src/phase_mechanical_plastic_kinehardening.f90 @@ -80,13 +80,14 @@ module function plastic_kinehardening_init() result(myPlasticity) mech, & pl - print'(/,a)', ' <<<+- phase:mechanics:plastic:kinehardening init -+>>>' - myPlasticity = plastic_active('kinehardening') Ninstances = count(myPlasticity) - print'(a,i2)', ' # instances: ',Ninstances; flush(IO_STDOUT) if(Ninstances == 0) return + print'(/,a)', ' <<<+- phase:mechanics:plastic:kinehardening init -+>>>' + print'(a,i0)', ' # phases: ',Ninstances; flush(IO_STDOUT) + + allocate(param(Ninstances)) allocate(state(Ninstances)) allocate(dotState(Ninstances)) diff --git a/src/phase_mechanical_plastic_none.f90 b/src/phase_mechanical_plastic_none.f90 index b95eaa039..053fe6507 100644 --- a/src/phase_mechanical_plastic_none.f90 +++ b/src/phase_mechanical_plastic_none.f90 @@ -16,32 +16,20 @@ module function plastic_none_init() result(myPlasticity) logical, dimension(:), allocatable :: myPlasticity integer :: & - Ninstances, & p, & Nconstituents class(tNode), pointer :: & - phases, & - phase, & - mech, & - pl + phases + + + myPlasticity = plastic_active('nonlocal') + if(count(myPlasticity) == 0) return print'(/,a)', ' <<<+- phase:mechanics:plastic:none init -+>>>' + print'(a,i0)', ' # phases: ',count(myPlasticity); flush(IO_STDOUT) phases => config_material%get('phase') - allocate(myPlasticity(phases%length), source = .false.) do p = 1, phases%length - phase => phases%get(p) - mech => phase%get('mechanics') - pl => mech%get ('plasticity') - if(pl%get_asString('type') == 'none') myPlasticity(p) = .true. - enddo - - Ninstances = count(myPlasticity) - print'(a,i2)', ' # instances: ',Ninstances; flush(IO_STDOUT) - if(Ninstances == 0) return - - do p = 1, phases%length - phase => phases%get(p) if(.not. myPlasticity(p)) cycle Nconstituents = count(material_phaseAt2 == p) call phase_allocateState(plasticState(p),Nconstituents,0,0,0) diff --git a/src/phase_mechanical_plastic_nonlocal.f90 b/src/phase_mechanical_plastic_nonlocal.f90 index 724086916..5cfd9cf4d 100644 --- a/src/phase_mechanical_plastic_nonlocal.f90 +++ b/src/phase_mechanical_plastic_nonlocal.f90 @@ -187,16 +187,16 @@ module function plastic_nonlocal_init() result(myPlasticity) mech, & pl - print'(/,a)', ' <<<+- phase:mechanics:plastic:nonlocal init -+>>>' - myPlasticity = plastic_active('nonlocal') Ninstances = count(myPlasticity) - print'(a,i2)', ' # instances: ',Ninstances; flush(IO_STDOUT) if(Ninstances == 0) then call geometry_plastic_nonlocal_disable return endif + print'(/,a)', ' <<<+- phase:mechanics:plastic:nonlocal init -+>>>' + print'(a,i0)', ' # phases: ',Ninstances; flush(IO_STDOUT) + print*, 'Reuber et al., Acta Materialia 71:333–348, 2014' print*, 'https://doi.org/10.1016/j.actamat.2014.03.012'//IO_EOL diff --git a/src/phase_mechanical_plastic_phenopowerlaw.f90 b/src/phase_mechanical_plastic_phenopowerlaw.f90 index f9791faa6..6a7d5ea7f 100644 --- a/src/phase_mechanical_plastic_phenopowerlaw.f90 +++ b/src/phase_mechanical_plastic_phenopowerlaw.f90 @@ -89,13 +89,15 @@ module function plastic_phenopowerlaw_init() result(myPlasticity) mech, & pl - print'(/,a)', ' <<<+- phase:mechanics:plastic:phenopowerlaw init -+>>>' myPlasticity = plastic_active('phenopowerlaw') Ninstances = count(myPlasticity) - print'(a,i2)', ' # instances: ',Ninstances; flush(IO_STDOUT) if(Ninstances == 0) return + print'(/,a)', ' <<<+- phase:mechanics:plastic:phenopowerlaw init -+>>>' + print'(a,i0)', ' # phases: ',Ninstances; flush(IO_STDOUT) + + allocate(param(Ninstances)) allocate(state(Ninstances)) allocate(dotState(Ninstances)) From 9ba185973fe539d530ab40276dc6624117661002 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 13 Feb 2021 13:33:58 +0100 Subject: [PATCH 17/25] more systematic directory structure --- PRIVATE | 2 +- examples/{SpectralMethod/Polycrystal => grid}/20grains.seeds | 0 .../{SpectralMethod/Polycrystal => grid}/20grains16x16x16.vtr | 0 .../{SpectralMethod/Polycrystal => grid}/20grains32x32x32.vtr | 0 .../{SpectralMethod/Polycrystal => grid}/20grains64x64x64.vtr | 0 examples/{SpectralMethod/Polycrystal => grid}/material.yaml | 0 examples/{SpectralMethod/Polycrystal => grid}/numerics.yaml | 0 examples/{SpectralMethod/Polycrystal => grid}/shearXY.yaml | 0 examples/{SpectralMethod/Polycrystal => grid}/shearZX.yaml | 0 examples/{SpectralMethod/Polycrystal => grid}/tensionX.yaml | 0 10 files changed, 1 insertion(+), 1 deletion(-) rename examples/{SpectralMethod/Polycrystal => grid}/20grains.seeds (100%) rename examples/{SpectralMethod/Polycrystal => grid}/20grains16x16x16.vtr (100%) rename examples/{SpectralMethod/Polycrystal => grid}/20grains32x32x32.vtr (100%) rename examples/{SpectralMethod/Polycrystal => grid}/20grains64x64x64.vtr (100%) rename examples/{SpectralMethod/Polycrystal => grid}/material.yaml (100%) rename examples/{SpectralMethod/Polycrystal => grid}/numerics.yaml (100%) rename examples/{SpectralMethod/Polycrystal => grid}/shearXY.yaml (100%) rename examples/{SpectralMethod/Polycrystal => grid}/shearZX.yaml (100%) rename examples/{SpectralMethod/Polycrystal => grid}/tensionX.yaml (100%) diff --git a/PRIVATE b/PRIVATE index 571b06ac4..d7f8460aa 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit 571b06ac4b42ed217eacfa5c71281e571bdab1df +Subproject commit d7f8460aa2afceb89f9b82054555d5dc6f9e5c43 diff --git a/examples/SpectralMethod/Polycrystal/20grains.seeds b/examples/grid/20grains.seeds similarity index 100% rename from examples/SpectralMethod/Polycrystal/20grains.seeds rename to examples/grid/20grains.seeds diff --git a/examples/SpectralMethod/Polycrystal/20grains16x16x16.vtr b/examples/grid/20grains16x16x16.vtr similarity index 100% rename from examples/SpectralMethod/Polycrystal/20grains16x16x16.vtr rename to examples/grid/20grains16x16x16.vtr diff --git a/examples/SpectralMethod/Polycrystal/20grains32x32x32.vtr b/examples/grid/20grains32x32x32.vtr similarity index 100% rename from examples/SpectralMethod/Polycrystal/20grains32x32x32.vtr rename to examples/grid/20grains32x32x32.vtr diff --git a/examples/SpectralMethod/Polycrystal/20grains64x64x64.vtr b/examples/grid/20grains64x64x64.vtr similarity index 100% rename from examples/SpectralMethod/Polycrystal/20grains64x64x64.vtr rename to examples/grid/20grains64x64x64.vtr diff --git a/examples/SpectralMethod/Polycrystal/material.yaml b/examples/grid/material.yaml similarity index 100% rename from examples/SpectralMethod/Polycrystal/material.yaml rename to examples/grid/material.yaml diff --git a/examples/SpectralMethod/Polycrystal/numerics.yaml b/examples/grid/numerics.yaml similarity index 100% rename from examples/SpectralMethod/Polycrystal/numerics.yaml rename to examples/grid/numerics.yaml diff --git a/examples/SpectralMethod/Polycrystal/shearXY.yaml b/examples/grid/shearXY.yaml similarity index 100% rename from examples/SpectralMethod/Polycrystal/shearXY.yaml rename to examples/grid/shearXY.yaml diff --git a/examples/SpectralMethod/Polycrystal/shearZX.yaml b/examples/grid/shearZX.yaml similarity index 100% rename from examples/SpectralMethod/Polycrystal/shearZX.yaml rename to examples/grid/shearZX.yaml diff --git a/examples/SpectralMethod/Polycrystal/tensionX.yaml b/examples/grid/tensionX.yaml similarity index 100% rename from examples/SpectralMethod/Polycrystal/tensionX.yaml rename to examples/grid/tensionX.yaml From b3231bf0a82989e3c90e422feba9fb389766b771 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 13 Feb 2021 14:15:41 +0100 Subject: [PATCH 18/25] avoid undefined return --- src/phase_damage.f90 | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/src/phase_damage.f90 b/src/phase_damage.f90 index ca08d8aea..2eda9540d 100644 --- a/src/phase_damage.f90 +++ b/src/phase_damage.f90 @@ -256,16 +256,20 @@ module function integrateDamageState(dt,co,ip,el) result(broken) converged_ ph = material_phaseAt(co,el) - if (phase_Nsources(ph) == 0) return me = material_phaseMemberAt(co,ip,el) + if (phase_Nsources(ph) == 0) then + broken = .false. + return + endif + converged_ = .true. broken = phase_damage_collectDotState(ph,me) if(broken) return size_so = damageState(ph)%sizeDotState damageState(ph)%state(1:size_so,me) = damageState(ph)%subState0(1:size_so,me) & - + damageState(ph)%dotState (1:size_so,me) * dt + + damageState(ph)%dotState (1:size_so,me) * dt source_dotState(1:size_so,2) = 0.0_pReal iteration: do NiterationState = 1, num%nState From 595ee7a35a3d4519d30b28e8a030a10764453362 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 13 Feb 2021 14:38:53 +0100 Subject: [PATCH 19/25] copy and paste error --- src/phase_mechanical_plastic_none.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/phase_mechanical_plastic_none.f90 b/src/phase_mechanical_plastic_none.f90 index 053fe6507..f16cf98fb 100644 --- a/src/phase_mechanical_plastic_none.f90 +++ b/src/phase_mechanical_plastic_none.f90 @@ -22,7 +22,7 @@ module function plastic_none_init() result(myPlasticity) phases - myPlasticity = plastic_active('nonlocal') + myPlasticity = plastic_active('none') if(count(myPlasticity) == 0) return print'(/,a)', ' <<<+- phase:mechanics:plastic:none init -+>>>' From 72c099dbbe9f6990d21458889bea5e355105dec9 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 13 Feb 2021 14:17:08 +0100 Subject: [PATCH 20/25] store data separetly --- PRIVATE | 2 +- src/lattice.f90 | 17 +++++++++++------ src/phase_mechanical_eigen.f90 | 15 +++++++-------- 3 files changed, 19 insertions(+), 15 deletions(-) diff --git a/PRIVATE b/PRIVATE index d7f8460aa..42ebe55f0 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit d7f8460aa2afceb89f9b82054555d5dc6f9e5c43 +Subproject commit 42ebe55f0ef1cd799115bd87e45d6025db42f7a7 diff --git a/src/lattice.f90 b/src/lattice.f90 index c9b6d99ef..e5f5453d4 100644 --- a/src/lattice.f90 +++ b/src/lattice.f90 @@ -459,7 +459,8 @@ subroutine lattice_init phase, & mech, & elasticity, & - thermal + thermal, & + damage print'(/,a)', ' <<<+- lattice init -+>>>'; flush(IO_STDOUT) @@ -535,13 +536,17 @@ subroutine lattice_init endif - lattice_D(1,1,ph) = phase%get_asFloat('D_11',defaultVal=0.0_pReal) - lattice_D(2,2,ph) = phase%get_asFloat('D_22',defaultVal=0.0_pReal) - lattice_D(3,3,ph) = phase%get_asFloat('D_33',defaultVal=0.0_pReal) - lattice_D(1:3,1:3,ph) = lattice_applyLatticeSymmetry33(lattice_D(1:3,1:3,ph), & + if (phase%contains('damage')) then + damage => phase%get('damage') + damage => damage%get(1) + lattice_D(1,1,ph) = damage%get_asFloat('D_11',defaultVal=0.0_pReal) + lattice_D(2,2,ph) = damage%get_asFloat('D_22',defaultVal=0.0_pReal) + lattice_D(3,3,ph) = damage%get_asFloat('D_33',defaultVal=0.0_pReal) + lattice_D(1:3,1:3,ph) = lattice_applyLatticeSymmetry33(lattice_D(1:3,1:3,ph), & phase%get_asString('lattice')) - lattice_M(ph) = phase%get_asFloat('M',defaultVal=0.0_pReal) + lattice_M(ph) = damage%get_asFloat('M',defaultVal=0.0_pReal) + endif ! SHOULD NOT BE PART OF LATTICE END call selfTest diff --git a/src/phase_mechanical_eigen.f90 b/src/phase_mechanical_eigen.f90 index 9eb3c5767..439f2bef9 100644 --- a/src/phase_mechanical_eigen.f90 +++ b/src/phase_mechanical_eigen.f90 @@ -74,8 +74,8 @@ module subroutine eigendeformation_init(phases) kinematics => phase%get('damage',defaultVal=emptyList) if(kinematics%length >0) then damage => kinematics%get(1) - if(damage%get_asString('type') == 'anisobrittle') phase_Nkinematics(ph) = phase_Nkinematics(ph) +1 - if(damage%get_asString('type') == 'isoductile') phase_Nkinematics(ph) = phase_Nkinematics(ph) +1 + if(damage%get_asString('type',defaultVal='n/a') == 'anisobrittle') phase_Nkinematics(ph) = phase_Nkinematics(ph) +1 + if(damage%get_asString('type',defaultVal='n/a') == 'isoductile' ) phase_Nkinematics(ph) = phase_Nkinematics(ph) +1 endif enddo @@ -113,7 +113,7 @@ function kinematics_active(kinematics_label,kinematics_length) result(active_ki kinematics => phase%get('kinematics',defaultVal=emptyList) do k = 1, kinematics%length kinematics_type => kinematics%get(k) - if(kinematics_type%get_asString('type') == kinematics_label) active_kinematics(k,p) = .true. + active_kinematics(k,p) = kinematics_type%get_asString('type') == kinematics_label enddo enddo @@ -136,17 +136,16 @@ function kinematics_active2(kinematics_label,kinematics_length) result(active_k phase, & kinematics, & kinematics_type - integer :: p,k + integer :: p phases => config_material%get('phase') allocate(active_kinematics(kinematics_length,phases%length), source = .false. ) do p = 1, phases%length phase => phases%get(p) kinematics => phase%get('damage',defaultVal=emptyList) - do k = 1, kinematics%length - kinematics_type => kinematics%get(k) - if(kinematics_type%get_asString('type') == kinematics_label) active_kinematics(k,p) = .true. - enddo + kinematics_type => kinematics%get(1) + if (.not. kinematics_type%contains('type')) continue + active_kinematics(1,p) = kinematics_type%get_asString('type',defaultVal='n/a') == kinematics_label enddo From 22a0aff488f69b07a757fdb8c1649072be2bbee8 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 13 Feb 2021 18:41:30 +0100 Subject: [PATCH 21/25] separting thermal and damage sources --- src/phase.f90 | 20 ++-- src/phase_damage.f90 | 16 +-- src/phase_mechanical.f90 | 9 +- src/phase_mechanical_eigen.f90 | 102 ++++++++++-------- ...phase_mechanical_eigen_cleavageopening.f90 | 36 +++---- ...hase_mechanical_eigen_slipplaneopening.f90 | 40 ++++--- ...hase_mechanical_eigen_thermalexpansion.f90 | 4 +- src/phase_thermal.f90 | 5 +- src/phase_thermal_dissipation.f90 | 9 +- src/phase_thermal_externalheat.f90 | 11 +- 10 files changed, 129 insertions(+), 123 deletions(-) diff --git a/src/phase.f90 b/src/phase.f90 index 8ea1fa5f7..eb6e61eb8 100644 --- a/src/phase.f90 +++ b/src/phase.f90 @@ -58,12 +58,9 @@ module phase type(tDebugOptions) :: debugCrystallite integer, dimension(:), allocatable, public :: & !< ToDo: should be protected (bug in Intel compiler) - thermal_Nsources, & - phase_Nsources, & !< number of source mechanisms active in each phase - phase_Nkinematics, & !< number of kinematic mechanisms active in each phase + phase_elasticityInstance, & phase_NstiffnessDegradations, & !< number of stiffness degradation mechanisms active in each phase - phase_plasticInstance, & !< instance of particular plasticity of each phase - phase_elasticityInstance !< instance of particular elasticity of each phase + phase_plasticInstance logical, dimension(:), allocatable, public :: & ! ToDo: should be protected (bug in Intel Compiler) phase_localPlasticity !< flags phases with local constitutive law @@ -351,7 +348,7 @@ subroutine phase_init !-------------------------------------------------------------------------------------------------- ! partition and initialize state plasticState(ph)%state = plasticState(ph)%state0 - if(phase_Nsources(ph) > 0) & + if(damageState(ph)%sizeState > 0) & damageState(ph)%state = damageState(ph)%state0 enddo PhaseLoop2 @@ -365,7 +362,7 @@ end subroutine phase_init !> @brief Allocate the components of the state structure for a given phase !-------------------------------------------------------------------------------------------------- subroutine phase_allocateState(state, & - Nconstituents,sizeState,sizeDotState,sizeDeltaState) + Nconstituents,sizeState,sizeDotState,sizeDeltaState) class(tState), intent(out) :: & state @@ -406,7 +403,7 @@ subroutine phase_restore(ce,includeL) do co = 1,homogenization_Nconstituents(material_homogenizationAt2(ce)) - if (phase_Nsources(material_phaseAt2(co,ce)) > 0) & + if (damageState(material_phaseAt2(co,ce))%sizeState > 0) & damageState(material_phaseAt2(co,ce))%state( :,material_phasememberAt2(co,ce)) = & damageState(material_phaseAt2(co,ce))%state0(:,material_phasememberAt2(co,ce)) enddo @@ -429,7 +426,7 @@ subroutine phase_forward() call thermal_forward() do ph = 1, size(damageState) - if (phase_Nsources(ph) > 0) & + if (damageState(ph)%sizeState > 0) & damageState(ph)%state0 = damageState(ph)%state enddo @@ -527,7 +524,7 @@ subroutine crystallite_init() phases => config_material%get('phase') do ph = 1, phases%length - if (phase_Nsources(ph) > 0) & + if (damageState(ph)%sizeState > 0) & allocate(damageState(ph)%subState0,source=damageState(ph)%state0) ! ToDo: hack enddo @@ -574,8 +571,7 @@ subroutine phase_windForward(ip,el) call mechanical_windForward(ph,me) - if (phase_Nsources(ph) > 0) & - damageState(ph)%state0(:,me) = damageState(ph)%state(:,me) + if(damageState(ph)%sizeState > 0) damageState(ph)%state0(:,me) = damageState(ph)%state(:,me) enddo diff --git a/src/phase_damage.f90 b/src/phase_damage.f90 index 2eda9540d..5e4e43017 100644 --- a/src/phase_damage.f90 +++ b/src/phase_damage.f90 @@ -18,6 +18,9 @@ submodule(phase) damagee integer(kind(DAMAGE_UNDEFINED_ID)), dimension(:), allocatable :: & phase_source !< active sources mechanisms of each phase + integer, dimension(:), allocatable :: & + phase_Nsources + type(tDataContainer), dimension(:), allocatable :: current interface @@ -156,6 +159,7 @@ module subroutine damage_init phase => phases%get(ph) sources => phase%get('damage',defaultVal=emptyList) if (sources%length > 1) error stop + phase_Nsources(ph) = sources%length enddo @@ -192,7 +196,6 @@ module subroutine phase_damage_getRateAndItsTangents(phiDot, dPhiDot_dPhi, phi, integer :: & ph, & co, & - so, & me phiDot = 0.0_pReal @@ -201,7 +204,7 @@ module subroutine phase_damage_getRateAndItsTangents(phiDot, dPhiDot_dPhi, phi, do co = 1, homogenization_Nconstituents(material_homogenizationAt(el)) ph = material_phaseAt(co,el) me = material_phasememberAt(co,ip,el) - do so = 1, phase_Nsources(ph) + select case(phase_source(ph)) case (DAMAGE_ISOBRITTLE_ID) call isobrittle_getRateAndItsTangent (localphiDot, dLocalphiDot_dPhi, phi, ph, me) @@ -222,7 +225,6 @@ module subroutine phase_damage_getRateAndItsTangents(phiDot, dPhiDot_dPhi, phi, end select phiDot = phiDot + localphiDot dPhiDot_dPhi = dPhiDot_dPhi + dLocalphiDot_dPhi - enddo enddo end subroutine phase_damage_getRateAndItsTangents @@ -258,7 +260,7 @@ module function integrateDamageState(dt,co,ip,el) result(broken) ph = material_phaseAt(co,el) me = material_phaseMemberAt(co,ip,el) - if (phase_Nsources(ph) == 0) then + if (damageState(ph)%sizeState == 0) then broken = .false. return endif @@ -377,7 +379,7 @@ function phase_damage_collectDotState(ph,me) result(broken) broken = .false. - if (phase_Nsources(ph)==1) then + if (damageState(ph)%sizeState > 0) then sourceType: select case (phase_source(ph)) @@ -420,7 +422,7 @@ function phase_damage_deltaState(Fe, ph, me) result(broken) broken = .false. - if (phase_Nsources(ph) == 0) return + if (damageState(ph)%sizeState == 0) return sourceType: select case (phase_source(ph)) @@ -461,7 +463,7 @@ function source_active(source_label) result(active_source) phase => phases%get(ph) sources => phase%get('damage',defaultVal=emptyList) src => sources%get(1) - active_source(ph) = src%get_asString('type') == source_label + active_source(ph) = src%get_asString('type',defaultVal = 'x') == source_label enddo diff --git a/src/phase_mechanical.f90 b/src/phase_mechanical.f90 index bc4b2cd47..ef962b7e6 100644 --- a/src/phase_mechanical.f90 +++ b/src/phase_mechanical.f90 @@ -3,6 +3,7 @@ !---------------------------------------------------------------------------------------------------- submodule(phase) mechanics + enum, bind(c); enumerator :: & ELASTICITY_UNDEFINED_ID, & ELASTICITY_HOOKE_ID, & @@ -22,8 +23,6 @@ submodule(phase) mechanics KINEMATICS_THERMAL_EXPANSION_ID end enum - integer(kind(KINEMATICS_UNDEFINED_ID)), dimension(:,:), allocatable :: & - phase_kinematics integer(kind(ELASTICITY_UNDEFINED_ID)), dimension(:), allocatable :: & phase_elasticity !< elasticity of each phase integer(kind(STIFFNESS_DEGRADATION_UNDEFINED_ID)), dimension(:,:), allocatable :: & @@ -1159,7 +1158,7 @@ module function crystallite_stress(dt,co,ip,el) result(converged_) subLp0 = phase_mechanical_Lp0(ph)%data(1:3,1:3,me) subState0 = plasticState(ph)%State0(:,me) - if (phase_Nsources(ph) > 0) & + if (damageState(ph)%sizeState > 0) & damageState(ph)%subState0(:,me) = damageState(ph)%state0(:,me) subFp0 = phase_mechanical_Fp0(ph)%data(1:3,1:3,me) @@ -1187,7 +1186,7 @@ module function crystallite_stress(dt,co,ip,el) result(converged_) subFp0 = phase_mechanical_Fp(ph)%data(1:3,1:3,me) subFi0 = phase_mechanical_Fi(ph)%data(1:3,1:3,me) subState0 = plasticState(ph)%state(:,me) - if (phase_Nsources(ph) > 0) & + if (damageState(ph)%sizeState > 0) & damageState(ph)%subState0(:,me) = damageState(ph)%state(:,me) endif @@ -1203,7 +1202,7 @@ module function crystallite_stress(dt,co,ip,el) result(converged_) phase_mechanical_Li(ph)%data(1:3,1:3,me) = subLi0 endif plasticState(ph)%state(:,me) = subState0 - if (phase_Nsources(ph) > 0) & + if (damageState(ph)%sizeState > 0) & damageState(ph)%state(:,me) = damageState(ph)%subState0(:,me) todo = subStep > num%subStepMinCryst ! still on track or already done (beyond repair) diff --git a/src/phase_mechanical_eigen.f90 b/src/phase_mechanical_eigen.f90 index 439f2bef9..d55ff8759 100644 --- a/src/phase_mechanical_eigen.f90 +++ b/src/phase_mechanical_eigen.f90 @@ -1,20 +1,26 @@ submodule(phase:mechanics) eigendeformation + integer, dimension(:), allocatable :: & + Nmodels + + integer(kind(KINEMATICS_UNDEFINED_ID)), dimension(:,:), allocatable :: & + model + integer(kind(KINEMATICS_UNDEFINED_ID)), dimension(:), allocatable :: & + model_damage + interface - module function kinematics_cleavage_opening_init(kinematics_length) result(myKinematics) - integer, intent(in) :: kinematics_length - logical, dimension(:,:), allocatable :: myKinematics + module function kinematics_cleavage_opening_init() result(myKinematics) + logical, dimension(:), allocatable :: myKinematics end function kinematics_cleavage_opening_init - module function kinematics_slipplane_opening_init(kinematics_length) result(myKinematics) - integer, intent(in) :: kinematics_length - logical, dimension(:,:), allocatable :: myKinematics + module function kinematics_slipplane_opening_init() result(myKinematics) + logical, dimension(:), allocatable :: myKinematics end function kinematics_slipplane_opening_init - module function kinematics_thermal_expansion_init(kinematics_length) result(myKinematics) + module function thermalexpansion_init(kinematics_length) result(myKinematics) integer, intent(in) :: kinematics_length logical, dimension(:,:), allocatable :: myKinematics - end function kinematics_thermal_expansion_init + end function thermalexpansion_init module subroutine kinematics_cleavage_opening_LiAndItsTangent(Ld, dLd_dTstar, S, ph,me) integer, intent(in) :: ph, me @@ -65,28 +71,27 @@ module subroutine eigendeformation_init(phases) print'(/,a)', ' <<<+- phase:mechanics:eigendeformation init -+>>>' !-------------------------------------------------------------------------------------------------- -! initialize kinematic mechanisms - allocate(phase_Nkinematics(phases%length),source = 0) +! explicit eigen mechanisms + allocate(Nmodels(phases%length),source = 0) + do ph = 1,phases%length phase => phases%get(ph) kinematics => phase%get('kinematics',defaultVal=emptyList) - phase_Nkinematics(ph) = kinematics%length - kinematics => phase%get('damage',defaultVal=emptyList) - if(kinematics%length >0) then - damage => kinematics%get(1) - if(damage%get_asString('type',defaultVal='n/a') == 'anisobrittle') phase_Nkinematics(ph) = phase_Nkinematics(ph) +1 - if(damage%get_asString('type',defaultVal='n/a') == 'isoductile' ) phase_Nkinematics(ph) = phase_Nkinematics(ph) +1 - endif + Nmodels(ph) = kinematics%length enddo - allocate(phase_kinematics(maxval(phase_Nkinematics),phases%length), source = KINEMATICS_undefined_ID) + allocate(model(maxval(Nmodels),phases%length), source = KINEMATICS_undefined_ID) - if(maxval(phase_Nkinematics) /= 0) then - where(kinematics_cleavage_opening_init(maxval(phase_Nkinematics))) phase_kinematics = KINEMATICS_cleavage_opening_ID - where(kinematics_slipplane_opening_init(maxval(phase_Nkinematics))) phase_kinematics = KINEMATICS_slipplane_opening_ID - where(kinematics_thermal_expansion_init(maxval(phase_Nkinematics))) phase_kinematics = KINEMATICS_thermal_expansion_ID + if(maxval(Nmodels) /= 0) then + where(thermalexpansion_init(maxval(Nmodels))) model = KINEMATICS_thermal_expansion_ID endif + allocate(model_damage(phases%length), source = KINEMATICS_UNDEFINED_ID) + + where(kinematics_cleavage_opening_init()) model_damage = KINEMATICS_cleavage_opening_ID + where(kinematics_slipplane_opening_init()) model_damage = KINEMATICS_slipplane_opening_ID + + end subroutine eigendeformation_init @@ -125,11 +130,10 @@ end function kinematics_active !-------------------------------------------------------------------------------------------------- !> @brief checks if a kinematic mechanism is active or not !-------------------------------------------------------------------------------------------------- -function kinematics_active2(kinematics_label,kinematics_length) result(active_kinematics) +function kinematics_active2(kinematics_label) result(active_kinematics) - character(len=*), intent(in) :: kinematics_label !< name of kinematic mechanism - integer, intent(in) :: kinematics_length !< max. number of kinematics in system - logical, dimension(:,:), allocatable :: active_kinematics + character(len=*), intent(in) :: kinematics_label !< name of kinematic mechanism + logical, dimension(:), allocatable :: active_kinematics class(tNode), pointer :: & phases, & @@ -139,13 +143,14 @@ function kinematics_active2(kinematics_label,kinematics_length) result(active_k integer :: p phases => config_material%get('phase') - allocate(active_kinematics(kinematics_length,phases%length), source = .false. ) + allocate(active_kinematics(phases%length), source = .false. ) do p = 1, phases%length phase => phases%get(p) kinematics => phase%get('damage',defaultVal=emptyList) + if(kinematics%length < 1) return kinematics_type => kinematics%get(1) if (.not. kinematics_type%contains('type')) continue - active_kinematics(1,p) = kinematics_type%get_asString('type',defaultVal='n/a') == kinematics_label + active_kinematics(p) = kinematics_type%get_asString('type',defaultVal='n/a') == kinematics_label enddo @@ -181,7 +186,9 @@ module subroutine phase_LiAndItsTangents(Li, dLi_dS, dLi_dFi, & detFi integer :: & k, i, j + logical :: active + active = .false. Li = 0.0_pReal dLi_dS = 0.0_pReal dLi_dFi = 0.0_pReal @@ -190,30 +197,37 @@ module subroutine phase_LiAndItsTangents(Li, dLi_dS, dLi_dFi, & plasticType: select case (phase_plasticity(ph)) case (PLASTICITY_isotropic_ID) plasticType call plastic_isotropic_LiAndItsTangent(my_Li, my_dLi_dS, S ,phase_plasticInstance(ph),me) - case default plasticType - my_Li = 0.0_pReal - my_dLi_dS = 0.0_pReal + Li = Li + my_Li + dLi_dS = dLi_dS + my_dLi_dS + active = .true. end select plasticType - Li = Li + my_Li - dLi_dS = dLi_dS + my_dLi_dS - KinematicsLoop: do k = 1, phase_Nkinematics(ph) - kinematicsType: select case (phase_kinematics(k,ph)) - case (KINEMATICS_cleavage_opening_ID) kinematicsType - call kinematics_cleavage_opening_LiAndItsTangent(my_Li, my_dLi_dS, S, ph, me) - case (KINEMATICS_slipplane_opening_ID) kinematicsType - call kinematics_slipplane_opening_LiAndItsTangent(my_Li, my_dLi_dS, S, ph, me) + KinematicsLoop: do k = 1, Nmodels(ph) + kinematicsType: select case (model(k,ph)) case (KINEMATICS_thermal_expansion_ID) kinematicsType call thermalexpansion_LiAndItsTangent(my_Li, my_dLi_dS, ph,me) - case default kinematicsType - my_Li = 0.0_pReal - my_dLi_dS = 0.0_pReal + Li = Li + my_Li + dLi_dS = dLi_dS + my_dLi_dS + active = .true. end select kinematicsType - Li = Li + my_Li - dLi_dS = dLi_dS + my_dLi_dS enddo KinematicsLoop + select case (model_damage(ph)) + case (KINEMATICS_cleavage_opening_ID) + call kinematics_cleavage_opening_LiAndItsTangent(my_Li, my_dLi_dS, S, ph, me) + Li = Li + my_Li + dLi_dS = dLi_dS + my_dLi_dS + active = .true. + case (KINEMATICS_slipplane_opening_ID) + call kinematics_slipplane_opening_LiAndItsTangent(my_Li, my_dLi_dS, S, ph, me) + Li = Li + my_Li + dLi_dS = dLi_dS + my_dLi_dS + active = .true. + end select + + if(.not. active) return + FiInv = math_inv33(Fi) detFi = math_det33(Fi) Li = matmul(matmul(Fi,Li),FiInv)*detFi !< push forward to intermediate configuration diff --git a/src/phase_mechanical_eigen_cleavageopening.f90 b/src/phase_mechanical_eigen_cleavageopening.f90 index 5bb9847f4..aaa1beb97 100644 --- a/src/phase_mechanical_eigen_cleavageopening.f90 +++ b/src/phase_mechanical_eigen_cleavageopening.f90 @@ -28,12 +28,11 @@ contains !> @brief module initialization !> @details reads in material parameters, allocates arrays, and does sanity checks !-------------------------------------------------------------------------------------------------- -module function kinematics_cleavage_opening_init(kinematics_length) result(myKinematics) +module function kinematics_cleavage_opening_init() result(myKinematics) - integer, intent(in) :: kinematics_length - logical, dimension(:,:), allocatable :: myKinematics + logical, dimension(:), allocatable :: myKinematics - integer :: Ninstances,p,k + integer :: p integer, dimension(:), allocatable :: N_cl !< active number of cleavage systems per family character(len=pStringLen) :: extmsg = '' class(tNode), pointer :: & @@ -42,24 +41,24 @@ module function kinematics_cleavage_opening_init(kinematics_length) result(myKin kinematics, & kinematic_type - print'(/,a)', ' <<<+- phase:mechanics:eigendeformation:cleavageopening init -+>>>' - myKinematics = kinematics_active2('anisobrittle',kinematics_length) - Ninstances = count(myKinematics) - print'(a,i2)', ' # instances: ',Ninstances; flush(IO_STDOUT) - if(Ninstances == 0) return + myKinematics = kinematics_active2('anisobrittle') + if(count(myKinematics) == 0) return + + print'(/,a)', ' <<<+- phase:mechanics:eigendeformation:cleavageopening init -+>>>' + print'(a,i2)', ' # phases: ',count(myKinematics); flush(IO_STDOUT) + phases => config_material%get('phase') allocate(param(phases%length)) do p = 1, phases%length - phase => phases%get(p) - if(count(myKinematics(:,p)) == 0) cycle - kinematics => phase%get('damage') - do k = 1, kinematics%length - if(myKinematics(k,p)) then - associate(prm => param(p)) - kinematic_type => kinematics%get(k) + if(myKinematics(p)) then + phase => phases%get(p) + kinematics => phase%get('damage') + + associate(prm => param(p)) + kinematic_type => kinematics%get(1) N_cl = kinematic_type%get_asInts('N_cl') prm%sum_N_cl = sum(abs(N_cl)) @@ -83,9 +82,8 @@ module function kinematics_cleavage_opening_init(kinematics_length) result(myKin !-------------------------------------------------------------------------------------------------- ! exit if any parameter is out of range if (extmsg /= '') call IO_error(211,ext_msg=trim(extmsg)//'(cleavage_opening)') - end associate - endif - enddo + end associate + endif enddo diff --git a/src/phase_mechanical_eigen_slipplaneopening.f90 b/src/phase_mechanical_eigen_slipplaneopening.f90 index 64988d87f..ccba98b14 100644 --- a/src/phase_mechanical_eigen_slipplaneopening.f90 +++ b/src/phase_mechanical_eigen_slipplaneopening.f90 @@ -32,12 +32,11 @@ contains !> @brief module initialization !> @details reads in material parameters, allocates arrays, and does sanity checks !-------------------------------------------------------------------------------------------------- -module function kinematics_slipplane_opening_init(kinematics_length) result(myKinematics) +module function kinematics_slipplane_opening_init() result(myKinematics) - integer, intent(in) :: kinematics_length - logical, dimension(:,:), allocatable :: myKinematics + logical, dimension(:), allocatable :: myKinematics - integer :: Ninstances,p,i,k + integer :: p,i character(len=pStringLen) :: extmsg = '' integer, dimension(:), allocatable :: N_sl real(pReal), dimension(:,:), allocatable :: d,n,t @@ -49,26 +48,26 @@ module function kinematics_slipplane_opening_init(kinematics_length) result(myKi kinematics, & kinematic_type - print'(/,a)', ' <<<+- phase:mechanics:eigendeformation:slipplaneopening init -+>>>' - myKinematics = kinematics_active2('isoductile',kinematics_length) - Ninstances = count(myKinematics) - print'(a,i2)', ' # instances: ',Ninstances; flush(IO_STDOUT) - if(Ninstances == 0) return + myKinematics = kinematics_active2('isoductile') + if(count(myKinematics) == 0) return + print'(/,a)', ' <<<+- phase:mechanics:eigendeformation:slipplaneopening init -+>>>' + print'(a,i2)', ' # phases: ',count(myKinematics); flush(IO_STDOUT) + phases => config_material%get('phase') allocate(param(phases%length)) do p = 1, phases%length - phase => phases%get(p) - mech => phase%get('mechanics') - pl => mech%get('plasticity') - if(count(myKinematics(:,p)) == 0) cycle - kinematics => phase%get('damage') - do k = 1, kinematics%length - if(myKinematics(k,p)) then - associate(prm => param(p)) - kinematic_type => kinematics%get(k) + if(myKinematics(p)) then + phase => phases%get(p) + mech => phase%get('mechanics') + pl => mech%get('plasticity') + + kinematics => phase%get('damage') + + associate(prm => param(p)) + kinematic_type => kinematics%get(1) prm%dot_o = kinematic_type%get_asFloat('dot_o') prm%q = kinematic_type%get_asFloat('q') @@ -103,9 +102,8 @@ module function kinematics_slipplane_opening_init(kinematics_length) result(myKi ! exit if any parameter is out of range if (extmsg /= '') call IO_error(211,ext_msg=trim(extmsg)//'(slipplane_opening)') - end associate - endif - enddo + end associate + endif enddo diff --git a/src/phase_mechanical_eigen_thermalexpansion.f90 b/src/phase_mechanical_eigen_thermalexpansion.f90 index 6630f0eb7..1041e0957 100644 --- a/src/phase_mechanical_eigen_thermalexpansion.f90 +++ b/src/phase_mechanical_eigen_thermalexpansion.f90 @@ -23,7 +23,7 @@ contains !> @brief module initialization !> @details reads in material parameters, allocates arrays, and does sanity checks !-------------------------------------------------------------------------------------------------- -module function kinematics_thermal_expansion_init(kinematics_length) result(myKinematics) +module function thermalexpansion_init(kinematics_length) result(myKinematics) integer, intent(in) :: kinematics_length logical, dimension(:,:), allocatable :: myKinematics @@ -77,7 +77,7 @@ module function kinematics_thermal_expansion_init(kinematics_length) result(myKi enddo -end function kinematics_thermal_expansion_init +end function thermalexpansion_init !-------------------------------------------------------------------------------------------------- diff --git a/src/phase_thermal.f90 b/src/phase_thermal.f90 index 5d3d89ccb..012554aef 100644 --- a/src/phase_thermal.f90 +++ b/src/phase_thermal.f90 @@ -3,6 +3,9 @@ !---------------------------------------------------------------------------------------------------- submodule(phase) thermal + integer, dimension(:), allocatable :: & + thermal_Nsources + type(tSourceState), allocatable, dimension(:) :: & thermalState @@ -36,8 +39,6 @@ submodule(phase) thermal end function externalheat_init - - module subroutine externalheat_dotState(ph, me) integer, intent(in) :: & ph, & diff --git a/src/phase_thermal_dissipation.f90 b/src/phase_thermal_dissipation.f90 index 3c0095401..9c75763dc 100644 --- a/src/phase_thermal_dissipation.f90 +++ b/src/phase_thermal_dissipation.f90 @@ -31,15 +31,14 @@ module function dissipation_init(source_length) result(mySources) phase, & sources, thermal, & src - integer :: Ninstances,so,Nconstituents,ph + integer :: so,Nconstituents,ph - print'(/,a)', ' <<<+- phase:thermal:dissipation init -+>>>' mySources = thermal_active('dissipation',source_length) + if(count(mySources) == 0) return + print'(/,a)', ' <<<+- phase:thermal:dissipation init -+>>>' + print'(a,i2)', ' # phases: ',count(mySources); flush(IO_STDOUT) - Ninstances = count(mySources) - print'(a,i2)', ' # instances: ',Ninstances; flush(IO_STDOUT) - if(Ninstances == 0) return phases => config_material%get('phase') allocate(param(phases%length)) diff --git a/src/phase_thermal_externalheat.f90 b/src/phase_thermal_externalheat.f90 index d2874ded0..c7f894215 100644 --- a/src/phase_thermal_externalheat.f90 +++ b/src/phase_thermal_externalheat.f90 @@ -8,7 +8,7 @@ submodule(phase:thermal) externalheat integer, dimension(:), allocatable :: & - source_thermal_externalheat_offset !< which source is my current thermal dissipation mechanism? + source_thermal_externalheat_offset !< which source is my current thermal dissipation mechanism? type :: tParameters !< container type for internal constitutive parameters real(pReal), dimension(:), allocatable :: & @@ -38,15 +38,14 @@ module function externalheat_init(source_length) result(mySources) phase, & sources, thermal, & src - integer :: Ninstances,so,Nconstituents,ph + integer :: so,Nconstituents,ph - print'(/,a)', ' <<<+- phase:thermal:externalheat init -+>>>' mySources = thermal_active('externalheat',source_length) + if(count(mySources) == 0) return + print'(/,a)', ' <<<+- phase:thermal:externalheat init -+>>>' + print'(a,i2)', ' # phases: ',count(mySources); flush(IO_STDOUT) - Ninstances = count(mySources) - print'(a,i2)', ' # instances: ',Ninstances; flush(IO_STDOUT) - if(Ninstances == 0) return phases => config_material%get('phase') allocate(param(phases%length)) From 8dc53344ecd169fa4885d3f8256b8588e101bbe2 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 13 Feb 2021 18:52:37 +0100 Subject: [PATCH 22/25] 'kinematics'=>'eigen', now part of 'mechanics' --- PRIVATE | 2 +- src/phase_mechanical_eigen.f90 | 12 ++++++++---- 2 files changed, 9 insertions(+), 5 deletions(-) diff --git a/PRIVATE b/PRIVATE index 42ebe55f0..fa5ef6139 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit 42ebe55f0ef1cd799115bd87e45d6025db42f7a7 +Subproject commit fa5ef61395381b103a4babe52be145f960af0bbe diff --git a/src/phase_mechanical_eigen.f90 b/src/phase_mechanical_eigen.f90 index d55ff8759..bee469385 100644 --- a/src/phase_mechanical_eigen.f90 +++ b/src/phase_mechanical_eigen.f90 @@ -66,7 +66,8 @@ module subroutine eigendeformation_init(phases) class(tNode), pointer :: & phase, & kinematics, & - damage + damage, & + mechanics print'(/,a)', ' <<<+- phase:mechanics:eigendeformation init -+>>>' @@ -76,7 +77,8 @@ module subroutine eigendeformation_init(phases) do ph = 1,phases%length phase => phases%get(ph) - kinematics => phase%get('kinematics',defaultVal=emptyList) + mechanics => phase%get('mechanics') + kinematics => mechanics%get('eigen',defaultVal=emptyList) Nmodels(ph) = kinematics%length enddo @@ -108,14 +110,16 @@ function kinematics_active(kinematics_label,kinematics_length) result(active_ki phases, & phase, & kinematics, & - kinematics_type + kinematics_type, & + mechanics integer :: p,k phases => config_material%get('phase') allocate(active_kinematics(kinematics_length,phases%length), source = .false. ) do p = 1, phases%length phase => phases%get(p) - kinematics => phase%get('kinematics',defaultVal=emptyList) + mechanics => phase%get('mechanics') + kinematics => mechanics%get('eigen',defaultVal=emptyList) do k = 1, kinematics%length kinematics_type => kinematics%get(k) active_kinematics(k,p) = kinematics_type%get_asString('type') == kinematics_label From 9481b168789de60187a0232789043c7ca58345b2 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 13 Feb 2021 18:57:41 +0100 Subject: [PATCH 23/25] missing renames --- src/homogenization_mechanical.f90 | 6 ++--- src/homogenization_mechanical_RGC.f90 | 4 +-- src/homogenization_mechanical_isostrain.f90 | 4 +-- src/homogenization_mechanical_pass.f90 | 4 +-- src/phase.f90 | 21 +++++++++++++++ src/phase_mechanical.f90 | 6 ++--- src/phase_mechanical_eigen.f90 | 26 +++---------------- ...phase_mechanical_eigen_cleavageopening.f90 | 4 +-- ...hase_mechanical_eigen_slipplaneopening.f90 | 4 +-- ...hase_mechanical_eigen_thermalexpansion.f90 | 4 +-- src/phase_mechanical_plastic.f90 | 4 +-- ...phase_mechanical_plastic_dislotungsten.f90 | 2 +- src/phase_mechanical_plastic_dislotwin.f90 | 2 +- src/phase_mechanical_plastic_isotropic.f90 | 2 +- ...phase_mechanical_plastic_kinehardening.f90 | 2 +- src/phase_mechanical_plastic_none.f90 | 2 +- src/phase_mechanical_plastic_nonlocal.f90 | 2 +- ...phase_mechanical_plastic_phenopowerlaw.f90 | 2 +- 18 files changed, 51 insertions(+), 50 deletions(-) diff --git a/src/homogenization_mechanical.f90 b/src/homogenization_mechanical.f90 index c05b576da..a19a61f72 100644 --- a/src/homogenization_mechanical.f90 +++ b/src/homogenization_mechanical.f90 @@ -2,7 +2,7 @@ !> @author Martin Diehl, KU Leuven !> @brief Partition F and homogenize P/dPdF !-------------------------------------------------------------------------------------------------- -submodule(homogenization) mechanics +submodule(homogenization) mechanical interface @@ -86,7 +86,7 @@ module subroutine mechanical_init(num_homog) class(tNode), pointer :: & num_homogMech - print'(/,a)', ' <<<+- homogenization:mechanics init -+>>>' + print'(/,a)', ' <<<+- homogenization:mechanical init -+>>>' allocate(homogenization_dPdF(3,3,3,3,discretization_nIPs*discretization_Nelems), source=0.0_pReal) homogenization_F0 = spread(math_I3,3,discretization_nIPs*discretization_Nelems) ! initialize to identity @@ -251,4 +251,4 @@ module subroutine mechanical_results(group_base,h) end subroutine mechanical_results -end submodule mechanics +end submodule mechanical diff --git a/src/homogenization_mechanical_RGC.f90 b/src/homogenization_mechanical_RGC.f90 index 546c2b65c..f1bcda380 100644 --- a/src/homogenization_mechanical_RGC.f90 +++ b/src/homogenization_mechanical_RGC.f90 @@ -6,7 +6,7 @@ !> @brief Relaxed grain cluster (RGC) homogenization scheme !> N_constituents is defined as p x q x r (cluster) !-------------------------------------------------------------------------------------------------- -submodule(homogenization:mechanics) RGC +submodule(homogenization:mechanical) RGC use rotations use lattice @@ -88,7 +88,7 @@ module subroutine mechanical_RGC_init(num_homogMech) homog, & homogMech - print'(/,a)', ' <<<+- homogenization:mechanics:RGC init -+>>>' + print'(/,a)', ' <<<+- homogenization:mechanical:RGC init -+>>>' Ninstances = count(homogenization_type == HOMOGENIZATION_RGC_ID) print'(a,i2)', ' # instances: ',Ninstances; flush(IO_STDOUT) diff --git a/src/homogenization_mechanical_isostrain.f90 b/src/homogenization_mechanical_isostrain.f90 index f9616de42..b492499d8 100644 --- a/src/homogenization_mechanical_isostrain.f90 +++ b/src/homogenization_mechanical_isostrain.f90 @@ -4,7 +4,7 @@ !> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH !> @brief Isostrain (full constraint Taylor assuption) homogenization scheme !-------------------------------------------------------------------------------------------------- -submodule(homogenization:mechanics) isostrain +submodule(homogenization:mechanical) isostrain enum, bind(c); enumerator :: & parallel_ID, & @@ -37,7 +37,7 @@ module subroutine mechanical_isostrain_init homog, & homogMech - print'(/,a)', ' <<<+- homogenization:mechanics:isostrain init -+>>>' + print'(/,a)', ' <<<+- homogenization:mechanical:isostrain init -+>>>' Ninstances = count(homogenization_type == HOMOGENIZATION_ISOSTRAIN_ID) print'(a,i2)', ' # instances: ',Ninstances; flush(IO_STDOUT) diff --git a/src/homogenization_mechanical_pass.f90 b/src/homogenization_mechanical_pass.f90 index 9e8f3e44c..6217e6836 100644 --- a/src/homogenization_mechanical_pass.f90 +++ b/src/homogenization_mechanical_pass.f90 @@ -4,7 +4,7 @@ !> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH !> @brief dummy homogenization homogenization scheme for 1 constituent per material point !-------------------------------------------------------------------------------------------------- -submodule(homogenization:mechanics) none +submodule(homogenization:mechanical) none contains @@ -18,7 +18,7 @@ module subroutine mechanical_pass_init h, & Nmaterialpoints - print'(/,a)', ' <<<+- homogenization:mechanics:pass init -+>>>' + print'(/,a)', ' <<<+- homogenization:mechanical:pass init -+>>>' Ninstances = count(homogenization_type == HOMOGENIZATION_NONE_ID) print'(a,i2)', ' # instances: ',Ninstances; flush(IO_STDOUT) diff --git a/src/phase.f90 b/src/phase.f90 index eb6e61eb8..f129ba92f 100644 --- a/src/phase.f90 +++ b/src/phase.f90 @@ -263,6 +263,27 @@ module phase el !< element end subroutine plastic_dependentState + + module subroutine kinematics_cleavage_opening_LiAndItsTangent(Ld, dLd_dTstar, S, ph,me) + integer, intent(in) :: ph, me + real(pReal), intent(in), dimension(3,3) :: & + S + real(pReal), intent(out), dimension(3,3) :: & + Ld !< damage velocity gradient + real(pReal), intent(out), dimension(3,3,3,3) :: & + dLd_dTstar !< derivative of Ld with respect to Tstar (4th-order tensor) + end subroutine kinematics_cleavage_opening_LiAndItsTangent + + module subroutine kinematics_slipplane_opening_LiAndItsTangent(Ld, dLd_dTstar, S, ph,me) + integer, intent(in) :: ph, me + real(pReal), intent(in), dimension(3,3) :: & + S + real(pReal), intent(out), dimension(3,3) :: & + Ld !< damage velocity gradient + real(pReal), intent(out), dimension(3,3,3,3) :: & + dLd_dTstar !< derivative of Ld with respect to Tstar (4th-order tensor) + end subroutine kinematics_slipplane_opening_LiAndItsTangent + end interface diff --git a/src/phase_mechanical.f90 b/src/phase_mechanical.f90 index ef962b7e6..c10a916da 100644 --- a/src/phase_mechanical.f90 +++ b/src/phase_mechanical.f90 @@ -1,7 +1,7 @@ !---------------------------------------------------------------------------------------------------- !> @brief internal microstructure state for all plasticity constitutive models !---------------------------------------------------------------------------------------------------- -submodule(phase) mechanics +submodule(phase) mechanical enum, bind(c); enumerator :: & @@ -202,7 +202,7 @@ module subroutine mechanical_init(phases) elastic, & stiffDegradation - print'(/,a)', ' <<<+- phase:mechanics init -+>>>' + print'(/,a)', ' <<<+- phase:mechanical init -+>>>' !------------------------------------------------------------------------------------------------- ! initialize elasticity (hooke) !ToDO: Maybe move to elastic submodule along with function homogenizedC? @@ -1499,4 +1499,4 @@ module subroutine phase_mechanical_setF(F,co,ip,el) end subroutine phase_mechanical_setF -end submodule mechanics +end submodule mechanical diff --git a/src/phase_mechanical_eigen.f90 b/src/phase_mechanical_eigen.f90 index bee469385..bee7dc2d2 100644 --- a/src/phase_mechanical_eigen.f90 +++ b/src/phase_mechanical_eigen.f90 @@ -1,4 +1,4 @@ -submodule(phase:mechanics) eigendeformation +submodule(phase:mechanical) eigen integer, dimension(:), allocatable :: & Nmodels @@ -22,26 +22,6 @@ submodule(phase:mechanics) eigendeformation logical, dimension(:,:), allocatable :: myKinematics end function thermalexpansion_init - module subroutine kinematics_cleavage_opening_LiAndItsTangent(Ld, dLd_dTstar, S, ph,me) - integer, intent(in) :: ph, me - real(pReal), intent(in), dimension(3,3) :: & - S - real(pReal), intent(out), dimension(3,3) :: & - Ld !< damage velocity gradient - real(pReal), intent(out), dimension(3,3,3,3) :: & - dLd_dTstar !< derivative of Ld with respect to Tstar (4th-order tensor) - end subroutine kinematics_cleavage_opening_LiAndItsTangent - - module subroutine kinematics_slipplane_opening_LiAndItsTangent(Ld, dLd_dTstar, S, ph,me) - integer, intent(in) :: ph, me - real(pReal), intent(in), dimension(3,3) :: & - S - real(pReal), intent(out), dimension(3,3) :: & - Ld !< damage velocity gradient - real(pReal), intent(out), dimension(3,3,3,3) :: & - dLd_dTstar !< derivative of Ld with respect to Tstar (4th-order tensor) - end subroutine kinematics_slipplane_opening_LiAndItsTangent - module subroutine thermalexpansion_LiAndItsTangent(Li, dLi_dTstar, ph,me) integer, intent(in) :: ph, me real(pReal), intent(out), dimension(3,3) :: & @@ -69,7 +49,7 @@ module subroutine eigendeformation_init(phases) damage, & mechanics - print'(/,a)', ' <<<+- phase:mechanics:eigendeformation init -+>>>' + print'(/,a)', ' <<<+- phase:mechanical:eigen init -+>>>' !-------------------------------------------------------------------------------------------------- ! explicit eigen mechanisms @@ -246,4 +226,4 @@ module subroutine phase_LiAndItsTangents(Li, dLi_dS, dLi_dFi, & end subroutine phase_LiAndItsTangents -end submodule eigendeformation +end submodule eigen diff --git a/src/phase_mechanical_eigen_cleavageopening.f90 b/src/phase_mechanical_eigen_cleavageopening.f90 index aaa1beb97..9f29b1634 100644 --- a/src/phase_mechanical_eigen_cleavageopening.f90 +++ b/src/phase_mechanical_eigen_cleavageopening.f90 @@ -4,7 +4,7 @@ !> @brief material subroutine incorporating kinematics resulting from opening of cleavage planes !> @details to be done !-------------------------------------------------------------------------------------------------- -submodule(phase:eigendeformation) cleavageopening +submodule(phase:eigen) cleavageopening type :: tParameters !< container type for internal constitutive parameters integer :: & @@ -45,7 +45,7 @@ module function kinematics_cleavage_opening_init() result(myKinematics) myKinematics = kinematics_active2('anisobrittle') if(count(myKinematics) == 0) return - print'(/,a)', ' <<<+- phase:mechanics:eigendeformation:cleavageopening init -+>>>' + print'(/,a)', ' <<<+- phase:mechanical:eigen:cleavageopening init -+>>>' print'(a,i2)', ' # phases: ',count(myKinematics); flush(IO_STDOUT) diff --git a/src/phase_mechanical_eigen_slipplaneopening.f90 b/src/phase_mechanical_eigen_slipplaneopening.f90 index ccba98b14..e8a7d65b9 100644 --- a/src/phase_mechanical_eigen_slipplaneopening.f90 +++ b/src/phase_mechanical_eigen_slipplaneopening.f90 @@ -4,7 +4,7 @@ !> @brief material subroutine incorporating kinematics resulting from opening of slip planes !> @details to be done !-------------------------------------------------------------------------------------------------- -submodule(phase:eigendeformation) slipplaneopening +submodule(phase:eigen) slipplaneopening integer, dimension(:), allocatable :: kinematics_slipplane_opening_instance @@ -51,7 +51,7 @@ module function kinematics_slipplane_opening_init() result(myKinematics) myKinematics = kinematics_active2('isoductile') if(count(myKinematics) == 0) return - print'(/,a)', ' <<<+- phase:mechanics:eigendeformation:slipplaneopening init -+>>>' + print'(/,a)', ' <<<+- phase:mechanical:eigen:slipplaneopening init -+>>>' print'(a,i2)', ' # phases: ',count(myKinematics); flush(IO_STDOUT) diff --git a/src/phase_mechanical_eigen_thermalexpansion.f90 b/src/phase_mechanical_eigen_thermalexpansion.f90 index 1041e0957..86e7fa907 100644 --- a/src/phase_mechanical_eigen_thermalexpansion.f90 +++ b/src/phase_mechanical_eigen_thermalexpansion.f90 @@ -3,7 +3,7 @@ !> @brief material subroutine incorporating kinematics resulting from thermal expansion !> @details to be done !-------------------------------------------------------------------------------------------------- -submodule(phase:eigendeformation) thermalexpansion +submodule(phase:eigen) thermalexpansion integer, dimension(:), allocatable :: kinematics_thermal_expansion_instance @@ -36,7 +36,7 @@ module function thermalexpansion_init(kinematics_length) result(myKinematics) kinematics, & kinematic_type - print'(/,a)', ' <<<+- phase:mechanics:eigendeformation:thermalexpansion init -+>>>' + print'(/,a)', ' <<<+- phase:mechanical:eigen:thermalexpansion init -+>>>' myKinematics = kinematics_active('thermal_expansion',kinematics_length) Ninstances = count(myKinematics) diff --git a/src/phase_mechanical_plastic.f90 b/src/phase_mechanical_plastic.f90 index ed1dc64fb..843273c72 100644 --- a/src/phase_mechanical_plastic.f90 +++ b/src/phase_mechanical_plastic.f90 @@ -1,4 +1,4 @@ -submodule(phase:mechanics) plastic +submodule(phase:mechanical) plastic interface @@ -226,7 +226,7 @@ contains module subroutine plastic_init - print'(/,a)', ' <<<+- phase:mechanics:plastic init -+>>>' + print'(/,a)', ' <<<+- phase:mechanical:plastic init -+>>>' where(plastic_none_init()) phase_plasticity = PLASTICITY_NONE_ID where(plastic_isotropic_init()) phase_plasticity = PLASTICITY_ISOTROPIC_ID diff --git a/src/phase_mechanical_plastic_dislotungsten.f90 b/src/phase_mechanical_plastic_dislotungsten.f90 index 4db546b67..6eb1e4abe 100644 --- a/src/phase_mechanical_plastic_dislotungsten.f90 +++ b/src/phase_mechanical_plastic_dislotungsten.f90 @@ -101,7 +101,7 @@ module function plastic_dislotungsten_init() result(myPlasticity) Ninstances = count(myPlasticity) if(Ninstances == 0) return - print'(/,a)', ' <<<+- phase:mechanics:plastic:dislotungsten init -+>>>' + print'(/,a)', ' <<<+- phase:mechanical:plastic:dislotungsten init -+>>>' print'(a,i0)', ' # phases: ',Ninstances; flush(IO_STDOUT) diff --git a/src/phase_mechanical_plastic_dislotwin.f90 b/src/phase_mechanical_plastic_dislotwin.f90 index 2a011dfda..ca385f417 100644 --- a/src/phase_mechanical_plastic_dislotwin.f90 +++ b/src/phase_mechanical_plastic_dislotwin.f90 @@ -148,7 +148,7 @@ module function plastic_dislotwin_init() result(myPlasticity) Ninstances = count(myPlasticity) if(Ninstances == 0) return - print'(/,a)', ' <<<+- phase:mechanics:plastic:dislotwin init -+>>>' + print'(/,a)', ' <<<+- phase:mechanical:plastic:dislotwin init -+>>>' print'(a,i0)', ' # phases: ',Ninstances; flush(IO_STDOUT) print*, 'Ma and Roters, Acta Materialia 52(12):3603–3612, 2004' diff --git a/src/phase_mechanical_plastic_isotropic.f90 b/src/phase_mechanical_plastic_isotropic.f90 index 4cb1bf28c..e315a38fb 100644 --- a/src/phase_mechanical_plastic_isotropic.f90 +++ b/src/phase_mechanical_plastic_isotropic.f90 @@ -72,7 +72,7 @@ module function plastic_isotropic_init() result(myPlasticity) Ninstances = count(myPlasticity) if(Ninstances == 0) return - print'(/,a)', ' <<<+- phase:mechanics:plastic:isotropic init -+>>>' + print'(/,a)', ' <<<+- phase:mechanical:plastic:isotropic init -+>>>' print'(a,i0)', ' # phses: ',Ninstances; flush(IO_STDOUT) diff --git a/src/phase_mechanical_plastic_kinehardening.f90 b/src/phase_mechanical_plastic_kinehardening.f90 index 0ef762ffe..1244988de 100644 --- a/src/phase_mechanical_plastic_kinehardening.f90 +++ b/src/phase_mechanical_plastic_kinehardening.f90 @@ -84,7 +84,7 @@ module function plastic_kinehardening_init() result(myPlasticity) Ninstances = count(myPlasticity) if(Ninstances == 0) return - print'(/,a)', ' <<<+- phase:mechanics:plastic:kinehardening init -+>>>' + print'(/,a)', ' <<<+- phase:mechanical:plastic:kinehardening init -+>>>' print'(a,i0)', ' # phases: ',Ninstances; flush(IO_STDOUT) diff --git a/src/phase_mechanical_plastic_none.f90 b/src/phase_mechanical_plastic_none.f90 index f16cf98fb..1510262cc 100644 --- a/src/phase_mechanical_plastic_none.f90 +++ b/src/phase_mechanical_plastic_none.f90 @@ -25,7 +25,7 @@ module function plastic_none_init() result(myPlasticity) myPlasticity = plastic_active('none') if(count(myPlasticity) == 0) return - print'(/,a)', ' <<<+- phase:mechanics:plastic:none init -+>>>' + print'(/,a)', ' <<<+- phase:mechanical:plastic:none init -+>>>' print'(a,i0)', ' # phases: ',count(myPlasticity); flush(IO_STDOUT) phases => config_material%get('phase') diff --git a/src/phase_mechanical_plastic_nonlocal.f90 b/src/phase_mechanical_plastic_nonlocal.f90 index 5cfd9cf4d..4e94f83ae 100644 --- a/src/phase_mechanical_plastic_nonlocal.f90 +++ b/src/phase_mechanical_plastic_nonlocal.f90 @@ -194,7 +194,7 @@ module function plastic_nonlocal_init() result(myPlasticity) return endif - print'(/,a)', ' <<<+- phase:mechanics:plastic:nonlocal init -+>>>' + print'(/,a)', ' <<<+- phase:mechanical:plastic:nonlocal init -+>>>' print'(a,i0)', ' # phases: ',Ninstances; flush(IO_STDOUT) print*, 'Reuber et al., Acta Materialia 71:333–348, 2014' diff --git a/src/phase_mechanical_plastic_phenopowerlaw.f90 b/src/phase_mechanical_plastic_phenopowerlaw.f90 index 6a7d5ea7f..bea7ef06b 100644 --- a/src/phase_mechanical_plastic_phenopowerlaw.f90 +++ b/src/phase_mechanical_plastic_phenopowerlaw.f90 @@ -94,7 +94,7 @@ module function plastic_phenopowerlaw_init() result(myPlasticity) Ninstances = count(myPlasticity) if(Ninstances == 0) return - print'(/,a)', ' <<<+- phase:mechanics:plastic:phenopowerlaw init -+>>>' + print'(/,a)', ' <<<+- phase:mechanical:plastic:phenopowerlaw init -+>>>' print'(a,i0)', ' # phases: ',Ninstances; flush(IO_STDOUT) From 5126370934dfcfc66dc808a05128ae2843a0be41 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 13 Feb 2021 19:20:52 +0100 Subject: [PATCH 24/25] cleavageopening+anisobrittle are strongly coupled --- src/phase_damage_anisobrittle.f90 | 62 +++++++++ ...phase_mechanical_eigen_cleavageopening.f90 | 122 ------------------ 2 files changed, 62 insertions(+), 122 deletions(-) diff --git a/src/phase_damage_anisobrittle.f90 b/src/phase_damage_anisobrittle.f90 index bfd069031..fd82b74c4 100644 --- a/src/phase_damage_anisobrittle.f90 +++ b/src/phase_damage_anisobrittle.f90 @@ -193,4 +193,66 @@ module subroutine anisobrittle_results(phase,group) end subroutine anisobrittle_results + +!-------------------------------------------------------------------------------------------------- +!> @brief contains the constitutive equation for calculating the velocity gradient +!-------------------------------------------------------------------------------------------------- +module subroutine kinematics_cleavage_opening_LiAndItsTangent(Ld, dLd_dTstar, S, ph,me) + + integer, intent(in) :: & + ph,me + real(pReal), intent(in), dimension(3,3) :: & + S + real(pReal), intent(out), dimension(3,3) :: & + Ld !< damage velocity gradient + real(pReal), intent(out), dimension(3,3,3,3) :: & + dLd_dTstar !< derivative of Ld with respect to Tstar (4th-order tensor) + + integer :: & + i, k, l, m, n + real(pReal) :: & + traction_d, traction_t, traction_n, traction_crit, & + udotd, dudotd_dt, udott, dudott_dt, udotn, dudotn_dt + + + Ld = 0.0_pReal + dLd_dTstar = 0.0_pReal + associate(prm => param(ph)) + do i = 1,prm%sum_N_cl + traction_crit = prm%g_crit(i)*damage_phi(ph,me)**2.0_pReal + + traction_d = math_tensordot(S,prm%cleavage_systems(1:3,1:3,1,i)) + if (abs(traction_d) > traction_crit + tol_math_check) then + udotd = sign(1.0_pReal,traction_d)* prm%dot_o * ((abs(traction_d) - traction_crit)/traction_crit)**prm%q + Ld = Ld + udotd*prm%cleavage_systems(1:3,1:3,1,i) + dudotd_dt = sign(1.0_pReal,traction_d)*udotd*prm%q / (abs(traction_d) - traction_crit) + forall (k=1:3,l=1:3,m=1:3,n=1:3) & + dLd_dTstar(k,l,m,n) = dLd_dTstar(k,l,m,n) & + + dudotd_dt*prm%cleavage_systems(k,l,1,i) * prm%cleavage_systems(m,n,1,i) + endif + + traction_t = math_tensordot(S,prm%cleavage_systems(1:3,1:3,2,i)) + if (abs(traction_t) > traction_crit + tol_math_check) then + udott = sign(1.0_pReal,traction_t)* prm%dot_o * ((abs(traction_t) - traction_crit)/traction_crit)**prm%q + Ld = Ld + udott*prm%cleavage_systems(1:3,1:3,2,i) + dudott_dt = sign(1.0_pReal,traction_t)*udott*prm%q / (abs(traction_t) - traction_crit) + forall (k=1:3,l=1:3,m=1:3,n=1:3) & + dLd_dTstar(k,l,m,n) = dLd_dTstar(k,l,m,n) & + + dudott_dt*prm%cleavage_systems(k,l,2,i) * prm%cleavage_systems(m,n,2,i) + endif + + traction_n = math_tensordot(S,prm%cleavage_systems(1:3,1:3,3,i)) + if (abs(traction_n) > traction_crit + tol_math_check) then + udotn = sign(1.0_pReal,traction_n)* prm%dot_o * ((abs(traction_n) - traction_crit)/traction_crit)**prm%q + Ld = Ld + udotn*prm%cleavage_systems(1:3,1:3,3,i) + dudotn_dt = sign(1.0_pReal,traction_n)*udotn*prm%q / (abs(traction_n) - traction_crit) + forall (k=1:3,l=1:3,m=1:3,n=1:3) & + dLd_dTstar(k,l,m,n) = dLd_dTstar(k,l,m,n) & + + dudotn_dt*prm%cleavage_systems(k,l,3,i) * prm%cleavage_systems(m,n,3,i) + endif + enddo + end associate + +end subroutine kinematics_cleavage_opening_LiAndItsTangent + end submodule anisobrittle diff --git a/src/phase_mechanical_eigen_cleavageopening.f90 b/src/phase_mechanical_eigen_cleavageopening.f90 index 9f29b1634..0f48be1a4 100644 --- a/src/phase_mechanical_eigen_cleavageopening.f90 +++ b/src/phase_mechanical_eigen_cleavageopening.f90 @@ -6,21 +6,6 @@ !-------------------------------------------------------------------------------------------------- submodule(phase:eigen) cleavageopening - type :: tParameters !< container type for internal constitutive parameters - integer :: & - sum_N_cl !< total number of cleavage planes - real(pReal) :: & - dot_o, & !< opening rate of cleavage planes - q !< damage rate sensitivity - real(pReal), dimension(:), allocatable :: & - g_crit - real(pReal), dimension(:,:,:,:), allocatable :: & - cleavage_systems - end type tParameters - - type(tParameters), dimension(:), allocatable :: param !< containers of constitutive parameters - - contains @@ -32,15 +17,6 @@ module function kinematics_cleavage_opening_init() result(myKinematics) logical, dimension(:), allocatable :: myKinematics - integer :: p - integer, dimension(:), allocatable :: N_cl !< active number of cleavage systems per family - character(len=pStringLen) :: extmsg = '' - class(tNode), pointer :: & - phases, & - phase, & - kinematics, & - kinematic_type - myKinematics = kinematics_active2('anisobrittle') if(count(myKinematics) == 0) return @@ -48,107 +24,9 @@ module function kinematics_cleavage_opening_init() result(myKinematics) print'(/,a)', ' <<<+- phase:mechanical:eigen:cleavageopening init -+>>>' print'(a,i2)', ' # phases: ',count(myKinematics); flush(IO_STDOUT) - - phases => config_material%get('phase') - allocate(param(phases%length)) - - do p = 1, phases%length - if(myKinematics(p)) then - phase => phases%get(p) - kinematics => phase%get('damage') - - associate(prm => param(p)) - kinematic_type => kinematics%get(1) - - N_cl = kinematic_type%get_asInts('N_cl') - prm%sum_N_cl = sum(abs(N_cl)) - - prm%q = kinematic_type%get_asFloat('q') - prm%dot_o = kinematic_type%get_asFloat('dot_o') - - prm%g_crit = kinematic_type%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%g_crit = math_expand(prm%g_crit,N_cl) - - ! sanity checks - if (prm%q <= 0.0_pReal) extmsg = trim(extmsg)//' q' - if (prm%dot_o <= 0.0_pReal) extmsg = trim(extmsg)//' dot_o' - if (any(prm%g_crit < 0.0_pReal)) extmsg = trim(extmsg)//' g_crit' - -!-------------------------------------------------------------------------------------------------- -! exit if any parameter is out of range - if (extmsg /= '') call IO_error(211,ext_msg=trim(extmsg)//'(cleavage_opening)') - end associate - endif - enddo - - end function kinematics_cleavage_opening_init -!-------------------------------------------------------------------------------------------------- -!> @brief contains the constitutive equation for calculating the velocity gradient -!-------------------------------------------------------------------------------------------------- -module subroutine kinematics_cleavage_opening_LiAndItsTangent(Ld, dLd_dTstar, S, ph,me) - integer, intent(in) :: & - ph,me - real(pReal), intent(in), dimension(3,3) :: & - S - real(pReal), intent(out), dimension(3,3) :: & - Ld !< damage velocity gradient - real(pReal), intent(out), dimension(3,3,3,3) :: & - dLd_dTstar !< derivative of Ld with respect to Tstar (4th-order tensor) - - integer :: & - i, k, l, m, n - real(pReal) :: & - traction_d, traction_t, traction_n, traction_crit, & - udotd, dudotd_dt, udott, dudott_dt, udotn, dudotn_dt - - - Ld = 0.0_pReal - dLd_dTstar = 0.0_pReal - associate(prm => param(ph)) - do i = 1,prm%sum_N_cl - traction_crit = prm%g_crit(i)*damage_phi(ph,me)**2.0_pReal - - traction_d = math_tensordot(S,prm%cleavage_systems(1:3,1:3,1,i)) - if (abs(traction_d) > traction_crit + tol_math_check) then - udotd = sign(1.0_pReal,traction_d)* prm%dot_o * ((abs(traction_d) - traction_crit)/traction_crit)**prm%q - Ld = Ld + udotd*prm%cleavage_systems(1:3,1:3,1,i) - dudotd_dt = sign(1.0_pReal,traction_d)*udotd*prm%q / (abs(traction_d) - traction_crit) - forall (k=1:3,l=1:3,m=1:3,n=1:3) & - dLd_dTstar(k,l,m,n) = dLd_dTstar(k,l,m,n) & - + dudotd_dt*prm%cleavage_systems(k,l,1,i) * prm%cleavage_systems(m,n,1,i) - endif - - traction_t = math_tensordot(S,prm%cleavage_systems(1:3,1:3,2,i)) - if (abs(traction_t) > traction_crit + tol_math_check) then - udott = sign(1.0_pReal,traction_t)* prm%dot_o * ((abs(traction_t) - traction_crit)/traction_crit)**prm%q - Ld = Ld + udott*prm%cleavage_systems(1:3,1:3,2,i) - dudott_dt = sign(1.0_pReal,traction_t)*udott*prm%q / (abs(traction_t) - traction_crit) - forall (k=1:3,l=1:3,m=1:3,n=1:3) & - dLd_dTstar(k,l,m,n) = dLd_dTstar(k,l,m,n) & - + dudott_dt*prm%cleavage_systems(k,l,2,i) * prm%cleavage_systems(m,n,2,i) - endif - - traction_n = math_tensordot(S,prm%cleavage_systems(1:3,1:3,3,i)) - if (abs(traction_n) > traction_crit + tol_math_check) then - udotn = sign(1.0_pReal,traction_n)* prm%dot_o * ((abs(traction_n) - traction_crit)/traction_crit)**prm%q - Ld = Ld + udotn*prm%cleavage_systems(1:3,1:3,3,i) - dudotn_dt = sign(1.0_pReal,traction_n)*udotn*prm%q / (abs(traction_n) - traction_crit) - forall (k=1:3,l=1:3,m=1:3,n=1:3) & - dLd_dTstar(k,l,m,n) = dLd_dTstar(k,l,m,n) & - + dudotn_dt*prm%cleavage_systems(k,l,3,i) * prm%cleavage_systems(m,n,3,i) - endif - enddo - end associate - -end subroutine kinematics_cleavage_opening_LiAndItsTangent end submodule cleavageopening From a34edda4d9c7098543980b9caca91ca15cf97d47 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 13 Feb 2021 22:06:44 +0100 Subject: [PATCH 25/25] fixed test for damage --- PRIVATE | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/PRIVATE b/PRIVATE index fa5ef6139..2ed5cd4ba 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit fa5ef61395381b103a4babe52be145f960af0bbe +Subproject commit 2ed5cd4ba97b44ad9c8b61ced94060aee57a2dd8