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

Merge branch 'separate-thermal-homog' into 'development' 

Separate thermal homog

See merge request damask/DAMASK!323
This commit is contained in:
Franz Roters 2021-01-19 11:38:17 +01:00
parent f2c233961b d2ae859b1c
commit 388e233486
9 changed files with 215 additions and 221 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
src/DAMASK_marc.f90
View File

@ -365,7 +365,7 @@ subroutine flux(f,ts,n,time)
f f
f(2) = 0.0_pReal f(2) = 0.0_pReal
call thermal_conduction_getSource(f(1), ts(3), n(3),mesh_FEM2DAMASK_elem(n(1))) call thermal_conduction_getSource(f(1), n(3),mesh_FEM2DAMASK_elem(n(1)))
end subroutine flux end subroutine flux

39
src/constitutive.f90
View File

@ -120,19 +120,15 @@ module constitutive
integer, intent(in) :: ph, me integer, intent(in) :: ph, me
end subroutine mech_initializeRestorationPoints end subroutine mech_initializeRestorationPoints
module subroutine thermal_initializeRestorationPoints(ph,me) module subroutine constitutive_thermal_initializeRestorationPoints(ph,me)
integer, intent(in) :: ph, me integer, intent(in) :: ph, me
end subroutine thermal_initializeRestorationPoints end subroutine constitutive_thermal_initializeRestorationPoints
module subroutine mech_windForward(ph,me) module subroutine mech_windForward(ph,me)
integer, intent(in) :: ph, me integer, intent(in) :: ph, me
end subroutine mech_windForward end subroutine mech_windForward
module subroutine thermal_windForward(ph,me)
integer, intent(in) :: ph, me
end subroutine thermal_windForward
module subroutine mech_forward() module subroutine mech_forward()
end subroutine mech_forward end subroutine mech_forward
@ -146,10 +142,6 @@ module constitutive
logical, intent(in) :: includeL logical, intent(in) :: includeL
end subroutine mech_restore 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) module function constitutive_mech_dPdF(dt,co,ip,el) result(dPdF)
real(pReal), intent(in) :: dt real(pReal), intent(in) :: dt
@ -207,21 +199,20 @@ module constitutive
integer, intent(in) :: co, ip, el integer, intent(in) :: co, ip, el
end subroutine constitutive_mech_setF end subroutine constitutive_mech_setF
module subroutine constitutive_thermal_setT(T,co,ip,el) module subroutine constitutive_thermal_setT(T,co,ce)
real(pReal), intent(in) :: T real(pReal), intent(in) :: T
integer, intent(in) :: co, ip, el integer, intent(in) :: co, ce
end subroutine constitutive_thermal_setT end subroutine constitutive_thermal_setT
! == cleaned:end =================================================================================== ! == 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 real(pReal), intent(in) :: Delta_t
integer, intent(in) :: & integer, intent(in) :: ph, me
el, & !< element index in element loop logical :: converged_
ip, & !< integration point index in ip loop
co !< grain index in grain loop end function thermal_stress
logical :: broken
end function integrateThermalState
module function integrateDamageState(dt,co,ip,el) result(broken) module function integrateDamageState(dt,co,ip,el) result(broken)
real(pReal), intent(in) :: dt real(pReal), intent(in) :: dt
@ -283,12 +274,10 @@ module constitutive
dPhiDot_dPhi dPhiDot_dPhi
end subroutine constitutive_damage_getRateAndItsTangents end subroutine constitutive_damage_getRateAndItsTangents
module subroutine constitutive_thermal_getRate(TDot, T,ip,el) module subroutine constitutive_thermal_getRate(TDot, ip,el)
integer, intent(in) :: & integer, intent(in) :: &
ip, & !< integration point number ip, & !< integration point number
el !< element number el !< element number
real(pReal), intent(in) :: &
T
real(pReal), intent(out) :: & real(pReal), intent(out) :: &
TDot TDot
end subroutine constitutive_thermal_getRate end subroutine constitutive_thermal_getRate
@ -394,18 +383,19 @@ module constitutive
converged, & converged, &
crystallite_init, & crystallite_init, &
crystallite_stress, & crystallite_stress, &
thermal_stress, &
constitutive_mech_dPdF, & constitutive_mech_dPdF, &
crystallite_orientations, & crystallite_orientations, &
crystallite_push33ToRef, & crystallite_push33ToRef, &
constitutive_restartWrite, & constitutive_restartWrite, &
constitutive_restartRead, & constitutive_restartRead, &
integrateThermalState, &
integrateDamageState, & integrateDamageState, &
constitutive_thermal_setT, & constitutive_thermal_setT, &
constitutive_mech_getP, & constitutive_mech_getP, &
constitutive_mech_setF, & constitutive_mech_setF, &
constitutive_mech_getF, & constitutive_mech_getF, &
constitutive_initializeRestorationPoints, & constitutive_initializeRestorationPoints, &
constitutive_thermal_initializeRestorationPoints, &
constitutive_windForward, & constitutive_windForward, &
KINEMATICS_UNDEFINED_ID ,& KINEMATICS_UNDEFINED_ID ,&
KINEMATICS_CLEAVAGE_OPENING_ID, & KINEMATICS_CLEAVAGE_OPENING_ID, &
@ -553,7 +543,6 @@ subroutine constitutive_restore(ip,el,includeL)
enddo enddo
call mech_restore(ip,el,includeL) call mech_restore(ip,el,includeL)
call thermal_restore(ip,el)
end subroutine constitutive_restore end subroutine constitutive_restore
@ -720,7 +709,6 @@ subroutine constitutive_initializeRestorationPoints(ip,el)
me = material_phaseMemberAt(co,ip,el) me = material_phaseMemberAt(co,ip,el)
call mech_initializeRestorationPoints(ph,me) call mech_initializeRestorationPoints(ph,me)
call thermal_initializeRestorationPoints(ph,me)
do so = 1, size(damageState(ph)%p) do so = 1, size(damageState(ph)%p)
damageState(ph)%p(so)%partitionedState0(:,me) = damageState(ph)%p(so)%state0(:,me) damageState(ph)%p(so)%partitionedState0(:,me) = damageState(ph)%p(so)%state0(:,me)
@ -750,7 +738,6 @@ subroutine constitutive_windForward(ip,el)
me = material_phaseMemberAt(co,ip,el) me = material_phaseMemberAt(co,ip,el)
call mech_windForward(ph,me) call mech_windForward(ph,me)
call thermal_windForward(ph,me)
do so = 1, phase_Nsources(material_phaseAt(co,el)) 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)%partitionedState0(:,me) = damageState(ph)%p(so)%state(:,me)

