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

separate state for thermal

This commit is contained in:
Martin Diehl 2021-01-07 22:15:18 +01:00
parent 2c64ae34e4
commit 27f4e4ce2a
9 changed files with 299 additions and 180 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

146
src/constitutive.f90
View File

@ -90,6 +90,7 @@ module constitutive
phase_kinematics !< active kinematic mechanisms of each phase phase_kinematics !< active kinematic mechanisms of each phase
integer, dimension(:), allocatable, public :: & !< ToDo: should be protected (bug in Intel compiler) integer, dimension(:), allocatable, public :: & !< ToDo: should be protected (bug in Intel compiler)
thermal_Nsources, &
phase_Nsources, & !< number of source mechanisms active in each phase phase_Nsources, & !< number of source mechanisms active in each phase
phase_Nkinematics, & !< number of kinematic mechanisms active in each phase phase_Nkinematics, & !< number of kinematic mechanisms active in each phase
phase_NstiffnessDegradations, & !< number of stiffness degradation mechanisms active in each phase phase_NstiffnessDegradations, & !< number of stiffness degradation mechanisms active in each phase
@ -233,6 +234,16 @@ module constitutive
! == cleaned:end =================================================================================== ! == cleaned:end ===================================================================================
module function integrateThermalState(dt,co,ip,el) result(broken)
real(pReal), intent(in) :: dt
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
module function crystallite_stress(dt,co,ip,el) result(converged_) module function crystallite_stress(dt,co,ip,el) result(converged_)
real(pReal), intent(in) :: dt real(pReal), intent(in) :: dt
integer, intent(in) :: co, ip, el integer, intent(in) :: co, ip, el
@ -665,31 +676,6 @@ function constitutive_damage_collectDotState(co,ip,el,ph,of) result(broken)
end function constitutive_damage_collectDotState end function constitutive_damage_collectDotState
!--------------------------------------------------------------------------------------------------
!> @brief contains the constitutive equation for calculating the rate of change of microstructure
!--------------------------------------------------------------------------------------------------
function constitutive_thermal_collectDotState(ph,me) result(broken)
integer, intent(in) :: ph, me
logical :: broken
integer :: i
broken = .false.
SourceLoop: do i = 1, phase_Nsources(ph)
if (phase_source(i,ph) == SOURCE_thermal_externalheat_ID) &
call source_thermal_externalheat_dotState(ph,me)
broken = broken .or. any(IEEE_is_NaN(sourceState(ph)%p(i)%dotState(:,me)))
enddo SourceLoop
end function constitutive_thermal_collectDotState
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief for constitutive models having an instantaneous change of state !> @brief for constitutive models having an instantaneous change of state
!> will return false if delta state is not needed/supported by the constitutive model !> will return false if delta state is not needed/supported by the constitutive model
@ -856,7 +842,7 @@ subroutine crystallite_init()
cMax, & !< maximum number of integration point components cMax, & !< maximum number of integration point components
iMax, & !< maximum number of integration points iMax, & !< maximum number of integration points
eMax !< maximum number of elements eMax !< maximum number of elements
class(tNode), pointer :: & class(tNode), pointer :: &
num_crystallite, & num_crystallite, &
@ -914,6 +900,9 @@ subroutine crystallite_init()
do so = 1, phase_Nsources(ph) do so = 1, phase_Nsources(ph)
allocate(sourceState(ph)%p(so)%subState0,source=sourceState(ph)%p(so)%state0) ! ToDo: hack allocate(sourceState(ph)%p(so)%subState0,source=sourceState(ph)%p(so)%state0) ! ToDo: hack
enddo enddo
do so = 1, thermal_Nsources(ph)
allocate(thermalState(ph)%p(so)%subState0,source=thermalState(ph)%p(so)%state0) ! ToDo: hack
enddo
enddo enddo
print'(a42,1x,i10)', ' # of elements: ', eMax print'(a42,1x,i10)', ' # of elements: ', eMax
@ -1144,111 +1133,6 @@ function integrateSourceState(dt,co,ip,el) result(broken)
end function integrateSourceState end function integrateSourceState
!--------------------------------------------------------------------------------------------------
!> @brief integrate stress, state with adaptive 1st order explicit Euler method
!> using Fixed Point Iteration to adapt the stepsize
!--------------------------------------------------------------------------------------------------
function integrateThermalState(dt,co,ip,el) result(broken)
real(pReal), intent(in) :: dt
integer, intent(in) :: &
el, & !< element index in element loop
ip, & !< integration point index in ip loop
co !< grain index in grain loop
integer :: &
NiterationState, & !< number of iterations in state loop
ph, &
me, &
so
integer, dimension(maxval(phase_Nsources)) :: &
size_so
real(pReal) :: &
zeta
real(pReal), dimension(constitutive_source_maxSizeDotState) :: &
r ! state residuum
real(pReal), dimension(constitutive_source_maxSizeDotState,2,maxval(phase_Nsources)) :: source_dotState
logical :: &
broken, converged_
ph = material_phaseAt(co,el)
me = material_phaseMemberAt(co,ip,el)
converged_ = .true.
broken = constitutive_thermal_collectDotState(ph,me)
if(broken) return
do so = 1, phase_Nsources(ph)
size_so(so) = thermalState(ph)%p(so)%sizeDotState
thermalState(ph)%p(so)%state(1:size_so(so),me) = thermalState(ph)%p(so)%subState0(1:size_so(so),me) &
+ thermalState(ph)%p(so)%dotState (1:size_so(so),me) * dt
source_dotState(1:size_so(so),2,so) = 0.0_pReal
enddo
iteration: do NiterationState = 1, num%nState
do so = 1, phase_Nsources(ph)
if(nIterationState > 1) source_dotState(1:size_so(so),2,so) = source_dotState(1:size_so(so),1,so)
source_dotState(1:size_so(so),1,so) = thermalState(ph)%p(so)%dotState(:,me)
enddo
broken = constitutive_thermal_collectDotState(ph,me)
broken = broken .or. constitutive_damage_collectDotState(co,ip,el,ph,me)
if(broken) exit iteration
do so = 1, phase_Nsources(ph)
zeta = damper(thermalState(ph)%p(so)%dotState(:,me), &
source_dotState(1:size_so(so),1,so),&
source_dotState(1:size_so(so),2,so))
thermalState(ph)%p(so)%dotState(:,me) = thermalState(ph)%p(so)%dotState(:,me) * zeta &
+ source_dotState(1:size_so(so),1,so)* (1.0_pReal - zeta)
r(1:size_so(so)) = thermalState(ph)%p(so)%state (1:size_so(so),me) &
- thermalState(ph)%p(so)%subState0(1:size_so(so),me) &
- thermalState(ph)%p(so)%dotState (1:size_so(so),me) * dt
thermalState(ph)%p(so)%state(1:size_so(so),me) = thermalState(ph)%p(so)%state(1:size_so(so),me) &
- r(1:size_so(so))
converged_ = converged_ .and. converged(r(1:size_so(so)), &
thermalState(ph)%p(so)%state(1:size_so(so),me), &
thermalState(ph)%p(so)%atol(1:size_so(so)))
enddo
if(converged_) then
broken = constitutive_damage_deltaState(mech_F_e(ph,me),co,ip,el,ph,me)
exit iteration
endif
enddo iteration
broken = broken .or. .not. converged_
contains
!--------------------------------------------------------------------------------------------------
!> @brief calculate the damping for correction of state and dot state
!--------------------------------------------------------------------------------------------------
real(pReal) pure function damper(current,previous,previous2)
real(pReal), dimension(:), intent(in) ::&
current, previous, previous2
real(pReal) :: dot_prod12, dot_prod22
dot_prod12 = dot_product(current - previous, previous - previous2)
dot_prod22 = dot_product(previous - previous2, previous - previous2)
if ((dot_product(current,previous) < 0.0_pReal .or. dot_prod12 < 0.0_pReal) .and. dot_prod22 > 0.0_pReal) then
damper = 0.75_pReal + 0.25_pReal * tanh(2.0_pReal + 4.0_pReal * dot_prod12 / dot_prod22)
else
damper = 1.0_pReal
endif
end function damper
end function integrateThermalState
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief determines whether a point is converged !> @brief determines whether a point is converged
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------

