takeover from 40_XX and 41_XX branch

easier to focus on thermal instead of doing all kinematics and sources
together
This commit is contained in:
Martin Diehl 2018-12-30 20:58:38 +01:00
parent 673eebd636
commit e0fa3e0b26
6 changed files with 76 additions and 134 deletions

View File

@ -188,7 +188,7 @@ subroutine constitutive_init()
call IO_checkAndRewind(FILEUNIT)
if (any(phase_kinematics == KINEMATICS_cleavage_opening_ID)) call kinematics_cleavage_opening_init(FILEUNIT)
if (any(phase_kinematics == KINEMATICS_slipplane_opening_ID)) call kinematics_slipplane_opening_init(FILEUNIT)
if (any(phase_kinematics == KINEMATICS_thermal_expansion_ID)) call kinematics_thermal_expansion_init(FILEUNIT)
if (any(phase_kinematics == KINEMATICS_thermal_expansion_ID)) call kinematics_thermal_expansion_init
if (any(phase_kinematics == KINEMATICS_vacancy_strain_ID)) call kinematics_vacancy_strain_init(FILEUNIT)
if (any(phase_kinematics == KINEMATICS_hydrogen_strain_ID)) call kinematics_hydrogen_strain_init(FILEUNIT)
close(FILEUNIT)

View File

