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

Merge branch 'New-Thermal' into 'development' 

New thermal

See merge request damask/DAMASK!61
This commit is contained in:
Franz Roters 2019-02-26 16:45:07 +01:00
parent c0757fe930 56e2c1264b
commit 07eff8eb7b
11 changed files with 460 additions and 903 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

76
src/constitutive.f90
View File

@ -38,11 +38,6 @@ contains
!> @brief allocates arrays pointing to array of the various constitutive modules !> @brief allocates arrays pointing to array of the various constitutive modules
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine constitutive_init() subroutine constitutive_init()
#if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800
use, intrinsic :: iso_fortran_env, only: &
compiler_version, &
compiler_options
#endif
use prec, only: & use prec, only: &
pReal pReal
use debug, only: & use debug, only: &
@ -53,12 +48,8 @@ subroutine constitutive_init()
use IO, only: & use IO, only: &
IO_error, & IO_error, &
IO_open_file, & IO_open_file, &
IO_checkAndRewind, &
IO_open_jobFile_stat, & IO_open_jobFile_stat, &
IO_write_jobFile, & IO_write_jobFile
IO_timeStamp
use config, only: &
config_phase
use config, only: & use config, only: &
material_Nphase, & material_Nphase, &
material_localFileExt, & material_localFileExt, &
@ -138,7 +129,7 @@ subroutine constitutive_init()
nonlocalConstitutionPresent = .false. nonlocalConstitutionPresent = .false.
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! parse plasticities from config file ! initialized plasticity
if (any(phase_plasticity == PLASTICITY_NONE_ID)) call plastic_none_init if (any(phase_plasticity == PLASTICITY_NONE_ID)) call plastic_none_init
if (any(phase_plasticity == PLASTICITY_ISOTROPIC_ID)) call plastic_isotropic_init if (any(phase_plasticity == PLASTICITY_ISOTROPIC_ID)) call plastic_isotropic_init
if (any(phase_plasticity == PLASTICITY_PHENOPOWERLAW_ID)) call plastic_phenopowerlaw_init if (any(phase_plasticity == PLASTICITY_PHENOPOWERLAW_ID)) call plastic_phenopowerlaw_init
@ -147,28 +138,20 @@ subroutine constitutive_init()
if (any(phase_plasticity == PLASTICITY_DISLOUCLA_ID)) call plastic_disloucla_init if (any(phase_plasticity == PLASTICITY_DISLOUCLA_ID)) call plastic_disloucla_init
if (any(phase_plasticity == PLASTICITY_NONLOCAL_ID)) call plastic_nonlocal_init if (any(phase_plasticity == PLASTICITY_NONLOCAL_ID)) call plastic_nonlocal_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! open material.config ! initialize source mechanisms
if (.not. IO_open_jobFile_stat(FILEUNIT,material_localFileExt)) & ! no local material configuration present... if (any(phase_source == SOURCE_thermal_dissipation_ID)) call source_thermal_dissipation_init
call IO_open_file(FILEUNIT,material_configFile) ! ... open material.config file if (any(phase_source == SOURCE_thermal_externalheat_ID)) call source_thermal_externalheat_init
!--------------------------------------------------------------------------------------------------
! parse source mechanisms from config file
if (any(phase_source == SOURCE_thermal_dissipation_ID)) call source_thermal_dissipation_init(FILEUNIT)
if (any(phase_source == SOURCE_thermal_externalheat_ID)) call source_thermal_externalheat_init(FILEUNIT)
if (any(phase_source == SOURCE_damage_isoBrittle_ID)) call source_damage_isoBrittle_init if (any(phase_source == SOURCE_damage_isoBrittle_ID)) call source_damage_isoBrittle_init
if (any(phase_source == SOURCE_damage_isoDuctile_ID)) call source_damage_isoDuctile_init if (any(phase_source == SOURCE_damage_isoDuctile_ID)) call source_damage_isoDuctile_init
if (any(phase_source == SOURCE_damage_anisoBrittle_ID)) call source_damage_anisoBrittle_init if (any(phase_source == SOURCE_damage_anisoBrittle_ID)) call source_damage_anisoBrittle_init
if (any(phase_source == SOURCE_damage_anisoDuctile_ID)) call source_damage_anisoDuctile_init if (any(phase_source == SOURCE_damage_anisoDuctile_ID)) call source_damage_anisoDuctile_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! parse kinematic mechanisms from config file ! initialize kinematic mechanisms
call IO_checkAndRewind(FILEUNIT)
if (any(phase_kinematics == KINEMATICS_cleavage_opening_ID)) call kinematics_cleavage_opening_init if (any(phase_kinematics == KINEMATICS_cleavage_opening_ID)) call kinematics_cleavage_opening_init
if (any(phase_kinematics == KINEMATICS_slipplane_opening_ID)) call kinematics_slipplane_opening_init if (any(phase_kinematics == KINEMATICS_slipplane_opening_ID)) call kinematics_slipplane_opening_init
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
close(FILEUNIT)
call config_deallocate('material.config/phase') call config_deallocate('material.config/phase')
@ -481,7 +464,7 @@ subroutine constitutive_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, &
case (PLASTICITY_KINEHARDENING_ID) plasticityType case (PLASTICITY_KINEHARDENING_ID) plasticityType
of = phasememberAt(ipc,ip,el) of = phasememberAt(ipc,ip,el)
instance = phase_plasticityInstance(material_phase(ipc,ip,el)) instance = phase_plasticityInstance(material_phase(ipc,ip,el))
call plastic_kinehardening_LpAndItsTangent (Lp,dLp_dMp,Mp,instance,of) call plastic_kinehardening_LpAndItsTangent (Lp,dLp_dMp, Mp,instance,of)
case (PLASTICITY_NONLOCAL_ID) plasticityType case (PLASTICITY_NONLOCAL_ID) plasticityType
call plastic_nonlocal_LpAndItsTangent (Lp,dLp_dMp,Mp, & call plastic_nonlocal_LpAndItsTangent (Lp,dLp_dMp,Mp, &
@ -528,8 +511,7 @@ subroutine constitutive_LiAndItsTangents(Li, dLi_dS, dLi_dFi, &
math_I3, & math_I3, &
math_inv33, & math_inv33, &
math_det33, & math_det33, &
math_mul33x33, & math_mul33x33
math_6toSym33
use material, only: & use material, only: &
phasememberAt, & phasememberAt, &
phase_plasticity, & phase_plasticity, &
@ -632,10 +614,11 @@ pure function constitutive_initialFi(ipc, ip, el)
use prec, only: & use prec, only: &
pReal pReal
use math, only: & use math, only: &
math_I3, & math_I3
math_inv33, &
math_mul33x33
use material, only: & use material, only: &
material_phase, &
material_homog, &
thermalMapping, &
phase_kinematics, & phase_kinematics, &
phase_Nkinematics, & phase_Nkinematics, &
material_phase, & material_phase, &
@ -652,14 +635,20 @@ pure function constitutive_initialFi(ipc, ip, el)
constitutive_initialFi !< composite initial intermediate deformation gradient constitutive_initialFi !< composite initial intermediate deformation gradient
integer(pInt) :: & integer(pInt) :: &
k !< counter in kinematics loop k !< counter in kinematics loop
integer(pInt) :: &
phase, &
homog, offset
constitutive_initialFi = math_I3 constitutive_initialFi = math_I3
phase = material_phase(ipc,ip,el)
KinematicsLoop: do k = 1_pInt, phase_Nkinematics(material_phase(ipc,ip,el)) !< Warning: small initial strain assumption KinematicsLoop: do k = 1_pInt, phase_Nkinematics(phase) !< Warning: small initial strain assumption
kinematicsType: select case (phase_kinematics(k,material_phase(ipc,ip,el))) kinematicsType: select case (phase_kinematics(k,phase))
case (KINEMATICS_thermal_expansion_ID) kinematicsType case (KINEMATICS_thermal_expansion_ID) kinematicsType
homog = material_homog(ip,el)
offset = thermalMapping(homog)%p(ip,el)
constitutive_initialFi = & constitutive_initialFi = &
constitutive_initialFi + kinematics_thermal_expansion_initialStrain(ipc, ip, el) constitutive_initialFi + kinematics_thermal_expansion_initialStrain(homog,phase,offset)
end select kinematicsType end select kinematicsType
enddo KinematicsLoop enddo KinematicsLoop
@ -764,7 +753,7 @@ end subroutine constitutive_hooke_SandItsTangents
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief contains the constitutive equation for calculating the rate of change of microstructure !> @brief contains the constitutive equation for calculating the rate of change of microstructure
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine constitutive_collectDotState(S, FeArray, Fi, FpArray, subdt, subfracArray,ipc, ip, el) subroutine constitutive_collectDotState(S, FeArray, Fi, FpArray, subdt, ipc, ip, el)
use prec, only: & use prec, only: &
pReal, & pReal, &
pLongInt pLongInt
@ -774,8 +763,6 @@ subroutine constitutive_collectDotState(S, FeArray, Fi, FpArray, subdt, subfracA
debug_levelBasic debug_levelBasic
use math, only: & use math, only: &
math_mul33x33, & math_mul33x33, &
math_6toSym33, &
math_sym33to6, &
math_mul33x33 math_mul33x33
use mesh, only: & use mesh, only: &
theMesh theMesh
@ -829,8 +816,6 @@ subroutine constitutive_collectDotState(S, FeArray, Fi, FpArray, subdt, subfracA
el !< element el !< element
real(pReal), intent(in) :: & real(pReal), intent(in) :: &
subdt !< timestep subdt !< timestep
real(pReal), intent(in), dimension(homogenization_maxNgrains,theMesh%elem%nIPs,theMesh%Nelems) :: &
subfracArray !< subfraction of timestep
real(pReal), intent(in), dimension(3,3,homogenization_maxNgrains,theMesh%elem%nIPs,theMesh%Nelems) :: & real(pReal), intent(in), dimension(3,3,homogenization_maxNgrains,theMesh%elem%nIPs,theMesh%Nelems) :: &
FeArray, & !< elastic deformation gradient FeArray, & !< elastic deformation gradient
FpArray !< plastic deformation gradient FpArray !< plastic deformation gradient
@ -891,13 +876,13 @@ subroutine constitutive_collectDotState(S, FeArray, Fi, FpArray, subdt, subfracA
call source_damage_anisoBrittle_dotState (S, ipc, ip, el) !< correct stress? call source_damage_anisoBrittle_dotState (S, ipc, ip, el) !< correct stress?
case (SOURCE_damage_isoDuctile_ID) sourceType case (SOURCE_damage_isoDuctile_ID) sourceType
call source_damage_isoDuctile_dotState ( ipc, ip, el) call source_damage_isoDuctile_dotState ( ipc, ip, el)
case (SOURCE_damage_anisoDuctile_ID) sourceType case (SOURCE_damage_anisoDuctile_ID) sourceType
call source_damage_anisoDuctile_dotState ( ipc, ip, el) call source_damage_anisoDuctile_dotState ( ipc, ip, el)
case (SOURCE_thermal_externalheat_ID) sourceType case (SOURCE_thermal_externalheat_ID) sourceType
call source_thermal_externalheat_dotState( ipc, ip, el) call source_thermal_externalheat_dotState( ipc, ip, el)
end select sourceType end select sourceType
@ -918,7 +903,6 @@ subroutine constitutive_collectDeltaState(S, Fe, Fi, ipc, ip, el)
debug_constitutive, & debug_constitutive, &
debug_levelBasic debug_levelBasic
use math, only: & use math, only: &
math_sym33to6, &
math_mul33x33 math_mul33x33
use material, only: & use material, only: &
phasememberAt, & phasememberAt, &
@ -984,14 +968,11 @@ end subroutine constitutive_collectDeltaState
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief returns array of constitutive results !> @brief returns array of constitutive results
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function constitutive_postResults(S, Fi, FeArray, ipc, ip, el) function constitutive_postResults(S, Fi, ipc, ip, el)
use prec, only: & use prec, only: &
pReal pReal
use math, only: & use math, only: &
math_6toSym33, &
math_mul33x33 math_mul33x33
use mesh, only: &
theMesh
use material, only: & use material, only: &
phasememberAt, & phasememberAt, &
phase_plasticityInstance, & phase_plasticityInstance, &
@ -1004,7 +985,6 @@ function constitutive_postResults(S, Fi, FeArray, ipc, ip, el)
material_homogenizationAt, & material_homogenizationAt, &
temperature, & temperature, &
thermalMapping, & thermalMapping, &
homogenization_maxNgrains, &
PLASTICITY_NONE_ID, & PLASTICITY_NONE_ID, &
PLASTICITY_ISOTROPIC_ID, & PLASTICITY_ISOTROPIC_ID, &
PLASTICITY_PHENOPOWERLAW_ID, & PLASTICITY_PHENOPOWERLAW_ID, &
@ -1047,8 +1027,6 @@ function constitutive_postResults(S, Fi, FeArray, ipc, ip, el)
constitutive_postResults constitutive_postResults
real(pReal), intent(in), dimension(3,3) :: & real(pReal), intent(in), dimension(3,3) :: &
Fi !< intermediate deformation gradient Fi !< intermediate deformation gradient
real(pReal), intent(in), dimension(3,3,homogenization_maxNgrains,theMesh%elem%nIPs,theMesh%Nelems) :: &
FeArray !< elastic deformation gradient
real(pReal), intent(in), dimension(3,3) :: & real(pReal), intent(in), dimension(3,3) :: &
S !< 2nd Piola Kirchhoff stress S !< 2nd Piola Kirchhoff stress
real(pReal), dimension(3,3) :: & real(pReal), dimension(3,3) :: &

6
src/crystallite.f90
View File

@ -247,7 +247,7 @@ subroutine crystallite_init
allocate(crystallite_sizePostResult(maxval(crystallite_Noutput), & allocate(crystallite_sizePostResult(maxval(crystallite_Noutput), &
size(config_crystallite)), source=0_pInt) size(config_crystallite)), source=0_pInt)
select case(numerics_integrator(1)) select case(numerics_integrator)
case(1_pInt) case(1_pInt)
integrateState => integrateStateFPI integrateState => integrateStateFPI
case(2_pInt) case(2_pInt)
@ -1068,7 +1068,7 @@ function crystallite_postResults(ipc, ip, el)
if (size(crystallite_postResults)-c > 0_pInt) & if (size(crystallite_postResults)-c > 0_pInt) &
crystallite_postResults(c+1:size(crystallite_postResults)) = & crystallite_postResults(c+1:size(crystallite_postResults)) = &
constitutive_postResults(math_6toSym33(crystallite_Tstar_v(1:6,ipc,ip,el)), crystallite_Fi(1:3,1:3,ipc,ip,el), & constitutive_postResults(math_6toSym33(crystallite_Tstar_v(1:6,ipc,ip,el)), crystallite_Fi(1:3,1:3,ipc,ip,el), &
crystallite_Fe, ipc, ip, el) ipc, ip, el)
end function crystallite_postResults end function crystallite_postResults
@ -2297,7 +2297,7 @@ subroutine update_dotState(timeFraction)
crystallite_Fe, & crystallite_Fe, &
crystallite_Fi(1:3,1:3,g,i,e), & crystallite_Fi(1:3,1:3,g,i,e), &
crystallite_Fp, & crystallite_Fp, &
crystallite_subdt(g,i,e)*timeFraction, crystallite_subFrac, g,i,e) crystallite_subdt(g,i,e)*timeFraction, g,i,e)
p = phaseAt(g,i,e); c = phasememberAt(g,i,e) p = phaseAt(g,i,e); c = phasememberAt(g,i,e)
NaN = any(IEEE_is_NaN(plasticState(p)%dotState(:,c))) NaN = any(IEEE_is_NaN(plasticState(p)%dotState(:,c)))
do s = 1_pInt, phase_Nsources(p) do s = 1_pInt, phase_Nsources(p)

36
src/homogenization.f90
View File

@ -108,30 +108,18 @@ subroutine homogenization_init
logical :: valid logical :: valid
!--------------------------------------------------------------------------------------------------
! open material.config
if (.not. IO_open_jobFile_stat(FILEUNIT,material_localFileExt)) & ! no local material configuration present...
call IO_open_file(FILEUNIT,material_configFile) ! ... open material.config file
!--------------------------------------------------------------------------------------------------
! parse homogenization from config file
if (any(homogenization_type == HOMOGENIZATION_NONE_ID)) call homogenization_none_init if (any(homogenization_type == HOMOGENIZATION_NONE_ID)) call homogenization_none_init
if (any(homogenization_type == HOMOGENIZATION_ISOSTRAIN_ID)) call homogenization_isostrain_init if (any(homogenization_type == HOMOGENIZATION_ISOSTRAIN_ID)) call homogenization_isostrain_init
if (any(homogenization_type == HOMOGENIZATION_RGC_ID)) call homogenization_RGC_init if (any(homogenization_type == HOMOGENIZATION_RGC_ID)) call homogenization_RGC_init
!-------------------------------------------------------------------------------------------------- if (any(thermal_type == THERMAL_isothermal_ID)) call thermal_isothermal_init
! parse thermal from config file if (any(thermal_type == THERMAL_adiabatic_ID)) call thermal_adiabatic_init
call IO_checkAndRewind(FILEUNIT) if (any(thermal_type == THERMAL_conduction_ID)) call thermal_conduction_init
if (any(thermal_type == THERMAL_isothermal_ID)) &
call thermal_isothermal_init()
if (any(thermal_type == THERMAL_adiabatic_ID)) &
call thermal_adiabatic_init(FILEUNIT)
if (any(thermal_type == THERMAL_conduction_ID)) &
call thermal_conduction_init(FILEUNIT)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! parse damage from config file ! open material.config
call IO_checkAndRewind(FILEUNIT) if (.not. IO_open_jobFile_stat(FILEUNIT,material_localFileExt)) & ! no local material configuration present...
call IO_open_file(FILEUNIT,material_configFile) ! ... open material.config file
if (any(damage_type == DAMAGE_none_ID)) & if (any(damage_type == DAMAGE_none_ID)) &
call damage_none_init() call damage_none_init()
if (any(damage_type == DAMAGE_local_ID)) & if (any(damage_type == DAMAGE_local_ID)) &
@ -895,6 +883,8 @@ function postResults(ip,el)
use mesh, only: & use mesh, only: &
mesh_element mesh_element
use material, only: & use material, only: &
thermalMapping, &
thermal_typeInstance, &
material_homogenizationAt, & material_homogenizationAt, &
homogenization_typeInstance,& homogenization_typeInstance,&
mappingHomogenization, & mappingHomogenization, &
@ -934,7 +924,7 @@ function postResults(ip,el)
postResults postResults
integer(pInt) :: & integer(pInt) :: &
startPos, endPos ,& startPos, endPos ,&
of, instance of, instance, homog
postResults = 0.0_pReal postResults = 0.0_pReal
@ -954,9 +944,13 @@ function postResults(ip,el)
chosenThermal: select case (thermal_type(mesh_element(3,el))) chosenThermal: select case (thermal_type(mesh_element(3,el)))
case (THERMAL_adiabatic_ID) chosenThermal case (THERMAL_adiabatic_ID) chosenThermal
postResults(startPos:endPos) = thermal_adiabatic_postResults(ip, el) homog = mappingHomogenization(2,ip,el)
postResults(startPos:endPos) = &
thermal_adiabatic_postResults(homog,thermal_typeInstance(homog),thermalMapping(homog)%p(ip,el))
case (THERMAL_conduction_ID) chosenThermal case (THERMAL_conduction_ID) chosenThermal
postResults(startPos:endPos) = thermal_conduction_postResults(ip, el) homog = mappingHomogenization(2,ip,el)
postResults(startPos:endPos) = &
thermal_conduction_postResults(homog,thermal_typeInstance(homog),thermalMapping(homog)%p(ip,el))
end select chosenThermal end select chosenThermal

300
src/kinematics_thermal_expansion.f90
View File

@ -4,34 +4,24 @@
!> @details to be done !> @details to be done
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module kinematics_thermal_expansion module kinematics_thermal_expansion
use prec, only: & use prec, only: &
pReal, & pReal, &
pInt pInt
implicit none implicit none
private 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 :: & type, private :: tParameters
kinematics_thermal_expansion_sizePostResult !< size of each post result output real(pReal), allocatable, dimension(:,:,:) :: &
expansion
end type tParameters
character(len=64), dimension(:,:), allocatable, target, public :: & type(tParameters), dimension(:), allocatable :: param
kinematics_thermal_expansion_output !< name of each post result output
integer(pInt), dimension(:), allocatable, target, public :: & public :: &
kinematics_thermal_expansion_Noutput !< number of outputs per instance of this damage kinematics_thermal_expansion_init, &
kinematics_thermal_expansion_initialStrain, &
enum, bind(c) ! ToDo kinematics need state machinery to deal with sizePostResult kinematics_thermal_expansion_LiAndItsTangent
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, &
kinematics_thermal_expansion_LiAndItsTangent
contains contains
@ -40,197 +30,129 @@ 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 kinematics_thermal_expansion_init(fileUnit) subroutine kinematics_thermal_expansion_init()
#if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800 use debug, only: &
use, intrinsic :: iso_fortran_env, only: & debug_level,&
compiler_version, & debug_constitutive,&
compiler_options debug_levelBasic
#endif use material, only: &
use debug, only: & phase_kinematics, &
debug_level,& KINEMATICS_thermal_expansion_label, &
debug_constitutive,& KINEMATICS_thermal_expansion_ID
debug_levelBasic use config, only: &
use IO, only: & config_phase
IO_read, &
IO_lc, &
IO_getTag, &
IO_isBlank, &
IO_stringPos, &
IO_stringValue, &
IO_floatValue, &
IO_intValue, &
IO_warning, &
IO_error, &
IO_timeStamp, &
IO_EOF
use material, only: &
phase_kinematics, &
phase_Nkinematics, &
phase_Noutput, &
KINEMATICS_thermal_expansion_label, &
KINEMATICS_thermal_expansion_ID
use config, only: &
material_Nphase, &
MATERIAL_partPhase
implicit none implicit none
integer(pInt), intent(in) :: fileUnit integer(pInt) :: &
Ninstance, &
p, i
real(pReal), dimension(:), allocatable :: &
temp
integer(pInt), allocatable, dimension(:) :: chunkPos write(6,'(/,a)') ' <<<+- kinematics_'//KINEMATICS_thermal_expansion_LABEL//' init -+>>>'
integer(pInt) :: maxNinstance,phase,instance,kinematics
character(len=65536) :: &
tag = '', &
line = ''
write(6,'(/,a)') ' <<<+- kinematics_'//KINEMATICS_thermal_expansion_LABEL//' init -+>>>' Ninstance = int(count(phase_kinematics == KINEMATICS_thermal_expansion_ID),pInt)
write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
#include "compilation_info.f90"
maxNinstance = int(count(phase_kinematics == KINEMATICS_thermal_expansion_ID),pInt) if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt) &
if (maxNinstance == 0_pInt) return write(6,'(a16,1x,i5,/)') '# instances:',Ninstance
if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt) & allocate(param(Ninstance))
write(6,'(a16,1x,i5,/)') '# instances:',maxNinstance
allocate(kinematics_thermal_expansion_offset(material_Nphase), source=0_pInt) do p = 1_pInt, size(phase_kinematics)
allocate(kinematics_thermal_expansion_instance(material_Nphase), source=0_pInt) if (all(phase_kinematics(:,p) /= KINEMATICS_thermal_expansion_ID)) cycle
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) ! ToDo: Here we need to decide how to extend the concept of instances to
allocate(kinematics_thermal_expansion_sizePostResult(maxval(phase_Noutput),maxNinstance),source=0_pInt) ! kinetics and sources. I would suggest that the same mechanism exists at maximum once per phase
allocate(kinematics_thermal_expansion_output(maxval(phase_Noutput),maxNinstance))
kinematics_thermal_expansion_output = ''
allocate(kinematics_thermal_expansion_Noutput(maxNinstance), source=0_pInt)
rewind(fileUnit) ! read up to three parameters (constant, linear, quadratic with T)
phase = 0_pInt temp = config_phase(p)%getFloats('thermal_expansion11')
do while (trim(line) /= IO_EOF .and. IO_lc(IO_getTag(line,'<','>')) /= MATERIAL_partPhase) ! wind forward to <phase> !lattice_thermalExpansion33(1,1,1:size(temp),p) = temp
line = IO_read(fileUnit) temp = config_phase(p)%getFloats('thermal_expansion22', &
enddo defaultVal=[(0.0_pReal, i=1,size(temp))],requiredSize=size(temp))
!lattice_thermalExpansion33(2,2,1:size(temp),p) = temp
parsingFile: do while (trim(line) /= IO_EOF) ! read through sections of phase part temp = config_phase(p)%getFloats('thermal_expansion33', &
line = IO_read(fileUnit) defaultVal=[(0.0_pReal, i=1,size(temp))],requiredSize=size(temp))
if (IO_isBlank(line)) cycle ! skip empty lines enddo
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
end subroutine kinematics_thermal_expansion_init end subroutine kinematics_thermal_expansion_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief report initial thermal strain based on current temperature deviation from reference !> @brief report initial thermal strain based on current temperature deviation from reference
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
pure function kinematics_thermal_expansion_initialStrain(ipc, ip, el) pure function kinematics_thermal_expansion_initialStrain(homog,phase,offset)
use material, only: & use material, only: &
material_phase, & temperature
material_homog, & use lattice, only: &
temperature, & lattice_thermalExpansion33, &
thermalMapping lattice_referenceTemperature
use lattice, only: &
lattice_thermalExpansion33, &
lattice_referenceTemperature
implicit none implicit none
integer(pInt), intent(in) :: & integer(pInt), intent(in) :: &
ipc, & !< grain number phase, &
ip, & !< integration point number homog, offset
el !< element number real(pReal), dimension(3,3) :: &
real(pReal), dimension(3,3) :: & kinematics_thermal_expansion_initialStrain !< initial thermal strain (should be small strain, though)
kinematics_thermal_expansion_initialStrain !< initial thermal strain (should be small strain, though)
integer(pInt) :: &
phase, &
homog, offset
phase = material_phase(ipc,ip,el)
homog = material_homog(ip,el)
offset = thermalMapping(homog)%p(ip,el)
kinematics_thermal_expansion_initialStrain = & kinematics_thermal_expansion_initialStrain = &
(temperature(homog)%p(offset) - lattice_referenceTemperature(phase))**1 / 1. * & (temperature(homog)%p(offset) - lattice_referenceTemperature(phase))**1 / 1. * &
lattice_thermalExpansion33(1:3,1:3,1,phase) + & ! constant coefficient lattice_thermalExpansion33(1:3,1:3,1,phase) + & ! constant coefficient
(temperature(homog)%p(offset) - lattice_referenceTemperature(phase))**2 / 2. * & (temperature(homog)%p(offset) - lattice_referenceTemperature(phase))**2 / 2. * &
lattice_thermalExpansion33(1:3,1:3,2,phase) + & ! linear coefficient lattice_thermalExpansion33(1:3,1:3,2,phase) + & ! linear coefficient
(temperature(homog)%p(offset) - lattice_referenceTemperature(phase))**3 / 3. * & (temperature(homog)%p(offset) - lattice_referenceTemperature(phase))**3 / 3. * &
lattice_thermalExpansion33(1:3,1:3,3,phase) ! quadratic coefficient lattice_thermalExpansion33(1:3,1:3,3,phase) ! quadratic coefficient
end function kinematics_thermal_expansion_initialStrain end function kinematics_thermal_expansion_initialStrain
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief contains the constitutive equation for calculating the velocity gradient !> @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: & use material, only: &
material_phase, & material_phase, &
material_homog, & material_homog, &
temperature, & temperature, &
temperatureRate, & temperatureRate, &
thermalMapping thermalMapping
use lattice, only: & use lattice, only: &
lattice_thermalExpansion33, & lattice_thermalExpansion33, &
lattice_referenceTemperature lattice_referenceTemperature
implicit none implicit none
integer(pInt), intent(in) :: & integer(pInt), intent(in) :: &
ipc, & !< grain number ipc, & !< grain number
ip, & !< integration point number ip, & !< integration point number
el !< element number el !< element number
real(pReal), intent(out), dimension(3,3) :: & real(pReal), intent(out), dimension(3,3) :: &
Li !< thermal velocity gradient Li !< thermal velocity gradient
real(pReal), intent(out), dimension(3,3,3,3) :: & 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) :: & integer(pInt) :: &
phase, & phase, &
homog, offset homog, offset
real(pReal) :: & real(pReal) :: &
T, TRef, TDot T, TRef, TDot
phase = material_phase(ipc,ip,el) phase = material_phase(ipc,ip,el)
homog = material_homog(ip,el) homog = material_homog(ip,el)
offset = thermalMapping(homog)%p(ip,el) offset = thermalMapping(homog)%p(ip,el)
T = temperature(homog)%p(offset) T = temperature(homog)%p(offset)
TDot = temperatureRate(homog)%p(offset) TDot = temperatureRate(homog)%p(offset)
TRef = lattice_referenceTemperature(phase) TRef = lattice_referenceTemperature(phase)
Li = TDot * ( & Li = TDot * ( &
lattice_thermalExpansion33(1:3,1:3,1,phase)*(T - TRef)**0 & ! constant coefficient lattice_thermalExpansion33(1:3,1:3,1,phase)*(T - TRef)**0 & ! constant coefficient
+ lattice_thermalExpansion33(1:3,1:3,2,phase)*(T - TRef)**1 & ! linear coefficient + lattice_thermalExpansion33(1:3,1:3,2,phase)*(T - TRef)**1 & ! linear coefficient
+ lattice_thermalExpansion33(1:3,1:3,3,phase)*(T - TRef)**2 & ! quadratic coefficient + lattice_thermalExpansion33(1:3,1:3,3,phase)*(T - TRef)**2 & ! quadratic coefficient
) / & ) / &
(1.0_pReal & (1.0_pReal &
+ lattice_thermalExpansion33(1:3,1:3,1,phase)*(T - TRef)**1 / 1. & + lattice_thermalExpansion33(1:3,1:3,1,phase)*(T - TRef)**1 / 1. &
+ lattice_thermalExpansion33(1:3,1:3,2,phase)*(T - TRef)**2 / 2. & + lattice_thermalExpansion33(1:3,1:3,2,phase)*(T - TRef)**2 / 2. &
+ lattice_thermalExpansion33(1:3,1:3,3,phase)*(T - TRef)**3 / 3. & + 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 end subroutine kinematics_thermal_expansion_LiAndItsTangent

10
src/material.f90
View File

@ -921,7 +921,7 @@ subroutine material_allocatePlasticState(phase,NofMyPhase,&
sizeState,sizeDotState,sizeDeltaState,& sizeState,sizeDotState,sizeDeltaState,&
Nslip,Ntwin,Ntrans) Nslip,Ntwin,Ntrans)
use numerics, only: & use numerics, only: &
numerics_integrator2 => numerics_integrator ! compatibility hack numerics_integrator
implicit none implicit none
integer(pInt), intent(in) :: & integer(pInt), intent(in) :: &
@ -933,8 +933,6 @@ subroutine material_allocatePlasticState(phase,NofMyPhase,&
Nslip, & Nslip, &
Ntwin, & Ntwin, &
Ntrans Ntrans
integer(pInt) :: numerics_integrator ! compatibility hack
numerics_integrator = numerics_integrator2(1) ! compatibility hack
plasticState(phase)%sizeState = sizeState plasticState(phase)%sizeState = sizeState
plasticState(phase)%sizeDotState = sizeDotState plasticState(phase)%sizeDotState = sizeDotState
@ -971,7 +969,7 @@ end subroutine material_allocatePlasticState
subroutine material_allocateSourceState(phase,of,NofMyPhase,& subroutine material_allocateSourceState(phase,of,NofMyPhase,&
sizeState,sizeDotState,sizeDeltaState) sizeState,sizeDotState,sizeDeltaState)
use numerics, only: & use numerics, only: &
numerics_integrator2 => numerics_integrator ! compatibility hack numerics_integrator
implicit none implicit none
integer(pInt), intent(in) :: & integer(pInt), intent(in) :: &
@ -979,8 +977,6 @@ subroutine material_allocateSourceState(phase,of,NofMyPhase,&
of, & of, &
NofMyPhase, & NofMyPhase, &
sizeState, sizeDotState,sizeDeltaState sizeState, sizeDotState,sizeDeltaState
integer(pInt) :: numerics_integrator ! compatibility hack
numerics_integrator = numerics_integrator2(1) ! compatibility hack
sourceState(phase)%p(of)%sizeState = sizeState sourceState(phase)%p(of)%sizeState = sizeState
sourceState(phase)%p(of)%sizeDotState = sizeDotState sourceState(phase)%p(of)%sizeDotState = sizeDotState
@ -1003,7 +999,7 @@ subroutine material_allocateSourceState(phase,of,NofMyPhase,&
if (numerics_integrator == 5_pInt) & if (numerics_integrator == 5_pInt) &
allocate(sourceState(phase)%p(of)%RKCK45dotState (6,sizeDotState,NofMyPhase), source=0.0_pReal) allocate(sourceState(phase)%p(of)%RKCK45dotState (6,sizeDotState,NofMyPhase), source=0.0_pReal)
allocate(plasticState(phase)%deltaState (sizeDeltaState,NofMyPhase), source=0.0_pReal) allocate(sourceState(phase)%p(of)%deltaState (sizeDeltaState,NofMyPhase), source=0.0_pReal)
end subroutine material_allocateSourceState end subroutine material_allocateSourceState

10
src/numerics.f90
View File

@ -23,12 +23,10 @@ module numerics
pert_method = 1_pInt, & !< method used in perturbation technique for tangent pert_method = 1_pInt, & !< method used in perturbation technique for tangent
randomSeed = 0_pInt, & !< fixed seeding for pseudo-random number generator, Default 0: use random seed randomSeed = 0_pInt, & !< fixed seeding for pseudo-random number generator, Default 0: use random seed
worldrank = 0_pInt, & !< MPI worldrank (/=0 for MPI simulations only) worldrank = 0_pInt, & !< MPI worldrank (/=0 for MPI simulations only)
worldsize = 0_pInt !< MPI worldsize (/=0 for MPI simulations only) worldsize = 0_pInt, & !< MPI worldsize (/=0 for MPI simulations only)
numerics_integrator = 1_pInt !< method used for state integration Default 1: fix-point iteration
integer(4), protected, public :: & integer(4), protected, public :: &
DAMASK_NumThreadsInt = 0 !< value stored in environment variable DAMASK_NUM_THREADS, set to zero if no OpenMP directive DAMASK_NumThreadsInt = 0 !< value stored in environment variable DAMASK_NUM_THREADS, set to zero if no OpenMP directive
!< ToDo: numerics_integrator is an array for historical reasons, only element 1 is used!
integer(pInt), dimension(2), protected, public :: &
numerics_integrator = 1_pInt !< method used for state integration (central & perturbed state), Default 1: fix-point iteration for both states
real(pReal), protected, public :: & real(pReal), protected, public :: &
relevantStrain = 1.0e-7_pReal, & !< strain increment considered significant (used by crystallite to determine whether strain inc is considered significant) relevantStrain = 1.0e-7_pReal, & !< strain increment considered significant (used by crystallite to determine whether strain inc is considered significant)
defgradTolerance = 1.0e-7_pReal, & !< deviation of deformation gradient that is still allowed (used by CPFEM to determine outdated ffn1) defgradTolerance = 1.0e-7_pReal, & !< deviation of deformation gradient that is still allowed (used by CPFEM to determine outdated ffn1)
@ -466,7 +464,7 @@ subroutine numerics_init
write(6,'(a24,1x,es8.1)') ' rTol_crystalliteState: ',rTol_crystalliteState write(6,'(a24,1x,es8.1)') ' rTol_crystalliteState: ',rTol_crystalliteState
write(6,'(a24,1x,es8.1)') ' rTol_crystalliteStress: ',rTol_crystalliteStress write(6,'(a24,1x,es8.1)') ' rTol_crystalliteStress: ',rTol_crystalliteStress
write(6,'(a24,1x,es8.1)') ' aTol_crystalliteStress: ',aTol_crystalliteStress write(6,'(a24,1x,es8.1)') ' aTol_crystalliteStress: ',aTol_crystalliteStress
write(6,'(a24,2(1x,i8))') ' integrator: ',numerics_integrator write(6,'(a24,1x,i8)') ' integrator: ',numerics_integrator
write(6,'(a24,1x,L8)') ' use ping pong scheme: ',usepingpong write(6,'(a24,1x,L8)') ' use ping pong scheme: ',usepingpong
write(6,'(a24,1x,es8.1,/)')' unitlength: ',numerics_unitlength write(6,'(a24,1x,es8.1,/)')' unitlength: ',numerics_unitlength
@ -589,7 +587,7 @@ subroutine numerics_init
if (rTol_crystalliteState <= 0.0_pReal) call IO_error(301_pInt,ext_msg='rTol_crystalliteState') if (rTol_crystalliteState <= 0.0_pReal) call IO_error(301_pInt,ext_msg='rTol_crystalliteState')
if (rTol_crystalliteStress <= 0.0_pReal) call IO_error(301_pInt,ext_msg='rTol_crystalliteStress') if (rTol_crystalliteStress <= 0.0_pReal) call IO_error(301_pInt,ext_msg='rTol_crystalliteStress')
if (aTol_crystalliteStress <= 0.0_pReal) call IO_error(301_pInt,ext_msg='aTol_crystalliteStress') if (aTol_crystalliteStress <= 0.0_pReal) call IO_error(301_pInt,ext_msg='aTol_crystalliteStress')
if (any(numerics_integrator <= 0_pInt) .or. any(numerics_integrator >= 6_pInt)) & if (numerics_integrator <= 0_pInt .or. numerics_integrator >= 6_pInt) &
call IO_error(301_pInt,ext_msg='integrator') call IO_error(301_pInt,ext_msg='integrator')
if (numerics_unitlength <= 0.0_pReal) call IO_error(301_pInt,ext_msg='unitlength') if (numerics_unitlength <= 0.0_pReal) call IO_error(301_pInt,ext_msg='unitlength')
if (absTol_RGC <= 0.0_pReal) call IO_error(301_pInt,ext_msg='absTol_RGC') if (absTol_RGC <= 0.0_pReal) call IO_error(301_pInt,ext_msg='absTol_RGC')

170
src/source_thermal_dissipation.f90
View File

@ -1,4 +1,5 @@
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
!> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH !> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH
!> @brief material subroutine for thermal source due to plastic dissipation !> @brief material subroutine for thermal source due to plastic dissipation
!> @details to be done !> @details to be done
@ -11,7 +12,6 @@ module source_thermal_dissipation
implicit none implicit none
private private
integer(pInt), dimension(:), allocatable, public, protected :: & integer(pInt), dimension(:), allocatable, public, protected :: &
source_thermal_dissipation_sizePostResults, & !< cumulative size of post results
source_thermal_dissipation_offset, & !< which source is my current thermal dissipation mechanism? source_thermal_dissipation_offset, & !< which source is my current thermal dissipation mechanism?
source_thermal_dissipation_instance !< instance of thermal dissipation source mechanism source_thermal_dissipation_instance !< instance of thermal dissipation source mechanism
@ -21,12 +21,18 @@ module source_thermal_dissipation
character(len=64), dimension(:,:), allocatable, target, public :: & character(len=64), dimension(:,:), allocatable, target, public :: &
source_thermal_dissipation_output !< name of each post result output source_thermal_dissipation_output !< name of each post result output
integer(pInt), dimension(:), allocatable, target, public :: &
source_thermal_dissipation_Noutput !< number of outputs per instance of this source
real(pReal), dimension(:), allocatable, private :: & real(pReal), dimension(:), allocatable, private :: &
source_thermal_dissipation_coldworkCoeff 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 :: & public :: &
source_thermal_dissipation_init, & source_thermal_dissipation_init, &
source_thermal_dissipation_getRateAndItsTangent source_thermal_dissipation_getRateAndItsTangent
@ -38,30 +44,13 @@ 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 source_thermal_dissipation_init(fileUnit) subroutine source_thermal_dissipation_init
#if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800
use, intrinsic :: iso_fortran_env, only: &
compiler_version, &
compiler_options
#endif
use debug, only: & use debug, only: &
debug_level,& debug_level,&
debug_constitutive,& debug_constitutive,&
debug_levelBasic 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
use material, only: & use material, only: &
material_allocateSourceState, &
phase_source, & phase_source, &
phase_Nsources, & phase_Nsources, &
phase_Noutput, & phase_Noutput, &
@ -70,144 +59,73 @@ subroutine source_thermal_dissipation_init(fileUnit)
material_phase, & material_phase, &
sourceState sourceState
use config, only: & use config, only: &
config_phase, &
material_Nphase, & material_Nphase, &
MATERIAL_partPhase MATERIAL_partPhase
use numerics,only: &
numerics_integrator
implicit none implicit none
integer(pInt), intent(in) :: fileUnit integer(pInt) :: Ninstance,instance,source,sourceOffset
integer(pInt) :: NofMyPhase,p
integer(pInt), allocatable, dimension(:) :: chunkPos
integer(pInt) :: maxNinstance,phase,instance,source,sourceOffset
integer(pInt) :: sizeState, sizeDotState, sizeDeltaState
integer(pInt) :: NofMyPhase
character(len=65536) :: &
tag = '', &
line = ''
write(6,'(/,a)') ' <<<+- source_'//SOURCE_thermal_dissipation_label//' init -+>>>' write(6,'(/,a)') ' <<<+- source_'//SOURCE_thermal_dissipation_label//' init -+>>>'
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) & 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_offset(material_Nphase), source=0_pInt)
allocate(source_thermal_dissipation_instance(material_Nphase), source=0_pInt) allocate(source_thermal_dissipation_instance(material_Nphase), source=0_pInt)
do phase = 1, material_Nphase do p = 1, material_Nphase
source_thermal_dissipation_instance(phase) = count(phase_source(:,1:phase) == SOURCE_thermal_dissipation_ID) source_thermal_dissipation_instance(p) = count(phase_source(:,1:p) == SOURCE_thermal_dissipation_ID)
do source = 1, phase_Nsources(phase) do source = 1, phase_Nsources(p)
if (phase_source(source,phase) == SOURCE_thermal_dissipation_ID) & if (phase_source(source,p) == SOURCE_thermal_dissipation_ID) &
source_thermal_dissipation_offset(phase) = source source_thermal_dissipation_offset(p) = source
enddo enddo
enddo enddo
allocate(source_thermal_dissipation_sizePostResults(maxNinstance), source=0_pInt) allocate(source_thermal_dissipation_sizePostResult(maxval(phase_Noutput),Ninstance),source=0_pInt)
allocate(source_thermal_dissipation_sizePostResult(maxval(phase_Noutput),maxNinstance),source=0_pInt) allocate(source_thermal_dissipation_output (maxval(phase_Noutput),Ninstance))
allocate(source_thermal_dissipation_output (maxval(phase_Noutput),maxNinstance))
source_thermal_dissipation_output = '' source_thermal_dissipation_output = ''
allocate(source_thermal_dissipation_Noutput(maxNinstance), source=0_pInt)
allocate(source_thermal_dissipation_coldworkCoeff(maxNinstance), source=0.0_pReal)
rewind(fileUnit) allocate(source_thermal_dissipation_coldworkCoeff(Ninstance), source=0.0_pReal)
phase = 0_pInt
do while (trim(line) /= IO_EOF .and. IO_lc(IO_getTag(line,'<','>')) /= MATERIAL_partPhase) ! wind forward to <phase> do p=1, size(config_phase)
line = IO_read(fileUnit) if (all(phase_source(:,p) /= SOURCE_THERMAL_DISSIPATION_ID)) cycle
instance = source_thermal_dissipation_instance(p)
source_thermal_dissipation_coldworkCoeff(instance) = config_phase(p)%getFloat('dissipation_coldworkcoeff')
NofMyPhase=count(material_phase==p)
sourceOffset = source_thermal_dissipation_offset(p)
call material_allocateSourceState(p,sourceOffset,NofMyPhase,0_pInt,0_pInt,0_pInt)
enddo 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_source(:,phase) == SOURCE_thermal_dissipation_ID)) then ! do not short-circuit here (.and. with next if statemen). It's not safe in Fortran
instance = source_thermal_dissipation_instance(phase) ! which instance of my source is present phase
chunkPos = IO_stringPos(line)
tag = IO_lc(IO_stringValue(line,chunkPos,1_pInt)) ! extract key
select case(tag)
case ('dissipation_coldworkcoeff')
source_thermal_dissipation_coldworkCoeff(instance) = IO_floatValue(line,chunkPos,2_pInt)
end select
endif; endif
enddo parsingFile
initializeInstances: do phase = 1_pInt, material_Nphase
if (any(phase_source(:,phase) == SOURCE_thermal_dissipation_ID)) then
NofMyPhase=count(material_phase==phase)
instance = source_thermal_dissipation_instance(phase)
sourceOffset = source_thermal_dissipation_offset(phase)
sizeDotState = 0_pInt
sizeDeltaState = 0_pInt
sizeState = 0_pInt
sourceState(phase)%p(sourceOffset)%sizeState = sizeState
sourceState(phase)%p(sourceOffset)%sizeDotState = sizeDotState
sourceState(phase)%p(sourceOffset)%sizeDeltaState = sizeDeltaState
sourceState(phase)%p(sourceOffset)%sizePostResults = source_thermal_dissipation_sizePostResults(instance)
allocate(sourceState(phase)%p(sourceOffset)%aTolState (sizeState), source=0.0_pReal)
allocate(sourceState(phase)%p(sourceOffset)%state0 (sizeState,NofMyPhase), source=0.0_pReal)
allocate(sourceState(phase)%p(sourceOffset)%partionedState0 (sizeState,NofMyPhase), source=0.0_pReal)
allocate(sourceState(phase)%p(sourceOffset)%subState0 (sizeState,NofMyPhase), source=0.0_pReal)
allocate(sourceState(phase)%p(sourceOffset)%state (sizeState,NofMyPhase), source=0.0_pReal)
allocate(sourceState(phase)%p(sourceOffset)%dotState (sizeDotState,NofMyPhase), source=0.0_pReal)
allocate(sourceState(phase)%p(sourceOffset)%deltaState (sizeDeltaState,NofMyPhase), source=0.0_pReal)
if (any(numerics_integrator == 1_pInt)) then
allocate(sourceState(phase)%p(sourceOffset)%previousDotState (sizeDotState,NofMyPhase), source=0.0_pReal)
allocate(sourceState(phase)%p(sourceOffset)%previousDotState2 (sizeDotState,NofMyPhase), source=0.0_pReal)
endif
if (any(numerics_integrator == 4_pInt)) &
allocate(sourceState(phase)%p(sourceOffset)%RK4dotState (sizeDotState,NofMyPhase), source=0.0_pReal)
if (any(numerics_integrator == 5_pInt)) &
allocate(sourceState(phase)%p(sourceOffset)%RKCK45dotState (6,sizeDotState,NofMyPhase),source=0.0_pReal)
endif
enddo initializeInstances
end subroutine source_thermal_dissipation_init end subroutine source_thermal_dissipation_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief returns local vacancy generation rate !> @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, Lp, phase)
use math, only: &
math_Mandel6to33
use material, only: &
phaseAt, phasememberAt
implicit none implicit none
integer(pInt), intent(in) :: & integer(pInt), intent(in) :: &
ipc, & !< grain number phase
ip, & !< integration point number real(pReal), intent(in), dimension(3,3) :: &
el !< element number Tstar
real(pReal), intent(in), dimension(6) :: &
Tstar_v !< 2nd Piola Kirchhoff stress tensor (Mandel)
real(pReal), intent(in), dimension(3,3) :: & real(pReal), intent(in), dimension(3,3) :: &
Lp Lp
real(pReal), intent(out) :: & real(pReal), intent(out) :: &
TDot, & TDot, &
dTDOT_dT dTDOT_dT
integer(pInt) :: & integer(pInt) :: &
instance, phase, constituent instance
phase = phaseAt(ipc,ip,el)
constituent = phasememberAt(ipc,ip,el)
instance = source_thermal_dissipation_instance(phase) instance = source_thermal_dissipation_instance(phase)
TDot = source_thermal_dissipation_coldworkCoeff(instance)* & TDot = source_thermal_dissipation_coldworkCoeff(instance)*sum(abs(Tstar*Lp))
sum(abs(math_Mandel6to33(Tstar_v)*Lp))
dTDOT_dT = 0.0_pReal dTDOT_dT = 0.0_pReal
end subroutine source_thermal_dissipation_getRateAndItsTangent end subroutine source_thermal_dissipation_getRateAndItsTangent