64
src/constitutive_mech.f90
View File

@ -599,20 +599,18 @@ module subroutine constitutive_plastic_dependentState(co, ip, el)
el !< element el !< element
integer :: & integer :: &
ho, & !< homogenization ph, &
tme, & !< thermal member position
instance, me instance, me
ho = material_homogenizationAt(el) ph = material_phaseAt(co,el)
tme = material_homogenizationMemberAt(ip,el) me = material_phasememberAt(co,ip,el)
me = material_phasememberAt(co,ip,el) instance = phase_plasticityInstance(ph)
instance = phase_plasticityInstance(material_phaseAt(co,el))
plasticityType: select case (phase_plasticity(material_phaseAt(co,el))) plasticityType: select case (phase_plasticity(material_phaseAt(co,el)))
case (PLASTICITY_DISLOTWIN_ID) plasticityType case (PLASTICITY_DISLOTWIN_ID) plasticityType
call plastic_dislotwin_dependentState(temperature(ho)%p(tme),instance,me) call plastic_dislotwin_dependentState(thermal_T(ph,me),instance,me)
case (PLASTICITY_DISLOTUNGSTEN_ID) plasticityType case (PLASTICITY_DISLOTUNGSTEN_ID) plasticityType
call plastic_dislotungsten_dependentState(instance,me) call plastic_dislotungsten_dependentState(instance,me)
@ -650,17 +648,13 @@ subroutine constitutive_plastic_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, &
real(pReal), dimension(3,3) :: & real(pReal), dimension(3,3) :: &
Mp !< Mandel stress work conjugate with Lp Mp !< Mandel stress work conjugate with Lp
integer :: & integer :: &
ho, & !< homogenization i, j, instance, me, ph
tme !< thermal member position
integer :: &
i, j, instance, me
ho = material_homogenizationAt(el)
tme = material_homogenizationMemberAt(ip,el)
Mp = matmul(matmul(transpose(Fi),Fi),S) Mp = matmul(matmul(transpose(Fi),Fi),S)
me = material_phasememberAt(co,ip,el) me = material_phasememberAt(co,ip,el)
instance = phase_plasticityInstance(material_phaseAt(co,el)) ph = material_phaseAt(co,el)
instance = phase_plasticityInstance(ph)
plasticityType: select case (phase_plasticity(material_phaseAt(co,el))) plasticityType: select case (phase_plasticity(material_phaseAt(co,el)))
@ -678,13 +672,13 @@ subroutine constitutive_plastic_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, &
call plastic_kinehardening_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,me) call plastic_kinehardening_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,me)
case (PLASTICITY_NONLOCAL_ID) plasticityType case (PLASTICITY_NONLOCAL_ID) plasticityType
call plastic_nonlocal_LpAndItsTangent(Lp,dLp_dMp,Mp, temperature(ho)%p(tme),instance,me,ip,el) call plastic_nonlocal_LpAndItsTangent(Lp,dLp_dMp,Mp, thermal_T(ph,me),instance,me,ip,el)
case (PLASTICITY_DISLOTWIN_ID) plasticityType case (PLASTICITY_DISLOTWIN_ID) plasticityType
call plastic_dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,temperature(ho)%p(tme),instance,me) call plastic_dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp, thermal_T(ph,me),instance,me)
case (PLASTICITY_DISLOTUNGSTEN_ID) plasticityType case (PLASTICITY_DISLOTUNGSTEN_ID) plasticityType
call plastic_dislotungsten_LpAndItsTangent(Lp,dLp_dMp,Mp,temperature(ho)%p(tme),instance,me) call plastic_dislotungsten_LpAndItsTangent(Lp,dLp_dMp,Mp, thermal_T(ph,me),instance,me)
end select plasticityType end select plasticityType
@ -700,52 +694,49 @@ end subroutine constitutive_plastic_LpAndItsTangents
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief contains the constitutive equation for calculating the rate of change of microstructure !> @brief contains the constitutive equation for calculating the rate of change of microstructure
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function mech_collectDotState(subdt,co,ip,el,ph,of) result(broken) function mech_collectDotState(subdt,co,ip,el,ph,me) result(broken)
integer, intent(in) :: & integer, intent(in) :: &
co, & !< component-ID of integration point co, & !< component-ID of integration point
ip, & !< integration point ip, & !< integration point
el, & !< element el, & !< element
ph, & ph, &
of me
real(pReal), intent(in) :: & real(pReal), intent(in) :: &
subdt !< timestep subdt !< timestep
real(pReal), dimension(3,3) :: & real(pReal), dimension(3,3) :: &
Mp Mp
integer :: & integer :: &
ho, & !< homogenization
tme, & !< thermal member position
instance instance
logical :: broken logical :: broken
ho = material_homogenizationAt(el)
tme = material_homogenizationMemberAt(ip,el)
instance = phase_plasticityInstance(ph) instance = phase_plasticityInstance(ph)
Mp = matmul(matmul(transpose(constitutive_mech_Fi(ph)%data(1:3,1:3,of)),& Mp = matmul(matmul(transpose(constitutive_mech_Fi(ph)%data(1:3,1:3,me)),&
constitutive_mech_Fi(ph)%data(1:3,1:3,of)),constitutive_mech_S(ph)%data(1:3,1:3,of)) constitutive_mech_Fi(ph)%data(1:3,1:3,me)),constitutive_mech_S(ph)%data(1:3,1:3,me))
plasticityType: select case (phase_plasticity(ph)) plasticityType: select case (phase_plasticity(ph))
case (PLASTICITY_ISOTROPIC_ID) plasticityType case (PLASTICITY_ISOTROPIC_ID) plasticityType
call plastic_isotropic_dotState(Mp,instance,of) call plastic_isotropic_dotState(Mp,instance,me)
case (PLASTICITY_PHENOPOWERLAW_ID) plasticityType case (PLASTICITY_PHENOPOWERLAW_ID) plasticityType
call plastic_phenopowerlaw_dotState(Mp,instance,of) call plastic_phenopowerlaw_dotState(Mp,instance,me)
case (PLASTICITY_KINEHARDENING_ID) plasticityType case (PLASTICITY_KINEHARDENING_ID) plasticityType
call plastic_kinehardening_dotState(Mp,instance,of) call plastic_kinehardening_dotState(Mp,instance,me)
case (PLASTICITY_DISLOTWIN_ID) plasticityType case (PLASTICITY_DISLOTWIN_ID) plasticityType
call plastic_dislotwin_dotState(Mp,temperature(ho)%p(tme),instance,of) call plastic_dislotwin_dotState(Mp,thermal_T(ph,me),instance,me)
case (PLASTICITY_DISLOTUNGSTEN_ID) plasticityType case (PLASTICITY_DISLOTUNGSTEN_ID) plasticityType
call plastic_disloTungsten_dotState(Mp,temperature(ho)%p(tme),instance,of) call plastic_disloTungsten_dotState(Mp,thermal_T(ph,me),instance,me)
case (PLASTICITY_NONLOCAL_ID) plasticityType case (PLASTICITY_NONLOCAL_ID) plasticityType
call plastic_nonlocal_dotState(Mp,temperature(ho)%p(tme),subdt,instance,of,ip,el) call plastic_nonlocal_dotState(Mp,thermal_T(ph,me),subdt,instance,me,ip,el)
end select plasticityType end select plasticityType
broken = any(IEEE_is_NaN(plasticState(ph)%dotState(:,of))) broken = any(IEEE_is_NaN(plasticState(ph)%dotState(:,me)))
end function mech_collectDotState end function mech_collectDotState
@ -1633,9 +1624,6 @@ module function crystallite_stress(dt,co,ip,el) result(converged_)
do so = 1, phase_Nsources(ph) do so = 1, phase_Nsources(ph)
damageState(ph)%p(so)%subState0(:,me) = damageState(ph)%p(so)%state(:,me) damageState(ph)%p(so)%subState0(:,me) = damageState(ph)%p(so)%state(:,me)
enddo enddo
do so = 1, thermal_Nsources(ph)
thermalState(ph)%p(so)%subState0(:,me) = thermalState(ph)%p(so)%state(:,me)
enddo
endif endif
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! cut back (reduced time and restore) ! cut back (reduced time and restore)
@ -1652,9 +1640,6 @@ module function crystallite_stress(dt,co,ip,el) result(converged_)
do so = 1, phase_Nsources(ph) do so = 1, phase_Nsources(ph)
damageState(ph)%p(so)%state(:,me) = damageState(ph)%p(so)%subState0(:,me) damageState(ph)%p(so)%state(:,me) = damageState(ph)%p(so)%subState0(:,me)
enddo 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) todo = subStep > num%subStepMinCryst ! still on track or already done (beyond repair)
endif endif
@ -1668,7 +1653,6 @@ module function crystallite_stress(dt,co,ip,el) result(converged_)
constitutive_mech_Fp(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) 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. integrateDamageState(subStep * dt,co,ip,el)
converged_ = converged_ .and. .not. integrateThermalState(subStep * dt,co,ip,el)
endif endif
enddo cutbackLooping enddo cutbackLooping

