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

Merge branch 'clean-thermal-2' into 'development' 

Clean thermal 2

See merge request damask/DAMASK!326
This commit is contained in:
Franz Roters 2021-01-26 10:14:31 +01:00
parent b783b66d3d 26c7969837
commit bd9bb5c7c3
14 changed files with 307 additions and 428 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

8
src/CPFEM.f90
View File

@ -180,10 +180,10 @@ subroutine CPFEM_general(mode, ffn, ffn1, temperature_inp, dt, elFE, ip, cauchyS
if (iand(mode, CPFEM_AGERESULTS) /= 0_pInt) call CPFEM_forward
chosenThermal1: select case (thermal_type(material_homogenizationAt(elCP)))
case (THERMAL_conduction_ID) chosenThermal1
temperature(material_homogenizationAt(elCP))%p(material_homogenizationMemberAt(ip,elCP)) = &
temperature_inp
end select chosenThermal1
!case (THERMAL_conduction_ID) chosenThermal1
! temperature(material_homogenizationAt(elCP))%p(material_homogenizationMemberAt(ip,elCP)) = &
! temperature_inp
end select chosenThermal1
homogenization_F0(1:3,1:3,ma) = ffn
homogenization_F(1:3,1:3,ma) = ffn1

2
src/DAMASK_marc.f90
View File

@ -351,7 +351,7 @@ end subroutine hypela2
!--------------------------------------------------------------------------------------------------
subroutine flux(f,ts,n,time)
use prec
use thermal_conduction
use homogenization
use discretization_marc
implicit none

2
src/commercialFEM_fileList.f90
View File

@ -42,8 +42,6 @@
#include "source_damage_anisoDuctile.f90"
#include "kinematics_cleavage_opening.f90"
#include "kinematics_slipplane_opening.f90"
#include "thermal_isothermal.f90"
#include "thermal_conduction.f90"
#include "damage_none.f90"
#include "damage_nonlocal.f90"
#include "homogenization.f90"

54
src/constitutive.f90
View File

@ -120,10 +120,6 @@ module constitutive
integer, intent(in) :: ph, me
end subroutine mech_initializeRestorationPoints
module subroutine constitutive_thermal_initializeRestorationPoints(ph,me)
integer, intent(in) :: ph, me
end subroutine constitutive_thermal_initializeRestorationPoints
module subroutine mech_windForward(ph,me)
integer, intent(in) :: ph, me
@ -137,8 +133,8 @@ module constitutive
end subroutine thermal_forward
module subroutine mech_restore(ip,el,includeL)
integer, intent(in) :: ip, el
module subroutine mech_restore(ce,includeL)
integer, intent(in) :: ce
logical, intent(in) :: includeL
end subroutine mech_restore
@ -204,11 +200,11 @@ module constitutive
integer, intent(in) :: co, ip, el
end subroutine constitutive_mech_setF
module subroutine constitutive_thermal_setT(T,co,ce)
real(pReal), intent(in) :: T
integer, intent(in) :: co, ce
end subroutine constitutive_thermal_setT
module subroutine constitutive_thermal_setField(T,dot_T, co,ce)
real(pReal), intent(in) :: T, dot_T
integer, intent(in) :: ce, co
end subroutine constitutive_thermal_setField
module subroutine constitutive_damage_set_phi(phi,co,ce)
real(pReal), intent(in) :: phi
integer, intent(in) :: co, ce
@ -245,9 +241,6 @@ module constitutive
end function constitutive_homogenizedC
module subroutine constitutive_damage_getRateAndItsTangents(phiDot, dPhiDot_dPhi, phi, ip, el)
integer, intent(in) :: &
ip, & !< integration point number
@ -317,11 +310,9 @@ module constitutive
dLd_dTstar !< derivative of Ld with respect to Tstar (4th-order tensor)
end subroutine kinematics_slipplane_opening_LiAndItsTangent
module subroutine kinematics_thermal_expansion_LiAndItsTangent(Li, dLi_dTstar, co, ip, el)
integer, intent(in) :: &
co, & !< grain number
ip, & !< integration point number
el !< element number
module subroutine kinematics_thermal_expansion_LiAndItsTangent(Li, dLi_dTstar, ph,me)
integer, intent(in) :: ph, me
!< element number
real(pReal), intent(out), dimension(3,3) :: &
Li !< thermal velocity gradient
real(pReal), intent(out), dimension(3,3,3,3) :: &
@ -362,14 +353,13 @@ module constitutive
constitutive_restartWrite, &
constitutive_restartRead, &
integrateDamageState, &
constitutive_thermal_setT, &
constitutive_thermal_setField, &
constitutive_damage_set_phi, &
constitutive_damage_get_phi, &
constitutive_mech_getP, &
constitutive_mech_setF, &
constitutive_mech_getF, &
constitutive_initializeRestorationPoints, &
constitutive_thermal_initializeRestorationPoints, &
constitutive_windForward, &
KINEMATICS_UNDEFINED_ID ,&
KINEMATICS_CLEAVAGE_OPENING_ID, &
@ -497,26 +487,24 @@ end subroutine constitutive_allocateState
!--------------------------------------------------------------------------------------------------
!> @brief Restore data after homog cutback.
!--------------------------------------------------------------------------------------------------
subroutine constitutive_restore(ip,el,includeL)
subroutine constitutive_restore(ce,includeL)
logical, intent(in) :: includeL
integer, intent(in) :: &
ip, & !< integration point number
el !< element number
integer, intent(in) :: ce
integer :: &
co, & !< constituent number
so
do co = 1,homogenization_Nconstituents(material_homogenizationAt(el))
do so = 1, phase_Nsources(material_phaseAt(co,el))
damageState(material_phaseAt(co,el))%p(so)%state( :,material_phasememberAt(co,ip,el)) = &
damageState(material_phaseAt(co,el))%p(so)%partitionedState0(:,material_phasememberAt(co,ip,el))
do co = 1,homogenization_Nconstituents(material_homogenizationAt2(ce))
do so = 1, phase_Nsources(material_phaseAt2(co,ce))
damageState(material_phaseAt2(co,ce))%p(so)%state( :,material_phasememberAt2(co,ce)) = &
damageState(material_phaseAt2(co,ce))%p(so)%partitionedState0(:,material_phasememberAt2(co,ce))
enddo
enddo
call mech_restore(ip,el,includeL)
call mech_restore(ce,includeL)
end subroutine constitutive_restore
@ -638,9 +626,6 @@ subroutine crystallite_init()
do so = 1, phase_Nsources(ph)
allocate(damageState(ph)%p(so)%subState0,source=damageState(ph)%p(so)%state0) ! ToDo: hack
enddo
do so = 1, thermal_Nsources(ph)
allocate(thermalState(ph)%p(so)%subState0,source=thermalState(ph)%p(so)%state0) ! ToDo: hack
enddo
enddo
print'(a42,1x,i10)', ' # of elements: ', eMax
@ -766,9 +751,6 @@ function crystallite_push33ToRef(co,ip,el, tensor33)
end function crystallite_push33ToRef
!--------------------------------------------------------------------------------------------------
!> @brief determines whether a point is converged
!--------------------------------------------------------------------------------------------------

21
src/constitutive_mech.f90
View File

@ -1586,9 +1586,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)%partitionedState0(:,me)
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)
subFi0 = constitutive_mech_partitionedFi0(ph)%data(1:3,1:3,me)
subF0 = constitutive_mech_partitionedF0(ph)%data(1:3,1:3,me)
@ -1656,11 +1653,9 @@ end function crystallite_stress
!--------------------------------------------------------------------------------------------------
!> @brief Restore data after homog cutback.
!--------------------------------------------------------------------------------------------------
module subroutine mech_restore(ip,el,includeL)
module subroutine mech_restore(ce,includeL)
integer, intent(in) :: &
ip, & !< integration point number
el !< element number
integer, intent(in) :: ce
logical, intent(in) :: &
includeL !< protect agains fake cutback
@ -1668,9 +1663,9 @@ module subroutine mech_restore(ip,el,includeL)
co, ph, me
do co = 1,homogenization_Nconstituents(material_homogenizationAt(el))
ph = material_phaseAt(co,el)
me = material_phaseMemberAt(co,ip,el)
do co = 1,homogenization_Nconstituents(material_homogenizationAt2(ce))
ph = material_phaseAt2(co,ce)
me = material_phaseMemberAt2(co,ce)
if (includeL) then
constitutive_mech_Lp(ph)%data(1:3,1:3,me) = constitutive_mech_partitionedLp0(ph)%data(1:3,1:3,me)
constitutive_mech_Li(ph)%data(1:3,1:3,me) = constitutive_mech_partitionedLi0(ph)%data(1:3,1:3,me)
@ -1961,7 +1956,7 @@ subroutine constitutive_LiAndItsTangents(Li, dLi_dS, dLi_dFi, &
detFi
integer :: &
k, i, j, &
instance, of
instance, of, me, ph
Li = 0.0_pReal
dLi_dS = 0.0_pReal
@ -1987,7 +1982,9 @@ subroutine constitutive_LiAndItsTangents(Li, dLi_dS, dLi_dFi, &
case (KINEMATICS_slipplane_opening_ID) kinematicsType
call kinematics_slipplane_opening_LiAndItsTangent(my_Li, my_dLi_dS, S, co, ip, el)
case (KINEMATICS_thermal_expansion_ID) kinematicsType
call kinematics_thermal_expansion_LiAndItsTangent(my_Li, my_dLi_dS, co, ip, el)
me = material_phaseMemberAt(co,ip,el)
ph = material_phaseAt(co,el)
call kinematics_thermal_expansion_LiAndItsTangent(my_Li, my_dLi_dS, ph,me)
case default kinematicsType
my_Li = 0.0_pReal
my_dLi_dS = 0.0_pReal

38
src/constitutive_thermal.f90
View File

@ -10,9 +10,9 @@ submodule(constitutive) constitutive_thermal
end enum
type :: tDataContainer
real(pReal), dimension(:), allocatable :: T
real(pReal), dimension(:), allocatable :: T, dot_T
end type tDataContainer
integer(kind(THERMAL_UNDEFINED_ID)), dimension(:,:), allocatable :: &
integer(kind(THERMAL_UNDEFINED_ID)), dimension(:,:), allocatable :: &
thermal_source
type(tDataContainer), dimension(:), allocatable :: current
@ -94,6 +94,7 @@ module subroutine thermal_init(phases)
Nconstituents = count(material_phaseAt == ph) * discretization_nIPs
allocate(current(ph)%T(Nconstituents),source=300.0_pReal)
allocate(current(ph)%dot_T(Nconstituents),source=0.0_pReal)
phase => phases%get(ph)
if(phase%contains('thermal')) then
thermal => phase%get('thermal')
@ -115,8 +116,8 @@ module subroutine thermal_init(phases)
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
deallocate(thermalState(ph)%p(so)%partitionedState0)
thermalState(ph)%p(so)%state = thermalState(ph)%p(so)%state0
enddo
thermal_source_maxSizeDotState = max(thermal_source_maxSizeDotState, &
@ -147,13 +148,11 @@ module subroutine constitutive_thermal_getRate(TDot, ip, el)
integer :: &
ph, &
homog, &
instance, &
me, &
so, &
co
homog = material_homogenizationAt(el)
instance = thermal_typeInstance(homog)
TDot = 0.0_pReal
do co = 1, homogenization_Nconstituents(homog)
@ -211,9 +210,6 @@ module function thermal_stress(Delta_t,ph,me) result(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
@ -238,28 +234,13 @@ function integrateThermalState(Delta_t, ph,me) result(broken)
do so = 1, thermal_Nsources(ph)
sizeDotState = thermalState(ph)%p(so)%sizeDotState
thermalState(ph)%p(so)%state(1:sizeDotState,me) = thermalState(ph)%p(so)%subState0(1:sizeDotState,me) &
thermalState(ph)%p(so)%state(1:sizeDotState,me) = thermalState(ph)%p(so)%state0(1:sizeDotState,me) &
+ thermalState(ph)%p(so)%dotState(1:sizeDotState,me) * Delta_t
enddo
end function integrateThermalState
module subroutine constitutive_thermal_initializeRestorationPoints(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)%state0(:,me)
enddo
end subroutine constitutive_thermal_initializeRestorationPoints
module subroutine thermal_forward()
integer :: ph, so
@ -291,15 +272,16 @@ end function thermal_T
!----------------------------------------------------------------------------------------------
!< @brief Set temperature
!----------------------------------------------------------------------------------------------
module subroutine constitutive_thermal_setT(T,co,ce)
module subroutine constitutive_thermal_setField(T,dot_T, co,ce)
real(pReal), intent(in) :: T
real(pReal), intent(in) :: T, dot_T
integer, intent(in) :: ce, co
current(material_phaseAt2(co,ce))%T(material_phaseMemberAt2(co,ce)) = T
current(material_phaseAt2(co,ce))%dot_T(material_phaseMemberAt2(co,ce)) = dot_T
end subroutine constitutive_thermal_setT
end subroutine constitutive_thermal_setField

2
src/grid/DAMASK_grid.f90
View File

@ -425,7 +425,7 @@ program DAMASK_grid
guess = .true. ! start guessing after first converged (sub)inc
if (worldrank == 0) then
write(statUnit,*) totalIncsCounter, time, cutBackLevel, &
solres%converged, solres%iterationsNeeded
solres(1)%converged, solres(1)%iterationsNeeded
flush(statUnit)
endif
elseif (cutBackLevel < maxCutBack) then ! further cutbacking tolerated?

22
src/grid/grid_thermal_spectral.f90
View File

@ -15,7 +15,6 @@ module grid_thermal_spectral
use IO
use spectral_utilities
use discretization_grid
use thermal_conduction
use homogenization
use YAML_types
use config
@ -132,7 +131,7 @@ subroutine grid_thermal_spectral_init
ce = 0
do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1)
ce = ce + 1
T_current(i,j,k) = temperature(material_homogenizationAt(ce))%p(material_homogenizationMemberAt(1,ce))
T_current(i,j,k) = homogenization_thermal_T(ce)
T_lastInc(i,j,k) = T_current(i,j,k)
T_stagInc(i,j,k) = T_current(i,j,k)
enddo; enddo; enddo
@ -188,10 +187,9 @@ function grid_thermal_spectral_solution(timeinc) result(solution)
ce = 0
do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1)
ce = ce + 1
call thermal_conduction_putTemperatureAndItsRate(T_current(i,j,k), &
call homogenization_thermal_setField(T_current(i,j,k), &
(T_current(i,j,k)-T_lastInc(i,j,k))/params%timeinc, &
1,ce)
homogenization_T(ce) = T_current(i,j,k)
ce)
enddo; enddo; enddo
call VecMin(solution_vec,devNull,T_min,ierr); CHKERRQ(ierr)
@ -229,11 +227,9 @@ subroutine grid_thermal_spectral_forward(cutBack)
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)
ce = ce + 1
call thermal_conduction_putTemperatureAndItsRate(T_current(i,j,k), &
(T_current(i,j,k) - &
T_lastInc(i,j,k))/params%timeinc, &
1,ce)
homogenization_T(ce) = T_current(i,j,k)
call homogenization_thermal_setField(T_current(i,j,k), &
(T_current(i,j,k)-T_lastInc(i,j,k))/params%timeinc, &
ce)
enddo; enddo; enddo
else
T_lastInc = T_current
@ -283,8 +279,8 @@ subroutine formResidual(in,x_scal,f_scal,dummy,ierr)
ce = ce + 1
call thermal_conduction_getSource(Tdot, 1,ce)
scalarField_real(i,j,k) = params%timeinc*(scalarField_real(i,j,k) + Tdot) &
+ thermal_conduction_getMassDensity (1,ce)* &
thermal_conduction_getSpecificHeat(1,ce)*(T_lastInc(i,j,k) - &
+ thermal_conduction_getMassDensity (ce)* &
thermal_conduction_getSpecificHeat(ce)*(T_lastInc(i,j,k) - &
T_current(i,j,k))&
+ mu_ref*T_current(i,j,k)
enddo; enddo; enddo
@ -315,7 +311,7 @@ subroutine updateReference
do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1)
ce = ce + 1
K_ref = K_ref + thermal_conduction_getConductivity(1,ce)
mu_ref = mu_ref + thermal_conduction_getMassDensity(1,ce)* thermal_conduction_getSpecificHeat(1,ce)
mu_ref = mu_ref + thermal_conduction_getMassDensity(ce)* thermal_conduction_getSpecificHeat(ce)
enddo; enddo; enddo
K_ref = K_ref*wgt
call MPI_Allreduce(MPI_IN_PLACE,K_ref,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr)

71
src/homogenization.f90
View File

@ -12,8 +12,6 @@ module homogenization
use material
use constitutive
use discretization
use thermal_isothermal
use thermal_conduction
use damage_none
use damage_nonlocal
use HDF5_utilities
@ -28,8 +26,6 @@ module homogenization
!--------------------------------------------------------------------------------------------------
! General variables for the homogenization at a material point
real(pReal), dimension(:), allocatable, public :: &
homogenization_T, &
homogenization_dot_T, &
homogenization_phi, &
homogenization_dot_phi
real(pReal), dimension(:,:,:), allocatable, public :: &
@ -75,8 +71,7 @@ module homogenization
el !< element number
end subroutine mech_partition
module subroutine thermal_partition(T,ce)
real(pReal), intent(in) :: T
module subroutine thermal_partition(ce)
integer, intent(in) :: ce
end subroutine thermal_partition
@ -112,11 +107,59 @@ module homogenization
logical, dimension(2) :: doneAndHappy
end function mech_updateState
module function thermal_conduction_getConductivity(ip,el) result(K)
integer, intent(in) :: &
ip, & !< integration point number
el !< element number
real(pReal), dimension(3,3) :: K
end function thermal_conduction_getConductivity
module function thermal_conduction_getSpecificHeat(ce) result(c_P)
integer, intent(in) :: ce
real(pReal) :: c_P
end function thermal_conduction_getSpecificHeat
module function thermal_conduction_getMassDensity(ce) result(rho)
integer, intent(in) :: ce
real(pReal) :: rho
end function thermal_conduction_getMassDensity
module subroutine homogenization_thermal_setField(T,dot_T, ce)
integer, intent(in) :: ce
real(pReal), intent(in) :: T, dot_T
end subroutine homogenization_thermal_setField
module subroutine thermal_conduction_results(ho,group)
integer, intent(in) :: ho
character(len=*), intent(in) :: group
end subroutine thermal_conduction_results
module function homogenization_thermal_T(ce) result(T)
integer, intent(in) :: ce
real(pReal) :: T
end function homogenization_thermal_T
module subroutine thermal_conduction_getSource(Tdot, ip,el)
integer, intent(in) :: &
ip, & !< integration point number
el !< element number
real(pReal), intent(out) :: Tdot
end subroutine thermal_conduction_getSource
end interface
public :: &
homogenization_init, &
materialpoint_stressAndItsTangent, &
thermal_conduction_getSpecificHeat, &
thermal_conduction_getConductivity, &
thermal_conduction_getMassDensity, &
thermal_conduction_getSource, &
homogenization_thermal_setfield, &
homogenization_thermal_T, &
homogenization_forward, &
homogenization_results, &
homogenization_restartRead, &
@ -154,9 +197,6 @@ subroutine homogenization_init()
call thermal_init()
call damage_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(homogenization_T)
if (any(damage_type == DAMAGE_none_ID)) call damage_none_init
if (any(damage_type == DAMAGE_nonlocal_ID)) call damage_nonlocal_init
@ -190,9 +230,9 @@ subroutine materialpoint_stressAndItsTangent(dt,FEsolving_execIP,FEsolving_execE
ho = material_homogenizationAt(el)
myNgrains = homogenization_Nconstituents(ho)
do ip = FEsolving_execIP(1),FEsolving_execIP(2)
me = material_homogenizationMemberAt(ip,el)
!--------------------------------------------------------------------------------------------------
! initialize restoration points
ce = (el-1)*discretization_nIPs + ip
me = material_homogenizationMemberAt2(ce)
call constitutive_initializeRestorationPoints(ip,el)
subFrac = 0.0_pReal
@ -226,7 +266,7 @@ subroutine materialpoint_stressAndItsTangent(dt,FEsolving_execIP,FEsolving_execE
else ! cutback makes sense
subStep = num%subStepSizeHomog * subStep ! crystallite had severe trouble, so do a significant cutback
call constitutive_restore(ip,el,subStep < 1.0_pReal)
call constitutive_restore(ce,subStep < 1.0_pReal)
if(homogState(ho)%sizeState > 0) homogState(ho)%State(:,me) = homogState(ho)%subState0(:,me)
if(damageState_h(ho)%sizeState > 0) damageState_h(ho)%State(:,me) = damageState_h(ho)%subState0(:,me)
@ -243,7 +283,6 @@ subroutine materialpoint_stressAndItsTangent(dt,FEsolving_execIP,FEsolving_execE
! deformation partitioning
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)
@ -255,7 +294,6 @@ subroutine materialpoint_stressAndItsTangent(dt,FEsolving_execIP,FEsolving_execE
if (.not. converged) then
doneAndHappy = [.true.,.false.]
else
ce = (el-1)*discretization_nIPs + ip
doneAndHappy = mech_updateState(dt*subStep, &
homogenization_F0(1:3,1:3,ce) &
+ (homogenization_F(1:3,1:3,ce)-homogenization_F0(1:3,1:3,ce)) &
@ -278,10 +316,9 @@ subroutine materialpoint_stressAndItsTangent(dt,FEsolving_execIP,FEsolving_execE
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)
call thermal_partition(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'

199
src/homogenization_thermal.f90
View File

@ -3,6 +3,22 @@
!--------------------------------------------------------------------------------------------------
submodule(homogenization) homogenization_thermal
use lattice
type :: tDataContainer
real(pReal), dimension(:), allocatable :: T, dot_T
end type tDataContainer
type(tDataContainer), dimension(:), allocatable :: current
type :: tParameters
character(len=pStringLen), allocatable, dimension(:) :: &
output
end type tParameters
type(tparameters), dimension(:), allocatable :: &
param
contains
@ -11,11 +27,37 @@ contains
!--------------------------------------------------------------------------------------------------
module subroutine thermal_init()
class(tNode), pointer :: &
configHomogenizations, &
configHomogenization, &
configHomogenizationThermal
integer :: ho
print'(/,a)', ' <<<+- homogenization_thermal init -+>>>'
allocate(homogenization_T(discretization_nIPs*discretization_Nelems))
allocate(homogenization_dot_T(discretization_nIPs*discretization_Nelems))
configHomogenizations => config_material%get('homogenization')
allocate(param(configHomogenizations%length))
allocate(current(configHomogenizations%length))
do ho = 1, configHomogenizations%length
allocate(current(ho)%T(count(material_homogenizationAt2==ho)), source=thermal_initialT(ho))
allocate(current(ho)%dot_T(count(material_homogenizationAt2==ho)), source=0.0_pReal)
configHomogenization => configHomogenizations%get(ho)
associate(prm => param(ho))
if (configHomogenization%contains('thermal')) then
configHomogenizationThermal => configHomogenization%get('thermal')
#if defined (__GFORTRAN__)
prm%output = output_asStrings(configHomogenizationThermal)
#else
prm%output = configHomogenizationThermal%get_asStrings('output',defaultVal=emptyStringArray)
#endif
else
prm%output = emptyStringArray
endif
end associate
enddo
end subroutine thermal_init
@ -23,15 +65,18 @@ end subroutine thermal_init
!--------------------------------------------------------------------------------------------------
!> @brief Partition temperature onto the individual constituents.
!--------------------------------------------------------------------------------------------------
module subroutine thermal_partition(T,ce)
module subroutine thermal_partition(ce)
real(pReal), intent(in) :: T
integer, intent(in) :: ce
real(pReal) :: T, dot_T
integer :: co
T = current(material_homogenizationAt2(ce))%T(material_homogenizationMemberAt2(ce))
dot_T = current(material_homogenizationAt2(ce))%dot_T(material_homogenizationMemberAt2(ce))
do co = 1, homogenization_Nconstituents(material_homogenizationAt2(ce))
call constitutive_thermal_setT(T,co,ce)
call constitutive_thermal_setField(T,dot_T,co,ce)
enddo
end subroutine thermal_partition
@ -44,9 +89,151 @@ 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)
!call constitutive_thermal_getRate(homogenization_dot_T((el-1)*discretization_nIPs+ip), ip,el)
end subroutine thermal_homogenize
!--------------------------------------------------------------------------------------------------
!> @brief return homogenized thermal conductivity in reference configuration
!--------------------------------------------------------------------------------------------------
module function thermal_conduction_getConductivity(ip,el) result(K)
integer, intent(in) :: &
ip, & !< integration point number
el !< element number
real(pReal), dimension(3,3) :: K
integer :: &
co
K = 0.0_pReal
do co = 1, homogenization_Nconstituents(material_homogenizationAt(el))
K = K + crystallite_push33ToRef(co,ip,el,lattice_K(:,:,material_phaseAt(co,el)))
enddo
K = K / real(homogenization_Nconstituents(material_homogenizationAt(el)),pReal)
end function thermal_conduction_getConductivity
!--------------------------------------------------------------------------------------------------
!> @brief returns homogenized specific heat capacity
!--------------------------------------------------------------------------------------------------
module function thermal_conduction_getSpecificHeat(ce) result(c_P)
integer, intent(in) :: ce
real(pReal) :: c_P
integer :: co
c_P = 0.0_pReal
do co = 1, homogenization_Nconstituents(material_homogenizationAt2(ce))
c_P = c_P + lattice_c_p(material_phaseAt2(co,ce))
enddo
c_P = c_P / real(homogenization_Nconstituents(material_homogenizationAt2(ce)),pReal)
end function thermal_conduction_getSpecificHeat
!--------------------------------------------------------------------------------------------------
!> @brief returns homogenized mass density
!--------------------------------------------------------------------------------------------------
module function thermal_conduction_getMassDensity(ce) result(rho)
integer, intent(in) :: ce
real(pReal) :: rho
integer :: co
rho = 0.0_pReal
do co = 1, homogenization_Nconstituents(material_homogenizationAt2(ce))
rho = rho + lattice_rho(material_phaseAt2(co,ce))
enddo
rho = rho / real(homogenization_Nconstituents(material_homogenizationAt2(ce)),pReal)
end function thermal_conduction_getMassDensity
!--------------------------------------------------------------------------------------------------
!> @brief Set thermal field and its rate (T and dot_T)
!--------------------------------------------------------------------------------------------------
module subroutine homogenization_thermal_setField(T,dot_T, ce)
integer, intent(in) :: ce
real(pReal), intent(in) :: T, dot_T
current(material_homogenizationAt2(ce))%T(material_homogenizationMemberAt2(ce)) = T
current(material_homogenizationAt2(ce))%dot_T(material_homogenizationMemberAt2(ce)) = dot_T
end subroutine homogenization_thermal_setField
!--------------------------------------------------------------------------------------------------
!> @brief writes results to HDF5 output file
!--------------------------------------------------------------------------------------------------
module subroutine thermal_conduction_results(ho,group)
integer, intent(in) :: ho
character(len=*), intent(in) :: group
integer :: o
associate(prm => param(ho))
outputsLoop: do o = 1,size(prm%output)
select case(trim(prm%output(o)))
case('T')
call results_writeDataset(group,current(ho)%T,'T','temperature','K')
end select
enddo outputsLoop
end associate
end subroutine thermal_conduction_results
module function homogenization_thermal_T(ce) result(T)
integer, intent(in) :: ce
real(pReal) :: T
T = current(material_homogenizationAt2(ce))%T(material_homogenizationMemberAt2(ce))
end function homogenization_thermal_T
!--------------------------------------------------------------------------------------------------
!> @brief return heat generation rate
!--------------------------------------------------------------------------------------------------
module subroutine thermal_conduction_getSource(Tdot, ip,el)
integer, intent(in) :: &
ip, & !< integration point number
el !< element number
real(pReal), intent(out) :: &
Tdot
integer :: &
homog
homog = material_homogenizationAt(el)
call constitutive_thermal_getRate(TDot, ip,el)
Tdot = Tdot/real(homogenization_Nconstituents(homog),pReal)
end subroutine thermal_conduction_getSource
end submodule homogenization_thermal

20
src/kinematics_thermal_expansion.f90
View File

@ -84,12 +84,9 @@ end function kinematics_thermal_expansion_init
!--------------------------------------------------------------------------------------------------
!> @brief constitutive equation for calculating the velocity gradient
!--------------------------------------------------------------------------------------------------
module subroutine kinematics_thermal_expansion_LiAndItsTangent(Li, dLi_dTstar, co, ip, el)
module subroutine kinematics_thermal_expansion_LiAndItsTangent(Li, dLi_dTstar, ph,me)
integer, intent(in) :: &
co, & !< grain number
ip, & !< integration point number
el !< element number
integer, intent(in) :: ph, me
real(pReal), intent(out), dimension(3,3) :: &
Li !< thermal velocity gradient
real(pReal), intent(out), dimension(3,3,3,3) :: &
@ -98,16 +95,13 @@ module subroutine kinematics_thermal_expansion_LiAndItsTangent(Li, dLi_dTstar, c
integer :: &
phase, &
homog
real(pReal) :: &
T, TDot
real(pReal) :: T, dot_T
phase = material_phaseAt(co,el)
homog = material_homogenizationAt(el)
T = temperature(homog)%p(material_homogenizationMemberAt(ip,el))
TDot = temperatureRate(homog)%p(material_homogenizationMemberAt(ip,el))
T = current(ph)%T(me)
dot_T = current(ph)%dot_T(me)
associate(prm => param(kinematics_thermal_expansion_instance(phase)))
Li = TDot * ( &
associate(prm => param(kinematics_thermal_expansion_instance(ph)))
Li = dot_T * ( &
prm%A(1:3,1:3,1)*(T - prm%T_ref)**0 & ! constant coefficient
+ prm%A(1:3,1:3,2)*(T - prm%T_ref)**1 & ! linear coefficient
+ prm%A(1:3,1:3,3)*(T - prm%T_ref)**2 & ! quadratic coefficient

6
src/material.f90
View File

@ -71,9 +71,7 @@ module material
material_orientation0 !< initial orientation of each grain,IP,element
type(group_float), allocatable, dimension(:), public :: &
temperature, & !< temperature field
damage, & !< damage field
temperatureRate !< temperature change rate field
damage !< damage field
public :: &
material_init, &
@ -107,9 +105,7 @@ subroutine material_init(restart)
allocate(homogState (size(material_name_homogenization)))
allocate(damageState_h (size(material_name_homogenization)))
allocate(temperature (size(material_name_homogenization)))
allocate(damage (size(material_name_homogenization)))
allocate(temperatureRate (size(material_name_homogenization)))
if (.not. restart) then

242
src/thermal_conduction.f90
View File

@ -1,242 +0,0 @@
!--------------------------------------------------------------------------------------------------
!> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH
!> @brief material subroutine for temperature evolution from heat conduction
!--------------------------------------------------------------------------------------------------
module thermal_conduction
use prec
use material
use config
use lattice
use results
use constitutive
use YAML_types
use discretization
implicit none
private
type :: tParameters
character(len=pStringLen), allocatable, dimension(:) :: &
output
end type tParameters
type(tparameters), dimension(:), allocatable :: &
param
public :: &
thermal_conduction_init, &
thermal_conduction_getSource, &
thermal_conduction_getConductivity, &
thermal_conduction_getSpecificHeat, &
thermal_conduction_getMassDensity, &
thermal_conduction_putTemperatureAndItsRate, &
thermal_conduction_results
contains
!--------------------------------------------------------------------------------------------------
!> @brief module initialization
!> @details reads in material parameters, allocates arrays, and does sanity checks
!--------------------------------------------------------------------------------------------------
subroutine thermal_conduction_init(T)
real(pReal), dimension(:), intent(inout) :: T
integer :: Ninstances,Nmaterialpoints,ho,ip,el,ce
class(tNode), pointer :: &
material_homogenization, &
homog, &
homogThermal
print'(/,a)', ' <<<+- thermal_conduction init -+>>>'; flush(6)
Ninstances = count(thermal_type == THERMAL_conduction_ID)
allocate(param(Ninstances))
material_homogenization => config_material%get('homogenization')
do ho = 1, size(material_name_homogenization)
if (thermal_type(ho) /= THERMAL_conduction_ID) cycle
homog => material_homogenization%get(ho)
homogThermal => homog%get('thermal')
associate(prm => param(thermal_typeInstance(ho)))
#if defined (__GFORTRAN__)
prm%output = output_asStrings(homogThermal)
#else
prm%output = homogThermal%get_asStrings('output',defaultVal=emptyStringArray)
#endif
Nmaterialpoints=count(material_homogenizationAt==ho)
allocate (temperature (ho)%p(Nmaterialpoints), source=thermal_initialT(ho))
allocate (temperatureRate(ho)%p(Nmaterialpoints), source=0.0_pReal)
end associate
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
!--------------------------------------------------------------------------------------------------
!> @brief return heat generation rate
!--------------------------------------------------------------------------------------------------
subroutine thermal_conduction_getSource(Tdot, ip,el)
integer, intent(in) :: &
ip, & !< integration point number
el !< element number
real(pReal), intent(out) :: &
Tdot
integer :: &
homog
homog = material_homogenizationAt(el)
call constitutive_thermal_getRate(TDot, ip,el)
Tdot = Tdot/real(homogenization_Nconstituents(homog),pReal)
end subroutine thermal_conduction_getSource
!--------------------------------------------------------------------------------------------------
!> @brief return homogenized thermal conductivity in reference configuration
!--------------------------------------------------------------------------------------------------
function thermal_conduction_getConductivity(ip,el)
integer, intent(in) :: &
ip, & !< integration point number
el !< element number
real(pReal), dimension(3,3) :: &
thermal_conduction_getConductivity
integer :: &
co
thermal_conduction_getConductivity = 0.0_pReal
do co = 1, homogenization_Nconstituents(material_homogenizationAt(el))
thermal_conduction_getConductivity = thermal_conduction_getConductivity + &
crystallite_push33ToRef(co,ip,el,lattice_K(:,:,material_phaseAt(co,el)))
enddo
thermal_conduction_getConductivity = thermal_conduction_getConductivity &
/ real(homogenization_Nconstituents(material_homogenizationAt(el)),pReal)
end function thermal_conduction_getConductivity
!--------------------------------------------------------------------------------------------------
!> @brief returns homogenized specific heat capacity
!--------------------------------------------------------------------------------------------------
function thermal_conduction_getSpecificHeat(ip,el)
integer, intent(in) :: &
ip, & !< integration point number
el !< element number
real(pReal) :: &
thermal_conduction_getSpecificHeat
integer :: &
co
thermal_conduction_getSpecificHeat = 0.0_pReal
do co = 1, homogenization_Nconstituents(material_homogenizationAt(el))
thermal_conduction_getSpecificHeat = thermal_conduction_getSpecificHeat &
+ lattice_c_p(material_phaseAt(co,el))
enddo
thermal_conduction_getSpecificHeat = thermal_conduction_getSpecificHeat &
/ real(homogenization_Nconstituents(material_homogenizationAt(el)),pReal)
end function thermal_conduction_getSpecificHeat
!--------------------------------------------------------------------------------------------------
!> @brief returns homogenized mass density
!--------------------------------------------------------------------------------------------------
function thermal_conduction_getMassDensity(ip,el)
integer, intent(in) :: &
ip, & !< integration point number
el !< element number
real(pReal) :: &
thermal_conduction_getMassDensity
integer :: &
co
thermal_conduction_getMassDensity = 0.0_pReal
do co = 1, homogenization_Nconstituents(material_homogenizationAt(el))
thermal_conduction_getMassDensity = thermal_conduction_getMassDensity &
+ lattice_rho(material_phaseAt(co,el))
enddo
thermal_conduction_getMassDensity = thermal_conduction_getMassDensity &
/ real(homogenization_Nconstituents(material_homogenizationAt(el)),pReal)
end function thermal_conduction_getMassDensity
!--------------------------------------------------------------------------------------------------
!> @brief updates thermal state with solution from heat conduction PDE
!--------------------------------------------------------------------------------------------------
subroutine thermal_conduction_putTemperatureAndItsRate(T,Tdot,ip,el)
integer, intent(in) :: &
ip, & !< integration point number
el !< element number
real(pReal), intent(in) :: &
T, &
Tdot
integer :: &
homog, &
offset
homog = material_homogenizationAt(el)
offset = material_homogenizationMemberAt(ip,el)
temperature (homog)%p(offset) = T
temperatureRate(homog)%p(offset) = Tdot
end subroutine thermal_conduction_putTemperatureAndItsRate
!--------------------------------------------------------------------------------------------------
!> @brief writes results to HDF5 output file
!--------------------------------------------------------------------------------------------------
subroutine thermal_conduction_results(homog,group)
integer, intent(in) :: homog
character(len=*), intent(in) :: group
integer :: o
associate(prm => param(damage_typeInstance(homog)))
outputsLoop: do o = 1,size(prm%output)
select case(trim(prm%output(o)))
case('T')
call results_writeDataset(group,temperature(homog)%p,'T',&
'temperature','K')
end select
enddo outputsLoop
end associate
end subroutine thermal_conduction_results
end module thermal_conduction

48
src/thermal_isothermal.f90
View File

@ -1,48 +0,0 @@
!--------------------------------------------------------------------------------------------------
!> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH
!> @brief material subroutine for isothermal temperature field
!--------------------------------------------------------------------------------------------------
module thermal_isothermal
use prec
use config
use material
use discretization
implicit none
public
contains
!--------------------------------------------------------------------------------------------------
!> @brief allocates fields, reads information from material configuration file
!--------------------------------------------------------------------------------------------------
subroutine thermal_isothermal_init(T)
real(pReal), dimension(:), intent(inout) :: T
integer :: Ninstances,Nmaterialpoints,ho,ip,el,ce
print'(/,a)', ' <<<+- thermal_isothermal init -+>>>'; flush(6)
do ho = 1, size(thermal_type)
if (thermal_type(ho) /= THERMAL_isothermal_ID) cycle
Nmaterialpoints = count(material_homogenizationAt == ho)
allocate(temperature (ho)%p(Nmaterialpoints),source=thermal_initialT(ho))
allocate(temperatureRate(ho)%p(Nmaterialpoints),source = 0.0_pReal)
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 module thermal_isothermal