From c729e7d53fcab996ad77dd3c3ff7756efc539059 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 24 Jan 2021 09:39:40 +0100 Subject: [PATCH 01/10] cell-based indexing --- src/constitutive.f90 | 20 +++++++++----------- src/constitutive_mech.f90 | 12 +++++------- src/homogenization.f90 | 10 ++++------ 3 files changed, 18 insertions(+), 24 deletions(-) diff --git a/src/constitutive.f90 b/src/constitutive.f90 index 16f3d3b59..16278957a 100644 --- a/src/constitutive.f90 +++ b/src/constitutive.f90 @@ -137,8 +137,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 @@ -497,26 +497,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 diff --git a/src/constitutive_mech.f90 b/src/constitutive_mech.f90 index cd1dedb61..d54367ad6 100644 --- a/src/constitutive_mech.f90 +++ b/src/constitutive_mech.f90 @@ -1656,11 +1656,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 +1666,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) diff --git a/src/homogenization.f90 b/src/homogenization.f90 index ed9b7fead..627c7c36e 100644 --- a/src/homogenization.f90 +++ b/src/homogenization.f90 @@ -190,9 +190,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 +226,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 +243,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 +254,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)) & From 4f059910ab286413190791745dbcb757befb55d7 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 24 Jan 2021 10:39:32 +0100 Subject: [PATCH 02/10] fix broken statistics reporting in case of multi-physics --- src/grid/DAMASK_grid.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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? From fede1dcd09e80e0a61baeb1b82f832b753668174 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 24 Jan 2021 11:16:17 +0100 Subject: [PATCH 03/10] dot_T is needed --- src/constitutive_thermal.f90 | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/src/constitutive_thermal.f90 b/src/constitutive_thermal.f90 index 2746ff889..610b0b4a6 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') @@ -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) From 983b59fe1efaf03395a5a0cf338fc9b340f3a341 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 24 Jan 2021 11:19:10 +0100 Subject: [PATCH 04/10] simpler access pattern --- src/constitutive.f90 | 14 +++----------- src/constitutive_mech.f90 | 6 ++++-- src/kinematics_thermal_expansion.f90 | 20 +++++++------------- 3 files changed, 14 insertions(+), 26 deletions(-) diff --git a/src/constitutive.f90 b/src/constitutive.f90 index 16278957a..e05336cc2 100644 --- a/src/constitutive.f90 +++ b/src/constitutive.f90 @@ -245,9 +245,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 +314,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) :: & @@ -764,9 +759,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 d54367ad6..fe07d2328 100644 --- a/src/constitutive_mech.f90 +++ b/src/constitutive_mech.f90 @@ -1959,7 +1959,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 @@ -1985,7 +1985,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/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 From 9f3fc6832531697307dc34eb5ac732b09a05a96a Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 24 Jan 2021 12:14:45 +0100 Subject: [PATCH 05/10] partitionedState not needed --- src/constitutive.f90 | 8 -------- src/constitutive_mech.f90 | 3 --- src/constitutive_thermal.f90 | 24 +++--------------------- src/homogenization.f90 | 1 - 4 files changed, 3 insertions(+), 33 deletions(-) diff --git a/src/constitutive.f90 b/src/constitutive.f90 index e05336cc2..6cf7b5e46 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 @@ -364,7 +360,6 @@ module constitutive constitutive_mech_setF, & constitutive_mech_getF, & constitutive_initializeRestorationPoints, & - constitutive_thermal_initializeRestorationPoints, & constitutive_windForward, & KINEMATICS_UNDEFINED_ID ,& KINEMATICS_CLEAVAGE_OPENING_ID, & @@ -631,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 diff --git a/src/constitutive_mech.f90 b/src/constitutive_mech.f90 index fe07d2328..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) diff --git a/src/constitutive_thermal.f90 b/src/constitutive_thermal.f90 index 610b0b4a6..a82262e81 100644 --- a/src/constitutive_thermal.f90 +++ b/src/constitutive_thermal.f90 @@ -116,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, & @@ -210,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 @@ -237,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 diff --git a/src/homogenization.f90 b/src/homogenization.f90 index 627c7c36e..86caafef3 100644 --- a/src/homogenization.f90 +++ b/src/homogenization.f90 @@ -279,7 +279,6 @@ subroutine materialpoint_stressAndItsTangent(dt,FEsolving_execIP,FEsolving_execE call thermal_partition(homogenization_T(ce),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' From a933fe57e70e27e5e596a1a6b0b7b1c3441e3f0c Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 24 Jan 2021 12:15:18 +0100 Subject: [PATCH 06/10] data structures for storing T at the homogenization level similar to thermal_conduction, better than homogenization_T --- src/homogenization_thermal.f90 | 44 ++++++++++++++++++++++++++++++++++ 1 file changed, 44 insertions(+) diff --git a/src/homogenization_thermal.f90 b/src/homogenization_thermal.f90 index de91a98d6..ce55d1573 100644 --- a/src/homogenization_thermal.f90 +++ b/src/homogenization_thermal.f90 @@ -4,6 +4,21 @@ submodule(homogenization) homogenization_thermal + 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,12 +26,41 @@ 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 From e22d76be9e309d28a7024dc4526b57ff5d486da1 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 24 Jan 2021 13:26:01 +0100 Subject: [PATCH 07/10] using new structure for thermal --- src/constitutive.f90 | 12 ++-- src/constitutive_thermal.f90 | 7 ++- src/grid/grid_thermal_spectral.f90 | 12 +++- src/homogenization.f90 | 34 ++++++++++- src/homogenization_thermal.f90 | 95 +++++++++++++++++++++++++++++- src/thermal_conduction.f90 | 87 --------------------------- 6 files changed, 142 insertions(+), 105 deletions(-) diff --git a/src/constitutive.f90 b/src/constitutive.f90 index 6cf7b5e46..149d09a49 100644 --- a/src/constitutive.f90 +++ b/src/constitutive.f90 @@ -200,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 @@ -353,7 +353,7 @@ module constitutive constitutive_restartWrite, & constitutive_restartRead, & integrateDamageState, & - constitutive_thermal_setT, & + constitutive_thermal_setField, & constitutive_damage_set_phi, & constitutive_damage_get_phi, & constitutive_mech_getP, & diff --git a/src/constitutive_thermal.f90 b/src/constitutive_thermal.f90 index a82262e81..57b8a3117 100644 --- a/src/constitutive_thermal.f90 +++ b/src/constitutive_thermal.f90 @@ -272,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/grid_thermal_spectral.f90 b/src/grid/grid_thermal_spectral.f90 index d164f11e7..44d9e22b1 100644 --- a/src/grid/grid_thermal_spectral.f90 +++ b/src/grid/grid_thermal_spectral.f90 @@ -191,6 +191,9 @@ function grid_thermal_spectral_solution(timeinc) result(solution) call thermal_conduction_putTemperatureAndItsRate(T_current(i,j,k), & (T_current(i,j,k)-T_lastInc(i,j,k))/params%timeinc, & 1,ce) + call homogenization_thermal_setField(T_current(i,j,k), & + (T_current(i,j,k)-T_lastInc(i,j,k))/params%timeinc, & + ce) homogenization_T(ce) = T_current(i,j,k) enddo; enddo; enddo @@ -233,6 +236,9 @@ subroutine grid_thermal_spectral_forward(cutBack) (T_current(i,j,k) - & T_lastInc(i,j,k))/params%timeinc, & 1,ce) + call homogenization_thermal_setField(T_current(i,j,k), & + (T_current(i,j,k)-T_lastInc(i,j,k))/params%timeinc, & + ce) homogenization_T(ce) = T_current(i,j,k) enddo; enddo; enddo else @@ -283,8 +289,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 +321,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 86caafef3..f316d3969 100644 --- a/src/homogenization.f90 +++ b/src/homogenization.f90 @@ -75,8 +75,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 +111,40 @@ 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 + end interface public :: & homogenization_init, & materialpoint_stressAndItsTangent, & + thermal_conduction_getSpecificHeat, & + thermal_conduction_getConductivity, & + thermal_conduction_getMassDensity, & + homogenization_thermal_setfield, & homogenization_forward, & homogenization_results, & homogenization_restartRead, & @@ -276,7 +304,7 @@ 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) if (.not. thermal_stress(dt,ph,material_phaseMemberAt(co,ip,el))) then diff --git a/src/homogenization_thermal.f90 b/src/homogenization_thermal.f90 index ce55d1573..80137f13a 100644 --- a/src/homogenization_thermal.f90 +++ b/src/homogenization_thermal.f90 @@ -3,6 +3,7 @@ !-------------------------------------------------------------------------------------------------- submodule(homogenization) homogenization_thermal + use lattice type :: tDataContainer real(pReal), dimension(:), allocatable :: T, dot_T @@ -67,15 +68,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 @@ -93,4 +97,89 @@ module subroutine thermal_homogenize(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 + end submodule homogenization_thermal diff --git a/src/thermal_conduction.f90 b/src/thermal_conduction.f90 index 1711c5512..ccca75026 100644 --- a/src/thermal_conduction.f90 +++ b/src/thermal_conduction.f90 @@ -26,9 +26,6 @@ module thermal_conduction public :: & thermal_conduction_init, & thermal_conduction_getSource, & - thermal_conduction_getConductivity, & - thermal_conduction_getSpecificHeat, & - thermal_conduction_getMassDensity, & thermal_conduction_putTemperatureAndItsRate, & thermal_conduction_results @@ -110,90 +107,6 @@ subroutine thermal_conduction_getSource(Tdot, ip,el) 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 !-------------------------------------------------------------------------------------------------- From c2ae2c919be1896d670b940ac2703a1c0f3c8490 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 24 Jan 2021 15:19:57 +0100 Subject: [PATCH 08/10] use new structure --- src/grid/grid_thermal_spectral.f90 | 4 +--- src/homogenization.f90 | 13 +++++++++---- src/homogenization_thermal.f90 | 16 ++++++++++++---- src/thermal_conduction.f90 | 13 +------------ src/thermal_isothermal.f90 | 13 +------------ 5 files changed, 24 insertions(+), 35 deletions(-) diff --git a/src/grid/grid_thermal_spectral.f90 b/src/grid/grid_thermal_spectral.f90 index 44d9e22b1..4b80c1d72 100644 --- a/src/grid/grid_thermal_spectral.f90 +++ b/src/grid/grid_thermal_spectral.f90 @@ -132,7 +132,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 @@ -194,7 +194,6 @@ function grid_thermal_spectral_solution(timeinc) result(solution) call homogenization_thermal_setField(T_current(i,j,k), & (T_current(i,j,k)-T_lastInc(i,j,k))/params%timeinc, & ce) - homogenization_T(ce) = T_current(i,j,k) enddo; enddo; enddo call VecMin(solution_vec,devNull,T_min,ierr); CHKERRQ(ierr) @@ -239,7 +238,6 @@ subroutine grid_thermal_spectral_forward(cutBack) call homogenization_thermal_setField(T_current(i,j,k), & (T_current(i,j,k)-T_lastInc(i,j,k))/params%timeinc, & ce) - homogenization_T(ce) = T_current(i,j,k) enddo; enddo; enddo else T_lastInc = T_current diff --git a/src/homogenization.f90 b/src/homogenization.f90 index f316d3969..51d339da9 100644 --- a/src/homogenization.f90 +++ b/src/homogenization.f90 @@ -28,8 +28,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 :: & @@ -136,6 +134,12 @@ module homogenization real(pReal), intent(in) :: T, dot_T end subroutine homogenization_thermal_setField + + module function homogenization_thermal_T(ce) result(T) + integer, intent(in) :: ce + real(pReal) :: T + end function homogenization_thermal_T + end interface public :: & @@ -145,6 +149,7 @@ module homogenization thermal_conduction_getConductivity, & thermal_conduction_getMassDensity, & homogenization_thermal_setfield, & + homogenization_thermal_T, & homogenization_forward, & homogenization_results, & homogenization_restartRead, & @@ -182,8 +187,8 @@ 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(thermal_type == THERMAL_isothermal_ID)) call thermal_isothermal_init() + if (any(thermal_type == THERMAL_conduction_ID)) call thermal_conduction_init() if (any(damage_type == DAMAGE_none_ID)) call damage_none_init if (any(damage_type == DAMAGE_nonlocal_ID)) call damage_nonlocal_init diff --git a/src/homogenization_thermal.f90 b/src/homogenization_thermal.f90 index 80137f13a..b9cbccc66 100644 --- a/src/homogenization_thermal.f90 +++ b/src/homogenization_thermal.f90 @@ -37,9 +37,6 @@ module subroutine thermal_init() 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)) @@ -92,7 +89,7 @@ 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 @@ -182,4 +179,15 @@ module subroutine homogenization_thermal_setField(T,dot_T, ce) end subroutine homogenization_thermal_setField + +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 + + end submodule homogenization_thermal diff --git a/src/thermal_conduction.f90 b/src/thermal_conduction.f90 index ccca75026..84521afde 100644 --- a/src/thermal_conduction.f90 +++ b/src/thermal_conduction.f90 @@ -36,9 +36,7 @@ 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 +subroutine thermal_conduction_init() integer :: Ninstances,Nmaterialpoints,ho,ip,el,ce class(tNode), pointer :: & @@ -73,15 +71,6 @@ subroutine thermal_conduction_init(T) 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 diff --git a/src/thermal_isothermal.f90 b/src/thermal_isothermal.f90 index 09e35931e..e3dee2901 100644 --- a/src/thermal_isothermal.f90 +++ b/src/thermal_isothermal.f90 @@ -16,9 +16,7 @@ contains !-------------------------------------------------------------------------------------------------- !> @brief allocates fields, reads information from material configuration file !-------------------------------------------------------------------------------------------------- -subroutine thermal_isothermal_init(T) - - real(pReal), dimension(:), intent(inout) :: T +subroutine thermal_isothermal_init() integer :: Ninstances,Nmaterialpoints,ho,ip,el,ce @@ -34,15 +32,6 @@ subroutine thermal_isothermal_init(T) 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 From 599dc2a2c67a033b9ac220a59eeb5d5480fb163b Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 24 Jan 2021 16:34:51 +0100 Subject: [PATCH 09/10] base HDF5 output on new data --- src/homogenization.f90 | 4 +++ src/homogenization_thermal.f90 | 23 +++++++++++++++ src/thermal_conduction.f90 | 52 ++-------------------------------- 3 files changed, 30 insertions(+), 49 deletions(-) diff --git a/src/homogenization.f90 b/src/homogenization.f90 index 51d339da9..67c5ab001 100644 --- a/src/homogenization.f90 +++ b/src/homogenization.f90 @@ -134,6 +134,10 @@ module homogenization 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 diff --git a/src/homogenization_thermal.f90 b/src/homogenization_thermal.f90 index b9cbccc66..cfbd2136d 100644 --- a/src/homogenization_thermal.f90 +++ b/src/homogenization_thermal.f90 @@ -180,6 +180,29 @@ module subroutine homogenization_thermal_setField(T,dot_T, ce) 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 diff --git a/src/thermal_conduction.f90 b/src/thermal_conduction.f90 index 84521afde..613612113 100644 --- a/src/thermal_conduction.f90 +++ b/src/thermal_conduction.f90 @@ -15,19 +15,10 @@ module thermal_conduction 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_putTemperatureAndItsRate, & - thermal_conduction_results + thermal_conduction_putTemperatureAndItsRate contains @@ -38,37 +29,22 @@ contains !-------------------------------------------------------------------------------------------------- subroutine thermal_conduction_init() - integer :: Ninstances,Nmaterialpoints,ho,ip,el,ce + integer :: Nmaterialpoints,ho class(tNode), pointer :: & - material_homogenization, & - homog, & - homogThermal + material_homogenization 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 end subroutine thermal_conduction_init @@ -119,26 +95,4 @@ subroutine thermal_conduction_putTemperatureAndItsRate(T,Tdot,ip,el) 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 From 26c7969837f27aa44a139912d127bfa6af81fd8b Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 24 Jan 2021 17:04:41 +0100 Subject: [PATCH 10/10] not needed anymore --- src/CPFEM.f90 | 8 +-- src/DAMASK_marc.f90 | 2 +- src/commercialFEM_fileList.f90 | 2 - src/grid/grid_thermal_spectral.f90 | 8 --- src/homogenization.f90 | 13 ++-- src/homogenization_thermal.f90 | 23 +++++++ src/material.f90 | 6 +- src/thermal_conduction.f90 | 98 ------------------------------ src/thermal_isothermal.f90 | 37 ----------- 9 files changed, 37 insertions(+), 160 deletions(-) delete mode 100644 src/thermal_conduction.f90 delete mode 100644 src/thermal_isothermal.f90 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/grid/grid_thermal_spectral.f90 b/src/grid/grid_thermal_spectral.f90 index 4b80c1d72..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 @@ -188,9 +187,6 @@ 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), & - (T_current(i,j,k)-T_lastInc(i,j,k))/params%timeinc, & - 1,ce) call homogenization_thermal_setField(T_current(i,j,k), & (T_current(i,j,k)-T_lastInc(i,j,k))/params%timeinc, & ce) @@ -231,10 +227,6 @@ 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) call homogenization_thermal_setField(T_current(i,j,k), & (T_current(i,j,k)-T_lastInc(i,j,k))/params%timeinc, & ce) diff --git a/src/homogenization.f90 b/src/homogenization.f90 index 67c5ab001..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 @@ -144,6 +142,13 @@ module homogenization 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 :: & @@ -152,6 +157,7 @@ module homogenization thermal_conduction_getSpecificHeat, & thermal_conduction_getConductivity, & thermal_conduction_getMassDensity, & + thermal_conduction_getSource, & homogenization_thermal_setfield, & homogenization_thermal_T, & homogenization_forward, & @@ -191,9 +197,6 @@ subroutine homogenization_init() call thermal_init() call damage_init() - if (any(thermal_type == THERMAL_isothermal_ID)) call thermal_isothermal_init() - if (any(thermal_type == THERMAL_conduction_ID)) call thermal_conduction_init() - if (any(damage_type == DAMAGE_none_ID)) call damage_none_init if (any(damage_type == DAMAGE_nonlocal_ID)) call damage_nonlocal_init diff --git a/src/homogenization_thermal.f90 b/src/homogenization_thermal.f90 index cfbd2136d..6a12730ac 100644 --- a/src/homogenization_thermal.f90 +++ b/src/homogenization_thermal.f90 @@ -213,4 +213,27 @@ module function homogenization_thermal_T(ce) result(T) 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/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 613612113..000000000 --- a/src/thermal_conduction.f90 +++ /dev/null @@ -1,98 +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 - - public :: & - thermal_conduction_init, & - thermal_conduction_getSource, & - thermal_conduction_putTemperatureAndItsRate - -contains - - -!-------------------------------------------------------------------------------------------------- -!> @brief module initialization -!> @details reads in material parameters, allocates arrays, and does sanity checks -!-------------------------------------------------------------------------------------------------- -subroutine thermal_conduction_init() - - integer :: Nmaterialpoints,ho - class(tNode), pointer :: & - material_homogenization - - - print'(/,a)', ' <<<+- thermal_conduction init -+>>>'; flush(6) - - material_homogenization => config_material%get('homogenization') - do ho = 1, size(material_name_homogenization) - if (thermal_type(ho) /= THERMAL_conduction_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 - -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 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 - - -end module thermal_conduction diff --git a/src/thermal_isothermal.f90 b/src/thermal_isothermal.f90 deleted file mode 100644 index e3dee2901..000000000 --- a/src/thermal_isothermal.f90 +++ /dev/null @@ -1,37 +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() - - 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 - -end subroutine thermal_isothermal_init - -end module thermal_isothermal