diff --git a/src/constitutive.f90 b/src/constitutive.f90 index 30c3edf74..b3b5ae875 100644 --- a/src/constitutive.f90 +++ b/src/constitutive.f90 @@ -38,11 +38,6 @@ contains !> @brief allocates arrays pointing to array of the various constitutive modules !-------------------------------------------------------------------------------------------------- subroutine constitutive_init() -#if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800 - use, intrinsic :: iso_fortran_env, only: & - compiler_version, & - compiler_options -#endif use prec, only: & pReal use debug, only: & @@ -53,12 +48,8 @@ subroutine constitutive_init() use IO, only: & IO_error, & IO_open_file, & - IO_checkAndRewind, & IO_open_jobFile_stat, & - IO_write_jobFile, & - IO_timeStamp - use config, only: & - config_phase + IO_write_jobFile use config, only: & material_Nphase, & material_localFileExt, & @@ -138,7 +129,7 @@ subroutine constitutive_init() 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_ISOTROPIC_ID)) call plastic_isotropic_init if (any(phase_plasticity == PLASTICITY_PHENOPOWERLAW_ID)) call plastic_phenopowerlaw_init @@ -146,29 +137,21 @@ subroutine constitutive_init() if (any(phase_plasticity == PLASTICITY_DISLOTWIN_ID)) call plastic_dislotwin_init if (any(phase_plasticity == PLASTICITY_DISLOUCLA_ID)) call plastic_disloucla_init if (any(phase_plasticity == PLASTICITY_NONLOCAL_ID)) call plastic_nonlocal_init - - -!-------------------------------------------------------------------------------------------------- -! 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 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) +! initialize source mechanisms + if (any(phase_source == SOURCE_thermal_dissipation_ID)) call source_thermal_dissipation_init + if (any(phase_source == SOURCE_thermal_externalheat_ID)) call source_thermal_externalheat_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_anisoBrittle_ID)) call source_damage_anisoBrittle_init if (any(phase_source == SOURCE_damage_anisoDuctile_ID)) call source_damage_anisoDuctile_init !-------------------------------------------------------------------------------------------------- -! parse kinematic mechanisms from config file - call IO_checkAndRewind(FILEUNIT) +! initialize kinematic mechanisms 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_thermal_expansion_ID)) call kinematics_thermal_expansion_init(FILEUNIT) - close(FILEUNIT) + if (any(phase_kinematics == KINEMATICS_thermal_expansion_ID)) call kinematics_thermal_expansion_init call config_deallocate('material.config/phase') @@ -481,7 +464,7 @@ subroutine constitutive_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, & case (PLASTICITY_KINEHARDENING_ID) plasticityType of = phasememberAt(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 call plastic_nonlocal_LpAndItsTangent (Lp,dLp_dMp,Mp, & @@ -528,8 +511,7 @@ subroutine constitutive_LiAndItsTangents(Li, dLi_dS, dLi_dFi, & math_I3, & math_inv33, & math_det33, & - math_mul33x33, & - math_6toSym33 + math_mul33x33 use material, only: & phasememberAt, & phase_plasticity, & @@ -632,10 +614,11 @@ pure function constitutive_initialFi(ipc, ip, el) use prec, only: & pReal use math, only: & - math_I3, & - math_inv33, & - math_mul33x33 + math_I3 use material, only: & + material_phase, & + material_homog, & + thermalMapping, & phase_kinematics, & phase_Nkinematics, & material_phase, & @@ -652,14 +635,20 @@ pure function constitutive_initialFi(ipc, ip, el) constitutive_initialFi !< composite initial intermediate deformation gradient integer(pInt) :: & k !< counter in kinematics loop + integer(pInt) :: & + phase, & + homog, offset 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 - kinematicsType: select case (phase_kinematics(k,material_phase(ipc,ip,el))) + KinematicsLoop: do k = 1_pInt, phase_Nkinematics(phase) !< Warning: small initial strain assumption + kinematicsType: select case (phase_kinematics(k,phase)) case (KINEMATICS_thermal_expansion_ID) kinematicsType + homog = material_homog(ip,el) + offset = thermalMapping(homog)%p(ip,el) constitutive_initialFi = & - constitutive_initialFi + kinematics_thermal_expansion_initialStrain(ipc, ip, el) + constitutive_initialFi + kinematics_thermal_expansion_initialStrain(homog,phase,offset) end select kinematicsType enddo KinematicsLoop @@ -764,7 +753,7 @@ end subroutine constitutive_hooke_SandItsTangents !-------------------------------------------------------------------------------------------------- !> @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: & pReal, & pLongInt @@ -774,8 +763,6 @@ subroutine constitutive_collectDotState(S, FeArray, Fi, FpArray, subdt, subfracA debug_levelBasic use math, only: & math_mul33x33, & - math_6toSym33, & - math_sym33to6, & math_mul33x33 use mesh, only: & theMesh @@ -829,8 +816,6 @@ subroutine constitutive_collectDotState(S, FeArray, Fi, FpArray, subdt, subfracA el !< element real(pReal), intent(in) :: & 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) :: & FeArray, & !< elastic 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? 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 - call source_damage_anisoDuctile_dotState ( ipc, ip, el) + call source_damage_anisoDuctile_dotState ( ipc, ip, el) 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 @@ -918,7 +903,6 @@ subroutine constitutive_collectDeltaState(S, Fe, Fi, ipc, ip, el) debug_constitutive, & debug_levelBasic use math, only: & - math_sym33to6, & math_mul33x33 use material, only: & phasememberAt, & @@ -984,14 +968,11 @@ end subroutine constitutive_collectDeltaState !-------------------------------------------------------------------------------------------------- !> @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: & pReal use math, only: & - math_6toSym33, & math_mul33x33 - use mesh, only: & - theMesh use material, only: & phasememberAt, & phase_plasticityInstance, & @@ -1004,7 +985,6 @@ function constitutive_postResults(S, Fi, FeArray, ipc, ip, el) material_homogenizationAt, & temperature, & thermalMapping, & - homogenization_maxNgrains, & PLASTICITY_NONE_ID, & PLASTICITY_ISOTROPIC_ID, & PLASTICITY_PHENOPOWERLAW_ID, & @@ -1047,8 +1027,6 @@ function constitutive_postResults(S, Fi, FeArray, ipc, ip, el) constitutive_postResults real(pReal), intent(in), dimension(3,3) :: & 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) :: & S !< 2nd Piola Kirchhoff stress real(pReal), dimension(3,3) :: & diff --git a/src/crystallite.f90 b/src/crystallite.f90 index d35adb323..5209f1867 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -247,7 +247,7 @@ subroutine crystallite_init allocate(crystallite_sizePostResult(maxval(crystallite_Noutput), & size(config_crystallite)), source=0_pInt) - select case(numerics_integrator(1)) + select case(numerics_integrator) case(1_pInt) integrateState => integrateStateFPI case(2_pInt) @@ -1068,7 +1068,7 @@ function crystallite_postResults(ipc, ip, el) if (size(crystallite_postResults)-c > 0_pInt) & 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), & - crystallite_Fe, ipc, ip, el) + ipc, ip, el) end function crystallite_postResults @@ -2297,7 +2297,7 @@ subroutine update_dotState(timeFraction) crystallite_Fe, & crystallite_Fi(1:3,1:3,g,i,e), & 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) NaN = any(IEEE_is_NaN(plasticState(p)%dotState(:,c))) do s = 1_pInt, phase_Nsources(p) diff --git a/src/homogenization.f90 b/src/homogenization.f90 index 8edecbc88..94acf8c82 100644 --- a/src/homogenization.f90 +++ b/src/homogenization.f90 @@ -108,30 +108,18 @@ subroutine homogenization_init 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_ISOSTRAIN_ID)) call homogenization_isostrain_init if (any(homogenization_type == HOMOGENIZATION_RGC_ID)) call homogenization_RGC_init -!-------------------------------------------------------------------------------------------------- -! parse thermal from config file - call IO_checkAndRewind(FILEUNIT) - 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) + if (any(thermal_type == THERMAL_isothermal_ID)) call thermal_isothermal_init + if (any(thermal_type == THERMAL_adiabatic_ID)) call thermal_adiabatic_init + if (any(thermal_type == THERMAL_conduction_ID)) call thermal_conduction_init !-------------------------------------------------------------------------------------------------- -! parse damage from config file - call IO_checkAndRewind(FILEUNIT) +! 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 if (any(damage_type == DAMAGE_none_ID)) & call damage_none_init() if (any(damage_type == DAMAGE_local_ID)) & @@ -895,6 +883,8 @@ function postResults(ip,el) use mesh, only: & mesh_element use material, only: & + thermalMapping, & + thermal_typeInstance, & material_homogenizationAt, & homogenization_typeInstance,& mappingHomogenization, & @@ -934,7 +924,7 @@ function postResults(ip,el) postResults integer(pInt) :: & startPos, endPos ,& - of, instance + of, instance, homog postResults = 0.0_pReal @@ -954,10 +944,14 @@ function postResults(ip,el) chosenThermal: select case (thermal_type(mesh_element(3,el))) 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 - 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 startPos = endPos + 1_pInt diff --git a/src/kinematics_thermal_expansion.f90 b/src/kinematics_thermal_expansion.f90 index 3696593ad..56caa6e4b 100644 --- a/src/kinematics_thermal_expansion.f90 +++ b/src/kinematics_thermal_expansion.f90 @@ -4,34 +4,24 @@ !> @details to be done !-------------------------------------------------------------------------------------------------- module kinematics_thermal_expansion - use prec, only: & - pReal, & - pInt + use prec, only: & + pReal, & + pInt - 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, & - kinematics_thermal_expansion_LiAndItsTangent + implicit none + private + + type, private :: tParameters + real(pReal), allocatable, dimension(:,:,:) :: & + expansion + end type tParameters + + type(tParameters), dimension(:), allocatable :: param + + public :: & + kinematics_thermal_expansion_init, & + kinematics_thermal_expansion_initialStrain, & + kinematics_thermal_expansion_LiAndItsTangent contains @@ -40,197 +30,129 @@ contains !> @brief module initialization !> @details reads in material parameters, allocates arrays, and does sanity checks !-------------------------------------------------------------------------------------------------- -subroutine kinematics_thermal_expansion_init(fileUnit) -#if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800 - use, intrinsic :: iso_fortran_env, only: & - compiler_version, & - compiler_options -#endif - 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_kinematics, & - phase_Nkinematics, & - phase_Noutput, & - KINEMATICS_thermal_expansion_label, & - KINEMATICS_thermal_expansion_ID - use config, only: & - material_Nphase, & - MATERIAL_partPhase - - implicit none - integer(pInt), intent(in) :: fileUnit - - integer(pInt), allocatable, dimension(:) :: chunkPos - integer(pInt) :: maxNinstance,phase,instance,kinematics - character(len=65536) :: & - tag = '', & - line = '' - - write(6,'(/,a)') ' <<<+- kinematics_'//KINEMATICS_thermal_expansion_LABEL//' init -+>>>' - write(6,'(a15,a)') ' Current time: ',IO_timeStamp() -#include "compilation_info.f90" - - maxNinstance = int(count(phase_kinematics == KINEMATICS_thermal_expansion_ID),pInt) - if (maxNinstance == 0_pInt) return +subroutine kinematics_thermal_expansion_init() + use debug, only: & + debug_level,& + debug_constitutive,& + debug_levelBasic + use material, only: & + phase_kinematics, & + KINEMATICS_thermal_expansion_label, & + KINEMATICS_thermal_expansion_ID + use config, only: & + config_phase - if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt) & - write(6,'(a16,1x,i5,/)') '# instances:',maxNinstance + implicit none + integer(pInt) :: & + Ninstance, & + p, i + real(pReal), dimension(:), allocatable :: & + temp + + write(6,'(/,a)') ' <<<+- kinematics_'//KINEMATICS_thermal_expansion_LABEL//' init -+>>>' - 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 + Ninstance = int(count(phase_kinematics == KINEMATICS_thermal_expansion_ID),pInt) + + if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt) & + write(6,'(a16,1x,i5,/)') '# instances:',Ninstance + + allocate(param(Ninstance)) + + do p = 1_pInt, size(phase_kinematics) + if (all(phase_kinematics(:,p) /= KINEMATICS_thermal_expansion_ID)) cycle + + ! ToDo: Here we need to decide how to extend the concept of instances to + ! kinetics and sources. I would suggest that the same mechanism exists at maximum once per phase + + ! read up to three parameters (constant, linear, quadratic with T) + temp = config_phase(p)%getFloats('thermal_expansion11') + !lattice_thermalExpansion33(1,1,1:size(temp),p) = temp + temp = config_phase(p)%getFloats('thermal_expansion22', & + defaultVal=[(0.0_pReal, i=1,size(temp))],requiredSize=size(temp)) + !lattice_thermalExpansion33(2,2,1:size(temp),p) = temp + temp = config_phase(p)%getFloats('thermal_expansion33', & + defaultVal=[(0.0_pReal, i=1,size(temp))],requiredSize=size(temp)) + enddo end subroutine kinematics_thermal_expansion_init + !-------------------------------------------------------------------------------------------------- !> @brief report initial thermal strain based on current temperature deviation from reference !-------------------------------------------------------------------------------------------------- -pure function kinematics_thermal_expansion_initialStrain(ipc, ip, el) - use material, only: & - material_phase, & - material_homog, & - temperature, & - thermalMapping - use lattice, only: & - lattice_thermalExpansion33, & - lattice_referenceTemperature +pure function kinematics_thermal_expansion_initialStrain(homog,phase,offset) + use material, only: & + temperature + use lattice, only: & + lattice_thermalExpansion33, & + lattice_referenceTemperature + + implicit none + integer(pInt), intent(in) :: & + phase, & + homog, offset + real(pReal), dimension(3,3) :: & + kinematics_thermal_expansion_initialStrain !< initial thermal strain (should be small strain, though) - implicit none - integer(pInt), intent(in) :: & - ipc, & !< grain number - ip, & !< integration point number - el !< element number - real(pReal), dimension(3,3) :: & - 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 = & - (temperature(homog)%p(offset) - lattice_referenceTemperature(phase))**1 / 1. * & - lattice_thermalExpansion33(1:3,1:3,1,phase) + & ! constant coefficient - (temperature(homog)%p(offset) - lattice_referenceTemperature(phase))**2 / 2. * & - lattice_thermalExpansion33(1:3,1:3,2,phase) + & ! linear coefficient - (temperature(homog)%p(offset) - lattice_referenceTemperature(phase))**3 / 3. * & - lattice_thermalExpansion33(1:3,1:3,3,phase) ! quadratic coefficient + + kinematics_thermal_expansion_initialStrain = & + (temperature(homog)%p(offset) - lattice_referenceTemperature(phase))**1 / 1. * & + lattice_thermalExpansion33(1:3,1:3,1,phase) + & ! constant coefficient + (temperature(homog)%p(offset) - lattice_referenceTemperature(phase))**2 / 2. * & + lattice_thermalExpansion33(1:3,1:3,2,phase) + & ! linear coefficient + (temperature(homog)%p(offset) - lattice_referenceTemperature(phase))**3 / 3. * & + lattice_thermalExpansion33(1:3,1:3,3,phase) ! quadratic coefficient 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) - use material, only: & - material_phase, & - material_homog, & - temperature, & - temperatureRate, & - thermalMapping - use lattice, only: & - lattice_thermalExpansion33, & - lattice_referenceTemperature - - implicit none - integer(pInt), intent(in) :: & - ipc, & !< grain number - ip, & !< integration point number - el !< element number - 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) - integer(pInt) :: & - phase, & - homog, offset - real(pReal) :: & - T, TRef, TDot - - phase = material_phase(ipc,ip,el) - homog = material_homog(ip,el) - offset = thermalMapping(homog)%p(ip,el) - T = temperature(homog)%p(offset) - TDot = temperatureRate(homog)%p(offset) - TRef = lattice_referenceTemperature(phase) - - Li = TDot * ( & - 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,3,phase)*(T - TRef)**2 & ! quadratic coefficient - ) / & - (1.0_pReal & - + 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,3,phase)*(T - TRef)**3 / 3. & - ) - dLi_dTstar3333 = 0.0_pReal +subroutine kinematics_thermal_expansion_LiAndItsTangent(Li, dLi_dTstar, ipc, ip, el) + use material, only: & + material_phase, & + material_homog, & + temperature, & + temperatureRate, & + thermalMapping + use lattice, only: & + lattice_thermalExpansion33, & + lattice_referenceTemperature + + implicit none + integer(pInt), intent(in) :: & + ipc, & !< grain number + ip, & !< integration point number + el !< element number + real(pReal), intent(out), dimension(3,3) :: & + Li !< thermal velocity gradient + real(pReal), intent(out), dimension(3,3,3,3) :: & + dLi_dTstar !< derivative of Li with respect to Tstar (4th-order tensor defined to be zero) + integer(pInt) :: & + phase, & + homog, offset + real(pReal) :: & + T, TRef, TDot + + phase = material_phase(ipc,ip,el) + homog = material_homog(ip,el) + offset = thermalMapping(homog)%p(ip,el) + T = temperature(homog)%p(offset) + TDot = temperatureRate(homog)%p(offset) + TRef = lattice_referenceTemperature(phase) + + Li = TDot * ( & + 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,3,phase)*(T - TRef)**2 & ! quadratic coefficient + ) / & + (1.0_pReal & + + 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,3,phase)*(T - TRef)**3 / 3. & + ) + dLi_dTstar = 0.0_pReal end subroutine kinematics_thermal_expansion_LiAndItsTangent diff --git a/src/material.f90 b/src/material.f90 index 8a8f36a55..49ee38ee3 100644 --- a/src/material.f90 +++ b/src/material.f90 @@ -921,7 +921,7 @@ subroutine material_allocatePlasticState(phase,NofMyPhase,& sizeState,sizeDotState,sizeDeltaState,& Nslip,Ntwin,Ntrans) use numerics, only: & - numerics_integrator2 => numerics_integrator ! compatibility hack + numerics_integrator implicit none integer(pInt), intent(in) :: & @@ -933,8 +933,6 @@ subroutine material_allocatePlasticState(phase,NofMyPhase,& Nslip, & Ntwin, & Ntrans - integer(pInt) :: numerics_integrator ! compatibility hack - numerics_integrator = numerics_integrator2(1) ! compatibility hack plasticState(phase)%sizeState = sizeState plasticState(phase)%sizeDotState = sizeDotState @@ -971,7 +969,7 @@ end subroutine material_allocatePlasticState subroutine material_allocateSourceState(phase,of,NofMyPhase,& sizeState,sizeDotState,sizeDeltaState) use numerics, only: & - numerics_integrator2 => numerics_integrator ! compatibility hack + numerics_integrator implicit none integer(pInt), intent(in) :: & @@ -979,8 +977,6 @@ subroutine material_allocateSourceState(phase,of,NofMyPhase,& of, & NofMyPhase, & 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)%sizeDotState = sizeDotState @@ -1003,7 +999,7 @@ subroutine material_allocateSourceState(phase,of,NofMyPhase,& if (numerics_integrator == 5_pInt) & 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 diff --git a/src/numerics.f90 b/src/numerics.f90 index 1678d0c48..9727a04a7 100644 --- a/src/numerics.f90 +++ b/src/numerics.f90 @@ -23,12 +23,10 @@ module numerics 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 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 :: & 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 :: & 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) @@ -466,7 +464,7 @@ subroutine numerics_init 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)') ' 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,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_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 (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') 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') diff --git a/src/source_thermal_dissipation.f90 b/src/source_thermal_dissipation.f90 index 994d26b41..db37c8286 100644 --- a/src/source_thermal_dissipation.f90 +++ b/src/source_thermal_dissipation.f90 @@ -1,4 +1,5 @@ !-------------------------------------------------------------------------------------------------- +!> @author Martin Diehl, 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 !> @details to be done @@ -11,7 +12,6 @@ module source_thermal_dissipation implicit none private 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_instance !< instance of thermal dissipation source mechanism @@ -20,13 +20,19 @@ module source_thermal_dissipation character(len=64), dimension(:,:), allocatable, target, public :: & 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 :: & 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 @@ -38,30 +44,13 @@ contains !> @brief module initialization !> @details reads in material parameters, allocates arrays, and does sanity checks !-------------------------------------------------------------------------------------------------- -subroutine source_thermal_dissipation_init(fileUnit) -#if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800 - use, intrinsic :: iso_fortran_env, only: & - compiler_version, & - compiler_options -#endif +subroutine source_thermal_dissipation_init 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: & + material_allocateSourceState, & phase_source, & phase_Nsources, & phase_Noutput, & @@ -70,144 +59,73 @@ subroutine source_thermal_dissipation_init(fileUnit) material_phase, & sourceState use config, only: & + config_phase, & material_Nphase, & MATERIAL_partPhase - use numerics,only: & - numerics_integrator implicit none - integer(pInt), intent(in) :: fileUnit - - 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 = '' + integer(pInt) :: Ninstance,instance,source,sourceOffset + integer(pInt) :: NofMyPhase,p 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) & - 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) - do phase = 1, material_Nphase - source_thermal_dissipation_instance(phase) = count(phase_source(:,1:phase) == SOURCE_thermal_dissipation_ID) - do source = 1, phase_Nsources(phase) - if (phase_source(source,phase) == SOURCE_thermal_dissipation_ID) & - source_thermal_dissipation_offset(phase) = source + do p = 1, material_Nphase + source_thermal_dissipation_instance(p) = count(phase_source(:,1:p) == SOURCE_thermal_dissipation_ID) + do source = 1, phase_Nsources(p) + if (phase_source(source,p) == SOURCE_thermal_dissipation_ID) & + source_thermal_dissipation_offset(p) = source 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_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) - 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) + 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 + 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 - - 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 + !-------------------------------------------------------------------------------------------------- !> @brief returns local vacancy generation rate !-------------------------------------------------------------------------------------------------- -subroutine source_thermal_dissipation_getRateAndItsTangent(TDot, dTDOT_dT, Tstar_v, Lp, ipc, ip, el) - use math, only: & - math_Mandel6to33 - use material, only: & - phaseAt, phasememberAt +subroutine source_thermal_dissipation_getRateAndItsTangent(TDot, dTDOT_dT, Tstar, Lp, phase) implicit none integer(pInt), intent(in) :: & - ipc, & !< grain number - ip, & !< integration point number - el !< element number - real(pReal), intent(in), dimension(6) :: & - Tstar_v !< 2nd Piola Kirchhoff stress tensor (Mandel) + phase + real(pReal), intent(in), dimension(3,3) :: & + Tstar real(pReal), intent(in), dimension(3,3) :: & Lp real(pReal), intent(out) :: & 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)* & - sum(abs(math_Mandel6to33(Tstar_v)*Lp)) + TDot = source_thermal_dissipation_coldworkCoeff(instance)*sum(abs(Tstar*Lp)) dTDOT_dT = 0.0_pReal end subroutine source_thermal_dissipation_getRateAndItsTangent diff --git a/src/source_thermal_externalheat.f90 b/src/source_thermal_externalheat.f90 index b7151aece..2bf4cac9c 100644 --- a/src/source_thermal_externalheat.f90 +++ b/src/source_thermal_externalheat.f90 @@ -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 Philip Eisenlohr, Michigan State University !> @brief material subroutine for variable heat source -!> @details to be done !-------------------------------------------------------------------------------------------------- module source_thermal_externalheat - use prec, only: & - pReal, & - pInt + 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 + implicit none + private + integer(pInt), dimension(:), allocatable, public, protected :: & + 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 + 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 + 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, target, public :: & + source_thermal_externalheat_Noutput !< number of outputs per instance of this source - integer(pInt), dimension(:), allocatable, private :: & - source_thermal_externalheat_nIntervals + type, private :: tParameters !< container type for internal constitutive parameters + real(pReal), dimension(:), allocatable :: & + time, & + heat_rate + integer(pInt) :: & + nIntervals + end type tParameters - real(pReal), dimension(:,:), allocatable, private :: & - source_thermal_externalheat_time, & - source_thermal_externalheat_rate + type(tParameters), dimension(:), allocatable, private :: param !< containers of constitutive parameters (len Ninstance) - public :: & - source_thermal_externalheat_init, & - source_thermal_externalheat_dotState, & - source_thermal_externalheat_getRateAndItsTangent + + public :: & + source_thermal_externalheat_init, & + source_thermal_externalheat_dotState, & + source_thermal_externalheat_getRateAndItsTangent contains @@ -44,170 +47,77 @@ contains !> @brief module initialization !> @details reads in material parameters, allocates arrays, and does sanity checks !-------------------------------------------------------------------------------------------------- -subroutine source_thermal_externalheat_init(fileUnit) -#if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800 - use, intrinsic :: iso_fortran_env, only: & - compiler_version, & - compiler_options -#endif - 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_phase, & - sourceState - use config, only: & - material_Nphase, & - MATERIAL_partPhase - use numerics,only: & - numerics_integrator - - implicit none - integer(pInt), intent(in) :: fileUnit - - integer(pInt), allocatable, dimension(:) :: chunkPos - 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 - - write(6,'(/,a)') ' <<<+- source_'//SOURCE_thermal_externalheat_label//' init -+>>>' - write(6,'(a15,a)') ' Current time: ',IO_timeStamp() -#include "compilation_info.f90" +subroutine source_thermal_externalheat_init + use debug, only: & + debug_level,& + debug_constitutive,& + debug_levelBasic + use material, only: & + material_allocateSourceState, & + material_phase, & + phase_source, & + phase_Nsources, & + phase_Noutput, & + SOURCE_thermal_externalheat_label, & + SOURCE_thermal_externalheat_ID + use config, only: & + config_phase, & + material_Nphase, & + MATERIAL_partPhase - 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 + implicit none - 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 + real(pReal), allocatable, dimension(:) :: tempVar + integer(pInt) :: maxNinstance,instance,source,sourceOffset + integer(pInt) :: NofMyPhase,p - 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 - chunkPos = IO_stringPos(line) - 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 - 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 + write(6,'(/,a)') ' <<<+- source_'//SOURCE_thermal_externalheat_label//' init -+>>>' - 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)%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 + + 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 p = 1, material_Nphase + 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 + + 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) - enddo initializeInstances + allocate(param(maxNinstance)) + + 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) + + tempVar = config_phase(p)%getFloats('externalheat_time') + param(instance)%nIntervals = size(tempVar) - 1_pInt + + param(instance)%time= tempVar + + tempVar = config_phase(p)%getFloats('externalheat_rate',requiredSize = size(tempVar)) + param(instance)%heat_rate = tempVar + + call material_allocateSourceState(p,sourceOffset,NofMyPhase,1_pInt,1_pInt,0_pInt) + + enddo + end subroutine source_thermal_externalheat_init + !-------------------------------------------------------------------------------------------------- !> @brief rate of change of state !> @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 !-------------------------------------------------------------------------------------------------- -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) - 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) - & - source_thermal_externalheat_time(instance,interval)) / & - (source_thermal_externalheat_time(instance,interval+1) - & - source_thermal_externalheat_time(instance,interval)) ! fractional time within segment + param(instance)%time(interval)) / & + (param(instance)%time(interval+1) - & + param(instance)%time(interval)) ! fractional time within segment 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) ) & - TDot = source_thermal_externalheat_rate(instance,interval ) * (1.0_pReal - frac_time) + & - source_thermal_externalheat_rate(instance,interval+1) * frac_time ! interpolate heat rate between segment boundaries... + TDot = param(instance)%heat_rate(interval ) * (1.0_pReal - frac_time) + & + param(instance)%heat_rate(interval+1) * frac_time ! interpolate heat rate between segment boundaries... ! ...or extrapolate if outside of bounds enddo dTDot_dT = 0.0 diff --git a/src/thermal_adiabatic.f90 b/src/thermal_adiabatic.f90 index 6a70ca7ee..c3290bdfe 100644 --- a/src/thermal_adiabatic.f90 +++ b/src/thermal_adiabatic.f90 @@ -10,12 +10,9 @@ module thermal_adiabatic implicit none private - integer(pInt), dimension(:), allocatable, public, protected :: & - thermal_adiabatic_sizePostResults !< cumulative size of post results integer(pInt), dimension(:,:), allocatable, target, public :: & thermal_adiabatic_sizePostResult !< size of each post result output - character(len=64), dimension(:,:), allocatable, target, public :: & thermal_adiabatic_output !< name of each post result output @@ -45,27 +42,7 @@ contains !> @brief module initialization !> @details reads in material parameters, allocates arrays, and does sanity checks !-------------------------------------------------------------------------------------------------- -subroutine thermal_adiabatic_init(fileUnit) -#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 +subroutine thermal_adiabatic_init use material, only: & thermal_type, & thermal_typeInstance, & @@ -79,106 +56,61 @@ subroutine thermal_adiabatic_init(fileUnit) thermal_initialT, & temperature, & temperatureRate + use config, only: & + material_partHomogenization, & + config_homogenization implicit none - integer(pInt), intent(in) :: fileUnit - - integer(pInt), allocatable, dimension(:) :: chunkPos - integer(pInt) :: maxNinstance,mySize=0_pInt,section,instance,o + integer(pInt) :: maxNinstance,section,instance,i integer(pInt) :: sizeState integer(pInt) :: NofMyHomog - character(len=65536) :: & - tag = '', & - line = '' + character(len=65536), dimension(0), parameter :: emptyStringArray = [character(len=65536)::] + character(len=65536), dimension(:), allocatable :: outputs 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) 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_output (maxval(homogenization_Noutput),maxNinstance)) thermal_adiabatic_output = '' allocate(thermal_adiabatic_outputID (maxval(homogenization_Noutput),maxNinstance),source=undefined_ID) 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 - 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) - if (thermal_type(section) == THERMAL_adiabatic_ID) then - NofMyHomog=count(material_homog==section) - instance = thermal_typeInstance(section) - -!-------------------------------------------------------------------------------------------------- -! Determine size of postResults array - outputsLoop: do o = 1_pInt,thermal_adiabatic_Noutput(instance) - select case(thermal_adiabatic_outputID(o,instance)) - case(temperature_ID) - mySize = 1_pInt - end select - - if (mySize > 0_pInt) then ! any meaningful output found - thermal_adiabatic_sizePostResult(o,instance) = mySize - thermal_adiabatic_sizePostResults(instance) = thermal_adiabatic_sizePostResults(instance) + mySize - endif - enddo outputsLoop + if (thermal_type(section) /= THERMAL_adiabatic_ID) cycle + NofMyHomog=count(material_homog==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_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) = outputs(i) + thermal_adiabatic_sizePostResult(thermal_adiabatic_Noutput(instance),instance) = 1_pInt + end select + enddo ! allocate state arrays - sizeState = 1_pInt - thermalState(section)%sizeState = sizeState - thermalState(section)%sizePostResults = thermal_adiabatic_sizePostResults(instance) - allocate(thermalState(section)%state0 (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)) + sizeState = 1_pInt + thermalState(section)%sizeState = sizeState + thermalState(section)%sizePostResults = sum(thermal_adiabatic_sizePostResult(:,instance)) + allocate(thermalState(section)%state0 (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)) - nullify(thermalMapping(section)%p) - thermalMapping(section)%p => mappingHomogenization(1,:,:) - deallocate(temperature(section)%p) - temperature(section)%p => thermalState(section)%state(1,:) - deallocate(temperatureRate(section)%p) - allocate (temperatureRate(section)%p(NofMyHomog), source=0.0_pReal) + nullify(thermalMapping(section)%p) + thermalMapping(section)%p => mappingHomogenization(1,:,:) + deallocate(temperature(section)%p) + temperature(section)%p => thermalState(section)%state(1,:) + deallocate(temperatureRate(section)%p) + allocate (temperatureRate(section)%p(NofMyHomog), source=0.0_pReal) - endif - enddo initializeInstances + end subroutine thermal_adiabatic_init !-------------------------------------------------------------------------------------------------- @@ -233,11 +165,12 @@ end function thermal_adiabatic_updateState !-------------------------------------------------------------------------------------------------- subroutine thermal_adiabatic_getSourceAndItsTangent(Tdot, dTdot_dT, T, ip, el) use math, only: & - math_Mandel6to33 + math_6toSym33 use material, only: & homogenization_Ngrains, & mappingHomogenization, & phaseAt, & + phasememberAt, & thermal_typeInstance, & phase_Nsources, & phase_source, & @@ -264,30 +197,30 @@ subroutine thermal_adiabatic_getSourceAndItsTangent(Tdot, dTdot_dT, T, ip, el) integer(pInt) :: & phase, & homog, & - offset, & instance, & grain, & - source + source, & + constituent homog = mappingHomogenization(2,ip,el) - offset = mappingHomogenization(1,ip,el) instance = thermal_typeInstance(homog) Tdot = 0.0_pReal 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), & + math_6toSym33(crystallite_Tstar_v(1:6,grain,ip,el)), & crystallite_Lp(1:3,1:3,grain,ip,el), & - grain, ip, el) + phase) 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 @@ -311,12 +244,9 @@ function thermal_adiabatic_getSpecificHeat(ip,el) lattice_specificHeat use material, only: & homogenization_Ngrains, & - mappingHomogenization, & material_phase use mesh, only: & mesh_element - use crystallite, only: & - crystallite_push33ToRef implicit none integer(pInt), intent(in) :: & @@ -325,11 +255,10 @@ function thermal_adiabatic_getSpecificHeat(ip,el) real(pReal) :: & thermal_adiabatic_getSpecificHeat integer(pInt) :: & - homog, grain + grain thermal_adiabatic_getSpecificHeat = 0.0_pReal - homog = mappingHomogenization(2,ip,el) do grain = 1, homogenization_Ngrains(mesh_element(3,el)) thermal_adiabatic_getSpecificHeat = thermal_adiabatic_getSpecificHeat + & @@ -341,6 +270,7 @@ function thermal_adiabatic_getSpecificHeat(ip,el) end function thermal_adiabatic_getSpecificHeat + !-------------------------------------------------------------------------------------------------- !> @brief returns homogenized mass density !-------------------------------------------------------------------------------------------------- @@ -353,9 +283,7 @@ function thermal_adiabatic_getMassDensity(ip,el) material_phase use mesh, only: & mesh_element - use crystallite, only: & - crystallite_push33ToRef - + implicit none integer(pInt), intent(in) :: & ip, & !< integration point number @@ -363,11 +291,10 @@ function thermal_adiabatic_getMassDensity(ip,el) real(pReal) :: & thermal_adiabatic_getMassDensity integer(pInt) :: & - homog, grain + grain thermal_adiabatic_getMassDensity = 0.0_pReal - - homog = mappingHomogenization(2,ip,el) + do grain = 1, homogenization_Ngrains(mesh_element(3,el)) thermal_adiabatic_getMassDensity = thermal_adiabatic_getMassDensity + & @@ -378,42 +305,38 @@ function thermal_adiabatic_getMassDensity(ip,el) thermal_adiabatic_getMassDensity/real(homogenization_Ngrains(mesh_element(3,el)),pReal) end function thermal_adiabatic_getMassDensity - + + !-------------------------------------------------------------------------------------------------- !> @brief return array of thermal results !-------------------------------------------------------------------------------------------------- -function thermal_adiabatic_postResults(ip,el) +function thermal_adiabatic_postResults(homog,instance,of) result(postResults) use material, only: & - mappingHomogenization, & - thermal_typeInstance, & - thermalMapping, & temperature implicit none - integer(pInt), intent(in) :: & - ip, & !< integration point - el !< element - real(pReal), dimension(thermal_adiabatic_sizePostResults(thermal_typeInstance(mappingHomogenization(2,ip,el)))) :: & - thermal_adiabatic_postResults + integer(pInt), intent(in) :: & + homog, & + instance, & + of + + real(pReal), dimension(sum(thermal_adiabatic_sizePostResult(:,instance))) :: & + postResults integer(pInt) :: & - instance, homog, offset, o, c - - homog = mappingHomogenization(2,ip,el) - offset = thermalMapping(homog)%p(ip,el) - instance = thermal_typeInstance(homog) + o, c c = 0_pInt - thermal_adiabatic_postResults = 0.0_pReal do o = 1_pInt,thermal_adiabatic_Noutput(instance) select case(thermal_adiabatic_outputID(o,instance)) case (temperature_ID) - thermal_adiabatic_postResults(c+1_pInt) = temperature(homog)%p(offset) + postResults(c+1_pInt) = temperature(homog)%p(of) c = c + 1 end select enddo + end function thermal_adiabatic_postResults end module thermal_adiabatic diff --git a/src/thermal_conduction.f90 b/src/thermal_conduction.f90 index 16497040b..88da0529b 100644 --- a/src/thermal_conduction.f90 +++ b/src/thermal_conduction.f90 @@ -10,12 +10,9 @@ module thermal_conduction implicit none private - integer(pInt), dimension(:), allocatable, public, protected :: & - thermal_conduction_sizePostResults !< cumulative size of post results integer(pInt), dimension(:,:), allocatable, target, public :: & thermal_conduction_sizePostResult !< size of each post result output - character(len=64), dimension(:,:), allocatable, target, public :: & thermal_conduction_output !< name of each post result output @@ -46,25 +43,7 @@ contains !> @brief module initialization !> @details reads in material parameters, allocates arrays, and does sanity checks !-------------------------------------------------------------------------------------------------- -subroutine thermal_conduction_init(fileUnit) -#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 +subroutine thermal_conduction_init use material, only: & thermal_type, & thermal_typeInstance, & @@ -79,107 +58,61 @@ subroutine thermal_conduction_init(fileUnit) temperature, & temperatureRate use config, only: & - material_partHomogenization + material_partHomogenization, & + config_homogenization implicit none - integer(pInt), intent(in) :: fileUnit - - integer(pInt), allocatable, dimension(:) :: chunkPos - integer(pInt) :: maxNinstance,mySize=0_pInt,section,instance,o + integer(pInt) :: maxNinstance,section,instance,i integer(pInt) :: sizeState integer(pInt) :: NofMyHomog - character(len=65536) :: & - tag = '', & - line = '' + character(len=65536), dimension(0), parameter :: emptyStringArray = [character(len=65536)::] + character(len=65536), dimension(:), allocatable :: outputs 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) 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_output (maxval(homogenization_Noutput),maxNinstance)) thermal_conduction_output = '' allocate(thermal_conduction_outputID (maxval(homogenization_Noutput),maxNinstance),source=undefined_ID) 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 - 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) - if (thermal_type(section) == THERMAL_conduction_ID) then - NofMyHomog=count(material_homog==section) - instance = thermal_typeInstance(section) + if (thermal_type(section) /= THERMAL_conduction_ID) cycle + NofMyHomog=count(material_homog==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 - sizeState = 0_pInt - thermalState(section)%sizeState = sizeState - thermalState(section)%sizePostResults = thermal_conduction_sizePostResults(instance) - allocate(thermalState(section)%state0 (sizeState,NofMyHomog)) - allocate(thermalState(section)%subState0(sizeState,NofMyHomog)) - allocate(thermalState(section)%state (sizeState,NofMyHomog)) + sizeState = 0_pInt + thermalState(section)%sizeState = sizeState + thermalState(section)%sizePostResults = sum(thermal_conduction_sizePostResult(:,instance)) + allocate(thermalState(section)%state0 (sizeState,NofMyHomog)) + allocate(thermalState(section)%subState0(sizeState,NofMyHomog)) + allocate(thermalState(section)%state (sizeState,NofMyHomog)) - nullify(thermalMapping(section)%p) - thermalMapping(section)%p => mappingHomogenization(1,:,:) - deallocate(temperature (section)%p) - allocate (temperature (section)%p(NofMyHomog), source=thermal_initialT(section)) - deallocate(temperatureRate(section)%p) - allocate (temperatureRate(section)%p(NofMyHomog), source=0.0_pReal) + nullify(thermalMapping(section)%p) + thermalMapping(section)%p => mappingHomogenization(1,:,:) + deallocate(temperature (section)%p) + allocate (temperature (section)%p(NofMyHomog), source=thermal_initialT(section)) + deallocate(temperatureRate(section)%p) + allocate (temperatureRate(section)%p(NofMyHomog), source=0.0_pReal) - endif - enddo initializeInstances + end subroutine thermal_conduction_init !-------------------------------------------------------------------------------------------------- @@ -187,11 +120,12 @@ end subroutine thermal_conduction_init !-------------------------------------------------------------------------------------------------- subroutine thermal_conduction_getSourceAndItsTangent(Tdot, dTdot_dT, T, ip, el) use math, only: & - math_Mandel6to33 + math_6toSym33 use material, only: & homogenization_Ngrains, & mappingHomogenization, & phaseAt, & + phasememberAt, & thermal_typeInstance, & phase_Nsources, & phase_source, & @@ -221,7 +155,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 +166,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), & + math_6toSym33(crystallite_Tstar_v(1:6,grain,ip,el)), & crystallite_Lp(1:3,1:3,grain,ip,el), & - grain, ip, el) + phase) 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 @@ -258,6 +194,7 @@ subroutine thermal_conduction_getSourceAndItsTangent(Tdot, dTdot_dT, T, ip, el) end subroutine thermal_conduction_getSourceAndItsTangent + !-------------------------------------------------------------------------------------------------- !> @brief returns homogenized thermal conductivity in reference configuration !-------------------------------------------------------------------------------------------------- @@ -266,7 +203,6 @@ function thermal_conduction_getConductivity33(ip,el) lattice_thermalConductivity33 use material, only: & homogenization_Ngrains, & - mappingHomogenization, & material_phase use mesh, only: & mesh_element @@ -280,10 +216,8 @@ function thermal_conduction_getConductivity33(ip,el) real(pReal), dimension(3,3) :: & thermal_conduction_getConductivity33 integer(pInt) :: & - homog, & grain - homog = mappingHomogenization(2,ip,el) thermal_conduction_getConductivity33 = 0.0_pReal do grain = 1, homogenization_Ngrains(mesh_element(3,el)) @@ -295,7 +229,8 @@ function thermal_conduction_getConductivity33(ip,el) thermal_conduction_getConductivity33/real(homogenization_Ngrains(mesh_element(3,el)),pReal) end function thermal_conduction_getConductivity33 - + + !-------------------------------------------------------------------------------------------------- !> @brief returns homogenized specific heat capacity !-------------------------------------------------------------------------------------------------- @@ -304,12 +239,9 @@ function thermal_conduction_getSpecificHeat(ip,el) lattice_specificHeat use material, only: & homogenization_Ngrains, & - mappingHomogenization, & material_phase use mesh, only: & mesh_element - use crystallite, only: & - crystallite_push33ToRef implicit none integer(pInt), intent(in) :: & @@ -318,11 +250,10 @@ function thermal_conduction_getSpecificHeat(ip,el) real(pReal) :: & thermal_conduction_getSpecificHeat integer(pInt) :: & - homog, grain + grain thermal_conduction_getSpecificHeat = 0.0_pReal - homog = mappingHomogenization(2,ip,el) do grain = 1, homogenization_Ngrains(mesh_element(3,el)) thermal_conduction_getSpecificHeat = thermal_conduction_getSpecificHeat + & @@ -342,12 +273,9 @@ function thermal_conduction_getMassDensity(ip,el) lattice_massDensity use material, only: & homogenization_Ngrains, & - mappingHomogenization, & material_phase use mesh, only: & mesh_element - use crystallite, only: & - crystallite_push33ToRef implicit none integer(pInt), intent(in) :: & @@ -356,22 +284,22 @@ function thermal_conduction_getMassDensity(ip,el) real(pReal) :: & thermal_conduction_getMassDensity integer(pInt) :: & - homog, grain + grain thermal_conduction_getMassDensity = 0.0_pReal - 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 = & thermal_conduction_getMassDensity/real(homogenization_Ngrains(mesh_element(3,el)),pReal) end function thermal_conduction_getMassDensity - + + !-------------------------------------------------------------------------------------------------- !> @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 + !-------------------------------------------------------------------------------------------------- !> @brief return array of thermal results !-------------------------------------------------------------------------------------------------- -function thermal_conduction_postResults(ip,el) +function thermal_conduction_postResults(homog,instance,of) result(postResults) use material, only: & - mappingHomogenization, & - thermal_typeInstance, & - temperature, & - thermalMapping + temperature implicit none integer(pInt), intent(in) :: & - ip, & !< integration point - el !< element - real(pReal), dimension(thermal_conduction_sizePostResults(thermal_typeInstance(mappingHomogenization(2,ip,el)))) :: & - thermal_conduction_postResults + homog, & + instance, & + of + + real(pReal), dimension(sum(thermal_conduction_sizePostResult(:,instance))) :: & + postResults integer(pInt) :: & - instance, homog, offset, o, c - - homog = mappingHomogenization(2,ip,el) - offset = thermalMapping(homog)%p(ip,el) - instance = thermal_typeInstance(homog) + o, c c = 0_pInt - thermal_conduction_postResults = 0.0_pReal do o = 1_pInt,thermal_conduction_Noutput(instance) select case(thermal_conduction_outputID(o,instance)) case (temperature_ID) - thermal_conduction_postResults(c+1_pInt) = temperature(homog)%p(offset) + postResults(c+1_pInt) = temperature(homog)%p(of) c = c + 1 end select enddo + end function thermal_conduction_postResults end module thermal_conduction diff --git a/src/thermal_isothermal.f90 b/src/thermal_isothermal.f90 index fb518fe24..7485cd34f 100644 --- a/src/thermal_isothermal.f90 +++ b/src/thermal_isothermal.f90 @@ -26,14 +26,14 @@ subroutine thermal_isothermal_init() pInt use IO, only: & IO_timeStamp + use config, only: & + material_Nhomogenization use material - use config implicit none integer(pInt) :: & homog, & - NofMyHomog, & - sizeState + NofMyHomog write(6,'(/,a)') ' <<<+- thermal_'//THERMAL_isothermal_label//' init -+>>>' write(6,'(a15,a)') ' Current time: ',IO_timeStamp() @@ -41,21 +41,19 @@ subroutine thermal_isothermal_init() initializeInstances: do homog = 1_pInt, material_Nhomogenization - myhomog: if (thermal_type(homog) == THERMAL_isothermal_ID) then - NofMyHomog = count(material_homog == homog) - sizeState = 0_pInt - thermalState(homog)%sizeState = sizeState - thermalState(homog)%sizePostResults = sizeState - allocate(thermalState(homog)%state0 (sizeState,NofMyHomog), source=0.0_pReal) - allocate(thermalState(homog)%subState0(sizeState,NofMyHomog), source=0.0_pReal) - allocate(thermalState(homog)%state (sizeState,NofMyHomog), source=0.0_pReal) + if (thermal_type(homog) /= THERMAL_isothermal_ID) cycle + NofMyHomog = count(material_homog == homog) + thermalState(homog)%sizeState = 0_pInt + thermalState(homog)%sizePostResults = 0_pInt + allocate(thermalState(homog)%state0 (0_pInt,NofMyHomog), source=0.0_pReal) + allocate(thermalState(homog)%subState0(0_pInt,NofMyHomog), source=0.0_pReal) + allocate(thermalState(homog)%state (0_pInt,NofMyHomog), source=0.0_pReal) - deallocate(temperature (homog)%p) - allocate (temperature (homog)%p(1), source=thermal_initialT(homog)) - deallocate(temperatureRate(homog)%p) - allocate (temperatureRate(homog)%p(1), source=0.0_pReal) + deallocate(temperature (homog)%p) + allocate (temperature (homog)%p(1), source=thermal_initialT(homog)) + deallocate(temperatureRate(homog)%p) + allocate (temperatureRate(homog)%p(1), source=0.0_pReal) - endif myhomog enddo initializeInstances