79
src/constitutive_thermal.f90
View File

@ -86,7 +86,7 @@ module subroutine thermal_init(phases)
Nconstituents = count(material_phaseAt == ph) * discretization_nIPs 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) phase => phases%get(ph)
if(phase%contains('thermal')) then if(phase%contains('thermal')) then
thermal => phase%get('thermal') thermal => phase%get('thermal')
@ -127,13 +127,11 @@ end subroutine thermal_init
!---------------------------------------------------------------------------------------------- !----------------------------------------------------------------------------------------------
!< @brief calculates thermal dissipation rate !< @brief calculates thermal dissipation rate
!---------------------------------------------------------------------------------------------- !----------------------------------------------------------------------------------------------
module subroutine constitutive_thermal_getRate(TDot, T, ip, el) module subroutine constitutive_thermal_getRate(TDot, ip, el)
integer, intent(in) :: & integer, intent(in) :: &
ip, & !< integration point number ip, & !< integration point number
el !< element number el !< element number
real(pReal), intent(in) :: &
T !< plastic velocity gradient
real(pReal), intent(out) :: & real(pReal), intent(out) :: &
TDot TDot
@ -197,30 +195,37 @@ function constitutive_thermal_collectDotState(ph,me) result(broken)
end function constitutive_thermal_collectDotState 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 !> @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 real(pReal), intent(in) :: Delta_t
integer, intent(in) :: & integer, intent(in) :: ph, me
el, & !< element index in element loop
ip, & !< integration point index in ip loop
co !< grain index in grain loop
logical :: & logical :: &
broken broken
integer :: & integer :: &
ph, &
me, &
so, & so, &
sizeDotState sizeDotState
ph = material_phaseAt(co,el)
me = material_phaseMemberAt(co,ip,el)
broken = constitutive_thermal_collectDotState(ph,me) broken = constitutive_thermal_collectDotState(ph,me)
if(broken) return if(broken) return
@ -233,7 +238,7 @@ module function integrateThermalState(Delta_t,co,ip,el) result(broken)
end function integrateThermalState end function integrateThermalState
module subroutine thermal_initializeRestorationPoints(ph,me) module subroutine constitutive_thermal_initializeRestorationPoints(ph,me)
integer, intent(in) :: ph, me integer, intent(in) :: ph, me
@ -244,24 +249,10 @@ module subroutine thermal_initializeRestorationPoints(ph,me)
thermalState(ph)%p(so)%partitionedState0(:,me) = thermalState(ph)%p(so)%state0(:,me) thermalState(ph)%p(so)%partitionedState0(:,me) = thermalState(ph)%p(so)%state0(:,me)
enddo 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() module subroutine thermal_forward()
integer :: ph, so integer :: ph, so
@ -276,26 +267,6 @@ module subroutine thermal_forward()
end 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) !< @brief Get temperature (for use by non-thermal physics)
!---------------------------------------------------------------------------------------------- !----------------------------------------------------------------------------------------------
@ -313,13 +284,13 @@ end function thermal_T
!---------------------------------------------------------------------------------------------- !----------------------------------------------------------------------------------------------
!< @brief Set temperature !< @brief Set temperature
!---------------------------------------------------------------------------------------------- !----------------------------------------------------------------------------------------------
module subroutine constitutive_thermal_setT(T,co,ip,el) module subroutine constitutive_thermal_setT(T,co,ce)
real(pReal), intent(in) :: T real(pReal), intent(in) :: T
integer, intent(in) :: co, ip, el integer, intent(in) :: ce, co
current(material_phaseAt(co,el))%T(material_phaseMemberAt(co,ip,el)) = T current(material_phaseAt2(co,ce))%T(material_phaseMemberAt2(co,ce)) = T
end subroutine constitutive_thermal_setT end subroutine constitutive_thermal_setT

