From e0fa3e0b2609c2429d5656306683b7caba1eabdf Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 30 Dec 2018 20:58:38 +0100 Subject: [PATCH 01/10] takeover from 40_XX and 41_XX branch easier to focus on thermal instead of doing all kinematics and sources together --- src/constitutive.f90 | 2 +- src/kinematics_thermal_expansion.f90 | 105 +++------------------------ src/source_thermal_dissipation.f90 | 48 +++++++----- src/source_thermal_externalheat.f90 | 33 ++++++--- src/thermal_adiabatic.f90 | 9 ++- src/thermal_conduction.f90 | 13 ++-- 6 files changed, 76 insertions(+), 134 deletions(-) diff --git a/src/constitutive.f90 b/src/constitutive.f90 index 8294047e7..1c32f73d9 100644 --- a/src/constitutive.f90 +++ b/src/constitutive.f90 @@ -188,7 +188,7 @@ subroutine constitutive_init() call IO_checkAndRewind(FILEUNIT) if (any(phase_kinematics == KINEMATICS_cleavage_opening_ID)) call kinematics_cleavage_opening_init(FILEUNIT) if (any(phase_kinematics == KINEMATICS_slipplane_opening_ID)) call kinematics_slipplane_opening_init(FILEUNIT) - if (any(phase_kinematics == KINEMATICS_thermal_expansion_ID)) call kinematics_thermal_expansion_init(FILEUNIT) + if (any(phase_kinematics == KINEMATICS_thermal_expansion_ID)) call kinematics_thermal_expansion_init if (any(phase_kinematics == KINEMATICS_vacancy_strain_ID)) call kinematics_vacancy_strain_init(FILEUNIT) if (any(phase_kinematics == KINEMATICS_hydrogen_strain_ID)) call kinematics_hydrogen_strain_init(FILEUNIT) close(FILEUNIT) diff --git a/src/kinematics_thermal_expansion.f90 b/src/kinematics_thermal_expansion.f90 index 3d1de3d0a..e8f0d71c7 100644 --- a/src/kinematics_thermal_expansion.f90 +++ b/src/kinematics_thermal_expansion.f90 @@ -10,24 +10,7 @@ module kinematics_thermal_expansion implicit none private - integer(pInt), dimension(:), allocatable, public, protected :: & - kinematics_thermal_expansion_sizePostResults, & !< cumulative size of post results - kinematics_thermal_expansion_offset, & !< which kinematics is my current damage mechanism? - kinematics_thermal_expansion_instance !< instance of damage kinematics mechanism - integer(pInt), dimension(:,:), allocatable, target, public :: & - kinematics_thermal_expansion_sizePostResult !< size of each post result output - - character(len=64), dimension(:,:), allocatable, target, public :: & - kinematics_thermal_expansion_output !< name of each post result output - - integer(pInt), dimension(:), allocatable, target, public :: & - kinematics_thermal_expansion_Noutput !< number of outputs per instance of this damage - -! enum, bind(c) ! ToDo kinematics need state machinery to deal with sizePostResult -! enumerator :: undefined_ID, & ! possible remedy is to decouple having state vars from having output -! thermalexpansionrate_ID ! which means to separate user-defined types tState + tOutput... -! end enum public :: & kinematics_thermal_expansion_init, & kinematics_thermal_expansion_initialStrain, & @@ -40,7 +23,7 @@ contains !> @brief module initialization !> @details reads in material parameters, allocates arrays, and does sanity checks !-------------------------------------------------------------------------------------------------- -subroutine kinematics_thermal_expansion_init(fileUnit) +subroutine kinematics_thermal_expansion_init() #if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800 use, intrinsic :: iso_fortran_env, only: & compiler_version, & @@ -51,37 +34,17 @@ subroutine kinematics_thermal_expansion_init(fileUnit) debug_constitutive,& debug_levelBasic use IO, only: & - IO_read, & - IO_lc, & - IO_getTag, & - IO_isBlank, & - IO_stringPos, & - IO_stringValue, & - IO_floatValue, & - IO_intValue, & - IO_warning, & - IO_error, & - IO_timeStamp, & - IO_EOF + IO_timeStamp use material, only: & phase_kinematics, & - phase_Nkinematics, & - phase_Noutput, & KINEMATICS_thermal_expansion_label, & KINEMATICS_thermal_expansion_ID use config, only: & - material_Nphase, & - MATERIAL_partPhase + config_phase implicit none - integer(pInt), intent(in) :: fileUnit - - integer(pInt), allocatable, dimension(:) :: chunkPos - integer(pInt) :: maxNinstance,phase,instance,kinematics - character(len=65536) :: & - tag = '', & - line = '' - + integer(pInt) maxNinstance + write(6,'(/,a)') ' <<<+- kinematics_'//KINEMATICS_thermal_expansion_LABEL//' init -+>>>' write(6,'(a15,a)') ' Current time: ',IO_timeStamp() #include "compilation_info.f90" @@ -91,58 +54,8 @@ subroutine kinematics_thermal_expansion_init(fileUnit) if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt) & write(6,'(a16,1x,i5,/)') '# instances:',maxNinstance - - allocate(kinematics_thermal_expansion_offset(material_Nphase), source=0_pInt) - allocate(kinematics_thermal_expansion_instance(material_Nphase), source=0_pInt) - do phase = 1, material_Nphase - kinematics_thermal_expansion_instance(phase) = count(phase_kinematics(:,1:phase) == kinematics_thermal_expansion_ID) - do kinematics = 1, phase_Nkinematics(phase) - if (phase_kinematics(kinematics,phase) == kinematics_thermal_expansion_ID) & - kinematics_thermal_expansion_offset(phase) = kinematics - enddo - enddo - - allocate(kinematics_thermal_expansion_sizePostResults(maxNinstance), source=0_pInt) - allocate(kinematics_thermal_expansion_sizePostResult(maxval(phase_Noutput),maxNinstance),source=0_pInt) - allocate(kinematics_thermal_expansion_output(maxval(phase_Noutput),maxNinstance)) - kinematics_thermal_expansion_output = '' - allocate(kinematics_thermal_expansion_Noutput(maxNinstance), source=0_pInt) - rewind(fileUnit) - phase = 0_pInt - do while (trim(line) /= IO_EOF .and. IO_lc(IO_getTag(line,'<','>')) /= MATERIAL_partPhase) ! wind forward to - line = IO_read(fileUnit) - enddo - - parsingFile: do while (trim(line) /= IO_EOF) ! read through sections of phase part - line = IO_read(fileUnit) - if (IO_isBlank(line)) cycle ! skip empty lines - if (IO_getTag(line,'<','>') /= '') then ! stop at next part - line = IO_read(fileUnit, .true.) ! reset IO_read - exit - endif - if (IO_getTag(line,'[',']') /= '') then ! next phase section - phase = phase + 1_pInt ! advance phase section counter - cycle ! skip to next line - endif - if (phase > 0_pInt ) then; if (any(phase_kinematics(:,phase) == KINEMATICS_thermal_expansion_ID)) then ! do not short-circuit here (.and. with next if statemen). It's not safe in Fortran - instance = kinematics_thermal_expansion_instance(phase) ! which instance of my damage is present phase - chunkPos = IO_stringPos(line) - tag = IO_lc(IO_stringValue(line,chunkPos,1_pInt)) ! extract key... - select case(tag) -! case ('(output)') -! output = IO_lc(IO_stringValue(line,chunkPos,2_pInt)) ! ...and corresponding output -! select case(output) -! case ('thermalexpansionrate') -! kinematics_thermal_expansion_Noutput(instance) = kinematics_thermal_expansion_Noutput(instance) + 1_pInt -! kinematics_thermal_expansion_outputID(kinematics_thermal_expansion_Noutput(instance),instance) = & -! thermalexpansionrate_ID -! kinematics_thermal_expansion_output(kinematics_thermal_expansion_Noutput(instance),instance) = output -! ToDo add sizePostResult loop afterwards... - - end select - endif; endif - enddo parsingFile +! ToDo: this subroutine should read in lattice_thermal_expansion. No need to make it a global array end subroutine kinematics_thermal_expansion_init @@ -187,7 +100,7 @@ end function kinematics_thermal_expansion_initialStrain !-------------------------------------------------------------------------------------------------- !> @brief contains the constitutive equation for calculating the velocity gradient !-------------------------------------------------------------------------------------------------- -subroutine kinematics_thermal_expansion_LiAndItsTangent(Li, dLi_dTstar3333, ipc, ip, el) +subroutine kinematics_thermal_expansion_LiAndItsTangent(Li, dLi_dTstar, ipc, ip, el) use material, only: & material_phase, & material_homog, & @@ -206,7 +119,7 @@ subroutine kinematics_thermal_expansion_LiAndItsTangent(Li, dLi_dTstar3333, ipc, real(pReal), intent(out), dimension(3,3) :: & Li !< thermal velocity gradient real(pReal), intent(out), dimension(3,3,3,3) :: & - dLi_dTstar3333 !< derivative of Li with respect to Tstar (4th-order tensor defined to be zero) + dLi_dTstar !< derivative of Li with respect to Tstar (4th-order tensor defined to be zero) integer(pInt) :: & phase, & homog, offset @@ -230,7 +143,7 @@ subroutine kinematics_thermal_expansion_LiAndItsTangent(Li, dLi_dTstar3333, ipc, + lattice_thermalExpansion33(1:3,1:3,2,phase)*(T - TRef)**2 / 2. & + lattice_thermalExpansion33(1:3,1:3,3,phase)*(T - TRef)**3 / 3. & ) - dLi_dTstar3333 = 0.0_pReal + dLi_dTstar = 0.0_pReal end subroutine kinematics_thermal_expansion_LiAndItsTangent diff --git a/src/source_thermal_dissipation.f90 b/src/source_thermal_dissipation.f90 index 994d26b41..290ad7efe 100644 --- a/src/source_thermal_dissipation.f90 +++ b/src/source_thermal_dissipation.f90 @@ -27,6 +27,15 @@ module source_thermal_dissipation real(pReal), dimension(:), allocatable, private :: & source_thermal_dissipation_coldworkCoeff + + type, private :: tParameters !< container type for internal constitutive parameters + real(pReal) :: & + coldworkCoeff + end type tParameters + + type(tParameters), dimension(:), allocatable, private :: param !< containers of constitutive parameters (len Ninstance) + + public :: & source_thermal_dissipation_init, & source_thermal_dissipation_getRateAndItsTangent @@ -70,6 +79,7 @@ subroutine source_thermal_dissipation_init(fileUnit) material_phase, & sourceState use config, only: & + config_phase, & material_Nphase, & MATERIAL_partPhase use numerics,only: & @@ -79,9 +89,9 @@ subroutine source_thermal_dissipation_init(fileUnit) integer(pInt), intent(in) :: fileUnit integer(pInt), allocatable, dimension(:) :: chunkPos - integer(pInt) :: maxNinstance,phase,instance,source,sourceOffset + integer(pInt) :: Ninstance,phase,instance,source,sourceOffset integer(pInt) :: sizeState, sizeDotState, sizeDeltaState - integer(pInt) :: NofMyPhase + integer(pInt) :: NofMyPhase,p character(len=65536) :: & tag = '', & line = '' @@ -90,10 +100,10 @@ subroutine source_thermal_dissipation_init(fileUnit) write(6,'(a15,a)') ' Current time: ',IO_timeStamp() #include "compilation_info.f90" - maxNinstance = int(count(phase_source == SOURCE_thermal_dissipation_ID),pInt) - if (maxNinstance == 0_pInt) return + Ninstance = int(count(phase_source == SOURCE_thermal_dissipation_ID),pInt) + if (Ninstance == 0_pInt) return if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt) & - write(6,'(a16,1x,i5,/)') '# instances:',maxNinstance + write(6,'(a16,1x,i5,/)') '# instances:',Ninstance allocate(source_thermal_dissipation_offset(material_Nphase), source=0_pInt) allocate(source_thermal_dissipation_instance(material_Nphase), source=0_pInt) @@ -105,12 +115,17 @@ subroutine source_thermal_dissipation_init(fileUnit) enddo enddo - allocate(source_thermal_dissipation_sizePostResults(maxNinstance), source=0_pInt) - allocate(source_thermal_dissipation_sizePostResult(maxval(phase_Noutput),maxNinstance),source=0_pInt) - allocate(source_thermal_dissipation_output (maxval(phase_Noutput),maxNinstance)) + allocate(source_thermal_dissipation_sizePostResults(Ninstance), source=0_pInt) + allocate(source_thermal_dissipation_sizePostResult(maxval(phase_Noutput),Ninstance),source=0_pInt) + allocate(source_thermal_dissipation_output (maxval(phase_Noutput),Ninstance)) source_thermal_dissipation_output = '' - allocate(source_thermal_dissipation_Noutput(maxNinstance), source=0_pInt) - allocate(source_thermal_dissipation_coldworkCoeff(maxNinstance), source=0.0_pReal) + allocate(source_thermal_dissipation_Noutput(Ninstance), source=0_pInt) + + allocate(source_thermal_dissipation_coldworkCoeff(Ninstance), source=0.0_pReal) + + do p=1, size(config_phase) + if (all(phase_source(:,p) /= SOURCE_thermal_dissipation_ID)) cycle + enddo rewind(fileUnit) phase = 0_pInt @@ -181,17 +196,14 @@ end subroutine source_thermal_dissipation_init !-------------------------------------------------------------------------------------------------- !> @brief returns local vacancy generation rate !-------------------------------------------------------------------------------------------------- -subroutine source_thermal_dissipation_getRateAndItsTangent(TDot, dTDOT_dT, Tstar_v, Lp, ipc, ip, el) +subroutine source_thermal_dissipation_getRateAndItsTangent(TDot, dTDOT_dT, Tstar_v, Lp, phase, constituent) use math, only: & math_Mandel6to33 - use material, only: & - phaseAt, phasememberAt implicit none integer(pInt), intent(in) :: & - ipc, & !< grain number - ip, & !< integration point number - el !< element number + phase, & + constituent real(pReal), intent(in), dimension(6) :: & Tstar_v !< 2nd Piola Kirchhoff stress tensor (Mandel) real(pReal), intent(in), dimension(3,3) :: & @@ -200,10 +212,8 @@ subroutine source_thermal_dissipation_getRateAndItsTangent(TDot, dTDOT_dT, Tstar TDot, & dTDOT_dT integer(pInt) :: & - instance, phase, constituent + instance - phase = phaseAt(ipc,ip,el) - constituent = phasememberAt(ipc,ip,el) instance = source_thermal_dissipation_instance(phase) TDot = source_thermal_dissipation_coldworkCoeff(instance)* & diff --git a/src/source_thermal_externalheat.f90 b/src/source_thermal_externalheat.f90 index b7151aece..eac1232f3 100644 --- a/src/source_thermal_externalheat.f90 +++ b/src/source_thermal_externalheat.f90 @@ -32,6 +32,18 @@ module source_thermal_externalheat source_thermal_externalheat_time, & source_thermal_externalheat_rate + + type, private :: tParameters !< container type for internal constitutive parameters + real(pReal), dimension(:), allocatable :: & + time, & + rate + integer(pInt) :: & + nInterval + end type tParameters + + type(tParameters), dimension(:), allocatable, private :: param !< containers of constitutive parameters (len Ninstance) + + public :: & source_thermal_externalheat_init, & source_thermal_externalheat_dotState, & @@ -76,6 +88,7 @@ subroutine source_thermal_externalheat_init(fileUnit) material_phase, & sourceState use config, only: & + config_phase, & material_Nphase, & MATERIAL_partPhase use numerics,only: & @@ -87,7 +100,7 @@ subroutine source_thermal_externalheat_init(fileUnit) integer(pInt), allocatable, dimension(:) :: chunkPos integer(pInt) :: maxNinstance,phase,instance,source,sourceOffset integer(pInt) :: sizeState, sizeDotState, sizeDeltaState - integer(pInt) :: NofMyPhase,interval + integer(pInt) :: NofMyPhase,interval,p character(len=65536) :: & tag = '', & line = '' @@ -117,11 +130,15 @@ subroutine source_thermal_externalheat_init(fileUnit) allocate(source_thermal_externalheat_output (maxval(phase_Noutput),maxNinstance)) source_thermal_externalheat_output = '' allocate(source_thermal_externalheat_Noutput(maxNinstance), source=0_pInt) - allocate(source_thermal_externalheat_nIntervals(maxNinstance), source=0_pInt) + allocate(source_thermal_externalheat_nIntervals(maxNinstance), source=0_pInt) allocate(temp_time(maxNinstance,1000), source=0.0_pReal) allocate(temp_rate(maxNinstance,1000), source=0.0_pReal) + do p=1, size(config_phase) + if (all(phase_source(:,p) /= SOURCE_thermal_externalheat_ID)) cycle + enddo + rewind(fileUnit) phase = 0_pInt do while (trim(line) /= IO_EOF .and. IO_lc(IO_getTag(line,'<','>')) /= MATERIAL_partPhase) ! wind forward to @@ -238,26 +255,22 @@ end subroutine source_thermal_externalheat_dotState !-------------------------------------------------------------------------------------------------- !> @brief returns local heat generation rate !-------------------------------------------------------------------------------------------------- -subroutine source_thermal_externalheat_getRateAndItsTangent(TDot, dTDot_dT, ipc, ip, el) +subroutine source_thermal_externalheat_getRateAndItsTangent(TDot, dTDot_dT, phase, constituent) use material, only: & - phaseAt, phasememberAt, & sourceState implicit none integer(pInt), intent(in) :: & - ipc, & !< grain number - ip, & !< integration point number - el !< element number + phase, & + constituent real(pReal), intent(out) :: & TDot, & dTDot_dT integer(pInt) :: & - instance, phase, constituent, sourceOffset, interval + instance, sourceOffset, interval real(pReal) :: & frac_time - phase = phaseAt(ipc,ip,el) - constituent = phasememberAt(ipc,ip,el) instance = source_thermal_externalheat_instance(phase) sourceOffset = source_thermal_externalheat_offset(phase) diff --git a/src/thermal_adiabatic.f90 b/src/thermal_adiabatic.f90 index 6a70ca7ee..dd52293eb 100644 --- a/src/thermal_adiabatic.f90 +++ b/src/thermal_adiabatic.f90 @@ -238,6 +238,7 @@ subroutine thermal_adiabatic_getSourceAndItsTangent(Tdot, dTdot_dT, T, ip, el) homogenization_Ngrains, & mappingHomogenization, & phaseAt, & + phasememberAt, & thermal_typeInstance, & phase_Nsources, & phase_source, & @@ -267,7 +268,8 @@ subroutine thermal_adiabatic_getSourceAndItsTangent(Tdot, dTdot_dT, T, ip, el) offset, & instance, & grain, & - source + source, & + constituent homog = mappingHomogenization(2,ip,el) offset = mappingHomogenization(1,ip,el) @@ -277,17 +279,18 @@ subroutine thermal_adiabatic_getSourceAndItsTangent(Tdot, dTdot_dT, T, ip, el) dTdot_dT = 0.0_pReal do grain = 1, homogenization_Ngrains(homog) phase = phaseAt(grain,ip,el) + constituent = phasememberAt(grain,ip,el) do source = 1, phase_Nsources(phase) select case(phase_source(source,phase)) case (SOURCE_thermal_dissipation_ID) call source_thermal_dissipation_getRateAndItsTangent(my_Tdot, my_dTdot_dT, & crystallite_Tstar_v(1:6,grain,ip,el), & crystallite_Lp(1:3,1:3,grain,ip,el), & - grain, ip, el) + phase, constituent) case (SOURCE_thermal_externalheat_ID) call source_thermal_externalheat_getRateAndItsTangent(my_Tdot, my_dTdot_dT, & - grain, ip, el) + phase, constituent) case default my_Tdot = 0.0_pReal diff --git a/src/thermal_conduction.f90 b/src/thermal_conduction.f90 index 16497040b..1a9014e89 100644 --- a/src/thermal_conduction.f90 +++ b/src/thermal_conduction.f90 @@ -192,6 +192,7 @@ subroutine thermal_conduction_getSourceAndItsTangent(Tdot, dTdot_dT, T, ip, el) homogenization_Ngrains, & mappingHomogenization, & phaseAt, & + phasememberAt, & thermal_typeInstance, & phase_Nsources, & phase_source, & @@ -221,7 +222,8 @@ subroutine thermal_conduction_getSourceAndItsTangent(Tdot, dTdot_dT, T, ip, el) offset, & instance, & grain, & - source + source, & + constituent homog = mappingHomogenization(2,ip,el) offset = mappingHomogenization(1,ip,el) @@ -231,17 +233,18 @@ subroutine thermal_conduction_getSourceAndItsTangent(Tdot, dTdot_dT, T, ip, el) dTdot_dT = 0.0_pReal do grain = 1, homogenization_Ngrains(homog) phase = phaseAt(grain,ip,el) + constituent = phasememberAt(grain,ip,el) do source = 1, phase_Nsources(phase) select case(phase_source(source,phase)) case (SOURCE_thermal_dissipation_ID) call source_thermal_dissipation_getRateAndItsTangent(my_Tdot, my_dTdot_dT, & crystallite_Tstar_v(1:6,grain,ip,el), & crystallite_Lp(1:3,1:3,grain,ip,el), & - grain, ip, el) + phase, constituent) case (SOURCE_thermal_externalheat_ID) call source_thermal_externalheat_getRateAndItsTangent(my_Tdot, my_dTdot_dT, & - grain, ip, el) + phase, constituent) case default my_Tdot = 0.0_pReal @@ -363,8 +366,8 @@ function thermal_conduction_getMassDensity(ip,el) homog = mappingHomogenization(2,ip,el) do grain = 1, homogenization_Ngrains(mesh_element(3,el)) - thermal_conduction_getMassDensity = thermal_conduction_getMassDensity + & - lattice_massDensity(material_phase(grain,ip,el)) + thermal_conduction_getMassDensity = thermal_conduction_getMassDensity & + + lattice_massDensity(material_phase(grain,ip,el)) enddo thermal_conduction_getMassDensity = & From ced7da4d62815f6c5663da612a0c4b46120fd45f Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 30 Dec 2018 21:54:50 +0100 Subject: [PATCH 02/10] avoid mappings in bottom end functions --- src/homogenization.f90 | 10 +++-- src/thermal_adiabatic.f90 | 35 ++++++--------- src/thermal_conduction.f90 | 90 ++++++++++++++++++-------------------- src/thermal_isothermal.f90 | 30 ++++++------- 4 files changed, 77 insertions(+), 88 deletions(-) diff --git a/src/homogenization.f90 b/src/homogenization.f90 index 82a97dc53..0bf0f80be 100644 --- a/src/homogenization.f90 +++ b/src/homogenization.f90 @@ -1091,6 +1091,8 @@ function homogenization_postResults(ip,el) use mesh, only: & mesh_element use material, only: & + thermalMapping, & + thermal_typeInstance, & mappingHomogenization, & homogState, & thermalState, & @@ -1153,7 +1155,7 @@ function homogenization_postResults(ip,el) + hydrogenfluxState(mappingHomogenization(2,ip,el))%sizePostResults) :: & homogenization_postResults integer(pInt) :: & - startPos, endPos + startPos, endPos, homog homogenization_postResults = 0.0_pReal @@ -1184,11 +1186,13 @@ function homogenization_postResults(ip,el) case (THERMAL_isothermal_ID) chosenThermal case (THERMAL_adiabatic_ID) chosenThermal + homog = mappingHomogenization(2,ip,el) homogenization_postResults(startPos:endPos) = & - thermal_adiabatic_postResults(ip, el) + thermal_adiabatic_postResults(homog,thermal_typeInstance(homog),thermalMapping(homog)%p(ip,el)) case (THERMAL_conduction_ID) chosenThermal + homog = mappingHomogenization(2,ip,el) homogenization_postResults(startPos:endPos) = & - thermal_conduction_postResults(ip, el) + thermal_conduction_postResults(homog,thermal_typeInstance(homog),thermalMapping(homog)%p(ip,el)) end select chosenThermal startPos = endPos + 1_pInt diff --git a/src/thermal_adiabatic.f90 b/src/thermal_adiabatic.f90 index dd52293eb..e44030e64 100644 --- a/src/thermal_adiabatic.f90 +++ b/src/thermal_adiabatic.f90 @@ -10,8 +10,6 @@ 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 @@ -98,7 +96,6 @@ subroutine thermal_adiabatic_init(fileUnit) 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 = '' @@ -157,14 +154,13 @@ subroutine thermal_adiabatic_init(fileUnit) 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 ! allocate state arrays sizeState = 1_pInt thermalState(section)%sizeState = sizeState - thermalState(section)%sizePostResults = thermal_adiabatic_sizePostResults(instance) + thermalState(section)%sizePostResults = sum(thermal_adiabatic_sizePostResult(:,instance)) allocate(thermalState(section)%state0 (sizeState,NofMyHomog), source=thermal_initialT(section)) allocate(thermalState(section)%subState0(sizeState,NofMyHomog), source=thermal_initialT(section)) allocate(thermalState(section)%state (sizeState,NofMyHomog), source=thermal_initialT(section)) @@ -344,6 +340,7 @@ function thermal_adiabatic_getSpecificHeat(ip,el) end function thermal_adiabatic_getSpecificHeat + !-------------------------------------------------------------------------------------------------- !> @brief returns homogenized mass density !-------------------------------------------------------------------------------------------------- @@ -381,42 +378,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 + 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 1a9014e89..d9e10eece 100644 --- a/src/thermal_conduction.f90 +++ b/src/thermal_conduction.f90 @@ -10,8 +10,6 @@ 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 @@ -99,7 +97,6 @@ subroutine thermal_conduction_init(fileUnit) 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 = '' @@ -144,42 +141,40 @@ subroutine thermal_conduction_init(fileUnit) 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) !-------------------------------------------------------------------------------------------------- ! 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 + 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 + if (mySize > 0_pInt) then ! any meaningful output found + thermal_conduction_sizePostResult(o,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 !-------------------------------------------------------------------------------------------------- @@ -261,6 +256,7 @@ subroutine thermal_conduction_getSourceAndItsTangent(Tdot, dTdot_dT, T, ip, el) end subroutine thermal_conduction_getSourceAndItsTangent + !-------------------------------------------------------------------------------------------------- !> @brief returns homogenized thermal conductivity in reference configuration !-------------------------------------------------------------------------------------------------- @@ -298,7 +294,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 !-------------------------------------------------------------------------------------------------- @@ -374,7 +371,8 @@ function thermal_conduction_getMassDensity(ip,el) 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 !-------------------------------------------------------------------------------------------------- @@ -403,41 +401,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 From 9d2c60e943e9471b47a7d016e76003bdf0037c8c Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 30 Dec 2018 22:30:21 +0100 Subject: [PATCH 03/10] don't read material.config during init --- src/homogenization.f90 | 6 +-- src/thermal_adiabatic.f90 | 84 ++++++++++---------------------------- src/thermal_conduction.f90 | 73 ++++++++------------------------- 3 files changed, 40 insertions(+), 123 deletions(-) diff --git a/src/homogenization.f90 b/src/homogenization.f90 index 0bf0f80be..3f20ef7b4 100644 --- a/src/homogenization.f90 +++ b/src/homogenization.f90 @@ -139,11 +139,11 @@ subroutine homogenization_init ! parse thermal from config file call IO_checkAndRewind(FILEUNIT) if (any(thermal_type == THERMAL_isothermal_ID)) & - call thermal_isothermal_init() + call thermal_isothermal_init if (any(thermal_type == THERMAL_adiabatic_ID)) & - call thermal_adiabatic_init(FILEUNIT) + call thermal_adiabatic_init if (any(thermal_type == THERMAL_conduction_ID)) & - call thermal_conduction_init(FILEUNIT) + call thermal_conduction_init !-------------------------------------------------------------------------------------------------- ! parse damage from config file diff --git a/src/thermal_adiabatic.f90 b/src/thermal_adiabatic.f90 index e44030e64..e0ad3214f 100644 --- a/src/thermal_adiabatic.f90 +++ b/src/thermal_adiabatic.f90 @@ -13,7 +13,6 @@ module thermal_adiabatic 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 @@ -43,27 +42,15 @@ contains !> @brief module initialization !> @details reads in material parameters, allocates arrays, and does sanity checks !-------------------------------------------------------------------------------------------------- -subroutine thermal_adiabatic_init(fileUnit) +subroutine thermal_adiabatic_init #if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800 use, intrinsic :: iso_fortran_env, only: & compiler_version, & compiler_options #endif use IO, only: & - IO_read, & - IO_lc, & - IO_getTag, & - IO_isBlank, & - IO_stringPos, & - IO_stringValue, & - IO_floatValue, & - IO_intValue, & - IO_warning, & IO_error, & - IO_timeStamp, & - IO_EOF - use config, only: & - material_partHomogenization + IO_timeStamp use material, only: & thermal_type, & thermal_typeInstance, & @@ -77,17 +64,16 @@ 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,mySize=0_pInt,section,instance,o,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() @@ -102,47 +88,20 @@ subroutine thermal_adiabatic_init(fileUnit) 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) + 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) + end select + enddo !-------------------------------------------------------------------------------------------------- ! Determine size of postResults array @@ -172,9 +131,8 @@ subroutine thermal_adiabatic_init(fileUnit) deallocate(temperatureRate(section)%p) allocate (temperatureRate(section)%p(NofMyHomog), source=0.0_pReal) - endif - enddo initializeInstances + end subroutine thermal_adiabatic_init !-------------------------------------------------------------------------------------------------- diff --git a/src/thermal_conduction.f90 b/src/thermal_conduction.f90 index d9e10eece..067871c59 100644 --- a/src/thermal_conduction.f90 +++ b/src/thermal_conduction.f90 @@ -13,7 +13,6 @@ module thermal_conduction 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 @@ -44,25 +43,15 @@ contains !> @brief module initialization !> @details reads in material parameters, allocates arrays, and does sanity checks !-------------------------------------------------------------------------------------------------- -subroutine thermal_conduction_init(fileUnit) +subroutine thermal_conduction_init #if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800 use, intrinsic :: iso_fortran_env, only: & compiler_version, & compiler_options #endif use IO, only: & - IO_read, & - IO_lc, & - IO_getTag, & - IO_isBlank, & - IO_stringPos, & - IO_stringValue, & - IO_floatValue, & - IO_intValue, & - IO_warning, & IO_error, & - IO_timeStamp, & - IO_EOF + IO_timeStamp use material, only: & thermal_type, & thermal_typeInstance, & @@ -77,18 +66,15 @@ 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,mySize=0_pInt,section,instance,o,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() @@ -103,47 +89,20 @@ subroutine thermal_conduction_init(fileUnit) 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) 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) + end select + enddo !-------------------------------------------------------------------------------------------------- ! Determine size of postResults array From 69d53ed869d55531f4cd2ef4b6b01104196529e8 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 12 Feb 2019 23:20:24 +0100 Subject: [PATCH 04/10] determining output size was overly complicated general cleaning --- src/thermal_adiabatic.f90 | 61 ++++++++++++-------------------------- src/thermal_conduction.f90 | 30 +++---------------- 2 files changed, 23 insertions(+), 68 deletions(-) diff --git a/src/thermal_adiabatic.f90 b/src/thermal_adiabatic.f90 index e0ad3214f..937c20275 100644 --- a/src/thermal_adiabatic.f90 +++ b/src/thermal_adiabatic.f90 @@ -69,7 +69,7 @@ subroutine thermal_adiabatic_init config_homogenization implicit none - integer(pInt) :: maxNinstance,mySize=0_pInt,section,instance,o,i + integer(pInt) :: maxNinstance,section,instance,i integer(pInt) :: sizeState integer(pInt) :: NofMyHomog character(len=65536), dimension(0), parameter :: emptyStringArray = [character(len=65536)::] @@ -100,36 +100,24 @@ subroutine thermal_adiabatic_init 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 -!-------------------------------------------------------------------------------------------------- -! 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 - endif - enddo outputsLoop - ! allocate state arrays - 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)) + 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) enddo initializeInstances @@ -186,8 +174,6 @@ end function thermal_adiabatic_updateState !> @brief returns heat generation rate !-------------------------------------------------------------------------------------------------- subroutine thermal_adiabatic_getSourceAndItsTangent(Tdot, dTdot_dT, T, ip, el) - use math, only: & - math_Mandel6to33 use material, only: & homogenization_Ngrains, & mappingHomogenization, & @@ -219,14 +205,12 @@ subroutine thermal_adiabatic_getSourceAndItsTangent(Tdot, dTdot_dT, T, ip, el) integer(pInt) :: & phase, & homog, & - offset, & instance, & grain, & source, & constituent homog = mappingHomogenization(2,ip,el) - offset = mappingHomogenization(1,ip,el) instance = thermal_typeInstance(homog) Tdot = 0.0_pReal @@ -268,12 +252,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) :: & @@ -282,11 +263,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 + & @@ -311,9 +291,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 @@ -321,11 +299,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 + & @@ -346,7 +323,7 @@ function thermal_adiabatic_postResults(homog,instance,of) result(postResults) temperature implicit none - integer(pInt), intent(in) :: & + integer(pInt), intent(in) :: & homog, & instance, & of diff --git a/src/thermal_conduction.f90 b/src/thermal_conduction.f90 index 067871c59..ab1b030c8 100644 --- a/src/thermal_conduction.f90 +++ b/src/thermal_conduction.f90 @@ -70,7 +70,7 @@ subroutine thermal_conduction_init config_homogenization implicit none - integer(pInt) :: maxNinstance,mySize=0_pInt,section,instance,o,i + integer(pInt) :: maxNinstance,section,instance,i integer(pInt) :: sizeState integer(pInt) :: NofMyHomog character(len=65536), dimension(0), parameter :: emptyStringArray = [character(len=65536)::] @@ -101,21 +101,10 @@ subroutine thermal_conduction_init 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 - endif - enddo outputsLoop ! allocate state arrays sizeState = 0_pInt @@ -224,7 +213,6 @@ function thermal_conduction_getConductivity33(ip,el) lattice_thermalConductivity33 use material, only: & homogenization_Ngrains, & - mappingHomogenization, & material_phase use mesh, only: & mesh_element @@ -238,10 +226,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)) @@ -263,12 +249,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) :: & @@ -277,11 +260,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 + & @@ -301,12 +283,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) :: & @@ -315,11 +294,10 @@ 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 & From c9fc7ea982ee776789b74c94b3f1fd60b636b665 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 12 Feb 2019 23:35:22 +0100 Subject: [PATCH 05/10] cleaning trying to find logic with less dependencies on the various mappings --- src/constitutive.f90 | 15 ++++++++--- src/kinematics_thermal_expansion.f90 | 40 +++++++++++++--------------- 2 files changed, 31 insertions(+), 24 deletions(-) diff --git a/src/constitutive.f90 b/src/constitutive.f90 index 1449c35e9..5dc59c2c3 100644 --- a/src/constitutive.f90 +++ b/src/constitutive.f90 @@ -649,6 +649,9 @@ pure function constitutive_initialFi(ipc, ip, el) math_inv33, & math_mul33x33 use material, only: & + material_phase, & + material_homog, & + thermalMapping, & phase_kinematics, & phase_Nkinematics, & material_phase, & @@ -665,14 +668,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 diff --git a/src/kinematics_thermal_expansion.f90 b/src/kinematics_thermal_expansion.f90 index e8f0d71c7..a44bc6902 100644 --- a/src/kinematics_thermal_expansion.f90 +++ b/src/kinematics_thermal_expansion.f90 @@ -10,7 +10,11 @@ module kinematics_thermal_expansion implicit none private - + + !type, private :: tParameters + ! real(pReal), allocatable, dimension(:) :: & + !end type tParameters + public :: & kinematics_thermal_expansion_init, & kinematics_thermal_expansion_initialStrain, & @@ -43,49 +47,43 @@ subroutine kinematics_thermal_expansion_init() config_phase implicit none - integer(pInt) maxNinstance + integer(pInt) :: & + Ninstance, & + p 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 + Ninstance = int(count(phase_kinematics == KINEMATICS_thermal_expansion_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 -! ToDo: this subroutine should read in lattice_thermal_expansion. No need to make it a global array + do p = 1_pInt, size(phase_kinematics) + if (all(phase_kinematics(:,p) /= KINEMATICS_thermal_expansion_ID)) cycle + 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) +pure function kinematics_thermal_expansion_initialStrain(homog,phase,offset) use material, only: & - material_phase, & - material_homog, & - temperature, & - thermalMapping + temperature 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), 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) + real(pReal), dimension(3,3) :: & + kinematics_thermal_expansion_initialStrain !< initial thermal strain (should be small strain, though) + kinematics_thermal_expansion_initialStrain = & (temperature(homog)%p(offset) - lattice_referenceTemperature(phase))**1 / 1. * & From af4ea76006177a23f1d58dc68b5ccabc4ce8a87a Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 22 Feb 2019 15:21:48 +0100 Subject: [PATCH 06/10] using central allocation facilities --- src/source_thermal_externalheat.f90 | 39 ++++------------------------- 1 file changed, 5 insertions(+), 34 deletions(-) diff --git a/src/source_thermal_externalheat.f90 b/src/source_thermal_externalheat.f90 index eac1232f3..724d36f17 100644 --- a/src/source_thermal_externalheat.f90 +++ b/src/source_thermal_externalheat.f90 @@ -57,11 +57,6 @@ contains !> @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,& @@ -77,9 +72,9 @@ subroutine source_thermal_externalheat_init(fileUnit) IO_intValue, & IO_warning, & IO_error, & - IO_timeStamp, & IO_EOF use material, only: & + material_allocateSourceState, & phase_source, & phase_Nsources, & phase_Noutput, & @@ -91,8 +86,6 @@ subroutine source_thermal_externalheat_init(fileUnit) config_phase, & material_Nphase, & MATERIAL_partPhase - use numerics,only: & - numerics_integrator implicit none integer(pInt), intent(in) :: fileUnit @@ -107,8 +100,7 @@ subroutine source_thermal_externalheat_init(fileUnit) 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" + maxNinstance = int(count(phase_source == SOURCE_thermal_externalheat_ID),pInt) if (maxNinstance == 0_pInt) return @@ -196,35 +188,14 @@ subroutine source_thermal_externalheat_init(fileUnit) 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) - + call material_allocateSourceState(phase,sourceOffset,NofMyPhase,1_pInt,1_pInt,0_pInt) + endif enddo initializeInstances end subroutine source_thermal_externalheat_init + !-------------------------------------------------------------------------------------------------- !> @brief rate of change of state !> @details state only contains current time to linearly interpolate given heat powers From e7268ce109da3bb5e149d05681b68291927fb899 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 22 Feb 2019 15:37:42 +0100 Subject: [PATCH 07/10] simpler structure: - do not read file - use function for allocation - do not constantly convert (3,3) <-> (6) --- src/constitutive.f90 | 12 +-- src/source_thermal_dissipation.f90 | 133 ++++++----------------------- src/thermal_adiabatic.f90 | 16 ++-- src/thermal_conduction.f90 | 16 +--- 4 files changed, 35 insertions(+), 142 deletions(-) diff --git a/src/constitutive.f90 b/src/constitutive.f90 index 2fba40328..48d4a3f8f 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: & @@ -55,8 +50,7 @@ subroutine constitutive_init() IO_open_file, & IO_checkAndRewind, & IO_open_jobFile_stat, & - IO_write_jobFile, & - IO_timeStamp + IO_write_jobFile use config, only: & config_phase use config, only: & @@ -158,7 +152,7 @@ subroutine constitutive_init() !-------------------------------------------------------------------------------------------------- ! parse source mechanisms from config file call IO_checkAndRewind(FILEUNIT) - if (any(phase_source == SOURCE_thermal_dissipation_ID)) call source_thermal_dissipation_init(FILEUNIT) + 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(FILEUNIT) 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 @@ -176,8 +170,6 @@ subroutine constitutive_init() call config_deallocate('material.config/phase') write(6,'(/,a)') ' <<<+- constitutive init -+>>>' - write(6,'(a15,a)') ' Current time: ',IO_timeStamp() -#include "compilation_info.f90" mainProcess: if (worldrank == 0) then !-------------------------------------------------------------------------------------------------- diff --git a/src/source_thermal_dissipation.f90 b/src/source_thermal_dissipation.f90 index 290ad7efe..ef0843239 100644 --- a/src/source_thermal_dissipation.f90 +++ b/src/source_thermal_dissipation.f90 @@ -47,30 +47,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, & @@ -82,23 +65,13 @@ subroutine source_thermal_dissipation_init(fileUnit) 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) :: Ninstance,phase,instance,source,sourceOffset - integer(pInt) :: sizeState, sizeDotState, sizeDeltaState + integer(pInt) :: Ninstance,instance,source,sourceOffset integer(pInt) :: NofMyPhase,p - character(len=65536) :: & - tag = '', & - line = '' write(6,'(/,a)') ' <<<+- source_'//SOURCE_thermal_dissipation_label//' init -+>>>' - write(6,'(a15,a)') ' Current time: ',IO_timeStamp() -#include "compilation_info.f90" + Ninstance = int(count(phase_source == SOURCE_thermal_dissipation_ID),pInt) if (Ninstance == 0_pInt) return @@ -107,11 +80,11 @@ subroutine source_thermal_dissipation_init(fileUnit) 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 @@ -124,88 +97,31 @@ subroutine source_thermal_dissipation_init(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 + 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 - - 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_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, phase, constituent) +subroutine source_thermal_dissipation_getRateAndItsTangent(TDot, dTDOT_dT, Tstar, Lp, phase) use math, only: & math_Mandel6to33 implicit none integer(pInt), intent(in) :: & - phase, & - constituent - real(pReal), intent(in), dimension(6) :: & - Tstar_v !< 2nd Piola Kirchhoff stress tensor (Mandel) + phase + real(pReal), intent(in), dimension(3,3) :: & + Tstar !< 2nd Piola Kirchhoff stress tensor (Mandel) real(pReal), intent(in), dimension(3,3) :: & Lp real(pReal), intent(out) :: & @@ -216,8 +132,7 @@ subroutine source_thermal_dissipation_getRateAndItsTangent(TDot, dTDOT_dT, Tstar 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/thermal_adiabatic.f90 b/src/thermal_adiabatic.f90 index 937c20275..947e28777 100644 --- a/src/thermal_adiabatic.f90 +++ b/src/thermal_adiabatic.f90 @@ -43,14 +43,8 @@ contains !> @details reads in material parameters, allocates arrays, and does sanity checks !-------------------------------------------------------------------------------------------------- subroutine thermal_adiabatic_init -#if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800 - use, intrinsic :: iso_fortran_env, only: & - compiler_version, & - compiler_options -#endif use IO, only: & - IO_error, & - IO_timeStamp + IO_error use material, only: & thermal_type, & thermal_typeInstance, & @@ -76,8 +70,6 @@ subroutine thermal_adiabatic_init 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 @@ -174,6 +166,8 @@ end function thermal_adiabatic_updateState !> @brief returns heat generation rate !-------------------------------------------------------------------------------------------------- subroutine thermal_adiabatic_getSourceAndItsTangent(Tdot, dTdot_dT, T, ip, el) + use math, only: & + math_6toSym33 use material, only: & homogenization_Ngrains, & mappingHomogenization, & @@ -222,9 +216,9 @@ subroutine thermal_adiabatic_getSourceAndItsTangent(Tdot, dTdot_dT, T, ip, el) 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), & - phase, constituent) + phase) case (SOURCE_thermal_externalheat_ID) call source_thermal_externalheat_getRateAndItsTangent(my_Tdot, my_dTdot_dT, & diff --git a/src/thermal_conduction.f90 b/src/thermal_conduction.f90 index ab1b030c8..3d5ca892e 100644 --- a/src/thermal_conduction.f90 +++ b/src/thermal_conduction.f90 @@ -44,14 +44,8 @@ contains !> @details reads in material parameters, allocates arrays, and does sanity checks !-------------------------------------------------------------------------------------------------- subroutine thermal_conduction_init -#if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800 - use, intrinsic :: iso_fortran_env, only: & - compiler_version, & - compiler_options -#endif use IO, only: & - IO_error, & - IO_timeStamp + IO_error use material, only: & thermal_type, & thermal_typeInstance, & @@ -77,8 +71,6 @@ subroutine thermal_conduction_init 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 @@ -130,7 +122,7 @@ 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, & @@ -181,9 +173,9 @@ subroutine thermal_conduction_getSourceAndItsTangent(Tdot, dTdot_dT, T, ip, el) 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), & - phase, constituent) + phase) case (SOURCE_thermal_externalheat_ID) call source_thermal_externalheat_getRateAndItsTangent(my_Tdot, my_dTdot_dT, & From 194824fd0fa25a581586428c9630797c36429c85 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 22 Feb 2019 20:37:41 +0100 Subject: [PATCH 08/10] WIP: cleaned no file reading getting rid of a number of obsolete dependencies --- src/constitutive.f90 | 2 +- src/crystallite.f90 | 2 +- src/kinematics_thermal_expansion.f90 | 229 ++++++++++++----------- src/material.f90 | 8 +- src/numerics.f90 | 10 +- src/source_thermal_dissipation.f90 | 1 + src/source_thermal_externalheat.f90 | 270 ++++++++++----------------- src/thermal_adiabatic.f90 | 2 - src/thermal_conduction.f90 | 2 - 9 files changed, 225 insertions(+), 301 deletions(-) diff --git a/src/constitutive.f90 b/src/constitutive.f90 index 48d4a3f8f..a340d0adf 100644 --- a/src/constitutive.f90 +++ b/src/constitutive.f90 @@ -153,7 +153,7 @@ subroutine constitutive_init() ! parse source mechanisms from config file call IO_checkAndRewind(FILEUNIT) 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(FILEUNIT) + 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 diff --git a/src/crystallite.f90 b/src/crystallite.f90 index c272abd07..ec9c782ea 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -248,7 +248,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) diff --git a/src/kinematics_thermal_expansion.f90 b/src/kinematics_thermal_expansion.f90 index a44bc6902..56caa6e4b 100644 --- a/src/kinematics_thermal_expansion.f90 +++ b/src/kinematics_thermal_expansion.f90 @@ -4,21 +4,24 @@ !> @details to be done !-------------------------------------------------------------------------------------------------- module kinematics_thermal_expansion - use prec, only: & - pReal, & - pInt + use prec, only: & + pReal, & + pInt - implicit none - private + implicit none + private - !type, private :: tParameters - ! real(pReal), allocatable, dimension(:) :: & - !end type tParameters + type, private :: tParameters + real(pReal), allocatable, dimension(:,:,:) :: & + expansion + end type tParameters - public :: & - kinematics_thermal_expansion_init, & - kinematics_thermal_expansion_initialStrain, & - kinematics_thermal_expansion_LiAndItsTangent + type(tParameters), dimension(:), allocatable :: param + + public :: & + kinematics_thermal_expansion_init, & + kinematics_thermal_expansion_initialStrain, & + kinematics_thermal_expansion_LiAndItsTangent contains @@ -28,120 +31,128 @@ contains !> @details reads in material parameters, allocates arrays, and does sanity checks !-------------------------------------------------------------------------------------------------- subroutine kinematics_thermal_expansion_init() -#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_timeStamp - use material, only: & - phase_kinematics, & - KINEMATICS_thermal_expansion_label, & - KINEMATICS_thermal_expansion_ID - use config, only: & - config_phase - - implicit none - integer(pInt) :: & - Ninstance, & - p + 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 - write(6,'(/,a)') ' <<<+- kinematics_'//KINEMATICS_thermal_expansion_LABEL//' init -+>>>' - write(6,'(a15,a)') ' Current time: ',IO_timeStamp() -#include "compilation_info.f90" - - Ninstance = int(count(phase_kinematics == KINEMATICS_thermal_expansion_ID),pInt) - if (Ninstance == 0_pInt) return + implicit none + integer(pInt) :: & + Ninstance, & + p, i + real(pReal), dimension(:), allocatable :: & + temp + + write(6,'(/,a)') ' <<<+- kinematics_'//KINEMATICS_thermal_expansion_LABEL//' init -+>>>' - if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt) & - write(6,'(a16,1x,i5,/)') '# instances:',Ninstance - - do p = 1_pInt, size(phase_kinematics) - if (all(phase_kinematics(:,p) /= KINEMATICS_thermal_expansion_ID)) cycle - enddo + 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(homog,phase,offset) - use material, only: & - temperature - use lattice, only: & - lattice_thermalExpansion33, & - lattice_referenceTemperature + 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) :: & - phase, & - homog, offset - real(pReal), dimension(3,3) :: & - kinematics_thermal_expansion_initialStrain !< initial thermal strain (should be small strain, though) - - - 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_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 + 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..291b73910 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 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 ef0843239..026a71726 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 diff --git a/src/source_thermal_externalheat.f90 b/src/source_thermal_externalheat.f90 index 724d36f17..2bf4cac9c 100644 --- a/src/source_thermal_externalheat.f90 +++ b/src/source_thermal_externalheat.f90 @@ -1,53 +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) - type, private :: tParameters !< container type for internal constitutive parameters - real(pReal), dimension(:), allocatable :: & - time, & - rate - integer(pInt) :: & - nInterval - end type tParameters - - type(tParameters), dimension(:), allocatable, private :: param !< containers of constitutive parameters (len Ninstance) - - - public :: & - source_thermal_externalheat_init, & - source_thermal_externalheat_dotState, & - source_thermal_externalheat_getRateAndItsTangent + public :: & + source_thermal_externalheat_init, & + source_thermal_externalheat_dotState, & + source_thermal_externalheat_getRateAndItsTangent contains @@ -56,143 +47,74 @@ contains !> @brief module initialization !> @details reads in material parameters, allocates arrays, and does sanity checks !-------------------------------------------------------------------------------------------------- -subroutine source_thermal_externalheat_init(fileUnit) - 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_EOF - use material, only: & - material_allocateSourceState, & - phase_source, & - phase_Nsources, & - phase_Noutput, & - SOURCE_thermal_externalheat_label, & - SOURCE_thermal_externalheat_ID, & - material_phase, & - sourceState - use config, only: & - config_phase, & - material_Nphase, & - MATERIAL_partPhase - - 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,p - character(len=65536) :: & - tag = '', & - line = '' - real(pReal), allocatable, dimension(:,:) :: temp_time, temp_rate - - write(6,'(/,a)') ' <<<+- source_'//SOURCE_thermal_externalheat_label//' init -+>>>' - +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) - - do p=1, size(config_phase) - if (all(phase_source(:,p) /= SOURCE_thermal_externalheat_ID)) cycle - enddo - - rewind(fileUnit) - phase = 0_pInt - do while (trim(line) /= IO_EOF .and. IO_lc(IO_getTag(line,'<','>')) /= MATERIAL_partPhase) ! wind forward to - 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) - - call material_allocateSourceState(phase,sourceOffset,NofMyPhase,1_pInt,1_pInt,0_pInt) + + 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 - endif + 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 @@ -245,16 +167,16 @@ subroutine source_thermal_externalheat_getRateAndItsTangent(TDot, dTDot_dT, phas 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 947e28777..c3290bdfe 100644 --- a/src/thermal_adiabatic.f90 +++ b/src/thermal_adiabatic.f90 @@ -43,8 +43,6 @@ contains !> @details reads in material parameters, allocates arrays, and does sanity checks !-------------------------------------------------------------------------------------------------- subroutine thermal_adiabatic_init - use IO, only: & - IO_error use material, only: & thermal_type, & thermal_typeInstance, & diff --git a/src/thermal_conduction.f90 b/src/thermal_conduction.f90 index 3d5ca892e..88da0529b 100644 --- a/src/thermal_conduction.f90 +++ b/src/thermal_conduction.f90 @@ -44,8 +44,6 @@ contains !> @details reads in material parameters, allocates arrays, and does sanity checks !-------------------------------------------------------------------------------------------------- subroutine thermal_conduction_init - use IO, only: & - IO_error use material, only: & thermal_type, & thermal_typeInstance, & From ad0ed4fdec46f7d96767269ee65a6210f8083fbd Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 22 Feb 2019 21:06:37 +0100 Subject: [PATCH 09/10] bugfix: wrong state was allocated --- src/material.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/material.f90 b/src/material.f90 index 291b73910..49ee38ee3 100644 --- a/src/material.f90 +++ b/src/material.f90 @@ -999,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 From b1bb68d52301e4942b857ed1f89b7057000bc3f3 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 22 Feb 2019 21:07:00 +0100 Subject: [PATCH 10/10] cleaning --- src/source_thermal_dissipation.f90 | 10 +--------- 1 file changed, 1 insertion(+), 9 deletions(-) diff --git a/src/source_thermal_dissipation.f90 b/src/source_thermal_dissipation.f90 index 026a71726..db37c8286 100644 --- a/src/source_thermal_dissipation.f90 +++ b/src/source_thermal_dissipation.f90 @@ -12,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 @@ -21,9 +20,6 @@ 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 @@ -89,11 +85,9 @@ subroutine source_thermal_dissipation_init enddo enddo - allocate(source_thermal_dissipation_sizePostResults(Ninstance), source=0_pInt) allocate(source_thermal_dissipation_sizePostResult(maxval(phase_Noutput),Ninstance),source=0_pInt) allocate(source_thermal_dissipation_output (maxval(phase_Noutput),Ninstance)) source_thermal_dissipation_output = '' - allocate(source_thermal_dissipation_Noutput(Ninstance), source=0_pInt) allocate(source_thermal_dissipation_coldworkCoeff(Ninstance), source=0.0_pReal) @@ -115,14 +109,12 @@ end subroutine source_thermal_dissipation_init !> @brief returns local vacancy generation rate !-------------------------------------------------------------------------------------------------- subroutine source_thermal_dissipation_getRateAndItsTangent(TDot, dTDOT_dT, Tstar, Lp, phase) - use math, only: & - math_Mandel6to33 implicit none integer(pInt), intent(in) :: & phase real(pReal), intent(in), dimension(3,3) :: & - Tstar !< 2nd Piola Kirchhoff stress tensor (Mandel) + Tstar real(pReal), intent(in), dimension(3,3) :: & Lp real(pReal), intent(out) :: &