@ -10,24 +10,7 @@ module kinematics_thermal_expansion
implicit none
private
integer(pInt), dimension(:), allocatable, public, protected :: &
kinematics_thermal_expansion_sizePostResults, & !< cumulative size of post results
kinematics_thermal_expansion_offset, & !< which kinematics is my current damage mechanism?
kinematics_thermal_expansion_instance !< instance of damage kinematics mechanism
integer(pInt), dimension(:,:), allocatable, target, public :: &
kinematics_thermal_expansion_sizePostResult !< size of each post result output
character(len=64), dimension(:,:), allocatable, target, public :: &
kinematics_thermal_expansion_output !< name of each post result output
integer(pInt), dimension(:), allocatable, target, public :: &
kinematics_thermal_expansion_Noutput !< number of outputs per instance of this damage
! enum, bind(c) ! ToDo kinematics need state machinery to deal with sizePostResult
! enumerator :: undefined_ID, & ! possible remedy is to decouple having state vars from having output
! thermalexpansionrate_ID ! which means to separate user-defined types tState + tOutput...
! end enum
public :: &
kinematics_thermal_expansion_init, &
kinematics_thermal_expansion_initialStrain, &
@ -40,7 +23,7 @@ contains
!> @brief module initialization
!> @details reads in material parameters, allocates arrays, and does sanity checks
!--------------------------------------------------------------------------------------------------
subroutine kinematics_thermal_expansion_init(fileUnit)
subroutine kinematics_thermal_expansion_init()
#if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800
use, intrinsic :: iso_fortran_env, only: &
compiler_version, &
@ -51,37 +34,17 @@ subroutine kinematics_thermal_expansion_init(fileUnit)
debug_constitutive,&
debug_levelBasic
use IO, only: &
IO_read, &
IO_lc, &
IO_getTag, &
IO_isBlank, &
IO_stringPos, &
IO_stringValue, &
IO_floatValue, &
IO_intValue, &
IO_warning, &
IO_error, &
IO_timeStamp, &
IO_EOF
IO_timeStamp
use material, only: &
phase_kinematics, &
phase_Nkinematics, &
phase_Noutput, &
KINEMATICS_thermal_expansion_label, &
KINEMATICS_thermal_expansion_ID
use config, only: &
material_Nphase, &
MATERIAL_partPhase
config_phase
implicit none
integer(pInt), intent(in) :: fileUnit
integer(pInt), allocatable, dimension(:) :: chunkPos
integer(pInt) :: maxNinstance,phase,instance,kinematics
character(len=65536) :: &
tag = '', &
line = ''
integer(pInt) maxNinstance
write(6,'(/,a)') ' <<<+- kinematics_'//KINEMATICS_thermal_expansion_LABEL//' init -+>>>'
write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
#include "compilation_info.f90"
@ -91,58 +54,8 @@ subroutine kinematics_thermal_expansion_init(fileUnit)
if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt) &
write(6,'(a16,1x,i5,/)') '# instances:',maxNinstance
allocate(kinematics_thermal_expansion_offset(material_Nphase), source=0_pInt)
allocate(kinematics_thermal_expansion_instance(material_Nphase), source=0_pInt)
do phase = 1, material_Nphase
kinematics_thermal_expansion_instance(phase) = count(phase_kinematics(:,1:phase) == kinematics_thermal_expansion_ID)
do kinematics = 1, phase_Nkinematics(phase)
if (phase_kinematics(kinematics,phase) == kinematics_thermal_expansion_ID) &
kinematics_thermal_expansion_offset(phase) = kinematics
enddo
enddo
allocate(kinematics_thermal_expansion_sizePostResults(maxNinstance), source=0_pInt)
allocate(kinematics_thermal_expansion_sizePostResult(maxval(phase_Noutput),maxNinstance),source=0_pInt)
allocate(kinematics_thermal_expansion_output(maxval(phase_Noutput),maxNinstance))
kinematics_thermal_expansion_output = ''
allocate(kinematics_thermal_expansion_Noutput(maxNinstance), source=0_pInt)
rewind(fileUnit)
phase = 0_pInt
do while (trim(line) /= IO_EOF .and. IO_lc(IO_getTag(line,'<','>')) /= MATERIAL_partPhase) ! wind forward to <phase>
line = IO_read(fileUnit)
enddo
parsingFile: do while (trim(line) /= IO_EOF) ! read through sections of phase part
line = IO_read(fileUnit)
if (IO_isBlank(line)) cycle ! skip empty lines
if (IO_getTag(line,'<','>') /= '') then ! stop at next part
line = IO_read(fileUnit, .true.) ! reset IO_read
exit
endif
if (IO_getTag(line,'[',']') /= '') then ! next phase section
phase = phase + 1_pInt ! advance phase section counter
cycle ! skip to next line
endif
if (phase > 0_pInt ) then; if (any(phase_kinematics(:,phase) == KINEMATICS_thermal_expansion_ID)) then ! do not short-circuit here (.and. with next if statemen). It's not safe in Fortran
instance = kinematics_thermal_expansion_instance(phase) ! which instance of my damage is present phase
chunkPos = IO_stringPos(line)
tag = IO_lc(IO_stringValue(line,chunkPos,1_pInt)) ! extract key...
select case(tag)
! case ('(output)')
! output = IO_lc(IO_stringValue(line,chunkPos,2_pInt)) ! ...and corresponding output
! select case(output)
! case ('thermalexpansionrate')
! kinematics_thermal_expansion_Noutput(instance) = kinematics_thermal_expansion_Noutput(instance) + 1_pInt
! kinematics_thermal_expansion_outputID(kinematics_thermal_expansion_Noutput(instance),instance) = &
! thermalexpansionrate_ID
! kinematics_thermal_expansion_output(kinematics_thermal_expansion_Noutput(instance),instance) = output
! ToDo add sizePostResult loop afterwards...
end select
endif; endif
enddo parsingFile
! ToDo: this subroutine should read in lattice_thermal_expansion. No need to make it a global array
end subroutine kinematics_thermal_expansion_init
@ -187,7 +100,7 @@ end function kinematics_thermal_expansion_initialStrain
!--------------------------------------------------------------------------------------------------
!> @brief contains the constitutive equation for calculating the velocity gradient
!--------------------------------------------------------------------------------------------------
subroutine kinematics_thermal_expansion_LiAndItsTangent(Li, dLi_dTstar3333, ipc, ip, el)
subroutine kinematics_thermal_expansion_LiAndItsTangent(Li, dLi_dTstar, ipc, ip, el)
use material, only: &
material_phase, &
material_homog, &
@ -206,7 +119,7 @@ subroutine kinematics_thermal_expansion_LiAndItsTangent(Li, dLi_dTstar3333, ipc,
real(pReal), intent(out), dimension(3,3) :: &
Li !< thermal velocity gradient
real(pReal), intent(out), dimension(3,3,3,3) :: &
dLi_dTstar3333 !< derivative of Li with respect to Tstar (4th-order tensor defined to be zero)
dLi_dTstar !< derivative of Li with respect to Tstar (4th-order tensor defined to be zero)
integer(pInt) :: &
phase, &
homog, offset
@ -230,7 +143,7 @@ subroutine kinematics_thermal_expansion_LiAndItsTangent(Li, dLi_dTstar3333, ipc,
+ lattice_thermalExpansion33(1:3,1:3,2,phase)*(T - TRef)**2 / 2. &
+ lattice_thermalExpansion33(1:3,1:3,3,phase)*(T - TRef)**3 / 3. &
)
dLi_dTstar3333 = 0.0_pReal
dLi_dTstar = 0.0_pReal
end subroutine kinematics_thermal_expansion_LiAndItsTangent

View File

@ -27,6 +27,15 @@ module source_thermal_dissipation
real(pReal), dimension(:), allocatable, private :: &
source_thermal_dissipation_coldworkCoeff
type, private :: tParameters !< container type for internal constitutive parameters
real(pReal) :: &
coldworkCoeff
end type tParameters
type(tParameters), dimension(:), allocatable, private :: param !< containers of constitutive parameters (len Ninstance)
public :: &
source_thermal_dissipation_init, &
source_thermal_dissipation_getRateAndItsTangent
@ -70,6 +79,7 @@ subroutine source_thermal_dissipation_init(fileUnit)
material_phase, &
sourceState
use config, only: &
config_phase, &
material_Nphase, &
MATERIAL_partPhase
use numerics,only: &
@ -79,9 +89,9 @@ subroutine source_thermal_dissipation_init(fileUnit)
integer(pInt), intent(in) :: fileUnit
integer(pInt), allocatable, dimension(:) :: chunkPos
integer(pInt) :: maxNinstance,phase,instance,source,sourceOffset
integer(pInt) :: Ninstance,phase,instance,source,sourceOffset
integer(pInt) :: sizeState, sizeDotState, sizeDeltaState
integer(pInt) :: NofMyPhase
integer(pInt) :: NofMyPhase,p
character(len=65536) :: &
tag = '', &
line = ''
@ -90,10 +100,10 @@ subroutine source_thermal_dissipation_init(fileUnit)
write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
#include "compilation_info.f90"
maxNinstance = int(count(phase_source == SOURCE_thermal_dissipation_ID),pInt)
if (maxNinstance == 0_pInt) return
Ninstance = int(count(phase_source == SOURCE_thermal_dissipation_ID),pInt)
if (Ninstance == 0_pInt) return
if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt) &
write(6,'(a16,1x,i5,/)') '# instances:',maxNinstance
write(6,'(a16,1x,i5,/)') '# instances:',Ninstance
allocate(source_thermal_dissipation_offset(material_Nphase), source=0_pInt)
allocate(source_thermal_dissipation_instance(material_Nphase), source=0_pInt)
@ -105,12 +115,17 @@ subroutine source_thermal_dissipation_init(fileUnit)
enddo
enddo
allocate(source_thermal_dissipation_sizePostResults(maxNinstance), source=0_pInt)
allocate(source_thermal_dissipation_sizePostResult(maxval(phase_Noutput),maxNinstance),source=0_pInt)
allocate(source_thermal_dissipation_output (maxval(phase_Noutput),maxNinstance))
allocate(source_thermal_dissipation_sizePostResults(Ninstance), source=0_pInt)
allocate(source_thermal_dissipation_sizePostResult(maxval(phase_Noutput),Ninstance),source=0_pInt)
allocate(source_thermal_dissipation_output (maxval(phase_Noutput),Ninstance))
source_thermal_dissipation_output = ''
allocate(source_thermal_dissipation_Noutput(maxNinstance), source=0_pInt)
allocate(source_thermal_dissipation_coldworkCoeff(maxNinstance), source=0.0_pReal)
allocate(source_thermal_dissipation_Noutput(Ninstance), source=0_pInt)
allocate(source_thermal_dissipation_coldworkCoeff(Ninstance), source=0.0_pReal)
do p=1, size(config_phase)
if (all(phase_source(:,p) /= SOURCE_thermal_dissipation_ID)) cycle
enddo
rewind(fileUnit)
phase = 0_pInt
@ -181,17 +196,14 @@ end subroutine source_thermal_dissipation_init
!--------------------------------------------------------------------------------------------------
!> @brief returns local vacancy generation rate
!--------------------------------------------------------------------------------------------------
subroutine source_thermal_dissipation_getRateAndItsTangent(TDot, dTDOT_dT, Tstar_v, Lp, ipc, ip, el)
subroutine source_thermal_dissipation_getRateAndItsTangent(TDot, dTDOT_dT, Tstar_v, Lp, phase, constituent)
use math, only: &
math_Mandel6to33
use material, only: &
phaseAt, phasememberAt
implicit none
integer(pInt), intent(in) :: &
ipc, & !< grain number
ip, & !< integration point number
el !< element number
phase, &
constituent
real(pReal), intent(in), dimension(6) :: &
Tstar_v !< 2nd Piola Kirchhoff stress tensor (Mandel)
real(pReal), intent(in), dimension(3,3) :: &
@ -200,10 +212,8 @@ subroutine source_thermal_dissipation_getRateAndItsTangent(TDot, dTDOT_dT, Tstar
TDot, &
dTDOT_dT
integer(pInt) :: &
instance, phase, constituent
instance
phase = phaseAt(ipc,ip,el)
constituent = phasememberAt(ipc,ip,el)
instance = source_thermal_dissipation_instance(phase)
TDot = source_thermal_dissipation_coldworkCoeff(instance)* &

View File

@ -32,6 +32,18 @@ module source_thermal_externalheat
source_thermal_externalheat_time, &
source_thermal_externalheat_rate
type, private :: tParameters !< container type for internal constitutive parameters
real(pReal), dimension(:), allocatable :: &
time, &
rate
integer(pInt) :: &
nInterval
end type tParameters
type(tParameters), dimension(:), allocatable, private :: param !< containers of constitutive parameters (len Ninstance)
public :: &
source_thermal_externalheat_init, &
source_thermal_externalheat_dotState, &
@ -76,6 +88,7 @@ subroutine source_thermal_externalheat_init(fileUnit)
material_phase, &
sourceState
use config, only: &
config_phase, &
material_Nphase, &
MATERIAL_partPhase
use numerics,only: &
@ -87,7 +100,7 @@ subroutine source_thermal_externalheat_init(fileUnit)
integer(pInt), allocatable, dimension(:) :: chunkPos
integer(pInt) :: maxNinstance,phase,instance,source,sourceOffset
integer(pInt) :: sizeState, sizeDotState, sizeDeltaState
integer(pInt) :: NofMyPhase,interval
integer(pInt) :: NofMyPhase,interval,p
character(len=65536) :: &
tag = '', &
line = ''
@ -117,11 +130,15 @@ subroutine source_thermal_externalheat_init(fileUnit)
allocate(source_thermal_externalheat_output (maxval(phase_Noutput),maxNinstance))
source_thermal_externalheat_output = ''
allocate(source_thermal_externalheat_Noutput(maxNinstance), source=0_pInt)
allocate(source_thermal_externalheat_nIntervals(maxNinstance), source=0_pInt)
allocate(source_thermal_externalheat_nIntervals(maxNinstance), source=0_pInt)
allocate(temp_time(maxNinstance,1000), source=0.0_pReal)
allocate(temp_rate(maxNinstance,1000), source=0.0_pReal)
do p=1, size(config_phase)
if (all(phase_source(:,p) /= SOURCE_thermal_externalheat_ID)) cycle
enddo
rewind(fileUnit)
phase = 0_pInt
do while (trim(line) /= IO_EOF .and. IO_lc(IO_getTag(line,'<','>')) /= MATERIAL_partPhase) ! wind forward to <phase>
@ -238,26 +255,22 @@ end subroutine source_thermal_externalheat_dotState
!--------------------------------------------------------------------------------------------------
!> @brief returns local heat generation rate
!--------------------------------------------------------------------------------------------------
subroutine source_thermal_externalheat_getRateAndItsTangent(TDot, dTDot_dT, ipc, ip, el)
subroutine source_thermal_externalheat_getRateAndItsTangent(TDot, dTDot_dT, phase, constituent)
use material, only: &
phaseAt, phasememberAt, &
sourceState
implicit none
integer(pInt), intent(in) :: &
ipc, & !< grain number
ip, & !< integration point number
el !< element number
phase, &
constituent
real(pReal), intent(out) :: &
TDot, &
dTDot_dT
integer(pInt) :: &
instance, phase, constituent, sourceOffset, interval
instance, sourceOffset, interval
real(pReal) :: &
frac_time
phase = phaseAt(ipc,ip,el)
constituent = phasememberAt(ipc,ip,el)
instance = source_thermal_externalheat_instance(phase)
sourceOffset = source_thermal_externalheat_offset(phase)

View File

@ -238,6 +238,7 @@ subroutine thermal_adiabatic_getSourceAndItsTangent(Tdot, dTdot_dT, T, ip, el)
homogenization_Ngrains, &
mappingHomogenization, &
phaseAt, &
phasememberAt, &
thermal_typeInstance, &
phase_Nsources, &
phase_source, &
@ -267,7 +268,8 @@ subroutine thermal_adiabatic_getSourceAndItsTangent(Tdot, dTdot_dT, T, ip, el)
offset, &
instance, &
grain, &
source
source, &
constituent
homog = mappingHomogenization(2,ip,el)
offset = mappingHomogenization(1,ip,el)
@ -277,17 +279,18 @@ subroutine thermal_adiabatic_getSourceAndItsTangent(Tdot, dTdot_dT, T, ip, el)
dTdot_dT = 0.0_pReal
do grain = 1, homogenization_Ngrains(homog)
phase = phaseAt(grain,ip,el)
constituent = phasememberAt(grain,ip,el)
do source = 1, phase_Nsources(phase)
select case(phase_source(source,phase))
case (SOURCE_thermal_dissipation_ID)
call source_thermal_dissipation_getRateAndItsTangent(my_Tdot, my_dTdot_dT, &
crystallite_Tstar_v(1:6,grain,ip,el), &
crystallite_Lp(1:3,1:3,grain,ip,el), &
grain, ip, el)
phase, constituent)
case (SOURCE_thermal_externalheat_ID)
call source_thermal_externalheat_getRateAndItsTangent(my_Tdot, my_dTdot_dT, &
grain, ip, el)
phase, constituent)
case default
my_Tdot = 0.0_pReal

