diff --git a/src/constitutive.f90 b/src/constitutive.f90 index c6415a883..2ec947f2b 100644 --- a/src/constitutive.f90 +++ b/src/constitutive.f90 @@ -120,19 +120,15 @@ module constitutive integer, intent(in) :: ph, me end subroutine mech_initializeRestorationPoints - module subroutine thermal_initializeRestorationPoints(ph,me) + module subroutine constitutive_thermal_initializeRestorationPoints(ph,me) integer, intent(in) :: ph, me - end subroutine thermal_initializeRestorationPoints + end subroutine constitutive_thermal_initializeRestorationPoints module subroutine mech_windForward(ph,me) integer, intent(in) :: ph, me end subroutine mech_windForward - module subroutine thermal_windForward(ph,me) - integer, intent(in) :: ph, me - end subroutine thermal_windForward - module subroutine mech_forward() end subroutine mech_forward @@ -146,10 +142,6 @@ module constitutive logical, intent(in) :: includeL end subroutine mech_restore - module subroutine thermal_restore(ip,el) - integer, intent(in) :: ip, el - end subroutine thermal_restore - module function constitutive_mech_dPdF(dt,co,ip,el) result(dPdF) real(pReal), intent(in) :: dt @@ -214,14 +206,13 @@ module constitutive ! == cleaned:end =================================================================================== - module function integrateThermalState(Delta_t,co,ip,el) result(broken) + module function thermal_stress(Delta_t,ph,me) result(converged_) + real(pReal), intent(in) :: Delta_t - integer, intent(in) :: & - el, & !< element index in element loop - ip, & !< integration point index in ip loop - co !< grain index in grain loop - logical :: broken - end function integrateThermalState + integer, intent(in) :: ph, me + logical :: converged_ + + end function thermal_stress module function integrateDamageState(dt,co,ip,el) result(broken) real(pReal), intent(in) :: dt @@ -394,18 +385,19 @@ module constitutive converged, & crystallite_init, & crystallite_stress, & + thermal_stress, & constitutive_mech_dPdF, & crystallite_orientations, & crystallite_push33ToRef, & constitutive_restartWrite, & constitutive_restartRead, & - integrateThermalState, & integrateDamageState, & constitutive_thermal_setT, & constitutive_mech_getP, & constitutive_mech_setF, & constitutive_mech_getF, & constitutive_initializeRestorationPoints, & + constitutive_thermal_initializeRestorationPoints, & constitutive_windForward, & KINEMATICS_UNDEFINED_ID ,& KINEMATICS_CLEAVAGE_OPENING_ID, & @@ -553,7 +545,6 @@ subroutine constitutive_restore(ip,el,includeL) enddo call mech_restore(ip,el,includeL) - call thermal_restore(ip,el) end subroutine constitutive_restore @@ -720,7 +711,6 @@ subroutine constitutive_initializeRestorationPoints(ip,el) me = material_phaseMemberAt(co,ip,el) call mech_initializeRestorationPoints(ph,me) - call thermal_initializeRestorationPoints(ph,me) do so = 1, size(damageState(ph)%p) damageState(ph)%p(so)%partitionedState0(:,me) = damageState(ph)%p(so)%state0(:,me) @@ -750,7 +740,6 @@ subroutine constitutive_windForward(ip,el) me = material_phaseMemberAt(co,ip,el) call mech_windForward(ph,me) - call thermal_windForward(ph,me) do so = 1, phase_Nsources(material_phaseAt(co,el)) damageState(ph)%p(so)%partitionedState0(:,me) = damageState(ph)%p(so)%state(:,me) diff --git a/src/constitutive_mech.f90 b/src/constitutive_mech.f90 index 8bc85354f..09aa797d8 100644 --- a/src/constitutive_mech.f90 +++ b/src/constitutive_mech.f90 @@ -1633,9 +1633,6 @@ module function crystallite_stress(dt,co,ip,el) result(converged_) do so = 1, phase_Nsources(ph) damageState(ph)%p(so)%subState0(:,me) = damageState(ph)%p(so)%state(:,me) enddo - do so = 1, thermal_Nsources(ph) - thermalState(ph)%p(so)%subState0(:,me) = thermalState(ph)%p(so)%state(:,me) - enddo endif !-------------------------------------------------------------------------------------------------- ! cut back (reduced time and restore) @@ -1652,9 +1649,6 @@ module function crystallite_stress(dt,co,ip,el) result(converged_) do so = 1, phase_Nsources(ph) damageState(ph)%p(so)%state(:,me) = damageState(ph)%p(so)%subState0(:,me) enddo - do so = 1, thermal_Nsources(ph) - thermalState(ph)%p(so)%state(:,me) = thermalState(ph)%p(so)%subState0(:,me) - enddo todo = subStep > num%subStepMinCryst ! still on track or already done (beyond repair) endif @@ -1668,7 +1662,6 @@ module function crystallite_stress(dt,co,ip,el) result(converged_) 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) converged_ = converged_ .and. .not. integrateDamageState(subStep * dt,co,ip,el) - converged_ = converged_ .and. .not. integrateThermalState(subStep * dt,co,ip,el) endif enddo cutbackLooping diff --git a/src/constitutive_thermal.f90 b/src/constitutive_thermal.f90 index 9787cb0e4..cabd03c1e 100644 --- a/src/constitutive_thermal.f90 +++ b/src/constitutive_thermal.f90 @@ -86,7 +86,7 @@ module subroutine thermal_init(phases) Nconstituents = count(material_phaseAt == ph) * discretization_nIPs - allocate(current(ph)%T(Nconstituents)) + allocate(current(ph)%T(Nconstituents),source=300.0_pReal) phase => phases%get(ph) if(phase%contains('thermal')) then thermal => phase%get('thermal') @@ -197,30 +197,37 @@ function constitutive_thermal_collectDotState(ph,me) result(broken) end function constitutive_thermal_collectDotState +module function thermal_stress(Delta_t,ph,me) result(converged_) + + real(pReal), intent(in) :: Delta_t + integer, intent(in) :: ph, me + logical :: converged_ + + integer :: so + + + do so = 1, thermal_Nsources(ph) + thermalState(ph)%p(so)%state(:,me) = thermalState(ph)%p(so)%subState0(:,me) + enddo + converged_ = .not. integrateThermalState(Delta_t,ph,me) + +end function thermal_stress + !-------------------------------------------------------------------------------------------------- !> @brief integrate state with 1st order explicit Euler method !-------------------------------------------------------------------------------------------------- -module function integrateThermalState(Delta_t,co,ip,el) result(broken) +function integrateThermalState(Delta_t, ph,me) result(broken) real(pReal), intent(in) :: Delta_t - integer, intent(in) :: & - el, & !< element index in element loop - ip, & !< integration point index in ip loop - co !< grain index in grain loop + integer, intent(in) :: ph, me logical :: & broken integer :: & - ph, & - me, & so, & sizeDotState - - ph = material_phaseAt(co,el) - me = material_phaseMemberAt(co,ip,el) - broken = constitutive_thermal_collectDotState(ph,me) if(broken) return @@ -233,7 +240,7 @@ module function integrateThermalState(Delta_t,co,ip,el) result(broken) end function integrateThermalState -module subroutine thermal_initializeRestorationPoints(ph,me) +module subroutine constitutive_thermal_initializeRestorationPoints(ph,me) integer, intent(in) :: ph, me @@ -244,24 +251,10 @@ module subroutine thermal_initializeRestorationPoints(ph,me) thermalState(ph)%p(so)%partitionedState0(:,me) = thermalState(ph)%p(so)%state0(:,me) enddo -end subroutine thermal_initializeRestorationPoints +end subroutine constitutive_thermal_initializeRestorationPoints -module subroutine thermal_windForward(ph,me) - - integer, intent(in) :: ph, me - - integer :: so - - - do so = 1, size(thermalState(ph)%p) - thermalState(ph)%p(so)%partitionedState0(:,me) = thermalState(ph)%p(so)%state(:,me) - enddo - -end subroutine thermal_windForward - - module subroutine thermal_forward() integer :: ph, so @@ -276,26 +269,6 @@ module subroutine thermal_forward() end subroutine thermal_forward -module subroutine thermal_restore(ip,el) - - integer, intent(in) :: ip, el - - integer :: co, ph, me, so - - - do co = 1, homogenization_Nconstituents(material_homogenizationAt(el)) - ph = material_phaseAt(co,el) - me = material_phaseMemberAt(co,ip,el) - - do so = 1, size(thermalState(ph)%p) - thermalState(ph)%p(so)%state(:,me) = thermalState(ph)%p(so)%partitionedState0(:,me) - enddo - - enddo - -end subroutine thermal_restore - - !---------------------------------------------------------------------------------------------- !< @brief Get temperature (for use by non-thermal physics) !---------------------------------------------------------------------------------------------- diff --git a/src/homogenization.f90 b/src/homogenization.f90 index 9112562b9..29d171062 100644 --- a/src/homogenization.f90 +++ b/src/homogenization.f90 @@ -161,7 +161,7 @@ subroutine materialpoint_stressAndItsTangent(dt,FEsolving_execIP,FEsolving_execE NiterationMPstate, & ip, & !< integration point number el, & !< element number - myNgrains, co, ce, ho, me + myNgrains, co, ce, ho, me, ph real(pReal) :: & subFrac, & subStep @@ -221,9 +221,8 @@ subroutine materialpoint_stressAndItsTangent(dt,FEsolving_execIP,FEsolving_execE if (subStep > num%subStepMinHomog) doneAndHappy = [.false.,.true.] NiterationMPstate = 0 - convergenceLooping: do while (.not. terminallyIll & - .and. .not. doneAndHappy(1) & - .and. NiterationMPstate < num%nMPstate) + convergenceLooping: do while (.not. (terminallyIll .or. doneAndHappy(1)) & + .and. NiterationMPstate < num%nMPstate) NiterationMPstate = NiterationMPstate + 1 !-------------------------------------------------------------------------------------------------- @@ -231,10 +230,9 @@ subroutine materialpoint_stressAndItsTangent(dt,FEsolving_execIP,FEsolving_execE if (.not. doneAndHappy(1)) then ce = (el-1)*discretization_nIPs + ip - 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) + 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) @@ -260,12 +258,29 @@ subroutine materialpoint_stressAndItsTangent(dt,FEsolving_execIP,FEsolving_execE !$OMP END PARALLEL DO if (.not. terminallyIll ) then - !$OMP PARALLEL DO PRIVATE(ho,myNgrains) + !$OMP PARALLEL DO PRIVATE(ho,ph) + do el = FEsolving_execElem(1),FEsolving_execElem(2) + if (terminallyIll) continue + ho = material_homogenizationAt(el) + do ip = FEsolving_execIP(1),FEsolving_execIP(2) + do co = 1, homogenization_Nconstituents(ho) + ph = material_phaseAt(co,el) + call constitutive_thermal_initializeRestorationPoints(ph,material_phaseMemberAt(co,ip,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 + enddo + enddo + enddo + !$OMP END PARALLEL DO + + !$OMP PARALLEL DO PRIVATE(ho) elementLooping3: do el = FEsolving_execElem(1),FEsolving_execElem(2) ho = material_homogenizationAt(el) - myNgrains = homogenization_Nconstituents(ho) IpLooping3: do ip = FEsolving_execIP(1),FEsolving_execIP(2) - do co = 1, myNgrains + do co = 1, homogenization_Nconstituents(ho) call crystallite_orientations(co,ip,el) enddo call mech_homogenize(dt,ip,el)