278
src/source_thermal_externalheat.f90
View File

@ -1,41 +1,44 @@
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
!> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH !> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH
!> @author Philip Eisenlohr, Michigan State University !> @author Philip Eisenlohr, Michigan State University
!> @brief material subroutine for variable heat source !> @brief material subroutine for variable heat source
!> @details to be done
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module source_thermal_externalheat module source_thermal_externalheat
use prec, only: & use prec, only: &
pReal, & pReal, &
pInt pInt
implicit none implicit none
private private
integer(pInt), dimension(:), allocatable, public, protected :: & integer(pInt), dimension(:), allocatable, public, protected :: &
source_thermal_externalheat_sizePostResults, & !< cumulative size of post results source_thermal_externalheat_offset, & !< which source is my current thermal dissipation mechanism?
source_thermal_externalheat_offset, & !< which source is my current thermal dissipation mechanism? source_thermal_externalheat_instance !< instance of thermal dissipation source mechanism
source_thermal_externalheat_instance !< instance of thermal dissipation source mechanism
integer(pInt), dimension(:,:), allocatable, target, public :: & integer(pInt), dimension(:,:), allocatable, target, public :: &
source_thermal_externalheat_sizePostResult !< size of each post result output source_thermal_externalheat_sizePostResult !< size of each post result output
character(len=64), dimension(:,:), allocatable, target, public :: & character(len=64), dimension(:,:), allocatable, target, public :: &
source_thermal_externalheat_output !< name of each post result output source_thermal_externalheat_output !< name of each post result output
integer(pInt), dimension(:), allocatable, target, public :: & integer(pInt), dimension(:), allocatable, target, public :: &
source_thermal_externalheat_Noutput !< number of outputs per instance of this source source_thermal_externalheat_Noutput !< number of outputs per instance of this source
integer(pInt), dimension(:), allocatable, private :: & type, private :: tParameters !< container type for internal constitutive parameters
source_thermal_externalheat_nIntervals real(pReal), dimension(:), allocatable :: &
time, &
heat_rate
integer(pInt) :: &
nIntervals
end type tParameters
real(pReal), dimension(:,:), allocatable, private :: & type(tParameters), dimension(:), allocatable, private :: param !< containers of constitutive parameters (len Ninstance)
source_thermal_externalheat_time, &
source_thermal_externalheat_rate
public :: &
source_thermal_externalheat_init, & public :: &
source_thermal_externalheat_dotState, & source_thermal_externalheat_init, &
source_thermal_externalheat_getRateAndItsTangent source_thermal_externalheat_dotState, &
source_thermal_externalheat_getRateAndItsTangent
contains contains
@ -44,170 +47,77 @@ 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 source_thermal_externalheat_init(fileUnit) subroutine source_thermal_externalheat_init
#if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800 use debug, only: &
use, intrinsic :: iso_fortran_env, only: & debug_level,&
compiler_version, & debug_constitutive,&
compiler_options debug_levelBasic
#endif use material, only: &
use debug, only: & material_allocateSourceState, &
debug_level,& material_phase, &
debug_constitutive,& phase_source, &
debug_levelBasic phase_Nsources, &
use IO, only: & phase_Noutput, &
IO_read, & SOURCE_thermal_externalheat_label, &
IO_lc, & SOURCE_thermal_externalheat_ID
IO_getTag, & use config, only: &
IO_isBlank, & config_phase, &
IO_stringPos, & material_Nphase, &
IO_stringValue, & MATERIAL_partPhase
IO_floatValue, &
IO_intValue, &
IO_warning, &
IO_error, &
IO_timeStamp, &
IO_EOF
use material, only: &
phase_source, &
phase_Nsources, &
phase_Noutput, &
SOURCE_thermal_externalheat_label, &
SOURCE_thermal_externalheat_ID, &
material_phase, &
sourceState
use config, only: &
material_Nphase, &
MATERIAL_partPhase
use numerics,only: &
numerics_integrator
implicit none implicit none
integer(pInt), intent(in) :: fileUnit
integer(pInt), allocatable, dimension(:) :: chunkPos real(pReal), allocatable, dimension(:) :: tempVar
integer(pInt) :: maxNinstance,phase,instance,source,sourceOffset integer(pInt) :: maxNinstance,instance,source,sourceOffset
integer(pInt) :: sizeState, sizeDotState, sizeDeltaState integer(pInt) :: NofMyPhase,p
integer(pInt) :: NofMyPhase,interval
character(len=65536) :: &
tag = '', &
line = ''
real(pReal), allocatable, dimension(:,:) :: temp_time, temp_rate
write(6,'(/,a)') ' <<<+- source_'//SOURCE_thermal_externalheat_label//' init -+>>>' write(6,'(/,a)') ' <<<+- source_'//SOURCE_thermal_externalheat_label//' init -+>>>'
write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
#include "compilation_info.f90"
maxNinstance = int(count(phase_source == SOURCE_thermal_externalheat_ID),pInt)
if (maxNinstance == 0_pInt) return
if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt) &
write(6,'(a16,1x,i5,/)') '# instances:',maxNinstance
allocate(source_thermal_externalheat_offset(material_Nphase), source=0_pInt) maxNinstance = int(count(phase_source == SOURCE_thermal_externalheat_ID),pInt)
allocate(source_thermal_externalheat_instance(material_Nphase), source=0_pInt) if (maxNinstance == 0_pInt) return
do phase = 1, material_Nphase if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt) &
source_thermal_externalheat_instance(phase) = count(phase_source(:,1:phase) == SOURCE_thermal_externalheat_ID) write(6,'(a16,1x,i5,/)') '# instances:',maxNinstance
do source = 1, phase_Nsources(phase)
if (phase_source(source,phase) == SOURCE_thermal_externalheat_ID) &
source_thermal_externalheat_offset(phase) = source
enddo
enddo
allocate(source_thermal_externalheat_sizePostResults(maxNinstance), source=0_pInt) allocate(source_thermal_externalheat_offset(material_Nphase), source=0_pInt)
allocate(source_thermal_externalheat_sizePostResult(maxval(phase_Noutput),maxNinstance),source=0_pInt) allocate(source_thermal_externalheat_instance(material_Nphase), source=0_pInt)
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(temp_time(maxNinstance,1000), source=0.0_pReal) do p = 1, material_Nphase
allocate(temp_rate(maxNinstance,1000), source=0.0_pReal) source_thermal_externalheat_instance(p) = count(phase_source(:,1:p) == SOURCE_thermal_externalheat_ID)
do source = 1, phase_Nsources(p)
if (phase_source(source,p) == SOURCE_thermal_externalheat_ID) &
source_thermal_externalheat_offset(p) = source
enddo
enddo
rewind(fileUnit) allocate(source_thermal_externalheat_sizePostResult(maxval(phase_Noutput),maxNinstance),source=0_pInt)
phase = 0_pInt allocate(source_thermal_externalheat_output (maxval(phase_Noutput),maxNinstance))
do while (trim(line) /= IO_EOF .and. IO_lc(IO_getTag(line,'<','>')) /= MATERIAL_partPhase) ! wind forward to <phase> source_thermal_externalheat_output = ''
line = IO_read(fileUnit) allocate(source_thermal_externalheat_Noutput(maxNinstance), source=0_pInt)
enddo
parsingFile: do while (trim(line) /= IO_EOF) ! read through sections of phase part allocate(param(maxNinstance))
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_source(:,phase) == SOURCE_thermal_externalheat_ID)) then ! do not short-circuit here (.and. with next if statemen). It's not safe in Fortran do p=1, size(config_phase)
if (all(phase_source(:,p) /= SOURCE_thermal_externalheat_ID)) cycle
instance = source_thermal_externalheat_instance(p)
sourceOffset = source_thermal_externalheat_offset(p)
NofMyPhase=count(material_phase==p)
instance = source_thermal_externalheat_instance(phase) ! which instance of my source is present phase tempVar = config_phase(p)%getFloats('externalheat_time')
chunkPos = IO_stringPos(line) param(instance)%nIntervals = size(tempVar) - 1_pInt
tag = IO_lc(IO_stringValue(line,chunkPos,1_pInt)) ! extract key
select case(tag)
case ('externalheat_time','externalheat_rate')
if (chunkPos(1) <= 2_pInt) &
call IO_error(150_pInt,ext_msg=trim(tag)//' ('//SOURCE_thermal_externalheat_label//')')
if ( source_thermal_externalheat_nIntervals(instance) > 0_pInt &
.and. source_thermal_externalheat_nIntervals(instance) /= chunkPos(1) - 2_pInt) &
call IO_error(150_pInt,ext_msg=trim(tag)//' ('//SOURCE_thermal_externalheat_label//')')
source_thermal_externalheat_nIntervals(instance) = chunkPos(1) - 2_pInt param(instance)%time= tempVar
do interval = 1, source_thermal_externalheat_nIntervals(instance) + 1_pInt
select case(tag)
case ('externalheat_time')
temp_time(instance, interval) = IO_floatValue(line,chunkPos,1_pInt + interval)
case ('externalheat_rate')
temp_rate(instance, interval) = IO_floatValue(line,chunkPos,1_pInt + interval)
end select
enddo
end select
endif; endif
enddo parsingFile
allocate(source_thermal_externalheat_time(maxNinstance,maxval(source_thermal_externalheat_nIntervals)+1_pInt), source=0.0_pReal) tempVar = config_phase(p)%getFloats('externalheat_rate',requiredSize = size(tempVar))
allocate(source_thermal_externalheat_rate(maxNinstance,maxval(source_thermal_externalheat_nIntervals)+1_pInt), source=0.0_pReal) param(instance)%heat_rate = tempVar
initializeInstances: do phase = 1_pInt, material_Nphase call material_allocateSourceState(p,sourceOffset,NofMyPhase,1_pInt,1_pInt,0_pInt)
if (any(phase_source(:,phase) == SOURCE_thermal_externalheat_ID)) then
NofMyPhase = count(material_phase==phase)
instance = source_thermal_externalheat_instance(phase)
sourceOffset = source_thermal_externalheat_offset(phase)
source_thermal_externalheat_time(instance,1:source_thermal_externalheat_nIntervals(instance)+1_pInt) = &
temp_time(instance,1:source_thermal_externalheat_nIntervals(instance)+1_pInt)
source_thermal_externalheat_rate(instance,1:source_thermal_externalheat_nIntervals(instance)+1_pInt) = &
temp_rate(instance,1:source_thermal_externalheat_nIntervals(instance)+1_pInt)
sizeDotState = 1_pInt enddo
sizeDeltaState = 0_pInt
sizeState = 1_pInt
sourceState(phase)%p(sourceOffset)%sizeState = sizeState
sourceState(phase)%p(sourceOffset)%sizeDotState = sizeDotState
sourceState(phase)%p(sourceOffset)%sizeDeltaState = sizeDeltaState
sourceState(phase)%p(sourceOffset)%sizePostResults = source_thermal_externalheat_sizePostResults(instance)
allocate(sourceState(phase)%p(sourceOffset)%aTolState (sizeState), source=0.00001_pReal)
allocate(sourceState(phase)%p(sourceOffset)%state0 (sizeState,NofMyPhase), source=0.0_pReal)
allocate(sourceState(phase)%p(sourceOffset)%partionedState0 (sizeState,NofMyPhase), source=0.0_pReal)
allocate(sourceState(phase)%p(sourceOffset)%subState0 (sizeState,NofMyPhase), source=0.0_pReal)
allocate(sourceState(phase)%p(sourceOffset)%state (sizeState,NofMyPhase), source=0.0_pReal)
allocate(sourceState(phase)%p(sourceOffset)%dotState (sizeDotState,NofMyPhase), source=0.0_pReal)
allocate(sourceState(phase)%p(sourceOffset)%deltaState (sizeDeltaState,NofMyPhase), source=0.0_pReal)
if (any(numerics_integrator == 1_pInt)) then
allocate(sourceState(phase)%p(sourceOffset)%previousDotState (sizeDotState,NofMyPhase), source=0.0_pReal)
allocate(sourceState(phase)%p(sourceOffset)%previousDotState2 (sizeDotState,NofMyPhase), source=0.0_pReal)
endif
if (any(numerics_integrator == 4_pInt)) &
allocate(sourceState(phase)%p(sourceOffset)%RK4dotState (sizeDotState,NofMyPhase), source=0.0_pReal)
if (any(numerics_integrator == 5_pInt)) &
allocate(sourceState(phase)%p(sourceOffset)%RKCK45dotState (6,sizeDotState,NofMyPhase),source=0.0_pReal)
endif
enddo initializeInstances
end subroutine source_thermal_externalheat_init end subroutine source_thermal_externalheat_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief rate of change of state !> @brief rate of change of state
!> @details state only contains current time to linearly interpolate given heat powers !> @details state only contains current time to linearly interpolate given heat powers
@ -238,39 +148,35 @@ end subroutine source_thermal_externalheat_dotState
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief returns local heat generation rate !> @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: & use material, only: &
phaseAt, phasememberAt, &
sourceState sourceState
implicit none implicit none
integer(pInt), intent(in) :: & integer(pInt), intent(in) :: &
ipc, & !< grain number phase, &
ip, & !< integration point number constituent
el !< element number
real(pReal), intent(out) :: & real(pReal), intent(out) :: &
TDot, & TDot, &
dTDot_dT dTDot_dT
integer(pInt) :: & integer(pInt) :: &
instance, phase, constituent, sourceOffset, interval instance, sourceOffset, interval
real(pReal) :: & real(pReal) :: &
frac_time frac_time
phase = phaseAt(ipc,ip,el)
constituent = phasememberAt(ipc,ip,el)
instance = source_thermal_externalheat_instance(phase) instance = source_thermal_externalheat_instance(phase)
sourceOffset = source_thermal_externalheat_offset(phase) sourceOffset = source_thermal_externalheat_offset(phase)
do interval = 1, source_thermal_externalheat_nIntervals(instance) ! scan through all rate segments do interval = 1, param(instance)%nIntervals ! scan through all rate segments
frac_time = (sourceState(phase)%p(sourceOffset)%state(1,constituent) - & frac_time = (sourceState(phase)%p(sourceOffset)%state(1,constituent) - &
source_thermal_externalheat_time(instance,interval)) / & param(instance)%time(interval)) / &
(source_thermal_externalheat_time(instance,interval+1) - & (param(instance)%time(interval+1) - &
source_thermal_externalheat_time(instance,interval)) ! fractional time within segment param(instance)%time(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 == source_thermal_externalheat_nIntervals(instance)) & .or. (frac_time >= 1.0_pReal .and. interval == param(instance)%nIntervals) &
.or. (frac_time >= 0.0_pReal .and. frac_time < 1.0_pReal) ) & .or. (frac_time >= 0.0_pReal .and. frac_time < 1.0_pReal) ) &
TDot = source_thermal_externalheat_rate(instance,interval ) * (1.0_pReal - frac_time) + & TDot = param(instance)%heat_rate(interval ) * (1.0_pReal - frac_time) + &
source_thermal_externalheat_rate(instance,interval+1) * frac_time ! interpolate heat rate between segment boundaries... param(instance)%heat_rate(interval+1) * frac_time ! interpolate heat rate between segment boundaries...
! ...or extrapolate if outside of bounds ! ...or extrapolate if outside of bounds
enddo enddo
dTDot_dT = 0.0 dTDot_dT = 0.0

189
src/thermal_adiabatic.f90
View File

@ -10,12 +10,9 @@ module thermal_adiabatic
implicit none implicit none
private private
integer(pInt), dimension(:), allocatable, public, protected :: &
thermal_adiabatic_sizePostResults !< cumulative size of post results
integer(pInt), dimension(:,:), allocatable, target, public :: & integer(pInt), dimension(:,:), allocatable, target, public :: &
thermal_adiabatic_sizePostResult !< size of each post result output thermal_adiabatic_sizePostResult !< size of each post result output
character(len=64), dimension(:,:), allocatable, target, public :: & character(len=64), dimension(:,:), allocatable, target, public :: &
thermal_adiabatic_output !< name of each post result output thermal_adiabatic_output !< name of each post result output
@ -45,27 +42,7 @@ 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_adiabatic_init(fileUnit) subroutine thermal_adiabatic_init
#if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800
use, intrinsic :: iso_fortran_env, only: &
compiler_version, &
compiler_options
#endif
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
use config, only: &
material_partHomogenization
use material, only: & use material, only: &
thermal_type, & thermal_type, &
thermal_typeInstance, & thermal_typeInstance, &
@ -79,106 +56,61 @@ subroutine thermal_adiabatic_init(fileUnit)
thermal_initialT, & thermal_initialT, &
temperature, & temperature, &
temperatureRate temperatureRate
use config, only: &
material_partHomogenization, &
config_homogenization
implicit none implicit none
integer(pInt), intent(in) :: fileUnit integer(pInt) :: maxNinstance,section,instance,i
integer(pInt), allocatable, dimension(:) :: chunkPos
integer(pInt) :: maxNinstance,mySize=0_pInt,section,instance,o
integer(pInt) :: sizeState integer(pInt) :: sizeState
integer(pInt) :: NofMyHomog integer(pInt) :: NofMyHomog
character(len=65536) :: & character(len=65536), dimension(0), parameter :: emptyStringArray = [character(len=65536)::]
tag = '', & character(len=65536), dimension(:), allocatable :: outputs
line = ''
write(6,'(/,a)') ' <<<+- thermal_'//THERMAL_ADIABATIC_label//' init -+>>>' write(6,'(/,a)') ' <<<+- thermal_'//THERMAL_ADIABATIC_label//' init -+>>>'
write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
#include "compilation_info.f90"
maxNinstance = int(count(thermal_type == THERMAL_adiabatic_ID),pInt) maxNinstance = int(count(thermal_type == THERMAL_adiabatic_ID),pInt)
if (maxNinstance == 0_pInt) return if (maxNinstance == 0_pInt) return
allocate(thermal_adiabatic_sizePostResults(maxNinstance), source=0_pInt)
allocate(thermal_adiabatic_sizePostResult (maxval(homogenization_Noutput),maxNinstance),source=0_pInt) allocate(thermal_adiabatic_sizePostResult (maxval(homogenization_Noutput),maxNinstance),source=0_pInt)
allocate(thermal_adiabatic_output (maxval(homogenization_Noutput),maxNinstance)) allocate(thermal_adiabatic_output (maxval(homogenization_Noutput),maxNinstance))
thermal_adiabatic_output = '' thermal_adiabatic_output = ''
allocate(thermal_adiabatic_outputID (maxval(homogenization_Noutput),maxNinstance),source=undefined_ID) allocate(thermal_adiabatic_outputID (maxval(homogenization_Noutput),maxNinstance),source=undefined_ID)
allocate(thermal_adiabatic_Noutput (maxNinstance), source=0_pInt) allocate(thermal_adiabatic_Noutput (maxNinstance), source=0_pInt)
rewind(fileUnit)
section = 0_pInt
do while (trim(line) /= IO_EOF .and. IO_lc(IO_getTag(line,'<','>')) /= material_partHomogenization)! wind forward to <homogenization>
line = IO_read(fileUnit)
enddo
parsingFile: do while (trim(line) /= IO_EOF) ! read through sections of homog 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 homog section
section = section + 1_pInt ! advance homog section counter
cycle ! skip to next line
endif
if (section > 0_pInt ) then; if (thermal_type(section) == THERMAL_adiabatic_ID) then ! do not short-circuit here (.and. with next if statemen). It's not safe in Fortran
instance = thermal_typeInstance(section) ! which instance of my thermal is present homog
chunkPos = IO_stringPos(line)
tag = IO_lc(IO_stringValue(line,chunkPos,1_pInt)) ! extract key
select case(tag)
case ('(output)')
select case(IO_lc(IO_stringValue(line,chunkPos,2_pInt)))
case ('temperature')
thermal_adiabatic_Noutput(instance) = thermal_adiabatic_Noutput(instance) + 1_pInt
thermal_adiabatic_outputID(thermal_adiabatic_Noutput(instance),instance) = temperature_ID
thermal_adiabatic_output(thermal_adiabatic_Noutput(instance),instance) = &
IO_lc(IO_stringValue(line,chunkPos,2_pInt))
end select
end select
endif; endif
enddo parsingFile
initializeInstances: do section = 1_pInt, size(thermal_type) initializeInstances: do section = 1_pInt, size(thermal_type)
if (thermal_type(section) == THERMAL_adiabatic_ID) then if (thermal_type(section) /= THERMAL_adiabatic_ID) cycle
NofMyHomog=count(material_homog==section) NofMyHomog=count(material_homog==section)
instance = thermal_typeInstance(section) instance = thermal_typeInstance(section)
outputs = config_homogenization(section)%getStrings('(output)',defaultVal=emptyStringArray)
!-------------------------------------------------------------------------------------------------- do i=1_pInt, size(outputs)
! Determine size of postResults array select case(outputs(i))
outputsLoop: do o = 1_pInt,thermal_adiabatic_Noutput(instance) case('temperature')
select case(thermal_adiabatic_outputID(o,instance)) thermal_adiabatic_Noutput(instance) = thermal_adiabatic_Noutput(instance) + 1_pInt
case(temperature_ID) thermal_adiabatic_outputID(thermal_adiabatic_Noutput(instance),instance) = temperature_ID
mySize = 1_pInt thermal_adiabatic_output(thermal_adiabatic_Noutput(instance),instance) = outputs(i)
end select thermal_adiabatic_sizePostResult(thermal_adiabatic_Noutput(instance),instance) = 1_pInt
end select
if (mySize > 0_pInt) then ! any meaningful output found enddo
thermal_adiabatic_sizePostResult(o,instance) = mySize
thermal_adiabatic_sizePostResults(instance) = thermal_adiabatic_sizePostResults(instance) + mySize
endif
enddo outputsLoop
! allocate state arrays ! allocate state arrays
sizeState = 1_pInt sizeState = 1_pInt
thermalState(section)%sizeState = sizeState thermalState(section)%sizeState = sizeState
thermalState(section)%sizePostResults = thermal_adiabatic_sizePostResults(instance) thermalState(section)%sizePostResults = sum(thermal_adiabatic_sizePostResult(:,instance))
allocate(thermalState(section)%state0 (sizeState,NofMyHomog), source=thermal_initialT(section)) allocate(thermalState(section)%state0 (sizeState,NofMyHomog), source=thermal_initialT(section))
allocate(thermalState(section)%subState0(sizeState,NofMyHomog), source=thermal_initialT(section)) allocate(thermalState(section)%subState0(sizeState,NofMyHomog), source=thermal_initialT(section))
allocate(thermalState(section)%state (sizeState,NofMyHomog), source=thermal_initialT(section)) allocate(thermalState(section)%state (sizeState,NofMyHomog), source=thermal_initialT(section))
nullify(thermalMapping(section)%p) nullify(thermalMapping(section)%p)
thermalMapping(section)%p => mappingHomogenization(1,:,:) thermalMapping(section)%p => mappingHomogenization(1,:,:)
deallocate(temperature(section)%p) deallocate(temperature(section)%p)
temperature(section)%p => thermalState(section)%state(1,:) temperature(section)%p => thermalState(section)%state(1,:)
deallocate(temperatureRate(section)%p) deallocate(temperatureRate(section)%p)
allocate (temperatureRate(section)%p(NofMyHomog), source=0.0_pReal) allocate (temperatureRate(section)%p(NofMyHomog), source=0.0_pReal)
endif
enddo initializeInstances enddo initializeInstances
end subroutine thermal_adiabatic_init end subroutine thermal_adiabatic_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -233,11 +165,12 @@ end function thermal_adiabatic_updateState
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine thermal_adiabatic_getSourceAndItsTangent(Tdot, dTdot_dT, T, ip, el) subroutine thermal_adiabatic_getSourceAndItsTangent(Tdot, dTdot_dT, T, ip, el)
use math, only: & use math, only: &
math_Mandel6to33 math_6toSym33
use material, only: & use material, only: &
homogenization_Ngrains, & homogenization_Ngrains, &
mappingHomogenization, & mappingHomogenization, &
phaseAt, & phaseAt, &
phasememberAt, &
thermal_typeInstance, & thermal_typeInstance, &
phase_Nsources, & phase_Nsources, &
phase_source, & phase_source, &
@ -264,30 +197,30 @@ subroutine thermal_adiabatic_getSourceAndItsTangent(Tdot, dTdot_dT, T, ip, el)
integer(pInt) :: & integer(pInt) :: &
phase, & phase, &
homog, & homog, &
offset, &
instance, & instance, &
grain, & grain, &
source source, &
constituent
homog = mappingHomogenization(2,ip,el) homog = mappingHomogenization(2,ip,el)
offset = mappingHomogenization(1,ip,el)
instance = thermal_typeInstance(homog) instance = thermal_typeInstance(homog)
Tdot = 0.0_pReal Tdot = 0.0_pReal
dTdot_dT = 0.0_pReal dTdot_dT = 0.0_pReal
do grain = 1, homogenization_Ngrains(homog) do grain = 1, homogenization_Ngrains(homog)
phase = phaseAt(grain,ip,el) phase = phaseAt(grain,ip,el)
constituent = phasememberAt(grain,ip,el)
do source = 1, phase_Nsources(phase) do source = 1, phase_Nsources(phase)
select case(phase_source(source,phase)) select case(phase_source(source,phase))
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, &
crystallite_Tstar_v(1:6,grain,ip,el), & math_6toSym33(crystallite_Tstar_v(1:6,grain,ip,el)), &
crystallite_Lp(1:3,1:3,grain,ip,el), & crystallite_Lp(1:3,1:3,grain,ip,el), &
grain, ip, el) phase)
case (SOURCE_thermal_externalheat_ID) case (SOURCE_thermal_externalheat_ID)
call source_thermal_externalheat_getRateAndItsTangent(my_Tdot, my_dTdot_dT, & call source_thermal_externalheat_getRateAndItsTangent(my_Tdot, my_dTdot_dT, &
grain, ip, el) phase, constituent)
case default case default
my_Tdot = 0.0_pReal my_Tdot = 0.0_pReal
@ -311,12 +244,9 @@ function thermal_adiabatic_getSpecificHeat(ip,el)
lattice_specificHeat lattice_specificHeat
use material, only: & use material, only: &
homogenization_Ngrains, & homogenization_Ngrains, &
mappingHomogenization, &
material_phase material_phase
use mesh, only: & use mesh, only: &
mesh_element mesh_element
use crystallite, only: &
crystallite_push33ToRef
implicit none implicit none
integer(pInt), intent(in) :: & integer(pInt), intent(in) :: &
@ -325,11 +255,10 @@ function thermal_adiabatic_getSpecificHeat(ip,el)
real(pReal) :: & real(pReal) :: &
thermal_adiabatic_getSpecificHeat thermal_adiabatic_getSpecificHeat
integer(pInt) :: & integer(pInt) :: &
homog, grain grain
thermal_adiabatic_getSpecificHeat = 0.0_pReal thermal_adiabatic_getSpecificHeat = 0.0_pReal
homog = mappingHomogenization(2,ip,el)
do grain = 1, homogenization_Ngrains(mesh_element(3,el)) do grain = 1, homogenization_Ngrains(mesh_element(3,el))
thermal_adiabatic_getSpecificHeat = thermal_adiabatic_getSpecificHeat + & thermal_adiabatic_getSpecificHeat = thermal_adiabatic_getSpecificHeat + &
@ -341,6 +270,7 @@ function thermal_adiabatic_getSpecificHeat(ip,el)
end function thermal_adiabatic_getSpecificHeat end function thermal_adiabatic_getSpecificHeat
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief returns homogenized mass density !> @brief returns homogenized mass density
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -353,8 +283,6 @@ function thermal_adiabatic_getMassDensity(ip,el)
material_phase material_phase
use mesh, only: & use mesh, only: &
mesh_element mesh_element
use crystallite, only: &
crystallite_push33ToRef
implicit none implicit none
integer(pInt), intent(in) :: & integer(pInt), intent(in) :: &
@ -363,11 +291,10 @@ function thermal_adiabatic_getMassDensity(ip,el)
real(pReal) :: & real(pReal) :: &
thermal_adiabatic_getMassDensity thermal_adiabatic_getMassDensity
integer(pInt) :: & integer(pInt) :: &
homog, grain grain
thermal_adiabatic_getMassDensity = 0.0_pReal thermal_adiabatic_getMassDensity = 0.0_pReal
homog = mappingHomogenization(2,ip,el)
do grain = 1, homogenization_Ngrains(mesh_element(3,el)) do grain = 1, homogenization_Ngrains(mesh_element(3,el))
thermal_adiabatic_getMassDensity = thermal_adiabatic_getMassDensity + & thermal_adiabatic_getMassDensity = thermal_adiabatic_getMassDensity + &
@ -379,41 +306,37 @@ function thermal_adiabatic_getMassDensity(ip,el)
end function thermal_adiabatic_getMassDensity end function thermal_adiabatic_getMassDensity
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief return array of thermal results !> @brief return array of thermal results
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function thermal_adiabatic_postResults(ip,el) function thermal_adiabatic_postResults(homog,instance,of) result(postResults)
use material, only: & use material, only: &
mappingHomogenization, &
thermal_typeInstance, &
thermalMapping, &
temperature temperature
implicit none implicit none
integer(pInt), intent(in) :: & integer(pInt), intent(in) :: &
ip, & !< integration point homog, &
el !< element instance, &
real(pReal), dimension(thermal_adiabatic_sizePostResults(thermal_typeInstance(mappingHomogenization(2,ip,el)))) :: & of
thermal_adiabatic_postResults
real(pReal), dimension(sum(thermal_adiabatic_sizePostResult(:,instance))) :: &
postResults
integer(pInt) :: & integer(pInt) :: &
instance, homog, offset, o, c o, c
homog = mappingHomogenization(2,ip,el)
offset = thermalMapping(homog)%p(ip,el)
instance = thermal_typeInstance(homog)
c = 0_pInt c = 0_pInt
thermal_adiabatic_postResults = 0.0_pReal
do o = 1_pInt,thermal_adiabatic_Noutput(instance) do o = 1_pInt,thermal_adiabatic_Noutput(instance)
select case(thermal_adiabatic_outputID(o,instance)) select case(thermal_adiabatic_outputID(o,instance))
case (temperature_ID) case (temperature_ID)
thermal_adiabatic_postResults(c+1_pInt) = temperature(homog)%p(offset) postResults(c+1_pInt) = temperature(homog)%p(of)
c = c + 1 c = c + 1
end select end select
enddo enddo
end function thermal_adiabatic_postResults end function thermal_adiabatic_postResults
end module thermal_adiabatic end module thermal_adiabatic