55
src/grid/grid_thermal_spectral.f90
View File

@ -16,6 +16,7 @@ module grid_thermal_spectral
use spectral_utilities use spectral_utilities
use discretization_grid use discretization_grid
use thermal_conduction use thermal_conduction
use homogenization
use YAML_types use YAML_types
use config use config
use material use material
@ -64,7 +65,7 @@ contains
subroutine grid_thermal_spectral_init subroutine grid_thermal_spectral_init
PetscInt, dimension(0:worldsize-1) :: localK PetscInt, dimension(0:worldsize-1) :: localK
integer :: i, j, k, cell integer :: i, j, k, ce
DM :: thermal_grid DM :: thermal_grid
PetscScalar, dimension(:,:,:), pointer :: x_scal PetscScalar, dimension(:,:,:), pointer :: x_scal
PetscErrorCode :: ierr PetscErrorCode :: ierr
@ -128,10 +129,10 @@ subroutine grid_thermal_spectral_init
allocate(T_current(grid(1),grid(2),grid3), source=0.0_pReal) allocate(T_current(grid(1),grid(2),grid3), source=0.0_pReal)
allocate(T_lastInc(grid(1),grid(2),grid3), source=0.0_pReal) allocate(T_lastInc(grid(1),grid(2),grid3), source=0.0_pReal)
allocate(T_stagInc(grid(1),grid(2),grid3), source=0.0_pReal) allocate(T_stagInc(grid(1),grid(2),grid3), source=0.0_pReal)
cell = 0 ce = 0
do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1)
cell = cell + 1 ce = ce + 1
T_current(i,j,k) = temperature(material_homogenizationAt(cell))%p(material_homogenizationMemberAt(1,cell)) T_current(i,j,k) = temperature(material_homogenizationAt(ce))%p(material_homogenizationMemberAt(1,ce))
T_lastInc(i,j,k) = T_current(i,j,k) T_lastInc(i,j,k) = T_current(i,j,k)
T_stagInc(i,j,k) = T_current(i,j,k) T_stagInc(i,j,k) = T_current(i,j,k)
enddo; enddo; enddo enddo; enddo; enddo
@ -151,7 +152,7 @@ function grid_thermal_spectral_solution(timeinc) result(solution)
real(pReal), intent(in) :: & real(pReal), intent(in) :: &
timeinc !< increment in time for current solution timeinc !< increment in time for current solution
integer :: i, j, k, cell integer :: i, j, k, ce
type(tSolutionState) :: solution type(tSolutionState) :: solution
PetscInt :: devNull PetscInt :: devNull
PetscReal :: T_min, T_max, stagNorm, solnNorm PetscReal :: T_min, T_max, stagNorm, solnNorm
@ -184,12 +185,13 @@ function grid_thermal_spectral_solution(timeinc) result(solution)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! updating thermal state ! updating thermal state
cell = 0 ce = 0
do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1)
cell = cell + 1 ce = ce + 1
call thermal_conduction_putTemperatureAndItsRate(T_current(i,j,k), & call thermal_conduction_putTemperatureAndItsRate(T_current(i,j,k), &
(T_current(i,j,k)-T_lastInc(i,j,k))/params%timeinc, & (T_current(i,j,k)-T_lastInc(i,j,k))/params%timeinc, &
1,cell) 1,ce)
homogenization_T(ce) = T_current(i,j,k)
enddo; enddo; enddo enddo; enddo; enddo
call VecMin(solution_vec,devNull,T_min,ierr); CHKERRQ(ierr) call VecMin(solution_vec,devNull,T_min,ierr); CHKERRQ(ierr)
@ -209,7 +211,7 @@ end function grid_thermal_spectral_solution
subroutine grid_thermal_spectral_forward(cutBack) subroutine grid_thermal_spectral_forward(cutBack)
logical, intent(in) :: cutBack logical, intent(in) :: cutBack
integer :: i, j, k, cell integer :: i, j, k, ce
DM :: dm_local DM :: dm_local
PetscScalar, dimension(:,:,:), pointer :: x_scal PetscScalar, dimension(:,:,:), pointer :: x_scal
PetscErrorCode :: ierr PetscErrorCode :: ierr
@ -220,17 +222,18 @@ subroutine grid_thermal_spectral_forward(cutBack)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! reverting thermal field state ! reverting thermal field state
cell = 0 ce = 0
call SNESGetDM(thermal_snes,dm_local,ierr); CHKERRQ(ierr) call SNESGetDM(thermal_snes,dm_local,ierr); CHKERRQ(ierr)
call DMDAVecGetArrayF90(dm_local,solution_vec,x_scal,ierr); CHKERRQ(ierr) !< get the data out of PETSc to work with call DMDAVecGetArrayF90(dm_local,solution_vec,x_scal,ierr); CHKERRQ(ierr) !< get the data out of PETSc to work with
x_scal(xstart:xend,ystart:yend,zstart:zend) = T_current x_scal(xstart:xend,ystart:yend,zstart:zend) = T_current
call DMDAVecRestoreArrayF90(dm_local,solution_vec,x_scal,ierr); CHKERRQ(ierr) call DMDAVecRestoreArrayF90(dm_local,solution_vec,x_scal,ierr); CHKERRQ(ierr)
do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1)
cell = cell + 1 ce = ce + 1
call thermal_conduction_putTemperatureAndItsRate(T_current(i,j,k), & call thermal_conduction_putTemperatureAndItsRate(T_current(i,j,k), &
(T_current(i,j,k) - & (T_current(i,j,k) - &
T_lastInc(i,j,k))/params%timeinc, & T_lastInc(i,j,k))/params%timeinc, &
1,cell) 1,ce)
homogenization_T(ce) = T_current(i,j,k)
enddo; enddo; enddo enddo; enddo; enddo
else else
T_lastInc = T_current T_lastInc = T_current
@ -255,7 +258,7 @@ subroutine formResidual(in,x_scal,f_scal,dummy,ierr)
f_scal f_scal
PetscObject :: dummy PetscObject :: dummy
PetscErrorCode :: ierr PetscErrorCode :: ierr
integer :: i, j, k, cell integer :: i, j, k, ce
real(pReal) :: Tdot real(pReal) :: Tdot
T_current = x_scal T_current = x_scal
@ -266,22 +269,22 @@ subroutine formResidual(in,x_scal,f_scal,dummy,ierr)
call utilities_FFTscalarForward call utilities_FFTscalarForward
call utilities_fourierScalarGradient !< calculate gradient of temperature field call utilities_fourierScalarGradient !< calculate gradient of temperature field
call utilities_FFTvectorBackward call utilities_FFTvectorBackward
cell = 0 ce = 0
do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1)
cell = cell + 1 ce = ce + 1
vectorField_real(1:3,i,j,k) = matmul(thermal_conduction_getConductivity(1,cell) - K_ref, & vectorField_real(1:3,i,j,k) = matmul(thermal_conduction_getConductivity(1,ce) - K_ref, &
vectorField_real(1:3,i,j,k)) vectorField_real(1:3,i,j,k))
enddo; enddo; enddo enddo; enddo; enddo
call utilities_FFTvectorForward call utilities_FFTvectorForward
call utilities_fourierVectorDivergence !< calculate temperature divergence in fourier field call utilities_fourierVectorDivergence !< calculate temperature divergence in fourier field
call utilities_FFTscalarBackward call utilities_FFTscalarBackward
cell = 0 ce = 0
do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1)
cell = cell + 1 ce = ce + 1
call thermal_conduction_getSource(Tdot, T_current(i,j,k), 1, cell) call thermal_conduction_getSource(Tdot, 1,ce)
scalarField_real(i,j,k) = params%timeinc*(scalarField_real(i,j,k) + Tdot) & scalarField_real(i,j,k) = params%timeinc*(scalarField_real(i,j,k) + Tdot) &
+ thermal_conduction_getMassDensity (1,cell)* & + thermal_conduction_getMassDensity (1,ce)* &
thermal_conduction_getSpecificHeat(1,cell)*(T_lastInc(i,j,k) - & thermal_conduction_getSpecificHeat(1,ce)*(T_lastInc(i,j,k) - &
T_current(i,j,k))& T_current(i,j,k))&
+ mu_ref*T_current(i,j,k) + mu_ref*T_current(i,j,k)
enddo; enddo; enddo enddo; enddo; enddo
@ -304,15 +307,15 @@ end subroutine formResidual
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine updateReference subroutine updateReference
integer :: i,j,k,cell,ierr integer :: i,j,k,ce,ierr
cell = 0 ce = 0
K_ref = 0.0_pReal K_ref = 0.0_pReal
mu_ref = 0.0_pReal mu_ref = 0.0_pReal
do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1)
cell = cell + 1 ce = ce + 1
K_ref = K_ref + thermal_conduction_getConductivity(1,cell) K_ref = K_ref + thermal_conduction_getConductivity(1,ce)
mu_ref = mu_ref + thermal_conduction_getMassDensity(1,cell)* thermal_conduction_getSpecificHeat(1,cell) mu_ref = mu_ref + thermal_conduction_getMassDensity(1,ce)* thermal_conduction_getSpecificHeat(1,ce)
enddo; enddo; enddo enddo; enddo; enddo
K_ref = K_ref*wgt K_ref = K_ref*wgt
call MPI_Allreduce(MPI_IN_PLACE,K_ref,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) call MPI_Allreduce(MPI_IN_PLACE,K_ref,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr)

