diff --git a/src/CPFEM2.f90 b/src/CPFEM2.f90 index cb8b1654a..51d9af152 100644 --- a/src/CPFEM2.f90 +++ b/src/CPFEM2.f90 @@ -22,6 +22,9 @@ module CPFEM2 use homogenization use constitutive use crystallite +#ifdef FEM + use FEM_Zoo +#endif implicit none private @@ -37,33 +40,29 @@ contains !-------------------------------------------------------------------------------------------------- !> @brief call (thread safe) all module initializations !-------------------------------------------------------------------------------------------------- -subroutine CPFEM_initAll() -#ifdef FEM - use FEM_Zoo, only: & - FEM_Zoo_init -#endif +subroutine CPFEM_initAll - call DAMASK_interface_init ! Spectral and FEM interface to commandline - call prec_init - call IO_init + call DAMASK_interface_init ! Spectral and FEM interface to commandline + call prec_init + call IO_init #ifdef FEM - call FEM_Zoo_init + call FEM_Zoo_init #endif - call numerics_init - call debug_init - call config_init - call math_init - call FE_init - call mesh_init - call lattice_init - call HDF5_utilities_init - call results_init - call material_init - call constitutive_init - call crystallite_init - call homogenization_init - call materialpoint_postResults - call CPFEM_init + call numerics_init + call debug_init + call config_init + call math_init + call FE_init + call mesh_init + call lattice_init + call HDF5_utilities_init + call results_init + call material_init + call constitutive_init + call crystallite_init + call homogenization_init + call materialpoint_postResults + call CPFEM_init end subroutine CPFEM_initAll @@ -72,50 +71,50 @@ end subroutine CPFEM_initAll !-------------------------------------------------------------------------------------------------- subroutine CPFEM_init - integer :: ph,homog - character(len=1024) :: rankStr, PlasticItem, HomogItem - integer(HID_T) :: fileHandle, groupPlasticID, groupHomogID + integer :: ph,homog + character(len=1024) :: rankStr, PlasticItem, HomogItem + integer(HID_T) :: fileHandle, groupPlasticID, groupHomogID - write(6,'(/,a)') ' <<<+- CPFEM init -+>>>' - flush(6) + write(6,'(/,a)') ' <<<+- CPFEM init -+>>>' + flush(6) - ! *** restore the last converged values of each essential variable from the binary file - if (restartRead) then - if (iand(debug_level(debug_CPFEM), debug_levelExtensive) /= 0) then - write(6,'(a)') '<< CPFEM >> restored state variables of last converged step from hdf5 file' - flush(6) - endif - - write(rankStr,'(a1,i0)')'_',worldrank + ! *** restore the last converged values of each essential variable from the binary file + if (restartRead) then + if (iand(debug_level(debug_CPFEM), debug_levelExtensive) /= 0) then + write(6,'(a)') '<< CPFEM >> restored state variables of last converged step from hdf5 file' + flush(6) + endif + + write(rankStr,'(a1,i0)')'_',worldrank - fileHandle = HDF5_openFile(trim(getSolverJobName())//trim(rankStr)//'.hdf5') - - call HDF5_read(fileHandle,material_phase, 'recordedPhase') - call HDF5_read(fileHandle,crystallite_F0, 'convergedF') - call HDF5_read(fileHandle,crystallite_Fp0, 'convergedFp') - call HDF5_read(fileHandle,crystallite_Fi0, 'convergedFi') - call HDF5_read(fileHandle,crystallite_Lp0, 'convergedLp') - call HDF5_read(fileHandle,crystallite_Li0, 'convergedLi') - call HDF5_read(fileHandle,crystallite_S0, 'convergedS') + fileHandle = HDF5_openFile(trim(getSolverJobName())//trim(rankStr)//'.hdf5') - groupPlasticID = HDF5_openGroup(fileHandle,'PlasticPhases') - do ph = 1,size(phase_plasticity) - write(PlasticItem,*) ph,'_' - call HDF5_read(groupPlasticID,plasticState(ph)%state0,trim(PlasticItem)//'convergedStateConst') - enddo - call HDF5_closeGroup(groupPlasticID) - - groupHomogID = HDF5_openGroup(fileHandle,'HomogStates') - do homog = 1, material_Nhomogenization - write(HomogItem,*) homog,'_' - call HDF5_read(groupHomogID,homogState(homog)%state0, trim(HomogItem)//'convergedStateHomog') - enddo - call HDF5_closeGroup(groupHomogID) + call HDF5_read(fileHandle,material_phase, 'recordedPhase') + call HDF5_read(fileHandle,crystallite_F0, 'convergedF') + call HDF5_read(fileHandle,crystallite_Fp0, 'convergedFp') + call HDF5_read(fileHandle,crystallite_Fi0, 'convergedFi') + call HDF5_read(fileHandle,crystallite_Lp0, 'convergedLp') + call HDF5_read(fileHandle,crystallite_Li0, 'convergedLi') + call HDF5_read(fileHandle,crystallite_S0, 'convergedS') + + groupPlasticID = HDF5_openGroup(fileHandle,'PlasticPhases') + do ph = 1,size(phase_plasticity) + write(PlasticItem,*) ph,'_' + call HDF5_read(groupPlasticID,plasticState(ph)%state0,trim(PlasticItem)//'convergedStateConst') + enddo + call HDF5_closeGroup(groupPlasticID) + + groupHomogID = HDF5_openGroup(fileHandle,'HomogStates') + do homog = 1, material_Nhomogenization + write(HomogItem,*) homog,'_' + call HDF5_read(groupHomogID,homogState(homog)%state0, trim(HomogItem)//'convergedStateHomog') + enddo + call HDF5_closeGroup(groupHomogID) - call HDF5_closeFile(fileHandle) - - restartRead = .false. - endif + call HDF5_closeFile(fileHandle) + + restartRead = .false. + endif end subroutine CPFEM_init @@ -125,68 +124,68 @@ end subroutine CPFEM_init !-------------------------------------------------------------------------------------------------- subroutine CPFEM_age - integer :: i, ph, homog, mySource - character(len=32) :: rankStr, PlasticItem, HomogItem - integer(HID_T) :: fileHandle, groupPlastic, groupHomog + integer :: i, ph, homog, mySource + character(len=32) :: rankStr, PlasticItem, HomogItem + integer(HID_T) :: fileHandle, groupPlastic, groupHomog - if (iand(debug_level(debug_CPFEM), debug_levelBasic) /= 0) & - write(6,'(a)') '<< CPFEM >> aging states' + if (iand(debug_level(debug_CPFEM), debug_levelBasic) /= 0) & + write(6,'(a)') '<< CPFEM >> aging states' - crystallite_F0 = crystallite_partionedF - crystallite_Fp0 = crystallite_Fp - crystallite_Lp0 = crystallite_Lp - crystallite_Fi0 = crystallite_Fi - crystallite_Li0 = crystallite_Li - crystallite_S0 = crystallite_S + crystallite_F0 = crystallite_partionedF + crystallite_Fp0 = crystallite_Fp + crystallite_Lp0 = crystallite_Lp + crystallite_Fi0 = crystallite_Fi + crystallite_Li0 = crystallite_Li + crystallite_S0 = crystallite_S - do i = 1, size(plasticState) - plasticState(i)%state0 = plasticState(i)%state - enddo - do i = 1, size(sourceState) - do mySource = 1,phase_Nsources(i) - sourceState(i)%p(mySource)%state0 = sourceState(i)%p(mySource)%state - enddo; enddo - do homog = 1, material_Nhomogenization - homogState (homog)%state0 = homogState (homog)%state - thermalState (homog)%state0 = thermalState (homog)%state - damageState (homog)%state0 = damageState (homog)%state - enddo + do i = 1, size(plasticState) + plasticState(i)%state0 = plasticState(i)%state + enddo + do i = 1, size(sourceState) + do mySource = 1,phase_Nsources(i) + sourceState(i)%p(mySource)%state0 = sourceState(i)%p(mySource)%state + enddo; enddo + do homog = 1, material_Nhomogenization + homogState (homog)%state0 = homogState (homog)%state + thermalState (homog)%state0 = thermalState (homog)%state + damageState (homog)%state0 = damageState (homog)%state + enddo - if (restartWrite) then - if (iand(debug_level(debug_CPFEM), debug_levelBasic) /= 0) & - write(6,'(a)') '<< CPFEM >> writing restart variables of last converged step to hdf5 file' - - write(rankStr,'(a1,i0)')'_',worldrank - fileHandle = HDF5_openFile(trim(getSolverJobName())//trim(rankStr)//'.hdf5','a') - - call HDF5_write(fileHandle,material_phase, 'recordedPhase') - call HDF5_write(fileHandle,crystallite_F0, 'convergedF') - call HDF5_write(fileHandle,crystallite_Fp0, 'convergedFp') - call HDF5_write(fileHandle,crystallite_Fi0, 'convergedFi') - call HDF5_write(fileHandle,crystallite_Lp0, 'convergedLp') - call HDF5_write(fileHandle,crystallite_Li0, 'convergedLi') - call HDF5_write(fileHandle,crystallite_S0, 'convergedS') - - groupPlastic = HDF5_addGroup(fileHandle,'PlasticPhases') - do ph = 1,size(phase_plasticity) - write(PlasticItem,*) ph,'_' - call HDF5_write(groupPlastic,plasticState(ph)%state0,trim(PlasticItem)//'convergedStateConst') - enddo - call HDF5_closeGroup(groupPlastic) + if (restartWrite) then + if (iand(debug_level(debug_CPFEM), debug_levelBasic) /= 0) & + write(6,'(a)') '<< CPFEM >> writing restart variables of last converged step to hdf5 file' + + write(rankStr,'(a1,i0)')'_',worldrank + fileHandle = HDF5_openFile(trim(getSolverJobName())//trim(rankStr)//'.hdf5','a') + + call HDF5_write(fileHandle,material_phase, 'recordedPhase') + call HDF5_write(fileHandle,crystallite_F0, 'convergedF') + call HDF5_write(fileHandle,crystallite_Fp0, 'convergedFp') + call HDF5_write(fileHandle,crystallite_Fi0, 'convergedFi') + call HDF5_write(fileHandle,crystallite_Lp0, 'convergedLp') + call HDF5_write(fileHandle,crystallite_Li0, 'convergedLi') + call HDF5_write(fileHandle,crystallite_S0, 'convergedS') + + groupPlastic = HDF5_addGroup(fileHandle,'PlasticPhases') + do ph = 1,size(phase_plasticity) + write(PlasticItem,*) ph,'_' + call HDF5_write(groupPlastic,plasticState(ph)%state0,trim(PlasticItem)//'convergedStateConst') + enddo + call HDF5_closeGroup(groupPlastic) - groupHomog = HDF5_addGroup(fileHandle,'HomogStates') - do homog = 1, material_Nhomogenization - write(HomogItem,*) homog,'_' - call HDF5_write(groupHomog,homogState(homog)%state0,trim(HomogItem)//'convergedStateHomog') - enddo - call HDF5_closeGroup(groupHomog) - - call HDF5_closeFile(fileHandle) - restartWrite = .false. - endif + groupHomog = HDF5_addGroup(fileHandle,'HomogStates') + do homog = 1, material_Nhomogenization + write(HomogItem,*) homog,'_' + call HDF5_write(groupHomog,homogState(homog)%state0,trim(HomogItem)//'convergedStateHomog') + enddo + call HDF5_closeGroup(groupHomog) + + call HDF5_closeFile(fileHandle) + restartWrite = .false. + endif - if (iand(debug_level(debug_CPFEM), debug_levelBasic) /= 0) & - write(6,'(a)') '<< CPFEM >> done aging states' + if (iand(debug_level(debug_CPFEM), debug_levelBasic) /= 0) & + write(6,'(a)') '<< CPFEM >> done aging states' end subroutine CPFEM_age @@ -196,17 +195,17 @@ end subroutine CPFEM_age !-------------------------------------------------------------------------------------------------- subroutine CPFEM_results(inc,time) - integer, intent(in) :: inc - real(pReal), intent(in) :: time - - call results_openJobFile - call results_addIncrement(inc,time) - call constitutive_results - call crystallite_results - call homogenization_results - call discretization_results - call results_removeLink('current') ! ToDo: put this into closeJobFile - call results_closeJobFile + integer, intent(in) :: inc + real(pReal), intent(in) :: time + + call results_openJobFile + call results_addIncrement(inc,time) + call constitutive_results + call crystallite_results + call homogenization_results + call discretization_results + call results_removeLink('current') ! ToDo: put this into closeJobFile + call results_closeJobFile end subroutine CPFEM_results diff --git a/src/FEsolving.f90 b/src/FEsolving.f90 index c188e66e2..7e4742139 100644 --- a/src/FEsolving.f90 +++ b/src/FEsolving.f90 @@ -5,32 +5,35 @@ !> @todo Descriptions for public variables needed !-------------------------------------------------------------------------------------------------- module FEsolving - use prec - - implicit none - private - integer, public :: & - restartInc = 1 !< needs description - - logical, public :: & - symmetricSolver = .false., & !< use a symmetric FEM solver - restartWrite = .false., & !< write current state to enable restart - restartRead = .false., & !< restart information to continue calculation from saved state - terminallyIll = .false. !< at least one material point is terminally ill - - integer, dimension(:,:), allocatable, public :: & - FEsolving_execIP !< for ping-pong scheme always range to max IP, otherwise one specific IP + use prec + use debug + use IO + use DAMASK_interface - integer, dimension(2), public :: & - FEsolving_execElem !< for ping-pong scheme always whole range, otherwise one specific element - - character(len=1024), public :: & - modelName !< needs description - - logical, dimension(:,:), allocatable, public :: & - calcMode !< do calculation or simply collect when using ping pong scheme + implicit none + private + integer, public :: & + restartInc = 1 !< needs description - public :: FE_init + logical, public :: & + symmetricSolver = .false., & !< use a symmetric FEM solver + restartWrite = .false., & !< write current state to enable restart + restartRead = .false., & !< restart information to continue calculation from saved state + terminallyIll = .false. !< at least one material point is terminally ill + + integer, dimension(:,:), allocatable, public :: & + FEsolving_execIP !< for ping-pong scheme always range to max IP, otherwise one specific IP + + integer, dimension(2), public :: & + FEsolving_execElem !< for ping-pong scheme always whole range, otherwise one specific element + + character(len=1024), public :: & + modelName !< needs description + + logical, dimension(:,:), allocatable, public :: & + calcMode !< do calculation or simply collect when using ping pong scheme + + public :: FE_init contains @@ -41,108 +44,93 @@ contains !> solver the information is provided by the interface module !-------------------------------------------------------------------------------------------------- subroutine FE_init - use debug, only: & - debug_level, & - debug_FEsolving, & - debug_levelBasic - use IO, only: & - IO_stringPos, & - IO_stringValue, & - IO_intValue, & - IO_lc, & -#if defined(Marc4DAMASK) || defined(Abaqus) - IO_open_inputFile, & - IO_open_logFile, & -#endif - IO_warning - use DAMASK_interface #if defined(Marc4DAMASK) || defined(Abaqus) - integer, parameter :: & - FILEUNIT = 222 - integer :: j - character(len=65536) :: tag, line - integer, allocatable, dimension(:) :: chunkPos + integer, parameter :: & + FILEUNIT = 222 + integer :: j + character(len=65536) :: tag, line + integer, allocatable, dimension(:) :: chunkPos #endif - write(6,'(/,a)') ' <<<+- FEsolving init -+>>>' + write(6,'(/,a)') ' <<<+- FEsolving init -+>>>' - modelName = getSolverJobName() + modelName = getSolverJobName() #if defined(Grid) || defined(FEM) - restartInc = interface_RestartInc + restartInc = interface_RestartInc - if(restartInc < 0) then - call IO_warning(warning_ID=34) - restartInc = 0 - endif - restartRead = restartInc > 0 ! only read in if "true" restart requested + if(restartInc < 0) then + call IO_warning(warning_ID=34) + restartInc = 0 + endif + restartRead = restartInc > 0 ! only read in if "true" restart requested #else - call IO_open_inputFile(FILEUNIT,modelName) - rewind(FILEUNIT) - do - read (FILEUNIT,'(a1024)',END=100) line - chunkPos = IO_stringPos(line) - tag = IO_lc(IO_stringValue(line,chunkPos,1)) ! extract key - select case(tag) - case ('solver') - read (FILEUNIT,'(a1024)',END=100) line ! next line - chunkPos = IO_stringPos(line) - symmetricSolver = (IO_intValue(line,chunkPos,2) /= 1) - case ('restart') - read (FILEUNIT,'(a1024)',END=100) line ! next line - chunkPos = IO_stringPos(line) - restartWrite = iand(IO_intValue(line,chunkPos,1),1) > 0 - restartRead = iand(IO_intValue(line,chunkPos,1),2) > 0 - case ('*restart') - do j=2,chunkPos(1) - restartWrite = (IO_lc(IO_StringValue(line,chunkPos,j)) == 'write') .or. restartWrite - restartRead = (IO_lc(IO_StringValue(line,chunkPos,j)) == 'read') .or. restartRead - enddo - if(restartWrite) then - do j=2,chunkPos(1) - restartWrite = (IO_lc(IO_StringValue(line,chunkPos,j)) /= 'frequency=0') .and. restartWrite - enddo - endif - end select - enddo - 100 close(FILEUNIT) + call IO_open_inputFile(FILEUNIT,modelName) + rewind(FILEUNIT) + do + read (FILEUNIT,'(a1024)',END=100) line + chunkPos = IO_stringPos(line) + tag = IO_lc(IO_stringValue(line,chunkPos,1)) ! extract key + select case(tag) + case ('solver') + read (FILEUNIT,'(a1024)',END=100) line ! next line + chunkPos = IO_stringPos(line) + symmetricSolver = (IO_intValue(line,chunkPos,2) /= 1) + case ('restart') + read (FILEUNIT,'(a1024)',END=100) line ! next line + chunkPos = IO_stringPos(line) + restartWrite = iand(IO_intValue(line,chunkPos,1),1) > 0 + restartRead = iand(IO_intValue(line,chunkPos,1),2) > 0 + case ('*restart') + do j=2,chunkPos(1) + restartWrite = (IO_lc(IO_StringValue(line,chunkPos,j)) == 'write') .or. restartWrite + restartRead = (IO_lc(IO_StringValue(line,chunkPos,j)) == 'read') .or. restartRead + enddo + if(restartWrite) then + do j=2,chunkPos(1) + restartWrite = (IO_lc(IO_StringValue(line,chunkPos,j)) /= 'frequency=0') .and. restartWrite + enddo + endif + end select + enddo + 100 close(FILEUNIT) - if (restartRead) then + if (restartRead) then #ifdef Marc4DAMASK - call IO_open_logFile(FILEUNIT) - rewind(FILEUNIT) - do - read (FILEUNIT,'(a1024)',END=200) line - chunkPos = IO_stringPos(line) - if ( IO_lc(IO_stringValue(line,chunkPos,1)) == 'restart' & - .and. IO_lc(IO_stringValue(line,chunkPos,2)) == 'file' & - .and. IO_lc(IO_stringValue(line,chunkPos,3)) == 'job' & - .and. IO_lc(IO_stringValue(line,chunkPos,4)) == 'id' ) & - modelName = IO_StringValue(line,chunkPos,6) - enddo -#else ! QUESTION: is this meaningful for the spectral/FEM case? - call IO_open_inputFile(FILEUNIT,modelName) - rewind(FILEUNIT) - do - read (FILEUNIT,'(a1024)',END=200) line - chunkPos = IO_stringPos(line) - if (IO_lc(IO_stringValue(line,chunkPos,1))=='*heading') then - read (FILEUNIT,'(a1024)',END=200) line - chunkPos = IO_stringPos(line) - modelName = IO_StringValue(line,chunkPos,1) - endif - enddo + call IO_open_logFile(FILEUNIT) + rewind(FILEUNIT) + do + read (FILEUNIT,'(a1024)',END=200) line + chunkPos = IO_stringPos(line) + if ( IO_lc(IO_stringValue(line,chunkPos,1)) == 'restart' & + .and. IO_lc(IO_stringValue(line,chunkPos,2)) == 'file' & + .and. IO_lc(IO_stringValue(line,chunkPos,3)) == 'job' & + .and. IO_lc(IO_stringValue(line,chunkPos,4)) == 'id' ) & + modelName = IO_StringValue(line,chunkPos,6) + enddo +#else + call IO_open_inputFile(FILEUNIT,modelName) + rewind(FILEUNIT) + do + read (FILEUNIT,'(a1024)',END=200) line + chunkPos = IO_stringPos(line) + if (IO_lc(IO_stringValue(line,chunkPos,1))=='*heading') then + read (FILEUNIT,'(a1024)',END=200) line + chunkPos = IO_stringPos(line) + modelName = IO_StringValue(line,chunkPos,1) + endif + enddo #endif - 200 close(FILEUNIT) - endif + 200 close(FILEUNIT) + endif #endif - if (iand(debug_level(debug_FEsolving),debug_levelBasic) /= 0) then - write(6,'(a21,l1)') ' restart writing: ', restartWrite - write(6,'(a21,l1)') ' restart reading: ', restartRead - if (restartRead) write(6,'(a,/)') ' restart Job: '//trim(modelName) - endif + if (iand(debug_level(debug_FEsolving),debug_levelBasic) /= 0) then + write(6,'(a21,l1)') ' restart writing: ', restartWrite + write(6,'(a21,l1)') ' restart reading: ', restartRead + if (restartRead) write(6,'(a,/)') ' restart Job: '//trim(modelName) + endif end subroutine FE_init diff --git a/src/element.f90 b/src/element.f90 index a7ebb9f5b..342ee3465 100644 --- a/src/element.f90 +++ b/src/element.f90 @@ -4,6 +4,7 @@ !-------------------------------------------------------------------------------------------------- module element use prec + use IO implicit none private @@ -805,8 +806,6 @@ module element contains subroutine tElement_init(self,elemType) - use IO, only: & - IO_error class(tElement) :: self integer, intent(in) :: elemType diff --git a/src/numerics.f90 b/src/numerics.f90 index 027296b92..0f63f87ab 100644 --- a/src/numerics.f90 +++ b/src/numerics.f90 @@ -5,6 +5,13 @@ !-------------------------------------------------------------------------------------------------- module numerics use prec + use IO + +#ifdef PETSc +#include + use petscsys +#endif +!$ use OMP_LIB implicit none private @@ -113,21 +120,6 @@ contains ! a sanity check !-------------------------------------------------------------------------------------------------- subroutine numerics_init - use IO, only: & - IO_read_ASCII, & - IO_error, & - IO_isBlank, & - IO_stringPos, & - IO_stringValue, & - IO_lc, & - IO_floatValue, & - IO_intValue, & - IO_warning -#ifdef PETSc -#include - use petscsys -#endif -!$ use OMP_LIB, only: omp_set_num_threads !$ integer :: gotDAMASK_NUM_THREADS = 1 integer :: i,j, ierr integer, allocatable, dimension(:) :: chunkPos diff --git a/src/thermal_isothermal.f90 b/src/thermal_isothermal.f90 index c15dfa8f1..0271813e2 100644 --- a/src/thermal_isothermal.f90 +++ b/src/thermal_isothermal.f90 @@ -3,48 +3,45 @@ !> @brief material subroutine for isothermal temperature field !-------------------------------------------------------------------------------------------------- module thermal_isothermal + use prec + use config + use material - implicit none - private - - public :: & - thermal_isothermal_init + implicit none + private + + public :: & + thermal_isothermal_init contains !-------------------------------------------------------------------------------------------------- !> @brief allocates all neccessary fields, reads information from material configuration file !-------------------------------------------------------------------------------------------------- -subroutine thermal_isothermal_init() - use prec, only: & - pReal - use config, only: & - material_Nhomogenization - use material +subroutine thermal_isothermal_init - integer :: & - homog, & - NofMyHomog + integer :: & + homog, & + NofMyHomog - write(6,'(/,a)') ' <<<+- thermal_'//THERMAL_isothermal_label//' init -+>>>' + write(6,'(/,a)') ' <<<+- thermal_'//THERMAL_isothermal_label//' init -+>>>' - initializeInstances: do homog = 1, material_Nhomogenization - - if (thermal_type(homog) /= THERMAL_isothermal_ID) cycle - NofMyHomog = count(material_homogenizationAt == homog) - thermalState(homog)%sizeState = 0 - thermalState(homog)%sizePostResults = 0 - allocate(thermalState(homog)%state0 (0,NofMyHomog), source=0.0_pReal) - allocate(thermalState(homog)%subState0(0,NofMyHomog), source=0.0_pReal) - allocate(thermalState(homog)%state (0,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) - - enddo initializeInstances + initializeInstances: do homog = 1, material_Nhomogenization + + if (thermal_type(homog) /= THERMAL_isothermal_ID) cycle + NofMyHomog = count(material_homogenizationAt == homog) + thermalState(homog)%sizeState = 0 + thermalState(homog)%sizePostResults = 0 + allocate(thermalState(homog)%state0 (0,NofMyHomog), source=0.0_pReal) + allocate(thermalState(homog)%subState0(0,NofMyHomog), source=0.0_pReal) + allocate(thermalState(homog)%state (0,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) + enddo initializeInstances end subroutine thermal_isothermal_init