diff --git a/src/constitutive.f90 b/src/constitutive.f90 index 149d09a49..1412f5259 100644 --- a/src/constitutive.f90 +++ b/src/constitutive.f90 @@ -115,12 +115,6 @@ module constitutive integer, intent(in) :: ph end subroutine damage_results - - module subroutine mech_initializeRestorationPoints(ph,me) - integer, intent(in) :: ph, me - end subroutine mech_initializeRestorationPoints - - module subroutine mech_windForward(ph,me) integer, intent(in) :: ph, me end subroutine mech_windForward @@ -359,7 +353,6 @@ module constitutive constitutive_mech_getP, & constitutive_mech_setF, & constitutive_mech_getF, & - constitutive_initializeRestorationPoints, & constitutive_windForward, & KINEMATICS_UNDEFINED_ID ,& KINEMATICS_CLEAVAGE_OPENING_ID, & @@ -404,11 +397,9 @@ subroutine constitutive_init PhaseLoop2:do ph = 1,phases%length !-------------------------------------------------------------------------------------------------- ! partition and initialize state - plasticState(ph)%partitionedState0 = plasticState(ph)%state0 - plasticState(ph)%state = plasticState(ph)%partitionedState0 + plasticState(ph)%state = plasticState(ph)%state0 forall(so = 1:phase_Nsources(ph)) - damageState(ph)%p(so)%partitionedState0 = damageState(ph)%p(so)%state0 - damageState(ph)%p(so)%state = damageState(ph)%p(so)%partitionedState0 + damageState(ph)%p(so)%state = damageState(ph)%p(so)%state0 end forall constitutive_source_maxSizeDotState = max(constitutive_source_maxSizeDotState, & @@ -473,7 +464,6 @@ subroutine constitutive_allocateState(state, & allocate(state%atol (sizeState), source=0.0_pReal) allocate(state%state0 (sizeState,Nconstituents), source=0.0_pReal) - allocate(state%partitionedState0(sizeState,Nconstituents), source=0.0_pReal) allocate(state%state (sizeState,Nconstituents), source=0.0_pReal) allocate(state%dotState (sizeDotState,Nconstituents), source=0.0_pReal) @@ -499,8 +489,8 @@ subroutine constitutive_restore(ce,includeL) do co = 1,homogenization_Nconstituents(material_homogenizationAt2(ce)) do so = 1, phase_Nsources(material_phaseAt2(co,ce)) - damageState(material_phaseAt2(co,ce))%p(so)%state( :,material_phasememberAt2(co,ce)) = & - damageState(material_phaseAt2(co,ce))%p(so)%partitionedState0(:,material_phasememberAt2(co,ce)) + damageState(material_phaseAt2(co,ce))%p(so)%state( :,material_phasememberAt2(co,ce)) = & + damageState(material_phaseAt2(co,ce))%p(so)%state0(:,material_phasememberAt2(co,ce)) enddo enddo @@ -569,7 +559,6 @@ subroutine crystallite_init() iMax, & !< maximum number of integration points eMax !< maximum number of elements - class(tNode), pointer :: & num_crystallite, & debug_crystallite, & ! pointer to debug options for crystallite @@ -651,33 +640,6 @@ subroutine crystallite_init() end subroutine crystallite_init -!-------------------------------------------------------------------------------------------------- -!> @brief Backup data for homog cutback. -!-------------------------------------------------------------------------------------------------- -subroutine constitutive_initializeRestorationPoints(ip,el) - - integer, intent(in) :: & - ip, & !< integration point number - el !< element number - integer :: & - co, & !< constituent number - so,ph, me - - do co = 1,homogenization_Nconstituents(material_homogenizationAt(el)) - ph = material_phaseAt(co,el) - me = material_phaseMemberAt(co,ip,el) - - call mech_initializeRestorationPoints(ph,me) - - do so = 1, size(damageState(ph)%p) - damageState(ph)%p(so)%partitionedState0(:,me) = damageState(ph)%p(so)%state0(:,me) - enddo - - enddo - -end subroutine constitutive_initializeRestorationPoints - - !-------------------------------------------------------------------------------------------------- !> @brief Wind homog inc forward. !-------------------------------------------------------------------------------------------------- @@ -699,7 +661,7 @@ subroutine constitutive_windForward(ip,el) call mech_windForward(ph,me) do so = 1, phase_Nsources(material_phaseAt(co,el)) - damageState(ph)%p(so)%partitionedState0(:,me) = damageState(ph)%p(so)%state(:,me) + damageState(ph)%p(so)%state0(:,me) = damageState(ph)%p(so)%state(:,me) enddo enddo diff --git a/src/constitutive_mech.f90 b/src/constitutive_mech.f90 index 56e03ac6f..be933ef4f 100644 --- a/src/constitutive_mech.f90 +++ b/src/constitutive_mech.f90 @@ -39,16 +39,7 @@ submodule(constitutive) constitutive_mech constitutive_mech_F0, & constitutive_mech_Li0, & constitutive_mech_Lp0, & - constitutive_mech_S0, & - ! converged value at end of last homogenization increment (RGC only) - constitutive_mech_partitionedFi0, & - constitutive_mech_partitionedFp0, & - constitutive_mech_partitionedF0, & - constitutive_mech_partitionedLi0, & - constitutive_mech_partitionedLp0, & - constitutive_mech_partitionedS0 - - + constitutive_mech_S0 integer(kind(PLASTICITY_undefined_ID)), dimension(:), allocatable :: & @@ -361,23 +352,17 @@ module subroutine mech_init(phases) allocate(constitutive_mech_Fe(phases%length)) allocate(constitutive_mech_Fi(phases%length)) allocate(constitutive_mech_Fi0(phases%length)) - allocate(constitutive_mech_partitionedFi0(phases%length)) allocate(constitutive_mech_Fp(phases%length)) allocate(constitutive_mech_Fp0(phases%length)) - allocate(constitutive_mech_partitionedFp0(phases%length)) allocate(constitutive_mech_F(phases%length)) allocate(constitutive_mech_F0(phases%length)) - allocate(constitutive_mech_partitionedF0(phases%length)) allocate(constitutive_mech_Li(phases%length)) allocate(constitutive_mech_Li0(phases%length)) - allocate(constitutive_mech_partitionedLi0(phases%length)) - allocate(constitutive_mech_partitionedLp0(phases%length)) allocate(constitutive_mech_Lp0(phases%length)) allocate(constitutive_mech_Lp(phases%length)) allocate(constitutive_mech_S(phases%length)) allocate(constitutive_mech_P(phases%length)) allocate(constitutive_mech_S0(phases%length)) - allocate(constitutive_mech_partitionedS0(phases%length)) do ph = 1, phases%length Nconstituents = count(material_phaseAt == ph) * discretization_nIPs @@ -385,23 +370,17 @@ module subroutine mech_init(phases) allocate(constitutive_mech_Fi(ph)%data(3,3,Nconstituents)) allocate(constitutive_mech_Fe(ph)%data(3,3,Nconstituents)) allocate(constitutive_mech_Fi0(ph)%data(3,3,Nconstituents)) - allocate(constitutive_mech_partitionedFi0(ph)%data(3,3,Nconstituents)) allocate(constitutive_mech_Fp(ph)%data(3,3,Nconstituents)) allocate(constitutive_mech_Fp0(ph)%data(3,3,Nconstituents)) - allocate(constitutive_mech_partitionedFp0(ph)%data(3,3,Nconstituents)) allocate(constitutive_mech_Li(ph)%data(3,3,Nconstituents)) allocate(constitutive_mech_Li0(ph)%data(3,3,Nconstituents)) - allocate(constitutive_mech_partitionedLi0(ph)%data(3,3,Nconstituents)) - allocate(constitutive_mech_partitionedLp0(ph)%data(3,3,Nconstituents)) allocate(constitutive_mech_Lp0(ph)%data(3,3,Nconstituents)) allocate(constitutive_mech_Lp(ph)%data(3,3,Nconstituents)) allocate(constitutive_mech_S(ph)%data(3,3,Nconstituents),source=0.0_pReal) allocate(constitutive_mech_P(ph)%data(3,3,Nconstituents),source=0.0_pReal) allocate(constitutive_mech_S0(ph)%data(3,3,Nconstituents),source=0.0_pReal) - allocate(constitutive_mech_partitionedS0(ph)%data(3,3,Nconstituents),source=0.0_pReal) allocate(constitutive_mech_F(ph)%data(3,3,Nconstituents)) allocate(constitutive_mech_F0(ph)%data(3,3,Nconstituents)) - allocate(constitutive_mech_partitionedF0(ph)%data(3,3,Nconstituents)) phase => phases%get(ph) mech => phase%get('mechanics') @@ -454,10 +433,6 @@ module subroutine mech_init(phases) constitutive_mech_Fi(ph)%data(1:3,1:3,me) = constitutive_mech_Fi0(ph)%data(1:3,1:3,me) constitutive_mech_F(ph)%data(1:3,1:3,me) = constitutive_mech_F0(ph)%data(1:3,1:3,me) - constitutive_mech_partitionedFi0(ph)%data(1:3,1:3,me) = constitutive_mech_Fi0(ph)%data(1:3,1:3,me) - constitutive_mech_partitionedFp0(ph)%data(1:3,1:3,me) = constitutive_mech_Fp0(ph)%data(1:3,1:3,me) - constitutive_mech_partitionedF0(ph)%data(1:3,1:3,me) = constitutive_mech_F0(ph)%data(1:3,1:3,me) - enddo enddo; enddo !$OMP END PARALLEL DO @@ -1464,26 +1439,6 @@ subroutine crystallite_results(group,ph) end subroutine crystallite_results -!-------------------------------------------------------------------------------------------------- -!> @brief Backup data for homog cutback. -!-------------------------------------------------------------------------------------------------- -module subroutine mech_initializeRestorationPoints(ph,me) - - integer, intent(in) :: ph, me - - - constitutive_mech_partitionedFi0(ph)%data(1:3,1:3,me) = constitutive_mech_Fi0(ph)%data(1:3,1:3,me) - constitutive_mech_partitionedFp0(ph)%data(1:3,1:3,me) = constitutive_mech_Fp0(ph)%data(1:3,1:3,me) - constitutive_mech_partitionedF0(ph)%data(1:3,1:3,me) = constitutive_mech_F0(ph)%data(1:3,1:3,me) - constitutive_mech_partitionedLi0(ph)%data(1:3,1:3,me) = constitutive_mech_Li0(ph)%data(1:3,1:3,me) - constitutive_mech_partitionedLp0(ph)%data(1:3,1:3,me) = constitutive_mech_Lp0(ph)%data(1:3,1:3,me) - constitutive_mech_partitionedS0(ph)%data(1:3,1:3,me) = constitutive_mech_S0(ph)%data(1:3,1:3,me) - - plasticState(ph)%partitionedState0(:,me) = plasticState(ph)%state0(:,me) - -end subroutine mech_initializeRestorationPoints - - !-------------------------------------------------------------------------------------------------- !> @brief Wind homog inc forward. !-------------------------------------------------------------------------------------------------- @@ -1492,14 +1447,14 @@ module subroutine mech_windForward(ph,me) integer, intent(in) :: ph, me - constitutive_mech_partitionedFp0(ph)%data(1:3,1:3,me) = constitutive_mech_Fp(ph)%data(1:3,1:3,me) - constitutive_mech_partitionedFi0(ph)%data(1:3,1:3,me) = constitutive_mech_Fi(ph)%data(1:3,1:3,me) - constitutive_mech_partitionedF0(ph)%data(1:3,1:3,me) = constitutive_mech_F(ph)%data(1:3,1:3,me) - constitutive_mech_partitionedLi0(ph)%data(1:3,1:3,me) = constitutive_mech_Li(ph)%data(1:3,1:3,me) - constitutive_mech_partitionedLp0(ph)%data(1:3,1:3,me) = constitutive_mech_Lp(ph)%data(1:3,1:3,me) - constitutive_mech_partitionedS0(ph)%data(1:3,1:3,me) = constitutive_mech_S(ph)%data(1:3,1:3,me) + constitutive_mech_Fp0(ph)%data(1:3,1:3,me) = constitutive_mech_Fp(ph)%data(1:3,1:3,me) + constitutive_mech_Fi0(ph)%data(1:3,1:3,me) = constitutive_mech_Fi(ph)%data(1:3,1:3,me) + constitutive_mech_F0(ph)%data(1:3,1:3,me) = constitutive_mech_F(ph)%data(1:3,1:3,me) + constitutive_mech_Li0(ph)%data(1:3,1:3,me) = constitutive_mech_Li(ph)%data(1:3,1:3,me) + constitutive_mech_Lp0(ph)%data(1:3,1:3,me) = constitutive_mech_Lp(ph)%data(1:3,1:3,me) + constitutive_mech_S0(ph)%data(1:3,1:3,me) = constitutive_mech_S(ph)%data(1:3,1:3,me) - plasticState(ph)%partitionedState0(:,me) = plasticState(ph)%state(:,me) + plasticState(ph)%State0(:,me) = plasticState(ph)%state(:,me) end subroutine mech_windForward @@ -1578,17 +1533,17 @@ module function crystallite_stress(dt,co,ip,el) result(converged_) me = material_phaseMemberAt(co,ip,el) sizeDotState = plasticState(ph)%sizeDotState - subLi0 = constitutive_mech_partitionedLi0(ph)%data(1:3,1:3,me) - subLp0 = constitutive_mech_partitionedLp0(ph)%data(1:3,1:3,me) - subState0 = plasticState(ph)%partitionedState0(:,me) + subLi0 = constitutive_mech_Li0(ph)%data(1:3,1:3,me) + subLp0 = constitutive_mech_Lp0(ph)%data(1:3,1:3,me) + subState0 = plasticState(ph)%State0(:,me) do so = 1, phase_Nsources(ph) - damageState(ph)%p(so)%subState0(:,me) = damageState(ph)%p(so)%partitionedState0(:,me) + damageState(ph)%p(so)%subState0(:,me) = damageState(ph)%p(so)%state0(:,me) enddo - subFp0 = constitutive_mech_partitionedFp0(ph)%data(1:3,1:3,me) - subFi0 = constitutive_mech_partitionedFi0(ph)%data(1:3,1:3,me) - subF0 = constitutive_mech_partitionedF0(ph)%data(1:3,1:3,me) + subFp0 = constitutive_mech_Fp0(ph)%data(1:3,1:3,me) + subFi0 = constitutive_mech_Fi0(ph)%data(1:3,1:3,me) + subF0 = constitutive_mech_F0(ph)%data(1:3,1:3,me) subFrac = 0.0_pReal subStep = 1.0_pReal/num%subStepSizeCryst todo = .true. @@ -1638,7 +1593,7 @@ module function crystallite_stress(dt,co,ip,el) result(converged_) ! prepare for integration if (todo) then subF = subF0 & - + subStep * (constitutive_mech_F(ph)%data(1:3,1:3,me) - constitutive_mech_partitionedF0(ph)%data(1:3,1:3,me)) + + subStep * (constitutive_mech_F(ph)%data(1:3,1:3,me) - constitutive_mech_F0(ph)%data(1:3,1:3,me)) constitutive_mech_Fe(ph)%data(1:3,1:3,me) = matmul(subF,math_inv33(matmul(constitutive_mech_Fi(ph)%data(1:3,1:3,me), & constitutive_mech_Fp(ph)%data(1:3,1:3,me)))) converged_ = .not. integrateState(subF0,subF,subFp0,subFi0,subState0(1:sizeDotState),subStep * dt,co,ip,el) @@ -1667,15 +1622,15 @@ module subroutine mech_restore(ce,includeL) ph = material_phaseAt2(co,ce) me = material_phaseMemberAt2(co,ce) if (includeL) then - constitutive_mech_Lp(ph)%data(1:3,1:3,me) = constitutive_mech_partitionedLp0(ph)%data(1:3,1:3,me) - constitutive_mech_Li(ph)%data(1:3,1:3,me) = constitutive_mech_partitionedLi0(ph)%data(1:3,1:3,me) + constitutive_mech_Lp(ph)%data(1:3,1:3,me) = constitutive_mech_Lp0(ph)%data(1:3,1:3,me) + constitutive_mech_Li(ph)%data(1:3,1:3,me) = constitutive_mech_Li0(ph)%data(1:3,1:3,me) endif ! maybe protecting everything from overwriting makes more sense - constitutive_mech_Fp(ph)%data(1:3,1:3,me) = constitutive_mech_partitionedFp0(ph)%data(1:3,1:3,me) - constitutive_mech_Fi(ph)%data(1:3,1:3,me) = constitutive_mech_partitionedFi0(ph)%data(1:3,1:3,me) - constitutive_mech_S(ph)%data(1:3,1:3,me) = constitutive_mech_partitionedS0(ph)%data(1:3,1:3,me) + constitutive_mech_Fp(ph)%data(1:3,1:3,me) = constitutive_mech_Fp0(ph)%data(1:3,1:3,me) + constitutive_mech_Fi(ph)%data(1:3,1:3,me) = constitutive_mech_Fi0(ph)%data(1:3,1:3,me) + constitutive_mech_S(ph)%data(1:3,1:3,me) = constitutive_mech_S0(ph)%data(1:3,1:3,me) - plasticState(ph)%state(:,me) = plasticState(ph)%partitionedState0(:,me) + plasticState(ph)%state(:,me) = plasticState(ph)%State0(:,me) enddo end subroutine mech_restore @@ -1727,8 +1682,8 @@ module function constitutive_mech_dPdF(dt,co,ip,el) result(dPdF) invFp = math_inv33(constitutive_mech_Fp(ph)%data(1:3,1:3,me)) invFi = math_inv33(constitutive_mech_Fi(ph)%data(1:3,1:3,me)) - invSubFp0 = math_inv33(constitutive_mech_partitionedFp0(ph)%data(1:3,1:3,me)) - invSubFi0 = math_inv33(constitutive_mech_partitionedFi0(ph)%data(1:3,1:3,me)) + invSubFp0 = math_inv33(constitutive_mech_Fp0(ph)%data(1:3,1:3,me)) + invSubFi0 = math_inv33(constitutive_mech_Fi0(ph)%data(1:3,1:3,me)) if (sum(abs(dLidS)) < tol_math_check) then dFidS = 0.0_pReal diff --git a/src/constitutive_thermal.f90 b/src/constitutive_thermal.f90 index 57b8a3117..831fd236d 100644 --- a/src/constitutive_thermal.f90 +++ b/src/constitutive_thermal.f90 @@ -116,12 +116,11 @@ module subroutine thermal_init(phases) PhaseLoop2:do ph = 1,phases%length do so = 1,thermal_Nsources(ph) - deallocate(thermalState(ph)%p(so)%partitionedState0) - thermalState(ph)%p(so)%state = thermalState(ph)%p(so)%state0 + thermalState(ph)%p(so)%state = thermalState(ph)%p(so)%state0 enddo - thermal_source_maxSizeDotState = max(thermal_source_maxSizeDotState, & - maxval(thermalState(ph)%p%sizeDotState)) + thermal_source_maxSizeDotState = max(thermal_source_maxSizeDotState, & + maxval(thermalState(ph)%p%sizeDotState)) enddo PhaseLoop2 !-------------------------------------------------------------------------------------------------- diff --git a/src/damage_none.f90 b/src/damage_none.f90 index 078d42af7..680110d55 100644 --- a/src/damage_none.f90 +++ b/src/damage_none.f90 @@ -27,7 +27,6 @@ subroutine damage_none_init Nmaterialpoints = count(material_homogenizationAt == h) damageState_h(h)%sizeState = 0 allocate(damageState_h(h)%state0 (0,Nmaterialpoints)) - allocate(damageState_h(h)%subState0(0,Nmaterialpoints)) allocate(damageState_h(h)%state (0,Nmaterialpoints)) allocate (damage(h)%p(Nmaterialpoints), source=1.0_pReal) diff --git a/src/damage_nonlocal.f90 b/src/damage_nonlocal.f90 index 4423c1e3a..f566ebbeb 100644 --- a/src/damage_nonlocal.f90 +++ b/src/damage_nonlocal.f90 @@ -14,28 +14,17 @@ module damage_nonlocal implicit none private - type :: tParameters - character(len=pStringLen), allocatable, dimension(:) :: & - output - end type tParameters - type, private :: tNumerics real(pReal) :: & charLength !< characteristic length scale for gradient problems end type tNumerics - - type(tparameters), dimension(:), allocatable :: & - param + type(tNumerics), private :: & num public :: & damage_nonlocal_init, & - damage_nonlocal_getSourceAndItsTangent, & - damage_nonlocal_getDiffusion, & - damage_nonlocal_getMobility, & - damage_nonlocal_putNonLocalDamage, & - damage_nonlocal_results + damage_nonlocal_getDiffusion contains @@ -48,9 +37,7 @@ subroutine damage_nonlocal_init integer :: Ninstances,Nmaterialpoints,h class(tNode), pointer :: & num_generic, & - material_homogenization, & - homog, & - homogDamage + material_homogenization print'(/,a)', ' <<<+- damage_nonlocal init -+>>>'; flush(6) @@ -60,58 +47,23 @@ subroutine damage_nonlocal_init num%charLength = num_generic%get_asFloat('charLength',defaultVal=1.0_pReal) Ninstances = count(damage_type == DAMAGE_nonlocal_ID) - allocate(param(Ninstances)) material_homogenization => config_material%get('homogenization') do h = 1, material_homogenization%length if (damage_type(h) /= DAMAGE_NONLOCAL_ID) cycle - homog => material_homogenization%get(h) - homogDamage => homog%get('damage') - associate(prm => param(damage_typeInstance(h))) - -#if defined (__GFORTRAN__) - prm%output = output_asStrings(homogDamage) -#else - prm%output = homogDamage%get_asStrings('output',defaultVal=emptyStringArray) -#endif Nmaterialpoints = count(material_homogenizationAt == h) damageState_h(h)%sizeState = 1 allocate(damageState_h(h)%state0 (1,Nmaterialpoints), source=1.0_pReal) - allocate(damageState_h(h)%subState0(1,Nmaterialpoints), source=1.0_pReal) allocate(damageState_h(h)%state (1,Nmaterialpoints), source=1.0_pReal) damage(h)%p => damageState_h(h)%state(1,:) - end associate enddo end subroutine damage_nonlocal_init -!-------------------------------------------------------------------------------------------------- -!> @brief calculates homogenized damage driving forces -!-------------------------------------------------------------------------------------------------- -subroutine damage_nonlocal_getSourceAndItsTangent(phiDot, dPhiDot_dPhi, phi, ip, el) - - integer, intent(in) :: & - ip, & !< integration point number - el !< element number - real(pReal), intent(in) :: & - phi - real(pReal) :: & - phiDot, dPhiDot_dPhi - - phiDot = 0.0_pReal - dPhiDot_dPhi = 0.0_pReal - - call constitutive_damage_getRateAndItsTangents(phiDot, dPhiDot_dPhi, phi, ip, el) - phiDot = phiDot/real(homogenization_Nconstituents(material_homogenizationAt(el)),pReal) - dPhiDot_dPhi = dPhiDot_dPhi/real(homogenization_Nconstituents(material_homogenizationAt(el)),pReal) - -end subroutine damage_nonlocal_getSourceAndItsTangent - - !-------------------------------------------------------------------------------------------------- !> @brief returns homogenized non local damage diffusion tensor in reference configuration !-------------------------------------------------------------------------------------------------- @@ -139,70 +91,4 @@ function damage_nonlocal_getDiffusion(ip,el) end function damage_nonlocal_getDiffusion -!-------------------------------------------------------------------------------------------------- -!> @brief Returns homogenized nonlocal damage mobility -!-------------------------------------------------------------------------------------------------- -real(pReal) function damage_nonlocal_getMobility(ip,el) - - integer, intent(in) :: & - ip, & !< integration point number - el !< element number - integer :: & - co - - damage_nonlocal_getMobility = 0.0_pReal - - do co = 1, homogenization_Nconstituents(material_homogenizationAt(el)) - damage_nonlocal_getMobility = damage_nonlocal_getMobility + lattice_M(material_phaseAt(co,el)) - enddo - - damage_nonlocal_getMobility = damage_nonlocal_getMobility/& - real(homogenization_Nconstituents(material_homogenizationAt(el)),pReal) - -end function damage_nonlocal_getMobility - - -!-------------------------------------------------------------------------------------------------- -!> @brief updated nonlocal damage field with solution from damage phase field PDE -!-------------------------------------------------------------------------------------------------- -subroutine damage_nonlocal_putNonLocalDamage(phi,ip,el) - - integer, intent(in) :: & - ip, & !< integration point number - el !< element number - real(pReal), intent(in) :: & - phi - integer :: & - homog, & - offset - - homog = material_homogenizationAt(el) - offset = material_homogenizationMemberAt(ip,el) - damage(homog)%p(offset) = phi - -end subroutine damage_nonlocal_putNonLocalDamage - - -!-------------------------------------------------------------------------------------------------- -!> @brief writes results to HDF5 output file -!-------------------------------------------------------------------------------------------------- -subroutine damage_nonlocal_results(homog,group) - - integer, intent(in) :: homog - character(len=*), intent(in) :: group - - integer :: o - - associate(prm => param(damage_typeInstance(homog))) - outputsLoop: do o = 1,size(prm%output) - select case(prm%output(o)) - case ('phi') - call results_writeDataset(group,damage(homog)%p,prm%output(o),& - 'damage indicator','-') - end select - enddo outputsLoop - end associate - -end subroutine damage_nonlocal_results - end module damage_nonlocal diff --git a/src/grid/grid_damage_spectral.f90 b/src/grid/grid_damage_spectral.f90 index a891f070d..dc78d9e44 100644 --- a/src/grid/grid_damage_spectral.f90 +++ b/src/grid/grid_damage_spectral.f90 @@ -198,7 +198,6 @@ function grid_damage_spectral_solution(timeinc) result(solution) do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) cell = cell + 1 call damage_nonlocal_putNonLocalDamage(phi_current(i,j,k),1,cell) - homogenization_phi(cell) = phi_current(i,j,k) enddo; enddo; enddo call VecMin(solution_vec,devNull,phi_min,ierr); CHKERRQ(ierr) @@ -236,7 +235,6 @@ subroutine grid_damage_spectral_forward(cutBack) do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) cell = cell + 1 call damage_nonlocal_putNonLocalDamage(phi_current(i,j,k),1,cell) - homogenization_phi(cell) = phi_current(i,j,k) enddo; enddo; enddo else phi_lastInc = phi_current diff --git a/src/homogenization.f90 b/src/homogenization.f90 index 64c2aade7..9b8c33c2e 100644 --- a/src/homogenization.f90 +++ b/src/homogenization.f90 @@ -25,9 +25,6 @@ module homogenization !-------------------------------------------------------------------------------------------------- ! General variables for the homogenization at a material point - real(pReal), dimension(:), allocatable, public :: & - homogenization_phi, & - homogenization_dot_phi real(pReal), dimension(:,:,:), allocatable, public :: & homogenization_F0, & !< def grad of IP at start of FE increment homogenization_F !< def grad of IP to be reached at end of FE increment @@ -41,10 +38,6 @@ module homogenization type :: tNumerics integer :: & nMPstate !< materialpoint state loop limit - real(pReal) :: & - subStepMinHomog, & !< minimum (relative) size of sub-step allowed during cutback in homogenization - subStepSizeHomog, & !< size of first substep when cutback in homogenization - stepIncreaseHomog !< increase of next substep size when previous substep converged in homogenization end type tNumerics type(tNumerics) :: num @@ -75,8 +68,7 @@ module homogenization integer, intent(in) :: ce end subroutine thermal_partition - module subroutine damage_partition(phi,ce) - real(pReal), intent(in) :: phi + module subroutine damage_partition(ce) integer, intent(in) :: ce end subroutine damage_partition @@ -149,6 +141,43 @@ module homogenization real(pReal), intent(out) :: Tdot end subroutine thermal_conduction_getSource + module function damage_nonlocal_getMobility(ip,el) result(M) + integer, intent(in) :: & + ip, & !< integration point number + el !< element number + integer :: & + co + real(pReal) :: M + end function damage_nonlocal_getMobility + + module subroutine damage_nonlocal_getSourceAndItsTangent(phiDot, dPhiDot_dPhi, phi, ip, el) + + integer, intent(in) :: & + ip, & !< integration point number + el !< element number + real(pReal), intent(in) :: & + phi + real(pReal) :: & + phiDot, dPhiDot_dPhi + end subroutine damage_nonlocal_getSourceAndItsTangent + + + module subroutine damage_nonlocal_putNonLocalDamage(phi,ip,el) + + integer, intent(in) :: & + ip, & !< integration point number + el !< element number + real(pReal), intent(in) :: & + phi + + end subroutine damage_nonlocal_putNonLocalDamage + + module subroutine damage_nonlocal_results(homog,group) + + integer, intent(in) :: homog + character(len=*), intent(in) :: group + + end subroutine damage_nonlocal_results end interface public :: & @@ -158,6 +187,9 @@ module homogenization thermal_conduction_getConductivity, & thermal_conduction_getMassDensity, & thermal_conduction_getSource, & + damage_nonlocal_getMobility, & + damage_nonlocal_getSourceAndItsTangent, & + damage_nonlocal_putNonLocalDamage, & homogenization_thermal_setfield, & homogenization_thermal_T, & homogenization_forward, & @@ -183,14 +215,7 @@ subroutine homogenization_init() num_homogGeneric => num_homog%get('generic',defaultVal=emptyDict) num%nMPstate = num_homogGeneric%get_asInt ('nMPstate', defaultVal=10) - num%subStepMinHomog = num_homogGeneric%get_asFloat('subStepMin', defaultVal=1.0e-3_pReal) - num%subStepSizeHomog = num_homogGeneric%get_asFloat('subStepSize', defaultVal=0.25_pReal) - num%stepIncreaseHomog = num_homogGeneric%get_asFloat('stepIncrease', defaultVal=1.5_pReal) - if (num%nMPstate < 1) call IO_error(301,ext_msg='nMPstate') - if (num%subStepMinHomog <= 0.0_pReal) call IO_error(301,ext_msg='subStepMinHomog') - if (num%subStepSizeHomog <= 0.0_pReal) call IO_error(301,ext_msg='subStepSizeHomog') - if (num%stepIncreaseHomog <= 0.0_pReal) call IO_error(301,ext_msg='stepIncreaseHomog') call mech_init(num_homog) @@ -216,16 +241,13 @@ subroutine materialpoint_stressAndItsTangent(dt,FEsolving_execIP,FEsolving_execE ip, & !< integration point number el, & !< element number myNgrains, co, ce, ho, me, ph - real(pReal) :: & - subFrac, & - subStep logical :: & converged logical, dimension(2) :: & doneAndHappy !$OMP PARALLEL - !$OMP DO PRIVATE(ce,me,ho,myNgrains,NiterationMPstate,subFrac,converged,subStep,doneAndHappy) + !$OMP DO PRIVATE(ce,me,ho,myNgrains,NiterationMPstate,converged,doneAndHappy) do el = FEsolving_execElem(1),FEsolving_execElem(2) ho = material_homogenizationAt(el) myNgrains = homogenization_Nconstituents(ho) @@ -233,78 +255,39 @@ subroutine materialpoint_stressAndItsTangent(dt,FEsolving_execIP,FEsolving_execE ce = (el-1)*discretization_nIPs + ip me = material_homogenizationMemberAt2(ce) - call constitutive_initializeRestorationPoints(ip,el) + call constitutive_restore(ce,.false.) ! wrong name (is more a forward function) - subFrac = 0.0_pReal - converged = .false. ! pretend failed step ... - subStep = 1.0_pReal/num%subStepSizeHomog ! ... larger then the requested calculation + if(homogState(ho)%sizeState > 0) homogState(ho)%State(:,me) = homogState(ho)%State0(:,me) + if(damageState_h(ho)%sizeState > 0) damageState_h(ho)%State(:,me) = damageState_h(ho)%State0(:,me) - if (homogState(ho)%sizeState > 0) homogState(ho)%subState0(:,me) = homogState(ho)%State0(:,me) - if (damageState_h(ho)%sizeState > 0) damageState_h(ho)%subState0(:,me) = damageState_h(ho)%State0(:,me) + doneAndHappy = [.false.,.true.] - cutBackLooping: do while (.not. terminallyIll .and. subStep > num%subStepMinHomog) + NiterationMPstate = 0 + convergenceLooping: do while (.not. (terminallyIll .or. doneAndHappy(1)) & + .and. NiterationMPstate < num%nMPstate) + NiterationMPstate = NiterationMPstate + 1 - if (converged) then - subFrac = subFrac + subStep - subStep = min(1.0_pReal-subFrac,num%stepIncreaseHomog*subStep) ! introduce flexibility for step increase/acceleration - steppingNeeded: if (subStep > num%subStepMinHomog) then + if (.not. doneAndHappy(1)) then + call mech_partition(homogenization_F(1:3,1:3,ce),ip,el) + converged = .true. + do co = 1, myNgrains + converged = converged .and. crystallite_stress(dt,co,ip,el) + enddo - ! wind forward grain starting point - call constitutive_windForward(ip,el) - - if(homogState(ho)%sizeState > 0) homogState(ho)%subState0(:,me) = homogState(ho)%State(:,me) - if(damageState_h(ho)%sizeState > 0) damageState_h(ho)%subState0(:,me) = damageState_h(ho)%State(:,me) - - endif steppingNeeded - elseif ( (myNgrains == 1 .and. subStep <= 1.0 ) .or. & ! single grain already tried internal subStepping in crystallite - num%subStepSizeHomog * subStep <= num%subStepMinHomog ) then ! would require too small subStep - ! cutback makes no sense - if (.not. terminallyIll) & ! so first signals terminally ill... - print*, ' Integration point ', ip,' at element ', el, ' terminally ill' - terminallyIll = .true. ! ...and kills all others - else ! cutback makes sense - subStep = num%subStepSizeHomog * subStep ! crystallite had severe trouble, so do a significant cutback - - call constitutive_restore(ce,subStep < 1.0_pReal) - - if(homogState(ho)%sizeState > 0) homogState(ho)%State(:,me) = homogState(ho)%subState0(:,me) - if(damageState_h(ho)%sizeState > 0) damageState_h(ho)%State(:,me) = damageState_h(ho)%subState0(:,me) + if (.not. converged) then + doneAndHappy = [.true.,.false.] + else + doneAndHappy = mech_updateState(dt,homogenization_F(1:3,1:3,ce),ip,el) + converged = all(doneAndHappy) + endif endif - if (subStep > num%subStepMinHomog) doneAndHappy = [.false.,.true.] - - NiterationMPstate = 0 - convergenceLooping: do while (.not. (terminallyIll .or. doneAndHappy(1)) & - .and. NiterationMPstate < num%nMPstate) - NiterationMPstate = NiterationMPstate + 1 - -!-------------------------------------------------------------------------------------------------- -! deformation partitioning - - if (.not. doneAndHappy(1)) then - call mech_partition( homogenization_F0(1:3,1:3,ce) & - + (homogenization_F(1:3,1:3,ce)-homogenization_F0(1:3,1:3,ce))*(subStep+subFrac), & - ip,el) - converged = .true. - do co = 1, myNgrains - converged = converged .and. crystallite_stress(dt*subStep,co,ip,el) - enddo - - if (.not. converged) then - doneAndHappy = [.true.,.false.] - else - doneAndHappy = mech_updateState(dt*subStep, & - homogenization_F0(1:3,1:3,ce) & - + (homogenization_F(1:3,1:3,ce)-homogenization_F0(1:3,1:3,ce)) & - *(subStep+subFrac), & - ip,el) - converged = all(doneAndHappy) - endif - endif - - enddo convergenceLooping - enddo cutBackLooping + enddo convergenceLooping + if (.not. converged) then + if (.not. terminallyIll) print*, ' Integration point ', ip,' at element ', el, ' terminally ill' + terminallyIll = .true. + endif enddo enddo !$OMP END DO @@ -330,6 +313,26 @@ subroutine materialpoint_stressAndItsTangent(dt,FEsolving_execIP,FEsolving_execE enddo !$OMP END DO +! !$OMP DO PRIVATE(ho,ph,ce) +! do el = FEsolving_execElem(1),FEsolving_execElem(2) +! if (terminallyIll) continue +! ho = material_homogenizationAt(el) +! do ip = FEsolving_execIP(1),FEsolving_execIP(2) +! ce = (el-1)*discretization_nIPs + ip +! call damage_partition(ce) +! do co = 1, homogenization_Nconstituents(ho) +! ph = material_phaseAt(co,el) +! if (.not. thermal_stress(dt,ph,material_phaseMemberAt(co,ip,el))) then +! if (.not. terminallyIll) & ! so first signals terminally ill... +! print*, ' Integration point ', ip,' at element ', el, ' terminally ill' +! terminallyIll = .true. ! ...and kills all others +! endif +! call thermal_homogenize(ip,el) +! enddo +! enddo +! enddo +! !$OMP END DO + !$OMP DO PRIVATE(ho) elementLooping3: do el = FEsolving_execElem(1),FEsolving_execElem(2) ho = material_homogenizationAt(el) diff --git a/src/homogenization_damage.f90 b/src/homogenization_damage.f90 index 479318340..3115e6a6f 100644 --- a/src/homogenization_damage.f90 +++ b/src/homogenization_damage.f90 @@ -3,19 +3,60 @@ !-------------------------------------------------------------------------------------------------- submodule(homogenization) homogenization_damage + use lattice + + type :: tDataContainer + real(pReal), dimension(:), allocatable :: phi + end type tDataContainer + + type(tDataContainer), dimension(:), allocatable :: current + + type :: tParameters + character(len=pStringLen), allocatable, dimension(:) :: & + output + end type tParameters + + type(tparameters), dimension(:), allocatable :: & + param contains + !-------------------------------------------------------------------------------------------------- !> @brief Allocate variables and set parameters. !-------------------------------------------------------------------------------------------------- module subroutine damage_init() + class(tNode), pointer :: & + configHomogenizations, & + configHomogenization, & + configHomogenizationDamage + integer :: ho + print'(/,a)', ' <<<+- homogenization_damage init -+>>>' - allocate(homogenization_phi(discretization_nIPs*discretization_Nelems)) - allocate(homogenization_dot_phi(discretization_nIPs*discretization_Nelems)) + + configHomogenizations => config_material%get('homogenization') + allocate(param(configHomogenizations%length)) + allocate(current(configHomogenizations%length)) + + do ho = 1, configHomogenizations%length + allocate(current(ho)%phi(count(material_homogenizationAt2==ho)), source=1.0_pReal) + configHomogenization => configHomogenizations%get(ho) + associate(prm => param(ho)) + if (configHomogenization%contains('damage')) then + configHomogenizationDamage => configHomogenization%get('damage') +#if defined (__GFORTRAN__) + prm%output = output_asStrings(configHomogenizationDamage) +#else + prm%output = configHomogenizationDamage%get_asStrings('output',defaultVal=emptyStringArray) +#endif + else + prm%output = emptyStringArray + endif + end associate + enddo end subroutine damage_init @@ -23,13 +64,15 @@ end subroutine damage_init !-------------------------------------------------------------------------------------------------- !> @brief Partition temperature onto the individual constituents. !-------------------------------------------------------------------------------------------------- -module subroutine damage_partition(phi,ce) +module subroutine damage_partition(ce) - real(pReal), intent(in) :: phi + real(pReal) :: phi integer, intent(in) :: ce integer :: co + + phi = current(material_homogenizationAt2(ce))%phi(material_homogenizationMemberAt2(ce)) do co = 1, homogenization_Nconstituents(material_homogenizationAt2(ce)) call constitutive_damage_set_phi(phi,co,ce) enddo @@ -37,4 +80,94 @@ module subroutine damage_partition(phi,ce) end subroutine damage_partition + +!-------------------------------------------------------------------------------------------------- +!> @brief Returns homogenized nonlocal damage mobility +!-------------------------------------------------------------------------------------------------- +module function damage_nonlocal_getMobility(ip,el) result(M) + + integer, intent(in) :: & + ip, & !< integration point number + el !< element number + integer :: & + co + real(pReal) :: M + + M = 0.0_pReal + + do co = 1, homogenization_Nconstituents(material_homogenizationAt(el)) + M = M + lattice_M(material_phaseAt(co,el)) + enddo + + M = M/real(homogenization_Nconstituents(material_homogenizationAt(el)),pReal) + +end function damage_nonlocal_getMobility + + +!-------------------------------------------------------------------------------------------------- +!> @brief calculates homogenized damage driving forces +!-------------------------------------------------------------------------------------------------- +module subroutine damage_nonlocal_getSourceAndItsTangent(phiDot, dPhiDot_dPhi, phi, ip, el) + + integer, intent(in) :: & + ip, & !< integration point number + el !< element number + real(pReal), intent(in) :: & + phi + real(pReal) :: & + phiDot, dPhiDot_dPhi + + phiDot = 0.0_pReal + dPhiDot_dPhi = 0.0_pReal + + call constitutive_damage_getRateAndItsTangents(phiDot, dPhiDot_dPhi, phi, ip, el) + phiDot = phiDot/real(homogenization_Nconstituents(material_homogenizationAt(el)),pReal) + dPhiDot_dPhi = dPhiDot_dPhi/real(homogenization_Nconstituents(material_homogenizationAt(el)),pReal) + +end subroutine damage_nonlocal_getSourceAndItsTangent + + +!-------------------------------------------------------------------------------------------------- +!> @brief updated nonlocal damage field with solution from damage phase field PDE +!-------------------------------------------------------------------------------------------------- +module subroutine damage_nonlocal_putNonLocalDamage(phi,ip,el) + + integer, intent(in) :: & + ip, & !< integration point number + el !< element number + real(pReal), intent(in) :: & + phi + integer :: & + homog, & + offset + + homog = material_homogenizationAt(el) + offset = material_homogenizationMemberAt(ip,el) + damage(homog)%p(offset) = phi + +end subroutine damage_nonlocal_putNonLocalDamage + + +!-------------------------------------------------------------------------------------------------- +!> @brief writes results to HDF5 output file +!-------------------------------------------------------------------------------------------------- +module subroutine damage_nonlocal_results(homog,group) + + integer, intent(in) :: homog + character(len=*), intent(in) :: group + + integer :: o + + associate(prm => param(damage_typeInstance(homog))) + outputsLoop: do o = 1,size(prm%output) + select case(prm%output(o)) + case ('phi') + call results_writeDataset(group,damage(homog)%p,prm%output(o),& + 'damage indicator','-') + end select + enddo outputsLoop + end associate + +end subroutine damage_nonlocal_results + end submodule homogenization_damage diff --git a/src/homogenization_mech_RGC.f90 b/src/homogenization_mech_RGC.f90 index 24d14e059..2c7a5a8cb 100644 --- a/src/homogenization_mech_RGC.f90 +++ b/src/homogenization_mech_RGC.f90 @@ -171,7 +171,6 @@ module subroutine mech_RGC_init(num_homogMech) homogState(h)%sizeState = sizeState allocate(homogState(h)%state0 (sizeState,Nmaterialpoints), source=0.0_pReal) - allocate(homogState(h)%subState0(sizeState,Nmaterialpoints), source=0.0_pReal) allocate(homogState(h)%state (sizeState,Nmaterialpoints), source=0.0_pReal) stt%relaxationVector => homogState(h)%state(1:nIntFaceTot,:) diff --git a/src/homogenization_mech_isostrain.f90 b/src/homogenization_mech_isostrain.f90 index a56104647..df8f5fc9d 100644 --- a/src/homogenization_mech_isostrain.f90 +++ b/src/homogenization_mech_isostrain.f90 @@ -64,7 +64,6 @@ module subroutine mech_isostrain_init Nmaterialpoints = count(material_homogenizationAt == h) homogState(h)%sizeState = 0 allocate(homogState(h)%state0 (0,Nmaterialpoints)) - allocate(homogState(h)%subState0(0,Nmaterialpoints)) allocate(homogState(h)%state (0,Nmaterialpoints)) end associate diff --git a/src/homogenization_mech_none.f90 b/src/homogenization_mech_none.f90 index d434d1ca0..767dbf557 100644 --- a/src/homogenization_mech_none.f90 +++ b/src/homogenization_mech_none.f90 @@ -32,7 +32,6 @@ module subroutine mech_none_init Nmaterialpoints = count(material_homogenizationAt == h) homogState(h)%sizeState = 0 allocate(homogState(h)%state0 (0,Nmaterialpoints)) - allocate(homogState(h)%subState0(0,Nmaterialpoints)) allocate(homogState(h)%state (0,Nmaterialpoints)) enddo diff --git a/src/prec.f90 b/src/prec.f90 index 4d73462c4..1a96c75a9 100644 --- a/src/prec.f90 +++ b/src/prec.f90 @@ -32,13 +32,13 @@ module prec real(pReal), dimension(:), pointer :: p end type group_float - ! http://stackoverflow.com/questions/3948210/can-i-have-a-pointer-to-an-item-in-an-allocatable-array type :: tState integer :: & sizeState = 0, & !< size of state sizeDotState = 0, & !< size of dot state, i.e. state(1:sizeDot) follows time evolution by dotState rates offsetDeltaState = 0, & !< index offset of delta state sizeDeltaState = 0 !< size of delta state, i.e. state(offset+1:offset+sizeDelta) follows time evolution by deltaState increments + ! http://stackoverflow.com/questions/3948210 real(pReal), pointer, dimension(:), contiguous :: & atol real(pReal), pointer, dimension(:,:), contiguous :: & ! a pointer is needed here because we might point to state/doState. However, they will never point to something, but are rather allocated and, hence, contiguous @@ -47,7 +47,6 @@ module prec dotState, & !< rate of state change deltaState !< increment of state change real(pReal), allocatable, dimension(:,:) :: & - partitionedState0, & subState0 end type