194
src/thermal_conduction.f90
View File

@ -10,12 +10,9 @@ module thermal_conduction
implicit none implicit none
private private
integer(pInt), dimension(:), allocatable, public, protected :: &
thermal_conduction_sizePostResults !< cumulative size of post results
integer(pInt), dimension(:,:), allocatable, target, public :: & integer(pInt), dimension(:,:), allocatable, target, public :: &
thermal_conduction_sizePostResult !< size of each post result output thermal_conduction_sizePostResult !< size of each post result output
character(len=64), dimension(:,:), allocatable, target, public :: & character(len=64), dimension(:,:), allocatable, target, public :: &
thermal_conduction_output !< name of each post result output thermal_conduction_output !< name of each post result output
@ -46,25 +43,7 @@ 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(fileUnit) subroutine thermal_conduction_init
#if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800
use, intrinsic :: iso_fortran_env, only: &
compiler_version, &
compiler_options
#endif
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
use material, only: & use material, only: &
thermal_type, & thermal_type, &
thermal_typeInstance, & thermal_typeInstance, &
@ -79,107 +58,61 @@ subroutine thermal_conduction_init(fileUnit)
temperature, & temperature, &
temperatureRate temperatureRate
use config, only: & use config, only: &
material_partHomogenization material_partHomogenization, &
config_homogenization
implicit none implicit none
integer(pInt), intent(in) :: fileUnit integer(pInt) :: maxNinstance,section,instance,i
integer(pInt), allocatable, dimension(:) :: chunkPos
integer(pInt) :: maxNinstance,mySize=0_pInt,section,instance,o
integer(pInt) :: sizeState integer(pInt) :: sizeState
integer(pInt) :: NofMyHomog integer(pInt) :: NofMyHomog
character(len=65536) :: & character(len=65536), dimension(0), parameter :: emptyStringArray = [character(len=65536)::]
tag = '', & character(len=65536), dimension(:), allocatable :: outputs
line = ''
write(6,'(/,a)') ' <<<+- thermal_'//THERMAL_CONDUCTION_label//' init -+>>>' write(6,'(/,a)') ' <<<+- thermal_'//THERMAL_CONDUCTION_label//' init -+>>>'
write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
#include "compilation_info.f90"
maxNinstance = int(count(thermal_type == THERMAL_conduction_ID),pInt) maxNinstance = int(count(thermal_type == THERMAL_conduction_ID),pInt)
if (maxNinstance == 0_pInt) return if (maxNinstance == 0_pInt) return
allocate(thermal_conduction_sizePostResults(maxNinstance), source=0_pInt)
allocate(thermal_conduction_sizePostResult (maxval(homogenization_Noutput),maxNinstance),source=0_pInt) allocate(thermal_conduction_sizePostResult (maxval(homogenization_Noutput),maxNinstance),source=0_pInt)
allocate(thermal_conduction_output (maxval(homogenization_Noutput),maxNinstance)) allocate(thermal_conduction_output (maxval(homogenization_Noutput),maxNinstance))
thermal_conduction_output = '' thermal_conduction_output = ''
allocate(thermal_conduction_outputID (maxval(homogenization_Noutput),maxNinstance),source=undefined_ID) allocate(thermal_conduction_outputID (maxval(homogenization_Noutput),maxNinstance),source=undefined_ID)
allocate(thermal_conduction_Noutput (maxNinstance), source=0_pInt) allocate(thermal_conduction_Noutput (maxNinstance), source=0_pInt)
rewind(fileUnit)
section = 0_pInt
do while (trim(line) /= IO_EOF .and. IO_lc(IO_getTag(line,'<','>')) /= material_partHomogenization)! wind forward to <homogenization>
line = IO_read(fileUnit)
enddo
parsingFile: do while (trim(line) /= IO_EOF) ! read through sections of homog 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 homog section
section = section + 1_pInt ! advance homog section counter
cycle ! skip to next line
endif
if (section > 0_pInt ) then; if (thermal_type(section) == THERMAL_conduction_ID) then ! do not short-circuit here (.and. with next if statemen). It's not safe in Fortran
instance = thermal_typeInstance(section) ! which instance of my thermal is present homog
chunkPos = IO_stringPos(line)
tag = IO_lc(IO_stringValue(line,chunkPos,1_pInt)) ! extract key
select case(tag)
case ('(output)')
select case(IO_lc(IO_stringValue(line,chunkPos,2_pInt)))
case ('temperature')
thermal_conduction_Noutput(instance) = thermal_conduction_Noutput(instance) + 1_pInt
thermal_conduction_outputID(thermal_conduction_Noutput(instance),instance) = temperature_ID
thermal_conduction_output(thermal_conduction_Noutput(instance),instance) = &
IO_lc(IO_stringValue(line,chunkPos,2_pInt))
end select
end select
endif; endif
enddo parsingFile
initializeInstances: do section = 1_pInt, size(thermal_type) initializeInstances: do section = 1_pInt, size(thermal_type)
if (thermal_type(section) == THERMAL_conduction_ID) then if (thermal_type(section) /= THERMAL_conduction_ID) cycle
NofMyHomog=count(material_homog==section) NofMyHomog=count(material_homog==section)
instance = thermal_typeInstance(section) instance = thermal_typeInstance(section)
outputs = config_homogenization(section)%getStrings('(output)',defaultVal=emptyStringArray)
do i=1_pInt, size(outputs)
select case(outputs(i))
case('temperature')
thermal_conduction_Noutput(instance) = thermal_conduction_Noutput(instance) + 1_pInt
thermal_conduction_outputID(thermal_conduction_Noutput(instance),instance) = temperature_ID
thermal_conduction_output(thermal_conduction_Noutput(instance),instance) = outputs(i)
thermal_conduction_sizePostResult(thermal_conduction_Noutput(instance),instance) = 1_pInt
end select
enddo
!--------------------------------------------------------------------------------------------------
! Determine size of postResults array
outputsLoop: do o = 1_pInt,thermal_conduction_Noutput(instance)
select case(thermal_conduction_outputID(o,instance))
case(temperature_ID)
mySize = 1_pInt
end select
if (mySize > 0_pInt) then ! any meaningful output found
thermal_conduction_sizePostResult(o,instance) = mySize
thermal_conduction_sizePostResults(instance) = thermal_conduction_sizePostResults(instance) + mySize
endif
enddo outputsLoop
! allocate state arrays ! allocate state arrays
sizeState = 0_pInt sizeState = 0_pInt
thermalState(section)%sizeState = sizeState thermalState(section)%sizeState = sizeState
thermalState(section)%sizePostResults = thermal_conduction_sizePostResults(instance) thermalState(section)%sizePostResults = sum(thermal_conduction_sizePostResult(:,instance))
allocate(thermalState(section)%state0 (sizeState,NofMyHomog)) allocate(thermalState(section)%state0 (sizeState,NofMyHomog))
allocate(thermalState(section)%subState0(sizeState,NofMyHomog)) allocate(thermalState(section)%subState0(sizeState,NofMyHomog))
allocate(thermalState(section)%state (sizeState,NofMyHomog)) allocate(thermalState(section)%state (sizeState,NofMyHomog))
nullify(thermalMapping(section)%p) nullify(thermalMapping(section)%p)
thermalMapping(section)%p => mappingHomogenization(1,:,:) thermalMapping(section)%p => mappingHomogenization(1,:,:)
deallocate(temperature (section)%p) deallocate(temperature (section)%p)
allocate (temperature (section)%p(NofMyHomog), source=thermal_initialT(section)) allocate (temperature (section)%p(NofMyHomog), source=thermal_initialT(section))
deallocate(temperatureRate(section)%p) deallocate(temperatureRate(section)%p)
allocate (temperatureRate(section)%p(NofMyHomog), source=0.0_pReal) allocate (temperatureRate(section)%p(NofMyHomog), source=0.0_pReal)
endif
enddo initializeInstances enddo initializeInstances
end subroutine thermal_conduction_init end subroutine thermal_conduction_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -187,11 +120,12 @@ end subroutine thermal_conduction_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine thermal_conduction_getSourceAndItsTangent(Tdot, dTdot_dT, T, ip, el) subroutine thermal_conduction_getSourceAndItsTangent(Tdot, dTdot_dT, T, ip, el)
use math, only: & use math, only: &
math_Mandel6to33 math_6toSym33
use material, only: & use material, only: &
homogenization_Ngrains, & homogenization_Ngrains, &
mappingHomogenization, & mappingHomogenization, &
phaseAt, & phaseAt, &
phasememberAt, &
thermal_typeInstance, & thermal_typeInstance, &
phase_Nsources, & phase_Nsources, &
phase_source, & phase_source, &
@ -221,7 +155,8 @@ subroutine thermal_conduction_getSourceAndItsTangent(Tdot, dTdot_dT, T, ip, el)
offset, & offset, &
instance, & instance, &
grain, & grain, &
source source, &
constituent
homog = mappingHomogenization(2,ip,el) homog = mappingHomogenization(2,ip,el)
offset = mappingHomogenization(1,ip,el) offset = mappingHomogenization(1,ip,el)
@ -231,17 +166,18 @@ subroutine thermal_conduction_getSourceAndItsTangent(Tdot, dTdot_dT, T, ip, el)
dTdot_dT = 0.0_pReal dTdot_dT = 0.0_pReal
do grain = 1, homogenization_Ngrains(homog) do grain = 1, homogenization_Ngrains(homog)
phase = phaseAt(grain,ip,el) phase = phaseAt(grain,ip,el)
constituent = phasememberAt(grain,ip,el)
do source = 1, phase_Nsources(phase) do source = 1, phase_Nsources(phase)
select case(phase_source(source,phase)) select case(phase_source(source,phase))
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, &
crystallite_Tstar_v(1:6,grain,ip,el), & math_6toSym33(crystallite_Tstar_v(1:6,grain,ip,el)), &
crystallite_Lp(1:3,1:3,grain,ip,el), & crystallite_Lp(1:3,1:3,grain,ip,el), &
grain, ip, el) phase)
case (SOURCE_thermal_externalheat_ID) case (SOURCE_thermal_externalheat_ID)
call source_thermal_externalheat_getRateAndItsTangent(my_Tdot, my_dTdot_dT, & call source_thermal_externalheat_getRateAndItsTangent(my_Tdot, my_dTdot_dT, &
grain, ip, el) phase, constituent)
case default case default
my_Tdot = 0.0_pReal my_Tdot = 0.0_pReal
@ -258,6 +194,7 @@ subroutine thermal_conduction_getSourceAndItsTangent(Tdot, dTdot_dT, T, ip, el)
end subroutine thermal_conduction_getSourceAndItsTangent end subroutine thermal_conduction_getSourceAndItsTangent
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief returns homogenized thermal conductivity in reference configuration !> @brief returns homogenized thermal conductivity in reference configuration
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -266,7 +203,6 @@ function thermal_conduction_getConductivity33(ip,el)
lattice_thermalConductivity33 lattice_thermalConductivity33
use material, only: & use material, only: &
homogenization_Ngrains, & homogenization_Ngrains, &
mappingHomogenization, &
material_phase material_phase
use mesh, only: & use mesh, only: &
mesh_element mesh_element
@ -280,10 +216,8 @@ function thermal_conduction_getConductivity33(ip,el)
real(pReal), dimension(3,3) :: & real(pReal), dimension(3,3) :: &
thermal_conduction_getConductivity33 thermal_conduction_getConductivity33
integer(pInt) :: & integer(pInt) :: &
homog, &
grain grain
homog = mappingHomogenization(2,ip,el)
thermal_conduction_getConductivity33 = 0.0_pReal thermal_conduction_getConductivity33 = 0.0_pReal
do grain = 1, homogenization_Ngrains(mesh_element(3,el)) do grain = 1, homogenization_Ngrains(mesh_element(3,el))
@ -296,6 +230,7 @@ function thermal_conduction_getConductivity33(ip,el)
end function thermal_conduction_getConductivity33 end function thermal_conduction_getConductivity33
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief returns homogenized specific heat capacity !> @brief returns homogenized specific heat capacity
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -304,12 +239,9 @@ function thermal_conduction_getSpecificHeat(ip,el)
lattice_specificHeat lattice_specificHeat
use material, only: & use material, only: &
homogenization_Ngrains, & homogenization_Ngrains, &
mappingHomogenization, &
material_phase material_phase
use mesh, only: & use mesh, only: &
mesh_element mesh_element
use crystallite, only: &
crystallite_push33ToRef
implicit none implicit none
integer(pInt), intent(in) :: & integer(pInt), intent(in) :: &
@ -318,11 +250,10 @@ function thermal_conduction_getSpecificHeat(ip,el)
real(pReal) :: & real(pReal) :: &
thermal_conduction_getSpecificHeat thermal_conduction_getSpecificHeat
integer(pInt) :: & integer(pInt) :: &
homog, grain grain
thermal_conduction_getSpecificHeat = 0.0_pReal thermal_conduction_getSpecificHeat = 0.0_pReal
homog = mappingHomogenization(2,ip,el)
do grain = 1, homogenization_Ngrains(mesh_element(3,el)) do grain = 1, homogenization_Ngrains(mesh_element(3,el))
thermal_conduction_getSpecificHeat = thermal_conduction_getSpecificHeat + & thermal_conduction_getSpecificHeat = thermal_conduction_getSpecificHeat + &
@ -342,12 +273,9 @@ function thermal_conduction_getMassDensity(ip,el)
lattice_massDensity lattice_massDensity
use material, only: & use material, only: &
homogenization_Ngrains, & homogenization_Ngrains, &
mappingHomogenization, &
material_phase material_phase
use mesh, only: & use mesh, only: &
mesh_element mesh_element
use crystallite, only: &
crystallite_push33ToRef
implicit none implicit none
integer(pInt), intent(in) :: & integer(pInt), intent(in) :: &
@ -356,15 +284,14 @@ function thermal_conduction_getMassDensity(ip,el)
real(pReal) :: & real(pReal) :: &
thermal_conduction_getMassDensity thermal_conduction_getMassDensity
integer(pInt) :: & integer(pInt) :: &
homog, grain grain
thermal_conduction_getMassDensity = 0.0_pReal thermal_conduction_getMassDensity = 0.0_pReal
homog = mappingHomogenization(2,ip,el)
do grain = 1, homogenization_Ngrains(mesh_element(3,el)) do grain = 1, homogenization_Ngrains(mesh_element(3,el))
thermal_conduction_getMassDensity = thermal_conduction_getMassDensity + & thermal_conduction_getMassDensity = thermal_conduction_getMassDensity &
lattice_massDensity(material_phase(grain,ip,el)) + lattice_massDensity(material_phase(grain,ip,el))
enddo enddo
thermal_conduction_getMassDensity = & thermal_conduction_getMassDensity = &
@ -372,6 +299,7 @@ function thermal_conduction_getMassDensity(ip,el)
end function thermal_conduction_getMassDensity end function thermal_conduction_getMassDensity
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief updates thermal state with solution from heat conduction PDE !> @brief updates thermal state with solution from heat conduction PDE
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -400,41 +328,37 @@ subroutine thermal_conduction_putTemperatureAndItsRate(T,Tdot,ip,el)
end subroutine thermal_conduction_putTemperatureAndItsRate end subroutine thermal_conduction_putTemperatureAndItsRate
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief return array of thermal results !> @brief return array of thermal results
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function thermal_conduction_postResults(ip,el) function thermal_conduction_postResults(homog,instance,of) result(postResults)
use material, only: & use material, only: &
mappingHomogenization, & temperature
thermal_typeInstance, &
temperature, &
thermalMapping
implicit none implicit none
integer(pInt), intent(in) :: & integer(pInt), intent(in) :: &
ip, & !< integration point homog, &
el !< element instance, &
real(pReal), dimension(thermal_conduction_sizePostResults(thermal_typeInstance(mappingHomogenization(2,ip,el)))) :: & of
thermal_conduction_postResults
real(pReal), dimension(sum(thermal_conduction_sizePostResult(:,instance))) :: &
postResults
integer(pInt) :: & integer(pInt) :: &
instance, homog, offset, o, c o, c
homog = mappingHomogenization(2,ip,el)
offset = thermalMapping(homog)%p(ip,el)
instance = thermal_typeInstance(homog)
c = 0_pInt c = 0_pInt
thermal_conduction_postResults = 0.0_pReal
do o = 1_pInt,thermal_conduction_Noutput(instance) do o = 1_pInt,thermal_conduction_Noutput(instance)
select case(thermal_conduction_outputID(o,instance)) select case(thermal_conduction_outputID(o,instance))
case (temperature_ID) case (temperature_ID)
thermal_conduction_postResults(c+1_pInt) = temperature(homog)%p(offset) postResults(c+1_pInt) = temperature(homog)%p(of)
c = c + 1 c = c + 1
end select end select
enddo enddo
end function thermal_conduction_postResults end function thermal_conduction_postResults
end module thermal_conduction end module thermal_conduction