View File

@ -192,6 +192,7 @@ subroutine thermal_conduction_getSourceAndItsTangent(Tdot, dTdot_dT, T, ip, el)
homogenization_Ngrains, &
mappingHomogenization, &
phaseAt, &
phasememberAt, &
thermal_typeInstance, &
phase_Nsources, &
phase_source, &
@ -221,7 +222,8 @@ subroutine thermal_conduction_getSourceAndItsTangent(Tdot, dTdot_dT, T, ip, el)
offset, &
instance, &
grain, &
source
source, &
constituent
homog = mappingHomogenization(2,ip,el)
offset = mappingHomogenization(1,ip,el)
@ -231,17 +233,18 @@ subroutine thermal_conduction_getSourceAndItsTangent(Tdot, dTdot_dT, T, ip, el)
dTdot_dT = 0.0_pReal
do grain = 1, homogenization_Ngrains(homog)
phase = phaseAt(grain,ip,el)
constituent = phasememberAt(grain,ip,el)
do source = 1, phase_Nsources(phase)
select case(phase_source(source,phase))
case (SOURCE_thermal_dissipation_ID)
call source_thermal_dissipation_getRateAndItsTangent(my_Tdot, my_dTdot_dT, &
crystallite_Tstar_v(1:6,grain,ip,el), &
crystallite_Lp(1:3,1:3,grain,ip,el), &
grain, ip, el)
phase, constituent)
case (SOURCE_thermal_externalheat_ID)
call source_thermal_externalheat_getRateAndItsTangent(my_Tdot, my_dTdot_dT, &
grain, ip, el)
phase, constituent)
case default
my_Tdot = 0.0_pReal
@ -363,8 +366,8 @@ function thermal_conduction_getMassDensity(ip,el)
homog = mappingHomogenization(2,ip,el)
do grain = 1, homogenization_Ngrains(mesh_element(3,el))
thermal_conduction_getMassDensity = thermal_conduction_getMassDensity + &
lattice_massDensity(material_phase(grain,ip,el))
thermal_conduction_getMassDensity = thermal_conduction_getMassDensity &
+ lattice_massDensity(material_phase(grain,ip,el))
enddo
thermal_conduction_getMassDensity = &