62
src/homogenization.f90
View File

@ -28,7 +28,8 @@ module homogenization
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! General variables for the homogenization at a material point ! General variables for the homogenization at a material point
real(pReal), dimension(:), allocatable, public :: & real(pReal), dimension(:), allocatable, public :: &
homogenization_T homogenization_T, &
homogenization_dot_T
real(pReal), dimension(:,:,:), allocatable, public :: & real(pReal), dimension(:,:,:), allocatable, public :: &
homogenization_F0, & !< def grad of IP at start of FE increment 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 homogenization_F !< def grad of IP to be reached at end of FE increment
@ -69,13 +70,15 @@ module homogenization
el !< element number el !< element number
end subroutine mech_partition end subroutine mech_partition
module subroutine thermal_partition(T,ip,el) module subroutine thermal_partition(T,ce)
real(pReal), intent(in) :: T real(pReal), intent(in) :: T
integer, intent(in) :: & integer, intent(in) :: ce
ip, & !< integration point
el !< element number
end subroutine thermal_partition end subroutine thermal_partition
module subroutine thermal_homogenize(ip,el)
integer, intent(in) :: ip,el
end subroutine thermal_homogenize
module subroutine mech_homogenize(dt,ip,el) module subroutine mech_homogenize(dt,ip,el)
real(pReal), intent(in) :: dt real(pReal), intent(in) :: dt
integer, intent(in) :: & integer, intent(in) :: &
@ -161,7 +164,7 @@ subroutine materialpoint_stressAndItsTangent(dt,FEsolving_execIP,FEsolving_execE
NiterationMPstate, & NiterationMPstate, &
ip, & !< integration point number ip, & !< integration point number
el, & !< element number el, & !< element number
myNgrains, co, ce, ho, me myNgrains, co, ce, ho, me, ph
real(pReal) :: & real(pReal) :: &
subFrac, & subFrac, &
subStep subStep
@ -170,8 +173,8 @@ subroutine materialpoint_stressAndItsTangent(dt,FEsolving_execIP,FEsolving_execE
logical, dimension(2) :: & logical, dimension(2) :: &
doneAndHappy doneAndHappy
!$OMP PARALLEL
!$OMP PARALLEL DO PRIVATE(ce,me,ho,myNgrains,NiterationMPstate,subFrac,converged,subStep,doneAndHappy) !$OMP DO PRIVATE(ce,me,ho,myNgrains,NiterationMPstate,subFrac,converged,subStep,doneAndHappy)
do el = FEsolving_execElem(1),FEsolving_execElem(2) do el = FEsolving_execElem(1),FEsolving_execElem(2)
ho = material_homogenizationAt(el) ho = material_homogenizationAt(el)
myNgrains = homogenization_Nconstituents(ho) myNgrains = homogenization_Nconstituents(ho)
@ -221,9 +224,8 @@ subroutine materialpoint_stressAndItsTangent(dt,FEsolving_execIP,FEsolving_execE
if (subStep > num%subStepMinHomog) doneAndHappy = [.false.,.true.] if (subStep > num%subStepMinHomog) doneAndHappy = [.false.,.true.]
NiterationMPstate = 0 NiterationMPstate = 0
convergenceLooping: do while (.not. terminallyIll & convergenceLooping: do while (.not. (terminallyIll .or. doneAndHappy(1)) &
.and. .not. doneAndHappy(1) & .and. NiterationMPstate < num%nMPstate)
.and. NiterationMPstate < num%nMPstate)
NiterationMPstate = NiterationMPstate + 1 NiterationMPstate = NiterationMPstate + 1
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -231,10 +233,9 @@ subroutine materialpoint_stressAndItsTangent(dt,FEsolving_execIP,FEsolving_execE
if (.not. doneAndHappy(1)) then if (.not. doneAndHappy(1)) then
ce = (el-1)*discretization_nIPs + ip ce = (el-1)*discretization_nIPs + ip
call mech_partition(homogenization_F0(1:3,1:3,ce) & call mech_partition( homogenization_F0(1:3,1:3,ce) &
+ (homogenization_F(1:3,1:3,ce)-homogenization_F0(1:3,1:3,ce))& + (homogenization_F(1:3,1:3,ce)-homogenization_F0(1:3,1:3,ce))*(subStep+subFrac), &
*(subStep+subFrac), & ip,el)
ip,el)
converged = .true. converged = .true.
do co = 1, myNgrains do co = 1, myNgrains
converged = converged .and. crystallite_stress(dt*subStep,co,ip,el) converged = converged .and. crystallite_stress(dt*subStep,co,ip,el)
@ -257,24 +258,45 @@ subroutine materialpoint_stressAndItsTangent(dt,FEsolving_execIP,FEsolving_execE
enddo cutBackLooping enddo cutBackLooping
enddo enddo
enddo enddo
!$OMP END PARALLEL DO !$OMP END DO
if (.not. terminallyIll ) then if (.not. terminallyIll ) then
!$OMP PARALLEL DO PRIVATE(ho,myNgrains) !$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 thermal_partition(homogenization_T(ce),ce)
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
call thermal_homogenize(ip,el)
enddo
enddo
enddo
!$OMP END DO
!$OMP DO PRIVATE(ho)
elementLooping3: do el = FEsolving_execElem(1),FEsolving_execElem(2) elementLooping3: do el = FEsolving_execElem(1),FEsolving_execElem(2)
ho = material_homogenizationAt(el) ho = material_homogenizationAt(el)
myNgrains = homogenization_Nconstituents(ho)
IpLooping3: do ip = FEsolving_execIP(1),FEsolving_execIP(2) 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) call crystallite_orientations(co,ip,el)
enddo enddo
call mech_homogenize(dt,ip,el) call mech_homogenize(dt,ip,el)
enddo IpLooping3 enddo IpLooping3
enddo elementLooping3 enddo elementLooping3
!$OMP END PARALLEL DO !$OMP END DO
else else
print'(/,a,/)', ' << HOMOG >> Material Point terminally ill' print'(/,a,/)', ' << HOMOG >> Material Point terminally ill'
endif endif
!$OMP END PARALLEL
end subroutine materialpoint_stressAndItsTangent end subroutine materialpoint_stressAndItsTangent

27
src/homogenization_thermal.f90
View File

@ -15,25 +15,38 @@ module subroutine thermal_init()
print'(/,a)', ' <<<+- homogenization_thermal init -+>>>' print'(/,a)', ' <<<+- homogenization_thermal init -+>>>'
allocate(homogenization_T(discretization_nIPs*discretization_Nelems)) allocate(homogenization_T(discretization_nIPs*discretization_Nelems))
allocate(homogenization_dot_T(discretization_nIPs*discretization_Nelems))
end subroutine thermal_init end subroutine thermal_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Partition T onto the individual constituents. !> @brief Partition temperature onto the individual constituents.
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module subroutine thermal_partition(T,ip,el) module subroutine thermal_partition(T,ce)
real(pReal), intent(in) :: T real(pReal), intent(in) :: T
integer, intent(in) :: & integer, intent(in) :: ce
ip, & !< integration point
el !< element number
integer :: co
call constitutive_thermal_setT(T,1,ip,el) do co = 1, homogenization_Nconstituents(material_homogenizationAt2(ce))
call constitutive_thermal_setT(T,co,ce)
enddo
end subroutine thermal_partition end subroutine thermal_partition
!--------------------------------------------------------------------------------------------------
!> @brief Homogenize temperature rates
!--------------------------------------------------------------------------------------------------
module subroutine thermal_homogenize(ip,el)
integer, intent(in) :: ip,el
call constitutive_thermal_getRate(homogenization_dot_T((el-1)*discretization_nIPs+ip), ip,el)
end subroutine thermal_homogenize
end submodule homogenization_thermal end submodule homogenization_thermal

50
src/material.f90
View File

@ -51,11 +51,15 @@ module material
thermal_initialT !< initial temperature per each homogenization thermal_initialT !< initial temperature per each homogenization
integer, dimension(:), allocatable, public, protected :: & ! (elem) integer, dimension(:), allocatable, public, protected :: & ! (elem)
material_homogenizationAt !< homogenization ID of each element material_homogenizationAt, & !< homogenization ID of each element
material_homogenizationAt2, & !< per cell
material_homogenizationMemberAt2 !< cell
integer, dimension(:,:), allocatable, public, protected :: & ! (ip,elem) integer, dimension(:,:), allocatable, public, protected :: & ! (ip,elem)
material_homogenizationMemberAt !< position of the element within its homogenization instance material_homogenizationMemberAt !< position of the element within its homogenization instance
integer, dimension(:,:), allocatable, public, protected :: & ! (constituent,elem) integer, dimension(:,:), allocatable, public, protected :: & ! (constituent,elem)
material_phaseAt !< phase ID of each element material_phaseAt, & !< phase ID of each element
material_phaseAt2, & !< per constituent,cell
material_phaseMemberAt2 !< per constituent, cell
integer, dimension(:,:,:), allocatable, public, protected :: & ! (constituent,IP,elem) integer, dimension(:,:,:), allocatable, public, protected :: & ! (constituent,IP,elem)
material_phaseMemberAt !< position of the element within its phase instance material_phaseMemberAt !< position of the element within its phase instance
@ -215,8 +219,8 @@ subroutine material_parseMaterial
real(pReal) :: & real(pReal) :: &
frac frac
integer :: & integer :: &
e, i, c, & el, ip, co, &
h h, ce
materials => config_material%get('material') materials => config_material%get('material')
phases => config_material%get('phase') phases => config_material%get('phase')
@ -241,29 +245,41 @@ subroutine material_parseMaterial
allocate(material_phaseAt(homogenization_maxNconstituents,discretization_Nelems),source=0) allocate(material_phaseAt(homogenization_maxNconstituents,discretization_Nelems),source=0)
allocate(material_phaseMemberAt(homogenization_maxNconstituents,discretization_nIPs,discretization_Nelems),source=0) allocate(material_phaseMemberAt(homogenization_maxNconstituents,discretization_nIPs,discretization_Nelems),source=0)
allocate(material_homogenizationAt2(discretization_nIPs*discretization_Nelems),source=0)
allocate(material_homogenizationMemberAt2(discretization_nIPs*discretization_Nelems),source=0)
allocate(material_phaseAt2(homogenization_maxNconstituents,discretization_nIPs*discretization_Nelems),source=0)
allocate(material_phaseMemberAt2(homogenization_maxNconstituents,discretization_nIPs*discretization_Nelems),source=0)
allocate(material_orientation0(homogenization_maxNconstituents,discretization_nIPs,discretization_Nelems)) allocate(material_orientation0(homogenization_maxNconstituents,discretization_nIPs,discretization_Nelems))
do e = 1, discretization_Nelems do el = 1, discretization_Nelems
material => materials%get(discretization_materialAt(e)) material => materials%get(discretization_materialAt(el))
constituents => material%get('constituents') constituents => material%get('constituents')
material_homogenizationAt(e) = homogenizations%getIndex(material%get_asString('homogenization')) material_homogenizationAt(el) = homogenizations%getIndex(material%get_asString('homogenization'))
do i = 1, discretization_nIPs do ip = 1, discretization_nIPs
counterHomogenization(material_homogenizationAt(e)) = counterHomogenization(material_homogenizationAt(e)) + 1 ce = (el-1)*discretization_nIPs + ip
material_homogenizationMemberAt(i,e) = counterHomogenization(material_homogenizationAt(e)) counterHomogenization(material_homogenizationAt(el)) = counterHomogenization(material_homogenizationAt(el)) + 1
material_homogenizationMemberAt(ip,el) = counterHomogenization(material_homogenizationAt(el))
material_homogenizationAt2(ce) = material_homogenizationAt(el)
material_homogenizationMemberAt2(ce) = material_homogenizationMemberAt(ip,el)
enddo enddo
frac = 0.0_pReal frac = 0.0_pReal
do c = 1, constituents%length do co = 1, constituents%length
constituent => constituents%get(c) constituent => constituents%get(co)
frac = frac + constituent%get_asFloat('fraction') frac = frac + constituent%get_asFloat('fraction')
material_phaseAt(c,e) = phases%getIndex(constituent%get_asString('phase')) material_phaseAt(co,el) = phases%getIndex(constituent%get_asString('phase'))
do i = 1, discretization_nIPs do ip = 1, discretization_nIPs
counterPhase(material_phaseAt(c,e)) = counterPhase(material_phaseAt(c,e)) + 1 ce = (el-1)*discretization_nIPs + ip
material_phaseMemberAt(c,i,e) = counterPhase(material_phaseAt(c,e)) counterPhase(material_phaseAt(co,el)) = counterPhase(material_phaseAt(co,el)) + 1
material_phaseMemberAt(co,ip,el) = counterPhase(material_phaseAt(co,el))
call material_orientation0(c,i,e)%fromQuaternion(constituent%get_asFloats('O',requiredSize=4)) ! should be done in crystallite material_phaseAt2(co,ce) = material_phaseAt(co,el)
material_phaseMemberAt2(co,ce) = material_phaseMemberAt(co,ip,el)
call material_orientation0(co,ip,el)%fromQuaternion(constituent%get_asFloats('O',requiredSize=4)) ! should be done in crystallite
enddo enddo
enddo enddo

6
src/thermal_conduction.f90
View File

@ -91,13 +91,11 @@ end subroutine thermal_conduction_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief return heat generation rate !> @brief return heat generation rate
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine thermal_conduction_getSource(Tdot, T,ip,el) subroutine thermal_conduction_getSource(Tdot, ip,el)
integer, intent(in) :: & integer, intent(in) :: &
ip, & !< integration point number ip, & !< integration point number
el !< element number el !< element number
real(pReal), intent(in) :: &
T
real(pReal), intent(out) :: & real(pReal), intent(out) :: &
Tdot Tdot
@ -105,7 +103,7 @@ subroutine thermal_conduction_getSource(Tdot, T,ip,el)
homog homog
homog = material_homogenizationAt(el) homog = material_homogenizationAt(el)
call constitutive_thermal_getRate(TDot, T,ip,el) call constitutive_thermal_getRate(TDot, ip,el)
Tdot = Tdot/real(homogenization_Nconstituents(homog),pReal) Tdot = Tdot/real(homogenization_Nconstituents(homog),pReal)