30
src/thermal_isothermal.f90
View File

@ -26,14 +26,14 @@ subroutine thermal_isothermal_init()
pInt pInt
use IO, only: & use IO, only: &
IO_timeStamp IO_timeStamp
use config, only: &
material_Nhomogenization
use material use material
use config
implicit none implicit none
integer(pInt) :: & integer(pInt) :: &
homog, & homog, &
NofMyHomog, & NofMyHomog
sizeState
write(6,'(/,a)') ' <<<+- thermal_'//THERMAL_isothermal_label//' init -+>>>' write(6,'(/,a)') ' <<<+- thermal_'//THERMAL_isothermal_label//' init -+>>>'
write(6,'(a15,a)') ' Current time: ',IO_timeStamp() write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
@ -41,21 +41,19 @@ subroutine thermal_isothermal_init()
initializeInstances: do homog = 1_pInt, material_Nhomogenization initializeInstances: do homog = 1_pInt, material_Nhomogenization
myhomog: if (thermal_type(homog) == THERMAL_isothermal_ID) then if (thermal_type(homog) /= THERMAL_isothermal_ID) cycle
NofMyHomog = count(material_homog == homog) NofMyHomog = count(material_homog == homog)
sizeState = 0_pInt thermalState(homog)%sizeState = 0_pInt
thermalState(homog)%sizeState = sizeState thermalState(homog)%sizePostResults = 0_pInt
thermalState(homog)%sizePostResults = sizeState allocate(thermalState(homog)%state0 (0_pInt,NofMyHomog), source=0.0_pReal)
allocate(thermalState(homog)%state0 (sizeState,NofMyHomog), source=0.0_pReal) allocate(thermalState(homog)%subState0(0_pInt,NofMyHomog), source=0.0_pReal)
allocate(thermalState(homog)%subState0(sizeState,NofMyHomog), source=0.0_pReal) allocate(thermalState(homog)%state (0_pInt,NofMyHomog), source=0.0_pReal)
allocate(thermalState(homog)%state (sizeState,NofMyHomog), source=0.0_pReal)
deallocate(temperature (homog)%p) deallocate(temperature (homog)%p)
allocate (temperature (homog)%p(1), source=thermal_initialT(homog)) allocate (temperature (homog)%p(1), source=thermal_initialT(homog))
deallocate(temperatureRate(homog)%p) deallocate(temperatureRate(homog)%p)
allocate (temperatureRate(homog)%p(1), source=0.0_pReal) allocate (temperatureRate(homog)%p(1), source=0.0_pReal)
endif myhomog
enddo initializeInstances enddo initializeInstances