From e0fa3e0b2609c2429d5656306683b7caba1eabdf Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 30 Dec 2018 20:58:38 +0100 Subject: [PATCH] takeover from 40_XX and 41_XX branch easier to focus on thermal instead of doing all kinematics and sources together --- src/constitutive.f90 | 2 +- src/kinematics_thermal_expansion.f90 | 105 +++------------------------ src/source_thermal_dissipation.f90 | 48 +++++++----- src/source_thermal_externalheat.f90 | 33 ++++++--- src/thermal_adiabatic.f90 | 9 ++- src/thermal_conduction.f90 | 13 ++-- 6 files changed, 76 insertions(+), 134 deletions(-) diff --git a/src/constitutive.f90 b/src/constitutive.f90 index 8294047e7..1c32f73d9 100644 --- a/src/constitutive.f90 +++ b/src/constitutive.f90 @@ -188,7 +188,7 @@ subroutine constitutive_init() call IO_checkAndRewind(FILEUNIT) if (any(phase_kinematics == KINEMATICS_cleavage_opening_ID)) call kinematics_cleavage_opening_init(FILEUNIT) if (any(phase_kinematics == KINEMATICS_slipplane_opening_ID)) call kinematics_slipplane_opening_init(FILEUNIT) - if (any(phase_kinematics == KINEMATICS_thermal_expansion_ID)) call kinematics_thermal_expansion_init(FILEUNIT) + if (any(phase_kinematics == KINEMATICS_thermal_expansion_ID)) call kinematics_thermal_expansion_init if (any(phase_kinematics == KINEMATICS_vacancy_strain_ID)) call kinematics_vacancy_strain_init(FILEUNIT) if (any(phase_kinematics == KINEMATICS_hydrogen_strain_ID)) call kinematics_hydrogen_strain_init(FILEUNIT) close(FILEUNIT) diff --git a/src/kinematics_thermal_expansion.f90 b/src/kinematics_thermal_expansion.f90 index 3d1de3d0a..e8f0d71c7 100644 --- a/src/kinematics_thermal_expansion.f90 +++ b/src/kinematics_thermal_expansion.f90 @@ -10,24 +10,7 @@ module kinematics_thermal_expansion implicit none private - integer(pInt), dimension(:), allocatable, public, protected :: & - kinematics_thermal_expansion_sizePostResults, & !< cumulative size of post results - kinematics_thermal_expansion_offset, & !< which kinematics is my current damage mechanism? - kinematics_thermal_expansion_instance !< instance of damage kinematics mechanism - integer(pInt), dimension(:,:), allocatable, target, public :: & - kinematics_thermal_expansion_sizePostResult !< size of each post result output - - character(len=64), dimension(:,:), allocatable, target, public :: & - kinematics_thermal_expansion_output !< name of each post result output - - integer(pInt), dimension(:), allocatable, target, public :: & - kinematics_thermal_expansion_Noutput !< number of outputs per instance of this damage - -! enum, bind(c) ! ToDo kinematics need state machinery to deal with sizePostResult -! enumerator :: undefined_ID, & ! possible remedy is to decouple having state vars from having output -! thermalexpansionrate_ID ! which means to separate user-defined types tState + tOutput... -! end enum public :: & kinematics_thermal_expansion_init, & kinematics_thermal_expansion_initialStrain, & @@ -40,7 +23,7 @@ contains !> @brief module initialization !> @details reads in material parameters, allocates arrays, and does sanity checks !-------------------------------------------------------------------------------------------------- -subroutine kinematics_thermal_expansion_init(fileUnit) +subroutine kinematics_thermal_expansion_init() #if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800 use, intrinsic :: iso_fortran_env, only: & compiler_version, & @@ -51,37 +34,17 @@ subroutine kinematics_thermal_expansion_init(fileUnit) debug_constitutive,& debug_levelBasic use IO, only: & - IO_read, & - IO_lc, & - IO_getTag, & - IO_isBlank, & - IO_stringPos, & - IO_stringValue, & - IO_floatValue, & - IO_intValue, & - IO_warning, & - IO_error, & - IO_timeStamp, & - IO_EOF + IO_timeStamp use material, only: & phase_kinematics, & - phase_Nkinematics, & - phase_Noutput, & KINEMATICS_thermal_expansion_label, & KINEMATICS_thermal_expansion_ID use config, only: & - material_Nphase, & - MATERIAL_partPhase + config_phase implicit none - integer(pInt), intent(in) :: fileUnit - - integer(pInt), allocatable, dimension(:) :: chunkPos - integer(pInt) :: maxNinstance,phase,instance,kinematics - character(len=65536) :: & - tag = '', & - line = '' - + integer(pInt) maxNinstance + write(6,'(/,a)') ' <<<+- kinematics_'//KINEMATICS_thermal_expansion_LABEL//' init -+>>>' write(6,'(a15,a)') ' Current time: ',IO_timeStamp() #include "compilation_info.f90" @@ -91,58 +54,8 @@ subroutine kinematics_thermal_expansion_init(fileUnit) if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt) & write(6,'(a16,1x,i5,/)') '# instances:',maxNinstance - - allocate(kinematics_thermal_expansion_offset(material_Nphase), source=0_pInt) - allocate(kinematics_thermal_expansion_instance(material_Nphase), source=0_pInt) - do phase = 1, material_Nphase - kinematics_thermal_expansion_instance(phase) = count(phase_kinematics(:,1:phase) == kinematics_thermal_expansion_ID) - do kinematics = 1, phase_Nkinematics(phase) - if (phase_kinematics(kinematics,phase) == kinematics_thermal_expansion_ID) & - kinematics_thermal_expansion_offset(phase) = kinematics - enddo - enddo - - allocate(kinematics_thermal_expansion_sizePostResults(maxNinstance), source=0_pInt) - allocate(kinematics_thermal_expansion_sizePostResult(maxval(phase_Noutput),maxNinstance),source=0_pInt) - allocate(kinematics_thermal_expansion_output(maxval(phase_Noutput),maxNinstance)) - kinematics_thermal_expansion_output = '' - allocate(kinematics_thermal_expansion_Noutput(maxNinstance), source=0_pInt) - rewind(fileUnit) - phase = 0_pInt - do while (trim(line) /= IO_EOF .and. IO_lc(IO_getTag(line,'<','>')) /= MATERIAL_partPhase) ! wind forward to - line = IO_read(fileUnit) - enddo - - parsingFile: do while (trim(line) /= IO_EOF) ! read through sections of phase part - line = IO_read(fileUnit) - if (IO_isBlank(line)) cycle ! skip empty lines - if (IO_getTag(line,'<','>') /= '') then ! stop at next part - line = IO_read(fileUnit, .true.) ! reset IO_read - exit - endif - if (IO_getTag(line,'[',']') /= '') then ! next phase section - phase = phase + 1_pInt ! advance phase section counter - cycle ! skip to next line - endif - if (phase > 0_pInt ) then; if (any(phase_kinematics(:,phase) == KINEMATICS_thermal_expansion_ID)) then ! do not short-circuit here (.and. with next if statemen). It's not safe in Fortran - instance = kinematics_thermal_expansion_instance(phase) ! which instance of my damage is present phase - chunkPos = IO_stringPos(line) - tag = IO_lc(IO_stringValue(line,chunkPos,1_pInt)) ! extract key... - select case(tag) -! case ('(output)') -! output = IO_lc(IO_stringValue(line,chunkPos,2_pInt)) ! ...and corresponding output -! select case(output) -! case ('thermalexpansionrate') -! kinematics_thermal_expansion_Noutput(instance) = kinematics_thermal_expansion_Noutput(instance) + 1_pInt -! kinematics_thermal_expansion_outputID(kinematics_thermal_expansion_Noutput(instance),instance) = & -! thermalexpansionrate_ID -! kinematics_thermal_expansion_output(kinematics_thermal_expansion_Noutput(instance),instance) = output -! ToDo add sizePostResult loop afterwards... - - end select - endif; endif - enddo parsingFile +! ToDo: this subroutine should read in lattice_thermal_expansion. No need to make it a global array end subroutine kinematics_thermal_expansion_init @@ -187,7 +100,7 @@ end function kinematics_thermal_expansion_initialStrain !-------------------------------------------------------------------------------------------------- !> @brief contains the constitutive equation for calculating the velocity gradient !-------------------------------------------------------------------------------------------------- -subroutine kinematics_thermal_expansion_LiAndItsTangent(Li, dLi_dTstar3333, ipc, ip, el) +subroutine kinematics_thermal_expansion_LiAndItsTangent(Li, dLi_dTstar, ipc, ip, el) use material, only: & material_phase, & material_homog, & @@ -206,7 +119,7 @@ subroutine kinematics_thermal_expansion_LiAndItsTangent(Li, dLi_dTstar3333, ipc, real(pReal), intent(out), dimension(3,3) :: & Li !< thermal velocity gradient real(pReal), intent(out), dimension(3,3,3,3) :: & - dLi_dTstar3333 !< derivative of Li with respect to Tstar (4th-order tensor defined to be zero) + dLi_dTstar !< derivative of Li with respect to Tstar (4th-order tensor defined to be zero) integer(pInt) :: & phase, & homog, offset @@ -230,7 +143,7 @@ subroutine kinematics_thermal_expansion_LiAndItsTangent(Li, dLi_dTstar3333, ipc, + lattice_thermalExpansion33(1:3,1:3,2,phase)*(T - TRef)**2 / 2. & + lattice_thermalExpansion33(1:3,1:3,3,phase)*(T - TRef)**3 / 3. & ) - dLi_dTstar3333 = 0.0_pReal + dLi_dTstar = 0.0_pReal end subroutine kinematics_thermal_expansion_LiAndItsTangent diff --git a/src/source_thermal_dissipation.f90 b/src/source_thermal_dissipation.f90 index 994d26b41..290ad7efe 100644 --- a/src/source_thermal_dissipation.f90 +++ b/src/source_thermal_dissipation.f90 @@ -27,6 +27,15 @@ module source_thermal_dissipation real(pReal), dimension(:), allocatable, private :: & source_thermal_dissipation_coldworkCoeff + + type, private :: tParameters !< container type for internal constitutive parameters + real(pReal) :: & + coldworkCoeff + end type tParameters + + type(tParameters), dimension(:), allocatable, private :: param !< containers of constitutive parameters (len Ninstance) + + public :: & source_thermal_dissipation_init, & source_thermal_dissipation_getRateAndItsTangent @@ -70,6 +79,7 @@ subroutine source_thermal_dissipation_init(fileUnit) material_phase, & sourceState use config, only: & + config_phase, & material_Nphase, & MATERIAL_partPhase use numerics,only: & @@ -79,9 +89,9 @@ subroutine source_thermal_dissipation_init(fileUnit) integer(pInt), intent(in) :: fileUnit integer(pInt), allocatable, dimension(:) :: chunkPos - integer(pInt) :: maxNinstance,phase,instance,source,sourceOffset + integer(pInt) :: Ninstance,phase,instance,source,sourceOffset integer(pInt) :: sizeState, sizeDotState, sizeDeltaState - integer(pInt) :: NofMyPhase + integer(pInt) :: NofMyPhase,p character(len=65536) :: & tag = '', & line = '' @@ -90,10 +100,10 @@ subroutine source_thermal_dissipation_init(fileUnit) write(6,'(a15,a)') ' Current time: ',IO_timeStamp() #include "compilation_info.f90" - maxNinstance = int(count(phase_source == SOURCE_thermal_dissipation_ID),pInt) - if (maxNinstance == 0_pInt) return + Ninstance = int(count(phase_source == SOURCE_thermal_dissipation_ID),pInt) + if (Ninstance == 0_pInt) return if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt) & - write(6,'(a16,1x,i5,/)') '# instances:',maxNinstance + write(6,'(a16,1x,i5,/)') '# instances:',Ninstance allocate(source_thermal_dissipation_offset(material_Nphase), source=0_pInt) allocate(source_thermal_dissipation_instance(material_Nphase), source=0_pInt) @@ -105,12 +115,17 @@ subroutine source_thermal_dissipation_init(fileUnit) enddo enddo - allocate(source_thermal_dissipation_sizePostResults(maxNinstance), source=0_pInt) - allocate(source_thermal_dissipation_sizePostResult(maxval(phase_Noutput),maxNinstance),source=0_pInt) - allocate(source_thermal_dissipation_output (maxval(phase_Noutput),maxNinstance)) + allocate(source_thermal_dissipation_sizePostResults(Ninstance), source=0_pInt) + allocate(source_thermal_dissipation_sizePostResult(maxval(phase_Noutput),Ninstance),source=0_pInt) + allocate(source_thermal_dissipation_output (maxval(phase_Noutput),Ninstance)) source_thermal_dissipation_output = '' - allocate(source_thermal_dissipation_Noutput(maxNinstance), source=0_pInt) - allocate(source_thermal_dissipation_coldworkCoeff(maxNinstance), source=0.0_pReal) + allocate(source_thermal_dissipation_Noutput(Ninstance), source=0_pInt) + + allocate(source_thermal_dissipation_coldworkCoeff(Ninstance), source=0.0_pReal) + + do p=1, size(config_phase) + if (all(phase_source(:,p) /= SOURCE_thermal_dissipation_ID)) cycle + enddo rewind(fileUnit) phase = 0_pInt @@ -181,17 +196,14 @@ end subroutine source_thermal_dissipation_init !-------------------------------------------------------------------------------------------------- !> @brief returns local vacancy generation rate !-------------------------------------------------------------------------------------------------- -subroutine source_thermal_dissipation_getRateAndItsTangent(TDot, dTDOT_dT, Tstar_v, Lp, ipc, ip, el) +subroutine source_thermal_dissipation_getRateAndItsTangent(TDot, dTDOT_dT, Tstar_v, Lp, phase, constituent) use math, only: & math_Mandel6to33 - use material, only: & - phaseAt, phasememberAt implicit none integer(pInt), intent(in) :: & - ipc, & !< grain number - ip, & !< integration point number - el !< element number + phase, & + constituent real(pReal), intent(in), dimension(6) :: & Tstar_v !< 2nd Piola Kirchhoff stress tensor (Mandel) real(pReal), intent(in), dimension(3,3) :: & @@ -200,10 +212,8 @@ subroutine source_thermal_dissipation_getRateAndItsTangent(TDot, dTDOT_dT, Tstar TDot, & dTDOT_dT integer(pInt) :: & - instance, phase, constituent + instance - phase = phaseAt(ipc,ip,el) - constituent = phasememberAt(ipc,ip,el) instance = source_thermal_dissipation_instance(phase) TDot = source_thermal_dissipation_coldworkCoeff(instance)* & diff --git a/src/source_thermal_externalheat.f90 b/src/source_thermal_externalheat.f90 index b7151aece..eac1232f3 100644 --- a/src/source_thermal_externalheat.f90 +++ b/src/source_thermal_externalheat.f90 @@ -32,6 +32,18 @@ module source_thermal_externalheat source_thermal_externalheat_time, & source_thermal_externalheat_rate + + type, private :: tParameters !< container type for internal constitutive parameters + real(pReal), dimension(:), allocatable :: & + time, & + rate + integer(pInt) :: & + nInterval + end type tParameters + + type(tParameters), dimension(:), allocatable, private :: param !< containers of constitutive parameters (len Ninstance) + + public :: & source_thermal_externalheat_init, & source_thermal_externalheat_dotState, & @@ -76,6 +88,7 @@ subroutine source_thermal_externalheat_init(fileUnit) material_phase, & sourceState use config, only: & + config_phase, & material_Nphase, & MATERIAL_partPhase use numerics,only: & @@ -87,7 +100,7 @@ subroutine source_thermal_externalheat_init(fileUnit) integer(pInt), allocatable, dimension(:) :: chunkPos integer(pInt) :: maxNinstance,phase,instance,source,sourceOffset integer(pInt) :: sizeState, sizeDotState, sizeDeltaState - integer(pInt) :: NofMyPhase,interval + integer(pInt) :: NofMyPhase,interval,p character(len=65536) :: & tag = '', & line = '' @@ -117,11 +130,15 @@ subroutine source_thermal_externalheat_init(fileUnit) allocate(source_thermal_externalheat_output (maxval(phase_Noutput),maxNinstance)) source_thermal_externalheat_output = '' allocate(source_thermal_externalheat_Noutput(maxNinstance), source=0_pInt) - allocate(source_thermal_externalheat_nIntervals(maxNinstance), source=0_pInt) + allocate(source_thermal_externalheat_nIntervals(maxNinstance), source=0_pInt) allocate(temp_time(maxNinstance,1000), source=0.0_pReal) allocate(temp_rate(maxNinstance,1000), source=0.0_pReal) + do p=1, size(config_phase) + if (all(phase_source(:,p) /= SOURCE_thermal_externalheat_ID)) cycle + enddo + rewind(fileUnit) phase = 0_pInt do while (trim(line) /= IO_EOF .and. IO_lc(IO_getTag(line,'<','>')) /= MATERIAL_partPhase) ! wind forward to @@ -238,26 +255,22 @@ end subroutine source_thermal_externalheat_dotState !-------------------------------------------------------------------------------------------------- !> @brief returns local heat generation rate !-------------------------------------------------------------------------------------------------- -subroutine source_thermal_externalheat_getRateAndItsTangent(TDot, dTDot_dT, ipc, ip, el) +subroutine source_thermal_externalheat_getRateAndItsTangent(TDot, dTDot_dT, phase, constituent) use material, only: & - phaseAt, phasememberAt, & sourceState implicit none integer(pInt), intent(in) :: & - ipc, & !< grain number - ip, & !< integration point number - el !< element number + phase, & + constituent real(pReal), intent(out) :: & TDot, & dTDot_dT integer(pInt) :: & - instance, phase, constituent, sourceOffset, interval + instance, sourceOffset, interval real(pReal) :: & frac_time - phase = phaseAt(ipc,ip,el) - constituent = phasememberAt(ipc,ip,el) instance = source_thermal_externalheat_instance(phase) sourceOffset = source_thermal_externalheat_offset(phase) diff --git a/src/thermal_adiabatic.f90 b/src/thermal_adiabatic.f90 index 6a70ca7ee..dd52293eb 100644 --- a/src/thermal_adiabatic.f90 +++ b/src/thermal_adiabatic.f90 @@ -238,6 +238,7 @@ subroutine thermal_adiabatic_getSourceAndItsTangent(Tdot, dTdot_dT, T, ip, el) homogenization_Ngrains, & mappingHomogenization, & phaseAt, & + phasememberAt, & thermal_typeInstance, & phase_Nsources, & phase_source, & @@ -267,7 +268,8 @@ subroutine thermal_adiabatic_getSourceAndItsTangent(Tdot, dTdot_dT, T, ip, el) offset, & instance, & grain, & - source + source, & + constituent homog = mappingHomogenization(2,ip,el) offset = mappingHomogenization(1,ip,el) @@ -277,17 +279,18 @@ subroutine thermal_adiabatic_getSourceAndItsTangent(Tdot, dTdot_dT, T, ip, el) dTdot_dT = 0.0_pReal do grain = 1, homogenization_Ngrains(homog) phase = phaseAt(grain,ip,el) + constituent = phasememberAt(grain,ip,el) do source = 1, phase_Nsources(phase) select case(phase_source(source,phase)) case (SOURCE_thermal_dissipation_ID) call source_thermal_dissipation_getRateAndItsTangent(my_Tdot, my_dTdot_dT, & crystallite_Tstar_v(1:6,grain,ip,el), & crystallite_Lp(1:3,1:3,grain,ip,el), & - grain, ip, el) + phase, constituent) case (SOURCE_thermal_externalheat_ID) call source_thermal_externalheat_getRateAndItsTangent(my_Tdot, my_dTdot_dT, & - grain, ip, el) + phase, constituent) case default my_Tdot = 0.0_pReal diff --git a/src/thermal_conduction.f90 b/src/thermal_conduction.f90 index 16497040b..1a9014e89 100644 --- a/src/thermal_conduction.f90 +++ b/src/thermal_conduction.f90 @@ -192,6 +192,7 @@ subroutine thermal_conduction_getSourceAndItsTangent(Tdot, dTdot_dT, T, ip, el) homogenization_Ngrains, & mappingHomogenization, & phaseAt, & + phasememberAt, & thermal_typeInstance, & phase_Nsources, & phase_source, & @@ -221,7 +222,8 @@ subroutine thermal_conduction_getSourceAndItsTangent(Tdot, dTdot_dT, T, ip, el) offset, & instance, & grain, & - source + source, & + constituent homog = mappingHomogenization(2,ip,el) offset = mappingHomogenization(1,ip,el) @@ -231,17 +233,18 @@ subroutine thermal_conduction_getSourceAndItsTangent(Tdot, dTdot_dT, T, ip, el) dTdot_dT = 0.0_pReal do grain = 1, homogenization_Ngrains(homog) phase = phaseAt(grain,ip,el) + constituent = phasememberAt(grain,ip,el) do source = 1, phase_Nsources(phase) select case(phase_source(source,phase)) case (SOURCE_thermal_dissipation_ID) call source_thermal_dissipation_getRateAndItsTangent(my_Tdot, my_dTdot_dT, & crystallite_Tstar_v(1:6,grain,ip,el), & crystallite_Lp(1:3,1:3,grain,ip,el), & - grain, ip, el) + phase, constituent) case (SOURCE_thermal_externalheat_ID) call source_thermal_externalheat_getRateAndItsTangent(my_Tdot, my_dTdot_dT, & - grain, ip, el) + phase, constituent) case default my_Tdot = 0.0_pReal @@ -363,8 +366,8 @@ function thermal_conduction_getMassDensity(ip,el) homog = mappingHomogenization(2,ip,el) do grain = 1, homogenization_Ngrains(mesh_element(3,el)) - thermal_conduction_getMassDensity = thermal_conduction_getMassDensity + & - lattice_massDensity(material_phase(grain,ip,el)) + thermal_conduction_getMassDensity = thermal_conduction_getMassDensity & + + lattice_massDensity(material_phase(grain,ip,el)) enddo thermal_conduction_getMassDensity = &