diff --git a/src/CPFEM.f90 b/src/CPFEM.f90 index 240688a8c..713aab5d7 100644 --- a/src/CPFEM.f90 +++ b/src/CPFEM.f90 @@ -180,10 +180,10 @@ subroutine CPFEM_general(mode, ffn, ffn1, temperature_inp, dt, elFE, ip, cauchyS if (iand(mode, CPFEM_AGERESULTS) /= 0_pInt) call CPFEM_forward chosenThermal1: select case (thermal_type(material_homogenizationAt(elCP))) - case (THERMAL_conduction_ID) chosenThermal1 - temperature(material_homogenizationAt(elCP))%p(material_homogenizationMemberAt(ip,elCP)) = & - temperature_inp - end select chosenThermal1 + !case (THERMAL_conduction_ID) chosenThermal1 + ! temperature(material_homogenizationAt(elCP))%p(material_homogenizationMemberAt(ip,elCP)) = & + ! temperature_inp + end select chosenThermal1 homogenization_F0(1:3,1:3,ma) = ffn homogenization_F(1:3,1:3,ma) = ffn1 diff --git a/src/DAMASK_marc.f90 b/src/DAMASK_marc.f90 index 9464595b1..360c911ff 100644 --- a/src/DAMASK_marc.f90 +++ b/src/DAMASK_marc.f90 @@ -351,7 +351,7 @@ end subroutine hypela2 !-------------------------------------------------------------------------------------------------- subroutine flux(f,ts,n,time) use prec - use thermal_conduction + use homogenization use discretization_marc implicit none diff --git a/src/commercialFEM_fileList.f90 b/src/commercialFEM_fileList.f90 index 963871c7a..a20fb1cee 100644 --- a/src/commercialFEM_fileList.f90 +++ b/src/commercialFEM_fileList.f90 @@ -42,8 +42,6 @@ #include "source_damage_anisoDuctile.f90" #include "kinematics_cleavage_opening.f90" #include "kinematics_slipplane_opening.f90" -#include "thermal_isothermal.f90" -#include "thermal_conduction.f90" #include "damage_none.f90" #include "damage_nonlocal.f90" #include "homogenization.f90" diff --git a/src/constitutive.f90 b/src/constitutive.f90 index 16f3d3b59..149d09a49 100644 --- a/src/constitutive.f90 +++ b/src/constitutive.f90 @@ -120,10 +120,6 @@ module constitutive integer, intent(in) :: ph, me end subroutine mech_initializeRestorationPoints - module subroutine constitutive_thermal_initializeRestorationPoints(ph,me) - integer, intent(in) :: ph, me - end subroutine constitutive_thermal_initializeRestorationPoints - module subroutine mech_windForward(ph,me) integer, intent(in) :: ph, me @@ -137,8 +133,8 @@ module constitutive end subroutine thermal_forward - module subroutine mech_restore(ip,el,includeL) - integer, intent(in) :: ip, el + module subroutine mech_restore(ce,includeL) + integer, intent(in) :: ce logical, intent(in) :: includeL end subroutine mech_restore @@ -204,11 +200,11 @@ module constitutive integer, intent(in) :: co, ip, el end subroutine constitutive_mech_setF - module subroutine constitutive_thermal_setT(T,co,ce) - real(pReal), intent(in) :: T - integer, intent(in) :: co, ce - end subroutine constitutive_thermal_setT - + module subroutine constitutive_thermal_setField(T,dot_T, co,ce) + real(pReal), intent(in) :: T, dot_T + integer, intent(in) :: ce, co + end subroutine constitutive_thermal_setField + module subroutine constitutive_damage_set_phi(phi,co,ce) real(pReal), intent(in) :: phi integer, intent(in) :: co, ce @@ -245,9 +241,6 @@ module constitutive end function constitutive_homogenizedC - - - module subroutine constitutive_damage_getRateAndItsTangents(phiDot, dPhiDot_dPhi, phi, ip, el) integer, intent(in) :: & ip, & !< integration point number @@ -317,11 +310,9 @@ module constitutive dLd_dTstar !< derivative of Ld with respect to Tstar (4th-order tensor) end subroutine kinematics_slipplane_opening_LiAndItsTangent - module subroutine kinematics_thermal_expansion_LiAndItsTangent(Li, dLi_dTstar, co, ip, el) - integer, intent(in) :: & - co, & !< grain number - ip, & !< integration point number - el !< element number + module subroutine kinematics_thermal_expansion_LiAndItsTangent(Li, dLi_dTstar, ph,me) + integer, intent(in) :: ph, me + !< element number real(pReal), intent(out), dimension(3,3) :: & Li !< thermal velocity gradient real(pReal), intent(out), dimension(3,3,3,3) :: & @@ -362,14 +353,13 @@ module constitutive constitutive_restartWrite, & constitutive_restartRead, & integrateDamageState, & - constitutive_thermal_setT, & + constitutive_thermal_setField, & constitutive_damage_set_phi, & constitutive_damage_get_phi, & constitutive_mech_getP, & constitutive_mech_setF, & constitutive_mech_getF, & constitutive_initializeRestorationPoints, & - constitutive_thermal_initializeRestorationPoints, & constitutive_windForward, & KINEMATICS_UNDEFINED_ID ,& KINEMATICS_CLEAVAGE_OPENING_ID, & @@ -497,26 +487,24 @@ end subroutine constitutive_allocateState !-------------------------------------------------------------------------------------------------- !> @brief Restore data after homog cutback. !-------------------------------------------------------------------------------------------------- -subroutine constitutive_restore(ip,el,includeL) +subroutine constitutive_restore(ce,includeL) logical, intent(in) :: includeL - integer, intent(in) :: & - ip, & !< integration point number - el !< element number + integer, intent(in) :: ce integer :: & co, & !< constituent number so - do co = 1,homogenization_Nconstituents(material_homogenizationAt(el)) - do so = 1, phase_Nsources(material_phaseAt(co,el)) - damageState(material_phaseAt(co,el))%p(so)%state( :,material_phasememberAt(co,ip,el)) = & - damageState(material_phaseAt(co,el))%p(so)%partitionedState0(:,material_phasememberAt(co,ip,el)) + do co = 1,homogenization_Nconstituents(material_homogenizationAt2(ce)) + do so = 1, phase_Nsources(material_phaseAt2(co,ce)) + damageState(material_phaseAt2(co,ce))%p(so)%state( :,material_phasememberAt2(co,ce)) = & + damageState(material_phaseAt2(co,ce))%p(so)%partitionedState0(:,material_phasememberAt2(co,ce)) enddo enddo - call mech_restore(ip,el,includeL) + call mech_restore(ce,includeL) end subroutine constitutive_restore @@ -638,9 +626,6 @@ subroutine crystallite_init() do so = 1, phase_Nsources(ph) allocate(damageState(ph)%p(so)%subState0,source=damageState(ph)%p(so)%state0) ! ToDo: hack enddo - do so = 1, thermal_Nsources(ph) - allocate(thermalState(ph)%p(so)%subState0,source=thermalState(ph)%p(so)%state0) ! ToDo: hack - enddo enddo print'(a42,1x,i10)', ' # of elements: ', eMax @@ -766,9 +751,6 @@ function crystallite_push33ToRef(co,ip,el, tensor33) end function crystallite_push33ToRef - - - !-------------------------------------------------------------------------------------------------- !> @brief determines whether a point is converged !-------------------------------------------------------------------------------------------------- diff --git a/src/constitutive_mech.f90 b/src/constitutive_mech.f90 index cd1dedb61..56e03ac6f 100644 --- a/src/constitutive_mech.f90 +++ b/src/constitutive_mech.f90 @@ -1586,9 +1586,6 @@ module function crystallite_stress(dt,co,ip,el) result(converged_) do so = 1, phase_Nsources(ph) damageState(ph)%p(so)%subState0(:,me) = damageState(ph)%p(so)%partitionedState0(:,me) enddo - do so = 1, thermal_Nsources(ph) - thermalState(ph)%p(so)%subState0(:,me) = thermalState(ph)%p(so)%partitionedState0(:,me) - enddo subFp0 = constitutive_mech_partitionedFp0(ph)%data(1:3,1:3,me) subFi0 = constitutive_mech_partitionedFi0(ph)%data(1:3,1:3,me) subF0 = constitutive_mech_partitionedF0(ph)%data(1:3,1:3,me) @@ -1656,11 +1653,9 @@ end function crystallite_stress !-------------------------------------------------------------------------------------------------- !> @brief Restore data after homog cutback. !-------------------------------------------------------------------------------------------------- -module subroutine mech_restore(ip,el,includeL) +module subroutine mech_restore(ce,includeL) - integer, intent(in) :: & - ip, & !< integration point number - el !< element number + integer, intent(in) :: ce logical, intent(in) :: & includeL !< protect agains fake cutback @@ -1668,9 +1663,9 @@ module subroutine mech_restore(ip,el,includeL) co, ph, me - do co = 1,homogenization_Nconstituents(material_homogenizationAt(el)) - ph = material_phaseAt(co,el) - me = material_phaseMemberAt(co,ip,el) + do co = 1,homogenization_Nconstituents(material_homogenizationAt2(ce)) + ph = material_phaseAt2(co,ce) + me = material_phaseMemberAt2(co,ce) if (includeL) then constitutive_mech_Lp(ph)%data(1:3,1:3,me) = constitutive_mech_partitionedLp0(ph)%data(1:3,1:3,me) constitutive_mech_Li(ph)%data(1:3,1:3,me) = constitutive_mech_partitionedLi0(ph)%data(1:3,1:3,me) @@ -1961,7 +1956,7 @@ subroutine constitutive_LiAndItsTangents(Li, dLi_dS, dLi_dFi, & detFi integer :: & k, i, j, & - instance, of + instance, of, me, ph Li = 0.0_pReal dLi_dS = 0.0_pReal @@ -1987,7 +1982,9 @@ subroutine constitutive_LiAndItsTangents(Li, dLi_dS, dLi_dFi, & case (KINEMATICS_slipplane_opening_ID) kinematicsType call kinematics_slipplane_opening_LiAndItsTangent(my_Li, my_dLi_dS, S, co, ip, el) case (KINEMATICS_thermal_expansion_ID) kinematicsType - call kinematics_thermal_expansion_LiAndItsTangent(my_Li, my_dLi_dS, co, ip, el) + me = material_phaseMemberAt(co,ip,el) + ph = material_phaseAt(co,el) + call kinematics_thermal_expansion_LiAndItsTangent(my_Li, my_dLi_dS, ph,me) case default kinematicsType my_Li = 0.0_pReal my_dLi_dS = 0.0_pReal diff --git a/src/constitutive_thermal.f90 b/src/constitutive_thermal.f90 index 2746ff889..57b8a3117 100644 --- a/src/constitutive_thermal.f90 +++ b/src/constitutive_thermal.f90 @@ -10,9 +10,9 @@ submodule(constitutive) constitutive_thermal end enum type :: tDataContainer - real(pReal), dimension(:), allocatable :: T + real(pReal), dimension(:), allocatable :: T, dot_T end type tDataContainer - integer(kind(THERMAL_UNDEFINED_ID)), dimension(:,:), allocatable :: & + integer(kind(THERMAL_UNDEFINED_ID)), dimension(:,:), allocatable :: & thermal_source type(tDataContainer), dimension(:), allocatable :: current @@ -94,6 +94,7 @@ module subroutine thermal_init(phases) Nconstituents = count(material_phaseAt == ph) * discretization_nIPs allocate(current(ph)%T(Nconstituents),source=300.0_pReal) + allocate(current(ph)%dot_T(Nconstituents),source=0.0_pReal) phase => phases%get(ph) if(phase%contains('thermal')) then thermal => phase%get('thermal') @@ -115,8 +116,8 @@ module subroutine thermal_init(phases) PhaseLoop2:do ph = 1,phases%length do so = 1,thermal_Nsources(ph) - thermalState(ph)%p(so)%partitionedState0 = thermalState(ph)%p(so)%state0 - thermalState(ph)%p(so)%state = thermalState(ph)%p(so)%partitionedState0 + deallocate(thermalState(ph)%p(so)%partitionedState0) + thermalState(ph)%p(so)%state = thermalState(ph)%p(so)%state0 enddo thermal_source_maxSizeDotState = max(thermal_source_maxSizeDotState, & @@ -147,13 +148,11 @@ module subroutine constitutive_thermal_getRate(TDot, ip, el) integer :: & ph, & homog, & - instance, & me, & so, & co homog = material_homogenizationAt(el) - instance = thermal_typeInstance(homog) TDot = 0.0_pReal do co = 1, homogenization_Nconstituents(homog) @@ -211,9 +210,6 @@ module function thermal_stress(Delta_t,ph,me) result(converged_) integer :: so - do so = 1, thermal_Nsources(ph) - thermalState(ph)%p(so)%state(:,me) = thermalState(ph)%p(so)%subState0(:,me) - enddo converged_ = .not. integrateThermalState(Delta_t,ph,me) end function thermal_stress @@ -238,28 +234,13 @@ function integrateThermalState(Delta_t, ph,me) result(broken) do so = 1, thermal_Nsources(ph) sizeDotState = thermalState(ph)%p(so)%sizeDotState - thermalState(ph)%p(so)%state(1:sizeDotState,me) = thermalState(ph)%p(so)%subState0(1:sizeDotState,me) & + thermalState(ph)%p(so)%state(1:sizeDotState,me) = thermalState(ph)%p(so)%state0(1:sizeDotState,me) & + thermalState(ph)%p(so)%dotState(1:sizeDotState,me) * Delta_t enddo end function integrateThermalState -module subroutine constitutive_thermal_initializeRestorationPoints(ph,me) - - integer, intent(in) :: ph, me - - integer :: so - - - do so = 1, size(thermalState(ph)%p) - thermalState(ph)%p(so)%partitionedState0(:,me) = thermalState(ph)%p(so)%state0(:,me) - enddo - -end subroutine constitutive_thermal_initializeRestorationPoints - - - module subroutine thermal_forward() integer :: ph, so @@ -291,15 +272,16 @@ end function thermal_T !---------------------------------------------------------------------------------------------- !< @brief Set temperature !---------------------------------------------------------------------------------------------- -module subroutine constitutive_thermal_setT(T,co,ce) +module subroutine constitutive_thermal_setField(T,dot_T, co,ce) - real(pReal), intent(in) :: T + real(pReal), intent(in) :: T, dot_T integer, intent(in) :: ce, co current(material_phaseAt2(co,ce))%T(material_phaseMemberAt2(co,ce)) = T + current(material_phaseAt2(co,ce))%dot_T(material_phaseMemberAt2(co,ce)) = dot_T -end subroutine constitutive_thermal_setT +end subroutine constitutive_thermal_setField diff --git a/src/grid/DAMASK_grid.f90 b/src/grid/DAMASK_grid.f90 index 02d7b4cc3..60c386c75 100644 --- a/src/grid/DAMASK_grid.f90 +++ b/src/grid/DAMASK_grid.f90 @@ -425,7 +425,7 @@ program DAMASK_grid guess = .true. ! start guessing after first converged (sub)inc if (worldrank == 0) then write(statUnit,*) totalIncsCounter, time, cutBackLevel, & - solres%converged, solres%iterationsNeeded + solres(1)%converged, solres(1)%iterationsNeeded flush(statUnit) endif elseif (cutBackLevel < maxCutBack) then ! further cutbacking tolerated? diff --git a/src/grid/grid_thermal_spectral.f90 b/src/grid/grid_thermal_spectral.f90 index d164f11e7..866914f5b 100644 --- a/src/grid/grid_thermal_spectral.f90 +++ b/src/grid/grid_thermal_spectral.f90 @@ -15,7 +15,6 @@ module grid_thermal_spectral use IO use spectral_utilities use discretization_grid - use thermal_conduction use homogenization use YAML_types use config @@ -132,7 +131,7 @@ subroutine grid_thermal_spectral_init ce = 0 do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) ce = ce + 1 - T_current(i,j,k) = temperature(material_homogenizationAt(ce))%p(material_homogenizationMemberAt(1,ce)) + T_current(i,j,k) = homogenization_thermal_T(ce) T_lastInc(i,j,k) = T_current(i,j,k) T_stagInc(i,j,k) = T_current(i,j,k) enddo; enddo; enddo @@ -188,10 +187,9 @@ function grid_thermal_spectral_solution(timeinc) result(solution) ce = 0 do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) ce = ce + 1 - call thermal_conduction_putTemperatureAndItsRate(T_current(i,j,k), & + call homogenization_thermal_setField(T_current(i,j,k), & (T_current(i,j,k)-T_lastInc(i,j,k))/params%timeinc, & - 1,ce) - homogenization_T(ce) = T_current(i,j,k) + ce) enddo; enddo; enddo call VecMin(solution_vec,devNull,T_min,ierr); CHKERRQ(ierr) @@ -229,11 +227,9 @@ subroutine grid_thermal_spectral_forward(cutBack) call DMDAVecRestoreArrayF90(dm_local,solution_vec,x_scal,ierr); CHKERRQ(ierr) do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) ce = ce + 1 - call thermal_conduction_putTemperatureAndItsRate(T_current(i,j,k), & - (T_current(i,j,k) - & - T_lastInc(i,j,k))/params%timeinc, & - 1,ce) - homogenization_T(ce) = T_current(i,j,k) + call homogenization_thermal_setField(T_current(i,j,k), & + (T_current(i,j,k)-T_lastInc(i,j,k))/params%timeinc, & + ce) enddo; enddo; enddo else T_lastInc = T_current @@ -283,8 +279,8 @@ subroutine formResidual(in,x_scal,f_scal,dummy,ierr) ce = ce + 1 call thermal_conduction_getSource(Tdot, 1,ce) scalarField_real(i,j,k) = params%timeinc*(scalarField_real(i,j,k) + Tdot) & - + thermal_conduction_getMassDensity (1,ce)* & - thermal_conduction_getSpecificHeat(1,ce)*(T_lastInc(i,j,k) - & + + thermal_conduction_getMassDensity (ce)* & + thermal_conduction_getSpecificHeat(ce)*(T_lastInc(i,j,k) - & T_current(i,j,k))& + mu_ref*T_current(i,j,k) enddo; enddo; enddo @@ -315,7 +311,7 @@ subroutine updateReference do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) ce = ce + 1 K_ref = K_ref + thermal_conduction_getConductivity(1,ce) - mu_ref = mu_ref + thermal_conduction_getMassDensity(1,ce)* thermal_conduction_getSpecificHeat(1,ce) + mu_ref = mu_ref + thermal_conduction_getMassDensity(ce)* thermal_conduction_getSpecificHeat(ce) enddo; enddo; enddo K_ref = K_ref*wgt call MPI_Allreduce(MPI_IN_PLACE,K_ref,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) diff --git a/src/homogenization.f90 b/src/homogenization.f90 index ed9b7fead..64c2aade7 100644 --- a/src/homogenization.f90 +++ b/src/homogenization.f90 @@ -12,8 +12,6 @@ module homogenization use material use constitutive use discretization - use thermal_isothermal - use thermal_conduction use damage_none use damage_nonlocal use HDF5_utilities @@ -28,8 +26,6 @@ module homogenization !-------------------------------------------------------------------------------------------------- ! General variables for the homogenization at a material point real(pReal), dimension(:), allocatable, public :: & - homogenization_T, & - homogenization_dot_T, & homogenization_phi, & homogenization_dot_phi real(pReal), dimension(:,:,:), allocatable, public :: & @@ -75,8 +71,7 @@ module homogenization el !< element number end subroutine mech_partition - module subroutine thermal_partition(T,ce) - real(pReal), intent(in) :: T + module subroutine thermal_partition(ce) integer, intent(in) :: ce end subroutine thermal_partition @@ -112,11 +107,59 @@ module homogenization logical, dimension(2) :: doneAndHappy end function mech_updateState + + module function thermal_conduction_getConductivity(ip,el) result(K) + + integer, intent(in) :: & + ip, & !< integration point number + el !< element number + real(pReal), dimension(3,3) :: K + + end function thermal_conduction_getConductivity + + module function thermal_conduction_getSpecificHeat(ce) result(c_P) + integer, intent(in) :: ce + real(pReal) :: c_P + end function thermal_conduction_getSpecificHeat + + module function thermal_conduction_getMassDensity(ce) result(rho) + integer, intent(in) :: ce + real(pReal) :: rho + end function thermal_conduction_getMassDensity + + module subroutine homogenization_thermal_setField(T,dot_T, ce) + integer, intent(in) :: ce + real(pReal), intent(in) :: T, dot_T + end subroutine homogenization_thermal_setField + + module subroutine thermal_conduction_results(ho,group) + integer, intent(in) :: ho + character(len=*), intent(in) :: group + end subroutine thermal_conduction_results + + module function homogenization_thermal_T(ce) result(T) + integer, intent(in) :: ce + real(pReal) :: T + end function homogenization_thermal_T + + module subroutine thermal_conduction_getSource(Tdot, ip,el) + integer, intent(in) :: & + ip, & !< integration point number + el !< element number + real(pReal), intent(out) :: Tdot + end subroutine thermal_conduction_getSource + end interface public :: & homogenization_init, & materialpoint_stressAndItsTangent, & + thermal_conduction_getSpecificHeat, & + thermal_conduction_getConductivity, & + thermal_conduction_getMassDensity, & + thermal_conduction_getSource, & + homogenization_thermal_setfield, & + homogenization_thermal_T, & homogenization_forward, & homogenization_results, & homogenization_restartRead, & @@ -154,9 +197,6 @@ subroutine homogenization_init() call thermal_init() call damage_init() - if (any(thermal_type == THERMAL_isothermal_ID)) call thermal_isothermal_init(homogenization_T) - if (any(thermal_type == THERMAL_conduction_ID)) call thermal_conduction_init(homogenization_T) - if (any(damage_type == DAMAGE_none_ID)) call damage_none_init if (any(damage_type == DAMAGE_nonlocal_ID)) call damage_nonlocal_init @@ -190,9 +230,9 @@ subroutine materialpoint_stressAndItsTangent(dt,FEsolving_execIP,FEsolving_execE ho = material_homogenizationAt(el) myNgrains = homogenization_Nconstituents(ho) do ip = FEsolving_execIP(1),FEsolving_execIP(2) - me = material_homogenizationMemberAt(ip,el) -!-------------------------------------------------------------------------------------------------- -! initialize restoration points + ce = (el-1)*discretization_nIPs + ip + me = material_homogenizationMemberAt2(ce) + call constitutive_initializeRestorationPoints(ip,el) subFrac = 0.0_pReal @@ -226,7 +266,7 @@ subroutine materialpoint_stressAndItsTangent(dt,FEsolving_execIP,FEsolving_execE else ! cutback makes sense subStep = num%subStepSizeHomog * subStep ! crystallite had severe trouble, so do a significant cutback - call constitutive_restore(ip,el,subStep < 1.0_pReal) + call constitutive_restore(ce,subStep < 1.0_pReal) if(homogState(ho)%sizeState > 0) homogState(ho)%State(:,me) = homogState(ho)%subState0(:,me) if(damageState_h(ho)%sizeState > 0) damageState_h(ho)%State(:,me) = damageState_h(ho)%subState0(:,me) @@ -243,7 +283,6 @@ subroutine materialpoint_stressAndItsTangent(dt,FEsolving_execIP,FEsolving_execE ! deformation partitioning if (.not. doneAndHappy(1)) then - ce = (el-1)*discretization_nIPs + ip call mech_partition( homogenization_F0(1:3,1:3,ce) & + (homogenization_F(1:3,1:3,ce)-homogenization_F0(1:3,1:3,ce))*(subStep+subFrac), & ip,el) @@ -255,7 +294,6 @@ subroutine materialpoint_stressAndItsTangent(dt,FEsolving_execIP,FEsolving_execE if (.not. converged) then doneAndHappy = [.true.,.false.] else - ce = (el-1)*discretization_nIPs + ip doneAndHappy = mech_updateState(dt*subStep, & homogenization_F0(1:3,1:3,ce) & + (homogenization_F(1:3,1:3,ce)-homogenization_F0(1:3,1:3,ce)) & @@ -278,10 +316,9 @@ subroutine materialpoint_stressAndItsTangent(dt,FEsolving_execIP,FEsolving_execE ho = material_homogenizationAt(el) do ip = FEsolving_execIP(1),FEsolving_execIP(2) ce = (el-1)*discretization_nIPs + ip - call thermal_partition(homogenization_T(ce),ce) + call thermal_partition(ce) do co = 1, homogenization_Nconstituents(ho) ph = material_phaseAt(co,el) - call constitutive_thermal_initializeRestorationPoints(ph,material_phaseMemberAt(co,ip,el)) if (.not. thermal_stress(dt,ph,material_phaseMemberAt(co,ip,el))) then if (.not. terminallyIll) & ! so first signals terminally ill... print*, ' Integration point ', ip,' at element ', el, ' terminally ill' diff --git a/src/homogenization_thermal.f90 b/src/homogenization_thermal.f90 index de91a98d6..6a12730ac 100644 --- a/src/homogenization_thermal.f90 +++ b/src/homogenization_thermal.f90 @@ -3,6 +3,22 @@ !-------------------------------------------------------------------------------------------------- submodule(homogenization) homogenization_thermal + use lattice + + type :: tDataContainer + real(pReal), dimension(:), allocatable :: T, dot_T + end type tDataContainer + + type(tDataContainer), dimension(:), allocatable :: current + + type :: tParameters + character(len=pStringLen), allocatable, dimension(:) :: & + output + end type tParameters + + type(tparameters), dimension(:), allocatable :: & + param + contains @@ -11,11 +27,37 @@ contains !-------------------------------------------------------------------------------------------------- module subroutine thermal_init() + class(tNode), pointer :: & + configHomogenizations, & + configHomogenization, & + configHomogenizationThermal + integer :: ho + print'(/,a)', ' <<<+- homogenization_thermal init -+>>>' - allocate(homogenization_T(discretization_nIPs*discretization_Nelems)) - allocate(homogenization_dot_T(discretization_nIPs*discretization_Nelems)) + + configHomogenizations => config_material%get('homogenization') + allocate(param(configHomogenizations%length)) + allocate(current(configHomogenizations%length)) + + do ho = 1, configHomogenizations%length + allocate(current(ho)%T(count(material_homogenizationAt2==ho)), source=thermal_initialT(ho)) + allocate(current(ho)%dot_T(count(material_homogenizationAt2==ho)), source=0.0_pReal) + configHomogenization => configHomogenizations%get(ho) + associate(prm => param(ho)) + if (configHomogenization%contains('thermal')) then + configHomogenizationThermal => configHomogenization%get('thermal') +#if defined (__GFORTRAN__) + prm%output = output_asStrings(configHomogenizationThermal) +#else + prm%output = configHomogenizationThermal%get_asStrings('output',defaultVal=emptyStringArray) +#endif + else + prm%output = emptyStringArray + endif + end associate + enddo end subroutine thermal_init @@ -23,15 +65,18 @@ end subroutine thermal_init !-------------------------------------------------------------------------------------------------- !> @brief Partition temperature onto the individual constituents. !-------------------------------------------------------------------------------------------------- -module subroutine thermal_partition(T,ce) +module subroutine thermal_partition(ce) - real(pReal), intent(in) :: T integer, intent(in) :: ce + real(pReal) :: T, dot_T integer :: co + + T = current(material_homogenizationAt2(ce))%T(material_homogenizationMemberAt2(ce)) + dot_T = current(material_homogenizationAt2(ce))%dot_T(material_homogenizationMemberAt2(ce)) do co = 1, homogenization_Nconstituents(material_homogenizationAt2(ce)) - call constitutive_thermal_setT(T,co,ce) + call constitutive_thermal_setField(T,dot_T,co,ce) enddo end subroutine thermal_partition @@ -44,9 +89,151 @@ module subroutine thermal_homogenize(ip,el) integer, intent(in) :: ip,el - call constitutive_thermal_getRate(homogenization_dot_T((el-1)*discretization_nIPs+ip), ip,el) + !call constitutive_thermal_getRate(homogenization_dot_T((el-1)*discretization_nIPs+ip), ip,el) end subroutine thermal_homogenize +!-------------------------------------------------------------------------------------------------- +!> @brief return homogenized thermal conductivity in reference configuration +!-------------------------------------------------------------------------------------------------- +module function thermal_conduction_getConductivity(ip,el) result(K) + + integer, intent(in) :: & + ip, & !< integration point number + el !< element number + real(pReal), dimension(3,3) :: K + + integer :: & + co + + + K = 0.0_pReal + + do co = 1, homogenization_Nconstituents(material_homogenizationAt(el)) + K = K + crystallite_push33ToRef(co,ip,el,lattice_K(:,:,material_phaseAt(co,el))) + enddo + + K = K / real(homogenization_Nconstituents(material_homogenizationAt(el)),pReal) + +end function thermal_conduction_getConductivity + + +!-------------------------------------------------------------------------------------------------- +!> @brief returns homogenized specific heat capacity +!-------------------------------------------------------------------------------------------------- +module function thermal_conduction_getSpecificHeat(ce) result(c_P) + + integer, intent(in) :: ce + real(pReal) :: c_P + + integer :: co + + + c_P = 0.0_pReal + + do co = 1, homogenization_Nconstituents(material_homogenizationAt2(ce)) + c_P = c_P + lattice_c_p(material_phaseAt2(co,ce)) + enddo + + c_P = c_P / real(homogenization_Nconstituents(material_homogenizationAt2(ce)),pReal) + +end function thermal_conduction_getSpecificHeat + + +!-------------------------------------------------------------------------------------------------- +!> @brief returns homogenized mass density +!-------------------------------------------------------------------------------------------------- +module function thermal_conduction_getMassDensity(ce) result(rho) + + integer, intent(in) :: ce + real(pReal) :: rho + + integer :: co + + + rho = 0.0_pReal + + do co = 1, homogenization_Nconstituents(material_homogenizationAt2(ce)) + rho = rho + lattice_rho(material_phaseAt2(co,ce)) + enddo + + rho = rho / real(homogenization_Nconstituents(material_homogenizationAt2(ce)),pReal) + +end function thermal_conduction_getMassDensity + + + +!-------------------------------------------------------------------------------------------------- +!> @brief Set thermal field and its rate (T and dot_T) +!-------------------------------------------------------------------------------------------------- +module subroutine homogenization_thermal_setField(T,dot_T, ce) + + integer, intent(in) :: ce + real(pReal), intent(in) :: T, dot_T + + + current(material_homogenizationAt2(ce))%T(material_homogenizationMemberAt2(ce)) = T + current(material_homogenizationAt2(ce))%dot_T(material_homogenizationMemberAt2(ce)) = dot_T + + +end subroutine homogenization_thermal_setField + + + +!-------------------------------------------------------------------------------------------------- +!> @brief writes results to HDF5 output file +!-------------------------------------------------------------------------------------------------- +module subroutine thermal_conduction_results(ho,group) + + integer, intent(in) :: ho + character(len=*), intent(in) :: group + + integer :: o + + associate(prm => param(ho)) + outputsLoop: do o = 1,size(prm%output) + select case(trim(prm%output(o))) + case('T') + call results_writeDataset(group,current(ho)%T,'T','temperature','K') + end select + enddo outputsLoop + end associate + +end subroutine thermal_conduction_results + + +module function homogenization_thermal_T(ce) result(T) + + integer, intent(in) :: ce + real(pReal) :: T + + T = current(material_homogenizationAt2(ce))%T(material_homogenizationMemberAt2(ce)) + +end function homogenization_thermal_T + + + +!-------------------------------------------------------------------------------------------------- +!> @brief return heat generation rate +!-------------------------------------------------------------------------------------------------- +module subroutine thermal_conduction_getSource(Tdot, ip,el) + + integer, intent(in) :: & + ip, & !< integration point number + el !< element number + real(pReal), intent(out) :: & + Tdot + + integer :: & + homog + + homog = material_homogenizationAt(el) + call constitutive_thermal_getRate(TDot, ip,el) + + Tdot = Tdot/real(homogenization_Nconstituents(homog),pReal) + +end subroutine thermal_conduction_getSource + + end submodule homogenization_thermal diff --git a/src/kinematics_thermal_expansion.f90 b/src/kinematics_thermal_expansion.f90 index 6d4a39632..6e4fa8263 100644 --- a/src/kinematics_thermal_expansion.f90 +++ b/src/kinematics_thermal_expansion.f90 @@ -84,12 +84,9 @@ end function kinematics_thermal_expansion_init !-------------------------------------------------------------------------------------------------- !> @brief constitutive equation for calculating the velocity gradient !-------------------------------------------------------------------------------------------------- -module subroutine kinematics_thermal_expansion_LiAndItsTangent(Li, dLi_dTstar, co, ip, el) +module subroutine kinematics_thermal_expansion_LiAndItsTangent(Li, dLi_dTstar, ph,me) - integer, intent(in) :: & - co, & !< grain number - ip, & !< integration point number - el !< element number + integer, intent(in) :: ph, me real(pReal), intent(out), dimension(3,3) :: & Li !< thermal velocity gradient real(pReal), intent(out), dimension(3,3,3,3) :: & @@ -98,16 +95,13 @@ module subroutine kinematics_thermal_expansion_LiAndItsTangent(Li, dLi_dTstar, c integer :: & phase, & homog - real(pReal) :: & - T, TDot + real(pReal) :: T, dot_T - phase = material_phaseAt(co,el) - homog = material_homogenizationAt(el) - T = temperature(homog)%p(material_homogenizationMemberAt(ip,el)) - TDot = temperatureRate(homog)%p(material_homogenizationMemberAt(ip,el)) + T = current(ph)%T(me) + dot_T = current(ph)%dot_T(me) - associate(prm => param(kinematics_thermal_expansion_instance(phase))) - Li = TDot * ( & + associate(prm => param(kinematics_thermal_expansion_instance(ph))) + Li = dot_T * ( & prm%A(1:3,1:3,1)*(T - prm%T_ref)**0 & ! constant coefficient + prm%A(1:3,1:3,2)*(T - prm%T_ref)**1 & ! linear coefficient + prm%A(1:3,1:3,3)*(T - prm%T_ref)**2 & ! quadratic coefficient diff --git a/src/material.f90 b/src/material.f90 index e165bd466..5f0f86fe5 100644 --- a/src/material.f90 +++ b/src/material.f90 @@ -71,9 +71,7 @@ module material material_orientation0 !< initial orientation of each grain,IP,element type(group_float), allocatable, dimension(:), public :: & - temperature, & !< temperature field - damage, & !< damage field - temperatureRate !< temperature change rate field + damage !< damage field public :: & material_init, & @@ -107,9 +105,7 @@ subroutine material_init(restart) allocate(homogState (size(material_name_homogenization))) allocate(damageState_h (size(material_name_homogenization))) - allocate(temperature (size(material_name_homogenization))) allocate(damage (size(material_name_homogenization))) - allocate(temperatureRate (size(material_name_homogenization))) if (.not. restart) then diff --git a/src/thermal_conduction.f90 b/src/thermal_conduction.f90 deleted file mode 100644 index 1711c5512..000000000 --- a/src/thermal_conduction.f90 +++ /dev/null @@ -1,242 +0,0 @@ -!-------------------------------------------------------------------------------------------------- -!> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH -!> @brief material subroutine for temperature evolution from heat conduction -!-------------------------------------------------------------------------------------------------- -module thermal_conduction - use prec - use material - use config - use lattice - use results - use constitutive - use YAML_types - use discretization - - implicit none - private - - type :: tParameters - character(len=pStringLen), allocatable, dimension(:) :: & - output - end type tParameters - - type(tparameters), dimension(:), allocatable :: & - param - - public :: & - thermal_conduction_init, & - thermal_conduction_getSource, & - thermal_conduction_getConductivity, & - thermal_conduction_getSpecificHeat, & - thermal_conduction_getMassDensity, & - thermal_conduction_putTemperatureAndItsRate, & - thermal_conduction_results - -contains - - -!-------------------------------------------------------------------------------------------------- -!> @brief module initialization -!> @details reads in material parameters, allocates arrays, and does sanity checks -!-------------------------------------------------------------------------------------------------- -subroutine thermal_conduction_init(T) - - real(pReal), dimension(:), intent(inout) :: T - - integer :: Ninstances,Nmaterialpoints,ho,ip,el,ce - class(tNode), pointer :: & - material_homogenization, & - homog, & - homogThermal - - - print'(/,a)', ' <<<+- thermal_conduction init -+>>>'; flush(6) - - Ninstances = count(thermal_type == THERMAL_conduction_ID) - allocate(param(Ninstances)) - - material_homogenization => config_material%get('homogenization') - do ho = 1, size(material_name_homogenization) - if (thermal_type(ho) /= THERMAL_conduction_ID) cycle - homog => material_homogenization%get(ho) - homogThermal => homog%get('thermal') - associate(prm => param(thermal_typeInstance(ho))) - -#if defined (__GFORTRAN__) - prm%output = output_asStrings(homogThermal) -#else - prm%output = homogThermal%get_asStrings('output',defaultVal=emptyStringArray) -#endif - - Nmaterialpoints=count(material_homogenizationAt==ho) - - allocate (temperature (ho)%p(Nmaterialpoints), source=thermal_initialT(ho)) - allocate (temperatureRate(ho)%p(Nmaterialpoints), source=0.0_pReal) - - end associate - enddo - - ce = 0 - do el = 1, discretization_Nelems - do ip = 1, discretization_nIPs - ce = ce + 1 - ho = material_homogenizationAt(el) - if (thermal_type(ho) == THERMAL_conduction_ID) T(ce) = thermal_initialT(ho) - enddo - enddo - -end subroutine thermal_conduction_init - - -!-------------------------------------------------------------------------------------------------- -!> @brief return heat generation rate -!-------------------------------------------------------------------------------------------------- -subroutine thermal_conduction_getSource(Tdot, ip,el) - - integer, intent(in) :: & - ip, & !< integration point number - el !< element number - real(pReal), intent(out) :: & - Tdot - - integer :: & - homog - - homog = material_homogenizationAt(el) - call constitutive_thermal_getRate(TDot, ip,el) - - Tdot = Tdot/real(homogenization_Nconstituents(homog),pReal) - -end subroutine thermal_conduction_getSource - - -!-------------------------------------------------------------------------------------------------- -!> @brief return homogenized thermal conductivity in reference configuration -!-------------------------------------------------------------------------------------------------- -function thermal_conduction_getConductivity(ip,el) - - integer, intent(in) :: & - ip, & !< integration point number - el !< element number - real(pReal), dimension(3,3) :: & - thermal_conduction_getConductivity - - integer :: & - co - - - thermal_conduction_getConductivity = 0.0_pReal - - do co = 1, homogenization_Nconstituents(material_homogenizationAt(el)) - thermal_conduction_getConductivity = thermal_conduction_getConductivity + & - crystallite_push33ToRef(co,ip,el,lattice_K(:,:,material_phaseAt(co,el))) - enddo - - thermal_conduction_getConductivity = thermal_conduction_getConductivity & - / real(homogenization_Nconstituents(material_homogenizationAt(el)),pReal) - -end function thermal_conduction_getConductivity - - -!-------------------------------------------------------------------------------------------------- -!> @brief returns homogenized specific heat capacity -!-------------------------------------------------------------------------------------------------- -function thermal_conduction_getSpecificHeat(ip,el) - - integer, intent(in) :: & - ip, & !< integration point number - el !< element number - real(pReal) :: & - thermal_conduction_getSpecificHeat - - integer :: & - co - - - thermal_conduction_getSpecificHeat = 0.0_pReal - - do co = 1, homogenization_Nconstituents(material_homogenizationAt(el)) - thermal_conduction_getSpecificHeat = thermal_conduction_getSpecificHeat & - + lattice_c_p(material_phaseAt(co,el)) - enddo - - thermal_conduction_getSpecificHeat = thermal_conduction_getSpecificHeat & - / real(homogenization_Nconstituents(material_homogenizationAt(el)),pReal) - -end function thermal_conduction_getSpecificHeat - - -!-------------------------------------------------------------------------------------------------- -!> @brief returns homogenized mass density -!-------------------------------------------------------------------------------------------------- -function thermal_conduction_getMassDensity(ip,el) - - integer, intent(in) :: & - ip, & !< integration point number - el !< element number - real(pReal) :: & - thermal_conduction_getMassDensity - - integer :: & - co - - - thermal_conduction_getMassDensity = 0.0_pReal - - do co = 1, homogenization_Nconstituents(material_homogenizationAt(el)) - thermal_conduction_getMassDensity = thermal_conduction_getMassDensity & - + lattice_rho(material_phaseAt(co,el)) - enddo - - thermal_conduction_getMassDensity = thermal_conduction_getMassDensity & - / real(homogenization_Nconstituents(material_homogenizationAt(el)),pReal) - -end function thermal_conduction_getMassDensity - - -!-------------------------------------------------------------------------------------------------- -!> @brief updates thermal state with solution from heat conduction PDE -!-------------------------------------------------------------------------------------------------- -subroutine thermal_conduction_putTemperatureAndItsRate(T,Tdot,ip,el) - - integer, intent(in) :: & - ip, & !< integration point number - el !< element number - real(pReal), intent(in) :: & - T, & - Tdot - integer :: & - homog, & - offset - - homog = material_homogenizationAt(el) - offset = material_homogenizationMemberAt(ip,el) - temperature (homog)%p(offset) = T - temperatureRate(homog)%p(offset) = Tdot - -end subroutine thermal_conduction_putTemperatureAndItsRate - - -!-------------------------------------------------------------------------------------------------- -!> @brief writes results to HDF5 output file -!-------------------------------------------------------------------------------------------------- -subroutine thermal_conduction_results(homog,group) - - integer, intent(in) :: homog - character(len=*), intent(in) :: group - - integer :: o - - associate(prm => param(damage_typeInstance(homog))) - outputsLoop: do o = 1,size(prm%output) - select case(trim(prm%output(o))) - case('T') - call results_writeDataset(group,temperature(homog)%p,'T',& - 'temperature','K') - end select - enddo outputsLoop - end associate - -end subroutine thermal_conduction_results - -end module thermal_conduction diff --git a/src/thermal_isothermal.f90 b/src/thermal_isothermal.f90 deleted file mode 100644 index 09e35931e..000000000 --- a/src/thermal_isothermal.f90 +++ /dev/null @@ -1,48 +0,0 @@ -!-------------------------------------------------------------------------------------------------- -!> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH -!> @brief material subroutine for isothermal temperature field -!-------------------------------------------------------------------------------------------------- -module thermal_isothermal - use prec - use config - use material - use discretization - - implicit none - public - -contains - -!-------------------------------------------------------------------------------------------------- -!> @brief allocates fields, reads information from material configuration file -!-------------------------------------------------------------------------------------------------- -subroutine thermal_isothermal_init(T) - - real(pReal), dimension(:), intent(inout) :: T - - integer :: Ninstances,Nmaterialpoints,ho,ip,el,ce - - print'(/,a)', ' <<<+- thermal_isothermal init -+>>>'; flush(6) - - do ho = 1, size(thermal_type) - if (thermal_type(ho) /= THERMAL_isothermal_ID) cycle - - Nmaterialpoints = count(material_homogenizationAt == ho) - - allocate(temperature (ho)%p(Nmaterialpoints),source=thermal_initialT(ho)) - allocate(temperatureRate(ho)%p(Nmaterialpoints),source = 0.0_pReal) - - enddo - - ce = 0 - do el = 1, discretization_Nelems - do ip = 1, discretization_nIPs - ce = ce + 1 - ho = material_homogenizationAt(el) - if (thermal_type(ho) == THERMAL_isothermal_ID) T(ce) = thermal_initialT(ho) - enddo - enddo - -end subroutine thermal_isothermal_init - -end module thermal_isothermal