separat handling of thermal constitutive response

This commit is contained in:
Martin Diehl 2021-01-17 09:30:42 +01:00
parent d1eb650f54
commit c8be2b321c
4 changed files with 57 additions and 87 deletions

View File

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

View File

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

View File

@ -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)
!----------------------------------------------------------------------------------------------

View File

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