From c97b83c9a13b6e2bb57c77e8bb227384a32b82b8 Mon Sep 17 00:00:00 2001 From: Pratheek Shanthraj Date: Mon, 27 Jul 2015 11:09:37 +0000 Subject: [PATCH] added user defined external heat source --- code/Makefile | 5 +- code/commercialFEM_fileList.f90 | 1 + code/constitutive.f90 | 25 ++- code/material.f90 | 5 + code/source_thermal_externalheat.f90 | 276 +++++++++++++++++++++++++++ code/thermal_adiabatic.f90 | 11 +- code/thermal_conduction.f90 | 9 +- 7 files changed, 322 insertions(+), 10 deletions(-) create mode 100644 code/source_thermal_externalheat.f90 diff --git a/code/Makefile b/code/Makefile index db0203c5c..dd5793a8c 100644 --- a/code/Makefile +++ b/code/Makefile @@ -307,7 +307,7 @@ COMPILE =$(OPENMP_FLAG_$(F90)) $(COMPILE_OPTIONS_$(F90)) $(STANDARD_CHEC COMPILE_MAXOPTI =$(OPENMP_FLAG_$(F90)) $(COMPILE_OPTIONS_$(F90)) $(STANDARD_CHECK_$(F90)) $(OPTIMIZATION_$(MAXOPTI)_$(F90)) $(INCLUDE_DIRS) $(PRECISION_$(F90)) ################################################################################################### SOURCE_FILES = \ - source_thermal_dissipation.o \ + source_thermal_dissipation.o source_thermal_externalheat.o \ source_damage_isoBrittle.o source_damage_isoDuctile.o source_damage_anisoBrittle.o source_damage_anisoDuctile.o \ source_vacancy_phenoplasticity.o source_vacancy_irradiation.o source_vacancy_thermalfluc.o @@ -517,6 +517,9 @@ constitutive.o: constitutive.f90 \ source_thermal_dissipation.o: source_thermal_dissipation.f90 \ lattice.o +source_thermal_externalheat.o: source_thermal_externalheat.f90 \ + lattice.o + source_damage_isoBrittle.o: source_damage_isoBrittle.f90 \ lattice.o diff --git a/code/commercialFEM_fileList.f90 b/code/commercialFEM_fileList.f90 index ffdb0015c..5d09b39ff 100644 --- a/code/commercialFEM_fileList.f90 +++ b/code/commercialFEM_fileList.f90 @@ -15,6 +15,7 @@ #include "material.f90" #include "lattice.f90" #include "source_thermal_dissipation.f90" +#include "source_thermal_externalheat.f90" #include "source_damage_isoBrittle.f90" #include "source_damage_isoDuctile.f90" #include "source_damage_anisoBrittle.f90" diff --git a/code/constitutive.f90 b/code/constitutive.f90 index 5c9015f02..353d41153 100644 --- a/code/constitutive.f90 +++ b/code/constitutive.f90 @@ -84,6 +84,7 @@ subroutine constitutive_init() PLASTICITY_titanmod_ID, & PLASTICITY_nonlocal_ID ,& SOURCE_thermal_dissipation_ID, & + SOURCE_thermal_externalheat_ID, & SOURCE_damage_isoBrittle_ID, & SOURCE_damage_isoDuctile_ID, & SOURCE_damage_anisoBrittle_ID, & @@ -106,6 +107,7 @@ subroutine constitutive_init() PLASTICITY_TITANMOD_label, & PLASTICITY_NONLOCAL_label, & SOURCE_thermal_dissipation_label, & + SOURCE_thermal_externalheat_label, & SOURCE_damage_isoBrittle_label, & SOURCE_damage_isoDuctile_label, & SOURCE_damage_anisoBrittle_label, & @@ -125,6 +127,7 @@ subroutine constitutive_init() use plastic_titanmod use plastic_nonlocal use source_thermal_dissipation + use source_thermal_externalheat use source_damage_isoBrittle use source_damage_isoDuctile use source_damage_anisoBrittle @@ -175,6 +178,7 @@ subroutine constitutive_init() 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(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(FILEUNIT) if (any(phase_source == SOURCE_damage_isoDuctile_ID)) call source_damage_isoDuctile_init(FILEUNIT) if (any(phase_source == SOURCE_damage_anisoBrittle_ID)) call source_damage_anisoBrittle_init(FILEUNIT) @@ -272,6 +276,12 @@ subroutine constitutive_init() thisNoutput => source_thermal_dissipation_Noutput thisOutput => source_thermal_dissipation_output thisSize => source_thermal_dissipation_sizePostResult + case (SOURCE_thermal_externalheat_ID) + instance = source_thermal_externalheat_instance(phase) + outputName = SOURCE_thermal_externalheat_label + thisNoutput => source_thermal_externalheat_Noutput + thisOutput => source_thermal_externalheat_output + thisSize => source_thermal_externalheat_sizePostResult case (SOURCE_damage_isoBrittle_ID) instance = source_damage_isoBrittle_instance(phase) outputName = SOURCE_damage_isoBrittle_label @@ -745,8 +755,6 @@ pure function constitutive_initialFi(ipc, ip, el) el !< element number real(pReal), dimension(3,3) :: & constitutive_initialFi !< composite initial intermediate deformation gradient - real(pReal), dimension(3,3) :: & - my_Fi !< individual intermediate deformation gradients integer(pInt) :: & kinematics @@ -911,7 +919,8 @@ subroutine constitutive_collectDotState(Tstar_v, FeArray, FpArray, subdt, subfra PLASTICITY_nonlocal_ID, & SOURCE_damage_isoDuctile_ID, & SOURCE_damage_anisoBrittle_ID, & - SOURCE_damage_anisoDuctile_ID + SOURCE_damage_anisoDuctile_ID, & + SOURCE_thermal_externalheat_ID use plastic_j2, only: & plastic_j2_dotState use plastic_phenopowerlaw, only: & @@ -932,6 +941,8 @@ subroutine constitutive_collectDotState(Tstar_v, FeArray, FpArray, subdt, subfra source_damage_anisoBrittle_dotState use source_damage_anisoDuctile, only: & source_damage_anisoDuctile_dotState + use source_thermal_externalheat, only: & + source_thermal_externalheat_dotState implicit none integer(pInt), intent(in) :: & @@ -985,11 +996,13 @@ subroutine constitutive_collectDotState(Tstar_v, FeArray, FpArray, subdt, subfra do mySource = 1_pInt, phase_Nsources(phase) select case (phase_source(mySource,phase)) case (SOURCE_damage_anisoBrittle_ID) - call source_damage_anisoBrittle_dotState(Tstar_v, ipc, ip, el) + call source_damage_anisoBrittle_dotState (Tstar_v, ipc, ip, el) case (SOURCE_damage_isoDuctile_ID) - call source_damage_isoDuctile_dotState ( ipc, ip, el) + call source_damage_isoDuctile_dotState ( ipc, ip, el) case (SOURCE_damage_anisoDuctile_ID) - call source_damage_anisoDuctile_dotState( ipc, ip, el) + call source_damage_anisoDuctile_dotState ( ipc, ip, el) + case (SOURCE_thermal_externalheat_ID) + call source_thermal_externalheat_dotState( ipc, ip, el) end select enddo diff --git a/code/material.f90 b/code/material.f90 index 8a117ed65..21850208e 100644 --- a/code/material.f90 +++ b/code/material.f90 @@ -33,6 +33,7 @@ module material PLASTICITY_titanmod_label = 'titanmod', & PLASTICITY_nonlocal_label = 'nonlocal', & SOURCE_thermal_dissipation_label = 'thermal_dissipation', & + SOURCE_thermal_externalheat_label = 'thermal_externalheat', & SOURCE_damage_isoBrittle_label = 'damage_isobrittle', & SOURCE_damage_isoDuctile_label = 'damage_isoductile', & SOURCE_damage_anisoBrittle_label = 'damage_anisobrittle', & @@ -85,6 +86,7 @@ module material enum, bind(c) enumerator :: SOURCE_undefined_ID, & SOURCE_thermal_dissipation_ID, & + SOURCE_thermal_externalheat_ID, & SOURCE_damage_isoBrittle_ID, & SOURCE_damage_isoDuctile_ID, & SOURCE_damage_anisoBrittle_ID, & @@ -317,6 +319,7 @@ module material PLASTICITY_titanmod_ID, & PLASTICITY_nonlocal_ID, & SOURCE_thermal_dissipation_ID, & + SOURCE_thermal_externalheat_ID, & SOURCE_damage_isoBrittle_ID, & SOURCE_damage_isoDuctile_ID, & SOURCE_damage_anisoBrittle_ID, & @@ -999,6 +1002,8 @@ subroutine material_parsePhase(fileUnit,myPart) select case (IO_lc(IO_stringValue(line,positions,2_pInt))) case (SOURCE_thermal_dissipation_label) phase_source(sourceCtr,section) = SOURCE_thermal_dissipation_ID + case (SOURCE_thermal_externalheat_label) + phase_source(sourceCtr,section) = SOURCE_thermal_externalheat_ID case (SOURCE_damage_isoBrittle_label) phase_source(sourceCtr,section) = SOURCE_damage_isoBrittle_ID case (SOURCE_damage_isoDuctile_label) diff --git a/code/source_thermal_externalheat.f90 b/code/source_thermal_externalheat.f90 new file mode 100644 index 000000000..af487e2a0 --- /dev/null +++ b/code/source_thermal_externalheat.f90 @@ -0,0 +1,276 @@ +!-------------------------------------------------------------------------------------------------- +! $Id$ +!-------------------------------------------------------------------------------------------------- +!> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH +!> @brief material subroutine for thermal source due to plastic dissipation +!> @details to be done +!-------------------------------------------------------------------------------------------------- +module source_thermal_externalheat + use prec, only: & + pReal, & + pInt + + implicit none + private + 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_instance !< instance of thermal dissipation source mechanism + + integer(pInt), dimension(:,:), allocatable, target, public :: & + source_thermal_externalheat_sizePostResult !< size of each post result output + + character(len=64), dimension(:,:), allocatable, target, public :: & + source_thermal_externalheat_output !< name of each post result output + + integer(pInt), dimension(:), allocatable, target, public :: & + source_thermal_externalheat_Noutput !< number of outputs per instance of this source + + integer(pInt), dimension(:), allocatable, private :: & + source_thermal_externalheat_nIntervals + + real(pReal), dimension(:,:), allocatable, private :: & + source_thermal_externalheat_time, & + source_thermal_externalheat_rate + + public :: & + source_thermal_externalheat_init, & + source_thermal_externalheat_dotState, & + source_thermal_externalheat_getRateAndItsTangent + +contains + + +!-------------------------------------------------------------------------------------------------- +!> @brief module initialization +!> @details reads in material parameters, allocates arrays, and does sanity checks +!-------------------------------------------------------------------------------------------------- +subroutine source_thermal_externalheat_init(fileUnit) + use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment) + use debug, only: & + debug_level,& + 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 + use material, only: & + phase_source, & + phase_Nsources, & + phase_Noutput, & + SOURCE_thermal_externalheat_label, & + SOURCE_thermal_externalheat_ID, & + material_Nphase, & + material_phase, & + sourceState, & + MATERIAL_partPhase + use numerics,only: & + worldrank, & + numerics_integrator + + implicit none + integer(pInt), intent(in) :: fileUnit + + integer(pInt), parameter :: MAXNCHUNKS = 7_pInt + integer(pInt), dimension(1+2*MAXNCHUNKS) :: positions + integer(pInt) :: maxNinstance,phase,instance,source,sourceOffset + integer(pInt) :: sizeState, sizeDotState, sizeDeltaState + integer(pInt) :: NofMyPhase,interval + character(len=65536) :: & + tag = '', & + line = '' + real(pReal), allocatable, dimension(:,:) :: temp_time, temp_rate + + mainProcess: if (worldrank == 0) then + write(6,'(/,a)') ' <<<+- source_'//SOURCE_thermal_externalheat_label//' init -+>>>' + write(6,'(a)') ' $Id$' + write(6,'(a15,a)') ' Current time: ',IO_timeStamp() +#include "compilation_info.f90" + endif mainProcess + + 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) + allocate(source_thermal_externalheat_instance(material_Nphase), source=0_pInt) + do phase = 1, material_Nphase + source_thermal_externalheat_instance(phase) = count(phase_source(:,1:phase) == SOURCE_thermal_externalheat_ID) + 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_sizePostResult(maxval(phase_Noutput),maxNinstance),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) + allocate(temp_rate(maxNinstance,1000), source=0.0_pReal) + + 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_source(:,phase) == SOURCE_thermal_externalheat_ID)) then ! do not short-circuit here (.and. with next if statemen). It's not safe in Fortran + + instance = source_thermal_externalheat_instance(phase) ! which instance of my source is present phase + positions = IO_stringPos(line,MAXNCHUNKS) + tag = IO_lc(IO_stringValue(line,positions,1_pInt)) ! extract key + select case(tag) + case ('externalheat_time') + if (positions(1) <= 2_pInt) & + call IO_error(150_pInt,ext_msg=trim(tag)//' ('//SOURCE_thermal_externalheat_label//')') + source_thermal_externalheat_nIntervals(instance) = positions(1) - 2_pInt + do interval = 1, source_thermal_externalheat_nIntervals(instance) + 1_pInt + temp_time(instance, interval) = IO_floatValue(line,positions,1_pInt + interval) + enddo + + case ('externalheat_rate') + do interval = 1, source_thermal_externalheat_nIntervals(instance) + 1_pInt + temp_rate(instance, interval) = IO_floatValue(line,positions,1_pInt + interval) + enddo + + end select + endif; endif + enddo parsingFile + + allocate(source_thermal_externalheat_time(maxNinstance,maxval(source_thermal_externalheat_nIntervals)+1_pInt), source=0.0_pReal) + allocate(source_thermal_externalheat_rate(maxNinstance,maxval(source_thermal_externalheat_nIntervals)+1_pInt), source=0.0_pReal) + + initializeInstances: do phase = 1_pInt, material_Nphase + 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 + 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)%state_backup (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) + allocate(sourceState(phase)%p(sourceOffset)%dotState_backup (sizeDotState,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 + +!-------------------------------------------------------------------------------------------------- +!> @brief calculates derived quantities from state +!-------------------------------------------------------------------------------------------------- +subroutine source_thermal_externalheat_dotState(ipc, ip, el) + use material, only: & + mappingConstitutive, & + sourceState + + implicit none + integer(pInt), intent(in) :: & + ipc, & !< component-ID of integration point + ip, & !< integration point + el !< element + integer(pInt) :: & + phase, & + constituent, & + sourceOffset + + phase = mappingConstitutive(2,ipc,ip,el) + constituent = mappingConstitutive(1,ipc,ip,el) + sourceOffset = source_thermal_externalheat_offset(phase) + + sourceState(phase)%p(sourceOffset)%dotState(1,constituent) = 1.0_pReal + +end subroutine source_thermal_externalheat_dotState + +!-------------------------------------------------------------------------------------------------- +!> @brief returns local vacancy generation rate +!-------------------------------------------------------------------------------------------------- +subroutine source_thermal_externalheat_getRateAndItsTangent(TDot, dTDot_dT, ipc, ip, el) + use material, only: & + mappingConstitutive, & + sourceState + + implicit none + integer(pInt), intent(in) :: & + ipc, & !< grain number + ip, & !< integration point number + el !< element number + real(pReal), intent(out) :: & + TDot, & + dTDot_dT + integer(pInt) :: & + instance, phase, constituent, sourceOffset, interval + real(pReal) :: & + norm_time + + phase = mappingConstitutive(2,ipc,ip,el) + constituent = mappingConstitutive(1,ipc,ip,el) + instance = source_thermal_externalheat_instance(phase) + sourceOffset = source_thermal_externalheat_offset(phase) + + do interval = 1, source_thermal_externalheat_nIntervals(instance) + norm_time = (sourceState(phase)%p(sourceOffset)%state(1,constituent) - & + source_thermal_externalheat_time(instance,interval)) / & + (source_thermal_externalheat_time(instance,interval+1) - & + source_thermal_externalheat_time(instance,interval)) + if (norm_time >= 0.0_pReal .and. norm_time < 1.0_pReal) & + TDot = source_thermal_externalheat_rate(instance,interval ) * (1.0_pReal - norm_time) + & + source_thermal_externalheat_rate(instance,interval+1) * norm_time + enddo + dTDot_dT = 0.0 + +end subroutine source_thermal_externalheat_getRateAndItsTangent + +end module source_thermal_externalheat diff --git a/code/thermal_adiabatic.f90 b/code/thermal_adiabatic.f90 index 62e0d58af..80debf2fc 100644 --- a/code/thermal_adiabatic.f90 +++ b/code/thermal_adiabatic.f90 @@ -217,7 +217,7 @@ function thermal_adiabatic_updateState(subdt, ip, el) T = thermalState(homog)%subState0(1,offset) call thermal_adiabatic_getSourceAndItsTangent(Tdot, dTdot_dT, T, ip, el) - T = T + subdt*thermal_adiabatic_getSpecificHeat(ip,el)*thermal_adiabatic_getMassDensity(ip,el)*Tdot + T = T + subdt*Tdot/(thermal_adiabatic_getSpecificHeat(ip,el)*thermal_adiabatic_getMassDensity(ip,el)) thermal_adiabatic_updateState = [ abs(T - thermalState(homog)%state(1,offset)) & <= err_thermal_tolAbs & @@ -244,9 +244,12 @@ subroutine thermal_adiabatic_getSourceAndItsTangent(Tdot, dTdot_dT, T, ip, el) thermal_typeInstance, & phase_Nsources, & phase_source, & - SOURCE_thermal_dissipation_ID + SOURCE_thermal_dissipation_ID, & + SOURCE_thermal_externalheat_ID use source_thermal_dissipation, only: & source_thermal_dissipation_getRateAndItsTangent + use source_thermal_externalheat, only: & + source_thermal_externalheat_getRateAndItsTangent use crystallite, only: & crystallite_Tstar_v, & crystallite_Lp @@ -285,6 +288,10 @@ subroutine thermal_adiabatic_getSourceAndItsTangent(Tdot, dTdot_dT, T, ip, el) crystallite_Lp(1:3,1:3,grain,ip,el), & grain, ip, el) + case (SOURCE_thermal_externalheat_ID) + call source_thermal_externalheat_getRateAndItsTangent(my_Tdot, my_dTdot_dT, & + grain, ip, el) + case default my_Tdot = 0.0_pReal my_dTdot_dT = 0.0_pReal diff --git a/code/thermal_conduction.f90 b/code/thermal_conduction.f90 index 2bbfb7b35..15ae8e99a 100644 --- a/code/thermal_conduction.f90 +++ b/code/thermal_conduction.f90 @@ -198,9 +198,12 @@ subroutine thermal_conduction_getSourceAndItsTangent(Tdot, dTdot_dT, T, ip, el) thermal_typeInstance, & phase_Nsources, & phase_source, & - SOURCE_thermal_dissipation_ID + SOURCE_thermal_dissipation_ID, & + SOURCE_thermal_externalheat_ID use source_thermal_dissipation, only: & source_thermal_dissipation_getRateAndItsTangent + use source_thermal_externalheat, only: & + source_thermal_externalheat_getRateAndItsTangent use crystallite, only: & crystallite_Tstar_v, & crystallite_Lp @@ -239,6 +242,10 @@ subroutine thermal_conduction_getSourceAndItsTangent(Tdot, dTdot_dT, T, ip, el) crystallite_Lp(1:3,1:3,grain,ip,el), & grain, ip, el) + case (SOURCE_thermal_externalheat_ID) + call source_thermal_externalheat_getRateAndItsTangent(my_Tdot, my_dTdot_dT, & + grain, ip, el) + case default my_Tdot = 0.0_pReal my_dTdot_dT = 0.0_pReal