10
src/constitutive_mech.f90
View File

@ -1588,6 +1588,9 @@ module function crystallite_stress(dt,co,ip,el) result(converged_)
do so = 1, phase_Nsources(ph) do so = 1, phase_Nsources(ph)
sourceState(ph)%p(so)%subState0(:,me) = sourceState(ph)%p(so)%partitionedState0(:,me) sourceState(ph)%p(so)%subState0(:,me) = sourceState(ph)%p(so)%partitionedState0(:,me)
enddo enddo
do so = 1, thermal_Nsources(ph)
thermalState(ph)%p(so)%subState0(:,me) = thermalState(ph)%p(so)%partitionedState0(:,me)
enddo
subFp0 = constitutive_mech_partitionedFp0(ph)%data(1:3,1:3,me) subFp0 = constitutive_mech_partitionedFp0(ph)%data(1:3,1:3,me)
subFi0 = constitutive_mech_partitionedFi0(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) subF0 = constitutive_mech_partitionedF0(ph)%data(1:3,1:3,me)
@ -1616,6 +1619,9 @@ module function crystallite_stress(dt,co,ip,el) result(converged_)
do so = 1, phase_Nsources(ph) do so = 1, phase_Nsources(ph)
sourceState(ph)%p(so)%subState0(:,me) = sourceState(ph)%p(so)%state(:,me) sourceState(ph)%p(so)%subState0(:,me) = sourceState(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)
@ -1632,6 +1638,9 @@ module function crystallite_stress(dt,co,ip,el) result(converged_)
do so = 1, phase_Nsources(ph) do so = 1, phase_Nsources(ph)
sourceState(ph)%p(so)%state(:,me) = sourceState(ph)%p(so)%subState0(:,me) sourceState(ph)%p(so)%state(:,me) = sourceState(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
@ -1645,6 +1654,7 @@ 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. integrateSourceState(subStep * dt,co,ip,el) converged_ = converged_ .and. .not. integrateSourceState(subStep * dt,co,ip,el)
converged_ = converged_ .and. .not. integrateThermalState(subStep * dt,co,ip,el)
endif endif
enddo cutbackLooping enddo cutbackLooping

212
src/constitutive_thermal.f90
View File

@ -6,9 +6,12 @@ submodule(constitutive) constitutive_thermal
type :: tDataContainer type :: tDataContainer
real(pReal), dimension(:), allocatable :: T real(pReal), dimension(:), allocatable :: T
end type tDataContainer end type tDataContainer
integer(kind(SOURCE_undefined_ID)), dimension(:,:), allocatable :: &
thermal_source
type(tDataContainer), dimension(:), allocatable :: current type(tDataContainer), dimension(:), allocatable :: current
integer :: thermal_source_maxSizeDotState
interface interface
module function source_thermal_dissipation_init(source_length) result(mySources) module function source_thermal_dissipation_init(source_length) result(mySources)
@ -60,30 +63,55 @@ module subroutine thermal_init(phases)
class(tNode), pointer :: & class(tNode), pointer :: &
phases phases
class(tNode), pointer :: &
phase, thermal, sources
integer :: & integer :: &
ph, & ph, so, &
Nconstituents Nconstituents
print'(/,a)', ' <<<+- constitutive_mech init -+>>>' print'(/,a)', ' <<<+- constitutive_thermal init -+>>>'
allocate(current(phases%length)) allocate(current(phases%length))
allocate(thermalState (phases%length))
allocate(thermal_Nsources(phases%length),source = 0)
do ph = 1, phases%length do ph = 1, phases%length
Nconstituents = count(material_phaseAt == ph) * discretization_nIPs Nconstituents = count(material_phaseAt == ph) * discretization_nIPs
allocate(current(ph)%T(Nconstituents)) allocate(current(ph)%T(Nconstituents))
phase => phases%get(ph)
if(phase%contains('thermal')) then
thermal => phase%get('thermal')
sources => thermal%get('source',defaultVal=emptyList)
thermal_Nsources(ph) = sources%length
endif
allocate(thermalstate(ph)%p(thermal_Nsources(ph)))
enddo enddo
! initialize source mechanisms allocate(thermal_source(maxval(thermal_Nsources),phases%length), source = SOURCE_undefined_ID)
if(maxval(phase_Nsources) /= 0) then
where(source_thermal_dissipation_init (maxval(phase_Nsources))) phase_source = SOURCE_thermal_dissipation_ID if(maxval(thermal_Nsources) /= 0) then
where(source_thermal_externalheat_init(maxval(phase_Nsources))) phase_source = SOURCE_thermal_externalheat_ID where(source_thermal_dissipation_init (maxval(thermal_Nsources))) thermal_source = SOURCE_thermal_dissipation_ID
where(source_thermal_externalheat_init(maxval(thermal_Nsources))) thermal_source = SOURCE_thermal_externalheat_ID
endif endif
thermal_source_maxSizeDotState = 0
PhaseLoop2:do ph = 1,phases%length
do so = 1,thermal_Nsources(ph)
thermalState(ph)%p(so)%partitionedState0 = thermalState(ph)%p(so)%state0
thermalState(ph)%p(so)%state = thermalState(ph)%p(so)%partitionedState0
enddo
thermal_source_maxSizeDotState = max(thermal_source_maxSizeDotState, &
maxval(thermalState(ph)%p%sizeDotState))
enddo PhaseLoop2
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!initialize kinematic mechanisms !initialize kinematic mechanisms
if(maxval(phase_Nkinematics) /= 0) where(kinematics_thermal_expansion_init(maxval(phase_Nkinematics))) & if(maxval(phase_Nkinematics) /= 0) where(kinematics_thermal_expansion_init(maxval(phase_Nkinematics))) &
@ -123,8 +151,8 @@ module subroutine constitutive_thermal_getRateAndItsTangents(TDot, dTDot_dT, T,
do co = 1, homogenization_Nconstituents(homog) do co = 1, homogenization_Nconstituents(homog)
ph = material_phaseAt(co,el) ph = material_phaseAt(co,el)
me = material_phasememberAt(co,ip,el) me = material_phasememberAt(co,ip,el)
do so = 1, phase_Nsources(ph) do so = 1, thermal_Nsources(ph)
select case(phase_source(so,ph)) select case(thermal_source(so,ph))
case (SOURCE_thermal_dissipation_ID) case (SOURCE_thermal_dissipation_ID)
call source_thermal_dissipation_getRateAndItsTangent(my_Tdot, my_dTdot_dT, & call source_thermal_dissipation_getRateAndItsTangent(my_Tdot, my_dTdot_dT, &
mech_S(ph,me),mech_L_p(ph,me), ph) mech_S(ph,me),mech_L_p(ph,me), ph)
@ -145,6 +173,131 @@ module subroutine constitutive_thermal_getRateAndItsTangents(TDot, dTDot_dT, T,
end subroutine constitutive_thermal_getRateAndItsTangents end subroutine constitutive_thermal_getRateAndItsTangents
!--------------------------------------------------------------------------------------------------
!> @brief contains the constitutive equation for calculating the rate of change of microstructure
!--------------------------------------------------------------------------------------------------
function constitutive_thermal_collectDotState(ph,me) result(broken)
integer, intent(in) :: ph, me
logical :: broken
integer :: i
broken = .false.
SourceLoop: do i = 1, thermal_Nsources(ph)
if (thermal_source(i,ph) == SOURCE_thermal_externalheat_ID) &
call source_thermal_externalheat_dotState(ph,me)
broken = broken .or. any(IEEE_is_NaN(thermalState(ph)%p(i)%dotState(:,me)))
enddo SourceLoop
end function constitutive_thermal_collectDotState
!--------------------------------------------------------------------------------------------------
!> @brief integrate stress, state with adaptive 1st order explicit Euler method
!> using Fixed Point Iteration to adapt the stepsize
!--------------------------------------------------------------------------------------------------
module function integrateThermalState(dt,co,ip,el) result(broken)
real(pReal), intent(in) :: dt
integer, intent(in) :: &
el, & !< element index in element loop
ip, & !< integration point index in ip loop
co !< grain index in grain loop
integer :: &
NiterationState, & !< number of iterations in state loop
ph, &
me, &
so
integer, dimension(maxval(thermal_Nsources)) :: &
size_so
real(pReal) :: &
zeta
real(pReal), dimension(thermal_source_maxSizeDotState) :: &
r ! state residuum
real(pReal), dimension(thermal_source_maxSizeDotState,2,maxval(thermal_Nsources)) :: source_dotState
logical :: &
broken, converged_
ph = material_phaseAt(co,el)
me = material_phaseMemberAt(co,ip,el)
converged_ = .true.
broken = constitutive_thermal_collectDotState(ph,me)
if(broken) return
do so = 1, thermal_Nsources(ph)
size_so(so) = thermalState(ph)%p(so)%sizeDotState
thermalState(ph)%p(so)%state(1:size_so(so),me) = thermalState(ph)%p(so)%subState0(1:size_so(so),me) &
+ thermalState(ph)%p(so)%dotState (1:size_so(so),me) * dt
source_dotState(1:size_so(so),2,so) = 0.0_pReal
enddo
iteration: do NiterationState = 1, num%nState
do so = 1, thermal_Nsources(ph)
if(nIterationState > 1) source_dotState(1:size_so(so),2,so) = source_dotState(1:size_so(so),1,so)
source_dotState(1:size_so(so),1,so) = thermalState(ph)%p(so)%dotState(:,me)
enddo
broken = constitutive_thermal_collectDotState(ph,me)
if(broken) exit iteration
do so = 1, thermal_Nsources(ph)
zeta = damper(thermalState(ph)%p(so)%dotState(:,me), &
source_dotState(1:size_so(so),1,so),&
source_dotState(1:size_so(so),2,so))
thermalState(ph)%p(so)%dotState(:,me) = thermalState(ph)%p(so)%dotState(:,me) * zeta &
+ source_dotState(1:size_so(so),1,so)* (1.0_pReal - zeta)
r(1:size_so(so)) = thermalState(ph)%p(so)%state (1:size_so(so),me) &
- thermalState(ph)%p(so)%subState0(1:size_so(so),me) &
- thermalState(ph)%p(so)%dotState (1:size_so(so),me) * dt
thermalState(ph)%p(so)%state(1:size_so(so),me) = thermalState(ph)%p(so)%state(1:size_so(so),me) &
- r(1:size_so(so))
converged_ = converged_ .and. converged(r(1:size_so(so)), &
thermalState(ph)%p(so)%state(1:size_so(so),me), &
thermalState(ph)%p(so)%atol(1:size_so(so)))
enddo
if(converged_) exit iteration
enddo iteration
broken = broken .or. .not. converged_
contains
!--------------------------------------------------------------------------------------------------
!> @brief calculate the damping for correction of state and dot state
!--------------------------------------------------------------------------------------------------
real(pReal) pure function damper(current,previous,previous2)
real(pReal), dimension(:), intent(in) ::&
current, previous, previous2
real(pReal) :: dot_prod12, dot_prod22
dot_prod12 = dot_product(current - previous, previous - previous2)
dot_prod22 = dot_product(previous - previous2, previous - previous2)
if ((dot_product(current,previous) < 0.0_pReal .or. dot_prod12 < 0.0_pReal) .and. dot_prod22 > 0.0_pReal) then
damper = 0.75_pReal + 0.25_pReal * tanh(2.0_pReal + 4.0_pReal * dot_prod12 / dot_prod22)
else
damper = 1.0_pReal
endif
end function damper
end function integrateThermalState
module subroutine thermal_initializeRestorationPoints(ph,me) module subroutine thermal_initializeRestorationPoints(ph,me)
@ -153,7 +306,7 @@ module subroutine thermal_initializeRestorationPoints(ph,me)
integer :: so integer :: so
do so = 1, size(sourceState(ph)%p) do so = 1, size(thermalState(ph)%p)
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
@ -168,7 +321,7 @@ module subroutine thermal_windForward(ph,me)
integer :: so integer :: so
do so = 1, size(sourceState(ph)%p) do so = 1, size(thermalState(ph)%p)
thermalState(ph)%p(so)%partitionedState0(:,me) = thermalState(ph)%p(so)%state(:,me) thermalState(ph)%p(so)%partitionedState0(:,me) = thermalState(ph)%p(so)%state(:,me)
enddo enddo
@ -181,7 +334,7 @@ module subroutine thermal_forward()
do ph = 1, size(thermalState) do ph = 1, size(thermalState)
do so = 1, size(sourceState(ph)%p) do so = 1, size(thermalState(ph)%p)
thermalState(ph)%p(so)%state0 = thermalState(ph)%p(so)%state thermalState(ph)%p(so)%state0 = thermalState(ph)%p(so)%state
enddo enddo
enddo enddo
@ -200,7 +353,7 @@ module subroutine thermal_restore(ip,el)
ph = material_phaseAt(co,el) ph = material_phaseAt(co,el)
me = material_phaseMemberAt(co,ip,el) me = material_phaseMemberAt(co,ip,el)
do so = 1, size(sourceState(ph)%p) do so = 1, size(thermalState(ph)%p)
thermalState(ph)%p(so)%state(:,me) = thermalState(ph)%p(so)%partitionedState0(:,me) thermalState(ph)%p(so)%state(:,me) = thermalState(ph)%p(so)%partitionedState0(:,me)
enddo enddo
@ -237,4 +390,39 @@ module subroutine constitutive_thermal_setT(T,co,ip,el)
end subroutine constitutive_thermal_setT end subroutine constitutive_thermal_setT
!--------------------------------------------------------------------------------------------------
!> @brief checks if a source mechanism is active or not
!--------------------------------------------------------------------------------------------------
function thermal_active(source_label,src_length) result(active_source)
character(len=*), intent(in) :: source_label !< name of source mechanism
integer, intent(in) :: src_length !< max. number of sources in system
logical, dimension(:,:), allocatable :: active_source
class(tNode), pointer :: &
phases, &
phase, &
sources, thermal, &
src
integer :: p,s
phases => config_material%get('phase')
allocate(active_source(src_length,phases%length), source = .false. )
do p = 1, phases%length
phase => phases%get(p)
if (phase%contains('thermal')) then
thermal => phase%get('thermal',defaultVal=emptyList)
sources => thermal%get('source',defaultVal=emptyList)
do s = 1, sources%length
src => sources%get(s)
if(src%get_asString('type') == source_label) active_source(s,p) = .true.
enddo
endif
enddo
end function thermal_active
end submodule constitutive_thermal end submodule constitutive_thermal

2
src/constitutive_thermal_dissipation.f90
View File

@ -62,7 +62,7 @@ module function source_thermal_dissipation_init(source_length) result(mySources)
src => sources%get(sourceOffset) src => sources%get(sourceOffset)
prm%kappa = src%get_asFloat('kappa') prm%kappa = src%get_asFloat('kappa')
Nconstituents = count(material_phaseAt==p) * discretization_nIPs Nconstituents = count(material_phaseAt==p) * discretization_nIPs
call constitutive_allocateState(sourceState(p)%p(sourceOffset),Nconstituents,0,0,0) call constitutive_allocateState(thermalState(p)%p(sourceOffset),Nconstituents,0,0,0)
end associate end associate
endif endif

24
src/constitutive_thermal_externalheat.f90
View File

@ -13,7 +13,7 @@ submodule(constitutive:constitutive_thermal) source_externalheat
type :: tParameters !< container type for internal constitutive parameters type :: tParameters !< container type for internal constitutive parameters
real(pReal), dimension(:), allocatable :: & real(pReal), dimension(:), allocatable :: &
t_n, & t_n, &
f_T f_T
integer :: & integer :: &
nIntervals nIntervals
@ -31,19 +31,20 @@ contains
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module function source_thermal_externalheat_init(source_length) result(mySources) module function source_thermal_externalheat_init(source_length) result(mySources)
integer, intent(in) :: source_length integer, intent(in) :: source_length
logical, dimension(:,:), allocatable :: mySources logical, dimension(:,:), allocatable :: mySources
class(tNode), pointer :: & class(tNode), pointer :: &
phases, & phases, &
phase, & phase, &
sources, & sources, thermal, &
src src
integer :: Ninstances,sourceOffset,Nconstituents,p integer :: Ninstances,sourceOffset,Nconstituents,p
print'(/,a)', ' <<<+- source_thermal_externalHeat init -+>>>' print'(/,a)', ' <<<+- source_thermal_externalHeat init -+>>>'
mySources = source_active('thermal_externalheat',source_length) mySources = thermal_active('externalheat',source_length)
Ninstances = count(mySources) Ninstances = count(mySources)
print'(a,i2)', ' # instances: ',Ninstances; flush(IO_STDOUT) print'(a,i2)', ' # instances: ',Ninstances; flush(IO_STDOUT)
if(Ninstances == 0) return if(Ninstances == 0) return
@ -54,15 +55,16 @@ module function source_thermal_externalheat_init(source_length) result(mySources
allocate(source_thermal_externalheat_instance(phases%length), source=0) allocate(source_thermal_externalheat_instance(phases%length), source=0)
do p = 1, phases%length do p = 1, phases%length
phase => phases%get(p) phase => phases%get(p)
if(any(mySources(:,p))) source_thermal_externalheat_instance(p) = count(mySources(:,1:p)) if(any(mySources(:,p))) source_thermal_externalheat_instance(p) = count(mySources(:,1:p))
if(count(mySources(:,p)) == 0) cycle if(count(mySources(:,p)) == 0) cycle
sources => phase%get('source') thermal => phase%get('thermal')
sources => thermal%get('source')
do sourceOffset = 1, sources%length do sourceOffset = 1, sources%length
if(mySources(sourceOffset,p)) then if(mySources(sourceOffset,p)) then
source_thermal_externalheat_offset(p) = sourceOffset source_thermal_externalheat_offset(p) = sourceOffset
associate(prm => param(source_thermal_externalheat_instance(p))) associate(prm => param(source_thermal_externalheat_instance(p)))
src => sources%get(sourceOffset) src => sources%get(sourceOffset)
prm%t_n = src%get_asFloats('t_n') prm%t_n = src%get_asFloats('t_n')
prm%nIntervals = size(prm%t_n) - 1 prm%nIntervals = size(prm%t_n) - 1
@ -70,7 +72,7 @@ module function source_thermal_externalheat_init(source_length) result(mySources
prm%f_T = src%get_asFloats('f_T',requiredSize = size(prm%t_n)) prm%f_T = src%get_asFloats('f_T',requiredSize = size(prm%t_n))
Nconstituents = count(material_phaseAt==p) * discretization_nIPs Nconstituents = count(material_phaseAt==p) * discretization_nIPs
call constitutive_allocateState(sourceState(p)%p(sourceOffset),Nconstituents,1,1,0) call constitutive_allocateState(thermalState(p)%p(sourceOffset),Nconstituents,1,1,0)
end associate end associate
endif endif
@ -95,7 +97,7 @@ module subroutine source_thermal_externalheat_dotState(phase, of)
sourceOffset = source_thermal_externalheat_offset(phase) sourceOffset = source_thermal_externalheat_offset(phase)
sourceState(phase)%p(sourceOffset)%dotState(1,of) = 1.0_pReal ! state is current time thermalState(phase)%p(sourceOffset)%dotState(1,of) = 1.0_pReal ! state is current time
end subroutine source_thermal_externalheat_dotState end subroutine source_thermal_externalheat_dotState
@ -121,7 +123,7 @@ module subroutine source_thermal_externalheat_getRateAndItsTangent(TDot, dTDot_d
associate(prm => param(source_thermal_externalheat_instance(phase))) associate(prm => param(source_thermal_externalheat_instance(phase)))
do interval = 1, prm%nIntervals ! scan through all rate segments do interval = 1, prm%nIntervals ! scan through all rate segments
frac_time = (sourceState(phase)%p(sourceOffset)%state(1,of) - prm%t_n(interval)) & frac_time = (thermalState(phase)%p(sourceOffset)%state(1,of) - prm%t_n(interval)) &
/ (prm%t_n(interval+1) - prm%t_n(interval)) ! fractional time within segment / (prm%t_n(interval+1) - prm%t_n(interval)) ! fractional time within segment
if ( (frac_time < 0.0_pReal .and. interval == 1) & if ( (frac_time < 0.0_pReal .and. interval == 1) &
.or. (frac_time >= 1.0_pReal .and. interval == prm%nIntervals) & .or. (frac_time >= 1.0_pReal .and. interval == prm%nIntervals) &

12
src/homogenization.f90
View File

@ -69,6 +69,13 @@ module homogenization
el !< element number el !< element number
end subroutine mech_partition end subroutine mech_partition
module subroutine thermal_partition(T,ip,el)
real(pReal), intent(in) :: T
integer, intent(in) :: &
ip, & !< integration point
el !< element number
end subroutine thermal_partition
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) :: &
@ -131,9 +138,10 @@ subroutine homogenization_init
call mech_init(num_homog) call mech_init(num_homog)
call thermal_init()
if (any(thermal_type == THERMAL_isothermal_ID)) call thermal_isothermal_init if (any(thermal_type == THERMAL_isothermal_ID)) call thermal_isothermal_init(homogenization_T)
if (any(thermal_type == THERMAL_conduction_ID)) call thermal_conduction_init if (any(thermal_type == THERMAL_conduction_ID)) call thermal_conduction_init(homogenization_T)
if (any(damage_type == DAMAGE_none_ID)) call damage_none_init if (any(damage_type == DAMAGE_none_ID)) call damage_none_init
if (any(damage_type == DAMAGE_nonlocal_ID)) call damage_nonlocal_init if (any(damage_type == DAMAGE_nonlocal_ID)) call damage_nonlocal_init

4
src/homogenization_thermal.f90
View File

@ -11,9 +11,11 @@ contains
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module subroutine thermal_init() module subroutine thermal_init()
print'(/,a)', ' <<<+- homogenization_thermal init -+>>>' print'(/,a)', ' <<<+- homogenization_thermal init -+>>>'
allocate(homogenization_T(discretization_nIPs*discretization_Nelems), source=0.0_pReal) allocate(homogenization_T(discretization_nIPs*discretization_Nelems))
end subroutine thermal_init end subroutine thermal_init

43
src/thermal_conduction.f90
View File

@ -10,6 +10,7 @@ module thermal_conduction
use results use results
use constitutive use constitutive
use YAML_types use YAML_types
use discretization
implicit none implicit none
private private
@ -38,25 +39,28 @@ contains
!> @brief module initialization !> @brief module initialization
!> @details reads in material parameters, allocates arrays, and does sanity checks !> @details reads in material parameters, allocates arrays, and does sanity checks
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine thermal_conduction_init subroutine thermal_conduction_init(T)
integer :: Ninstances,Nmaterialpoints,h real(pReal), dimension(:), intent(inout) :: T
integer :: Ninstances,Nmaterialpoints,ho,ip,el,ce
class(tNode), pointer :: & class(tNode), pointer :: &
material_homogenization, & material_homogenization, &
homog, & homog, &
homogThermal homogThermal
print'(/,a)', ' <<<+- thermal_conduction init -+>>>'; flush(6) print'(/,a)', ' <<<+- thermal_conduction init -+>>>'; flush(6)
Ninstances = count(thermal_type == THERMAL_conduction_ID) Ninstances = count(thermal_type == THERMAL_conduction_ID)
allocate(param(Ninstances)) allocate(param(Ninstances))
material_homogenization => config_material%get('homogenization') material_homogenization => config_material%get('homogenization')
do h = 1, size(material_name_homogenization) do ho = 1, size(material_name_homogenization)
if (thermal_type(h) /= THERMAL_conduction_ID) cycle if (thermal_type(ho) /= THERMAL_conduction_ID) cycle
homog => material_homogenization%get(h) homog => material_homogenization%get(ho)
homogThermal => homog%get('thermal') homogThermal => homog%get('thermal')
associate(prm => param(thermal_typeInstance(h))) associate(prm => param(thermal_typeInstance(ho)))
#if defined (__GFORTRAN__) #if defined (__GFORTRAN__)
prm%output = output_asStrings(homogThermal) prm%output = output_asStrings(homogThermal)
@ -64,14 +68,23 @@ subroutine thermal_conduction_init
prm%output = homogThermal%get_asStrings('output',defaultVal=emptyStringArray) prm%output = homogThermal%get_asStrings('output',defaultVal=emptyStringArray)
#endif #endif
Nmaterialpoints=count(material_homogenizationAt==h) Nmaterialpoints=count(material_homogenizationAt==ho)
allocate (temperature (h)%p(Nmaterialpoints), source=thermal_initialT(h)) allocate (temperature (ho)%p(Nmaterialpoints), source=thermal_initialT(ho))
allocate (temperatureRate(h)%p(Nmaterialpoints), source=0.0_pReal) allocate (temperatureRate(ho)%p(Nmaterialpoints), source=0.0_pReal)
end associate end associate
enddo enddo
ce = 0
do el = 1, discretization_Nelems
do ip = 1, discretization_nIPs
ce = ce + 1
ho = material_homogenizationAt(el)
if (thermal_type(ho) == THERMAL_conduction_ID) T(ce) = thermal_initialT(ho)
enddo
enddo
end subroutine thermal_conduction_init end subroutine thermal_conduction_init
@ -89,12 +102,12 @@ subroutine thermal_conduction_getSourceAndItsTangent(Tdot, dTdot_dT, T, ip, el)
Tdot, dTdot_dT Tdot, dTdot_dT
integer :: & integer :: &
homog homog
Tdot = 0.0_pReal Tdot = 0.0_pReal
dTdot_dT = 0.0_pReal dTdot_dT = 0.0_pReal
homog = material_homogenizationAt(el) homog = material_homogenizationAt(el)
call constitutive_thermal_getRateAndItsTangents(TDot, dTDot_dT, T, ip, el) call constitutive_thermal_getRateAndItsTangents(TDot, dTDot_dT, T, ip, el)
Tdot = Tdot/real(homogenization_Nconstituents(homog),pReal) Tdot = Tdot/real(homogenization_Nconstituents(homog),pReal)
dTdot_dT = dTdot_dT/real(homogenization_Nconstituents(homog),pReal) dTdot_dT = dTdot_dT/real(homogenization_Nconstituents(homog),pReal)
@ -112,13 +125,13 @@ function thermal_conduction_getConductivity(ip,el)
el !< element number el !< element number
real(pReal), dimension(3,3) :: & real(pReal), dimension(3,3) :: &
thermal_conduction_getConductivity thermal_conduction_getConductivity
integer :: & integer :: &
co co
thermal_conduction_getConductivity = 0.0_pReal thermal_conduction_getConductivity = 0.0_pReal
do co = 1, homogenization_Nconstituents(material_homogenizationAt(el)) do co = 1, homogenization_Nconstituents(material_homogenizationAt(el))
thermal_conduction_getConductivity = thermal_conduction_getConductivity + & thermal_conduction_getConductivity = thermal_conduction_getConductivity + &
crystallite_push33ToRef(co,ip,el,lattice_K(:,:,material_phaseAt(co,el))) crystallite_push33ToRef(co,ip,el,lattice_K(:,:,material_phaseAt(co,el)))
@ -168,7 +181,7 @@ function thermal_conduction_getMassDensity(ip,el)
el !< element number el !< element number
real(pReal) :: & real(pReal) :: &
thermal_conduction_getMassDensity thermal_conduction_getMassDensity
integer :: & integer :: &
co co

26
src/thermal_isothermal.f90
View File

@ -6,6 +6,7 @@ module thermal_isothermal
use prec use prec
use config use config
use material use material
use discretization
implicit none implicit none
public public
@ -15,22 +16,33 @@ contains
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief allocates fields, reads information from material configuration file !> @brief allocates fields, reads information from material configuration file
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine thermal_isothermal_init subroutine thermal_isothermal_init(T)
integer :: h,Nmaterialpoints real(pReal), dimension(:), intent(inout) :: T
integer :: Ninstances,Nmaterialpoints,ho,ip,el,ce
print'(/,a)', ' <<<+- thermal_isothermal init -+>>>'; flush(6) print'(/,a)', ' <<<+- thermal_isothermal init -+>>>'; flush(6)
do h = 1, size(material_name_homogenization) do ho = 1, size(thermal_type)
if (thermal_type(h) /= THERMAL_isothermal_ID) cycle if (thermal_type(ho) /= THERMAL_isothermal_ID) cycle
Nmaterialpoints = count(material_homogenizationAt == h) Nmaterialpoints = count(material_homogenizationAt == ho)
allocate(temperature (h)%p(Nmaterialpoints),source=thermal_initialT(h)) allocate(temperature (ho)%p(Nmaterialpoints),source=thermal_initialT(ho))
allocate(temperatureRate(h)%p(Nmaterialpoints),source = 0.0_pReal) allocate(temperatureRate(ho)%p(Nmaterialpoints),source = 0.0_pReal)
enddo enddo
ce = 0
do el = 1, discretization_Nelems
do ip = 1, discretization_nIPs
ce = ce + 1
ho = material_homogenizationAt(el)
if (thermal_type(ho) == THERMAL_isothermal_ID) T(ce) = thermal_initialT(ho)
enddo
enddo
end subroutine thermal_isothermal_init end subroutine thermal_isothermal_init
end module thermal_isothermal end module thermal_isothermal