diff --git a/src/constitutive.f90 b/src/constitutive.f90 index f9537a136..6a023022f 100644 --- a/src/constitutive.f90 +++ b/src/constitutive.f90 @@ -90,6 +90,7 @@ module constitutive phase_kinematics !< active kinematic mechanisms of each phase integer, dimension(:), allocatable, public :: & !< ToDo: should be protected (bug in Intel compiler) + thermal_Nsources, & phase_Nsources, & !< number of source mechanisms active in each phase phase_Nkinematics, & !< number of kinematic mechanisms active in each phase phase_NstiffnessDegradations, & !< number of stiffness degradation mechanisms active in each phase @@ -233,6 +234,16 @@ module constitutive ! == cleaned:end =================================================================================== + module function integrateThermalState(dt,co,ip,el) result(broken) + + real(pReal), intent(in) :: dt + integer, intent(in) :: & + el, & !< element index in element loop + ip, & !< integration point index in ip loop + co !< grain index in grain loop + logical :: broken + end function + module function crystallite_stress(dt,co,ip,el) result(converged_) real(pReal), intent(in) :: dt integer, intent(in) :: co, ip, el @@ -665,31 +676,6 @@ function constitutive_damage_collectDotState(co,ip,el,ph,of) result(broken) end function constitutive_damage_collectDotState -!-------------------------------------------------------------------------------------------------- -!> @brief contains the constitutive equation for calculating the rate of change of microstructure -!-------------------------------------------------------------------------------------------------- -function constitutive_thermal_collectDotState(ph,me) result(broken) - - integer, intent(in) :: ph, me - logical :: broken - - integer :: i - - - broken = .false. - - SourceLoop: do i = 1, phase_Nsources(ph) - - if (phase_source(i,ph) == SOURCE_thermal_externalheat_ID) & - call source_thermal_externalheat_dotState(ph,me) - - broken = broken .or. any(IEEE_is_NaN(sourceState(ph)%p(i)%dotState(:,me))) - - enddo SourceLoop - -end function constitutive_thermal_collectDotState - - !-------------------------------------------------------------------------------------------------- !> @brief for constitutive models having an instantaneous change of state !> will return false if delta state is not needed/supported by the constitutive model @@ -856,7 +842,7 @@ subroutine crystallite_init() cMax, & !< maximum number of integration point components iMax, & !< maximum number of integration points eMax !< maximum number of elements - + class(tNode), pointer :: & num_crystallite, & @@ -914,6 +900,9 @@ subroutine crystallite_init() do so = 1, phase_Nsources(ph) allocate(sourceState(ph)%p(so)%subState0,source=sourceState(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 @@ -1144,111 +1133,6 @@ function integrateSourceState(dt,co,ip,el) result(broken) end function integrateSourceState - -!-------------------------------------------------------------------------------------------------- -!> @brief integrate stress, state with adaptive 1st order explicit Euler method -!> using Fixed Point Iteration to adapt the stepsize -!-------------------------------------------------------------------------------------------------- -function integrateThermalState(dt,co,ip,el) result(broken) - - real(pReal), intent(in) :: dt - integer, intent(in) :: & - el, & !< element index in element loop - ip, & !< integration point index in ip loop - co !< grain index in grain loop - - integer :: & - NiterationState, & !< number of iterations in state loop - ph, & - me, & - so - integer, dimension(maxval(phase_Nsources)) :: & - size_so - real(pReal) :: & - zeta - real(pReal), dimension(constitutive_source_maxSizeDotState) :: & - r ! state residuum - real(pReal), dimension(constitutive_source_maxSizeDotState,2,maxval(phase_Nsources)) :: source_dotState - logical :: & - broken, converged_ - - - ph = material_phaseAt(co,el) - me = material_phaseMemberAt(co,ip,el) - - converged_ = .true. - broken = constitutive_thermal_collectDotState(ph,me) - if(broken) return - - do so = 1, phase_Nsources(ph) - size_so(so) = thermalState(ph)%p(so)%sizeDotState - thermalState(ph)%p(so)%state(1:size_so(so),me) = thermalState(ph)%p(so)%subState0(1:size_so(so),me) & - + thermalState(ph)%p(so)%dotState (1:size_so(so),me) * dt - source_dotState(1:size_so(so),2,so) = 0.0_pReal - enddo - - iteration: do NiterationState = 1, num%nState - - do so = 1, phase_Nsources(ph) - if(nIterationState > 1) source_dotState(1:size_so(so),2,so) = source_dotState(1:size_so(so),1,so) - source_dotState(1:size_so(so),1,so) = thermalState(ph)%p(so)%dotState(:,me) - enddo - - broken = constitutive_thermal_collectDotState(ph,me) - broken = broken .or. constitutive_damage_collectDotState(co,ip,el,ph,me) - if(broken) exit iteration - - do so = 1, phase_Nsources(ph) - zeta = damper(thermalState(ph)%p(so)%dotState(:,me), & - source_dotState(1:size_so(so),1,so),& - source_dotState(1:size_so(so),2,so)) - thermalState(ph)%p(so)%dotState(:,me) = thermalState(ph)%p(so)%dotState(:,me) * zeta & - + source_dotState(1:size_so(so),1,so)* (1.0_pReal - zeta) - r(1:size_so(so)) = thermalState(ph)%p(so)%state (1:size_so(so),me) & - - thermalState(ph)%p(so)%subState0(1:size_so(so),me) & - - thermalState(ph)%p(so)%dotState (1:size_so(so),me) * dt - thermalState(ph)%p(so)%state(1:size_so(so),me) = thermalState(ph)%p(so)%state(1:size_so(so),me) & - - r(1:size_so(so)) - converged_ = converged_ .and. converged(r(1:size_so(so)), & - thermalState(ph)%p(so)%state(1:size_so(so),me), & - thermalState(ph)%p(so)%atol(1:size_so(so))) - enddo - - if(converged_) then - broken = constitutive_damage_deltaState(mech_F_e(ph,me),co,ip,el,ph,me) - exit iteration - endif - - enddo iteration - - broken = broken .or. .not. converged_ - - - contains - - !-------------------------------------------------------------------------------------------------- - !> @brief calculate the damping for correction of state and dot state - !-------------------------------------------------------------------------------------------------- - real(pReal) pure function damper(current,previous,previous2) - - real(pReal), dimension(:), intent(in) ::& - current, previous, previous2 - - real(pReal) :: dot_prod12, dot_prod22 - - dot_prod12 = dot_product(current - previous, previous - previous2) - dot_prod22 = dot_product(previous - previous2, previous - previous2) - if ((dot_product(current,previous) < 0.0_pReal .or. dot_prod12 < 0.0_pReal) .and. dot_prod22 > 0.0_pReal) then - damper = 0.75_pReal + 0.25_pReal * tanh(2.0_pReal + 4.0_pReal * dot_prod12 / dot_prod22) - else - damper = 1.0_pReal - endif - - end function damper - -end function integrateThermalState - - !-------------------------------------------------------------------------------------------------- !> @brief determines whether a point is converged !-------------------------------------------------------------------------------------------------- diff --git a/src/constitutive_mech.f90 b/src/constitutive_mech.f90 index fccccf00c..9a065a829 100644 --- a/src/constitutive_mech.f90 +++ b/src/constitutive_mech.f90 @@ -1588,6 +1588,9 @@ module function crystallite_stress(dt,co,ip,el) result(converged_) do so = 1, phase_Nsources(ph) sourceState(ph)%p(so)%subState0(:,me) = sourceState(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) @@ -1616,6 +1619,9 @@ module function crystallite_stress(dt,co,ip,el) result(converged_) do so = 1, phase_Nsources(ph) sourceState(ph)%p(so)%subState0(:,me) = sourceState(ph)%p(so)%state(:,me) enddo + do so = 1, thermal_Nsources(ph) + thermalState(ph)%p(so)%subState0(:,me) = thermalState(ph)%p(so)%state(:,me) + enddo endif !-------------------------------------------------------------------------------------------------- ! cut back (reduced time and restore) @@ -1632,6 +1638,9 @@ module function crystallite_stress(dt,co,ip,el) result(converged_) do so = 1, phase_Nsources(ph) sourceState(ph)%p(so)%state(:,me) = sourceState(ph)%p(so)%subState0(:,me) enddo + do so = 1, thermal_Nsources(ph) + thermalState(ph)%p(so)%state(:,me) = thermalState(ph)%p(so)%subState0(:,me) + enddo todo = subStep > num%subStepMinCryst ! still on track or already done (beyond repair) endif @@ -1645,6 +1654,7 @@ module function crystallite_stress(dt,co,ip,el) result(converged_) constitutive_mech_Fp(ph)%data(1:3,1:3,me)))) converged_ = .not. integrateState(subF0,subF,subFp0,subFi0,subState0(1:sizeDotState),subStep * dt,co,ip,el) converged_ = converged_ .and. .not. integrateSourceState(subStep * dt,co,ip,el) + converged_ = converged_ .and. .not. integrateThermalState(subStep * dt,co,ip,el) endif enddo cutbackLooping diff --git a/src/constitutive_thermal.f90 b/src/constitutive_thermal.f90 index f1675f0a1..c86a286f9 100644 --- a/src/constitutive_thermal.f90 +++ b/src/constitutive_thermal.f90 @@ -6,9 +6,12 @@ submodule(constitutive) constitutive_thermal type :: tDataContainer real(pReal), dimension(:), allocatable :: T end type tDataContainer + integer(kind(SOURCE_undefined_ID)), dimension(:,:), allocatable :: & + thermal_source type(tDataContainer), dimension(:), allocatable :: current + integer :: thermal_source_maxSizeDotState interface module function source_thermal_dissipation_init(source_length) result(mySources) @@ -60,30 +63,55 @@ module subroutine thermal_init(phases) class(tNode), pointer :: & phases + class(tNode), pointer :: & + phase, thermal, sources + integer :: & - ph, & + ph, so, & Nconstituents - print'(/,a)', ' <<<+- constitutive_mech init -+>>>' + print'(/,a)', ' <<<+- constitutive_thermal init -+>>>' allocate(current(phases%length)) + allocate(thermalState (phases%length)) + allocate(thermal_Nsources(phases%length),source = 0) do ph = 1, phases%length Nconstituents = count(material_phaseAt == ph) * discretization_nIPs allocate(current(ph)%T(Nconstituents)) + phase => phases%get(ph) + if(phase%contains('thermal')) then + thermal => phase%get('thermal') + sources => thermal%get('source',defaultVal=emptyList) + thermal_Nsources(ph) = sources%length + endif + allocate(thermalstate(ph)%p(thermal_Nsources(ph))) enddo -! initialize source mechanisms - if(maxval(phase_Nsources) /= 0) then - where(source_thermal_dissipation_init (maxval(phase_Nsources))) phase_source = SOURCE_thermal_dissipation_ID - where(source_thermal_externalheat_init(maxval(phase_Nsources))) phase_source = SOURCE_thermal_externalheat_ID + allocate(thermal_source(maxval(thermal_Nsources),phases%length), source = SOURCE_undefined_ID) + + if(maxval(thermal_Nsources) /= 0) then + where(source_thermal_dissipation_init (maxval(thermal_Nsources))) thermal_source = SOURCE_thermal_dissipation_ID + where(source_thermal_externalheat_init(maxval(thermal_Nsources))) thermal_source = SOURCE_thermal_externalheat_ID endif + thermal_source_maxSizeDotState = 0 + 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 + enddo + + thermal_source_maxSizeDotState = max(thermal_source_maxSizeDotState, & + maxval(thermalState(ph)%p%sizeDotState)) + enddo PhaseLoop2 + !-------------------------------------------------------------------------------------------------- !initialize kinematic mechanisms if(maxval(phase_Nkinematics) /= 0) where(kinematics_thermal_expansion_init(maxval(phase_Nkinematics))) & @@ -123,8 +151,8 @@ module subroutine constitutive_thermal_getRateAndItsTangents(TDot, dTDot_dT, T, do co = 1, homogenization_Nconstituents(homog) ph = material_phaseAt(co,el) me = material_phasememberAt(co,ip,el) - do so = 1, phase_Nsources(ph) - select case(phase_source(so,ph)) + do so = 1, thermal_Nsources(ph) + select case(thermal_source(so,ph)) case (SOURCE_thermal_dissipation_ID) call source_thermal_dissipation_getRateAndItsTangent(my_Tdot, my_dTdot_dT, & mech_S(ph,me),mech_L_p(ph,me), ph) @@ -145,6 +173,131 @@ module subroutine constitutive_thermal_getRateAndItsTangents(TDot, dTDot_dT, T, end subroutine constitutive_thermal_getRateAndItsTangents +!-------------------------------------------------------------------------------------------------- +!> @brief contains the constitutive equation for calculating the rate of change of microstructure +!-------------------------------------------------------------------------------------------------- +function constitutive_thermal_collectDotState(ph,me) result(broken) + + integer, intent(in) :: ph, me + logical :: broken + + integer :: i + + + broken = .false. + + SourceLoop: do i = 1, thermal_Nsources(ph) + + if (thermal_source(i,ph) == SOURCE_thermal_externalheat_ID) & + call source_thermal_externalheat_dotState(ph,me) + + broken = broken .or. any(IEEE_is_NaN(thermalState(ph)%p(i)%dotState(:,me))) + + enddo SourceLoop + +end function constitutive_thermal_collectDotState + + +!-------------------------------------------------------------------------------------------------- +!> @brief integrate stress, state with adaptive 1st order explicit Euler method +!> using Fixed Point Iteration to adapt the stepsize +!-------------------------------------------------------------------------------------------------- +module function integrateThermalState(dt,co,ip,el) result(broken) + + real(pReal), intent(in) :: dt + integer, intent(in) :: & + el, & !< element index in element loop + ip, & !< integration point index in ip loop + co !< grain index in grain loop + + integer :: & + NiterationState, & !< number of iterations in state loop + ph, & + me, & + so + integer, dimension(maxval(thermal_Nsources)) :: & + size_so + real(pReal) :: & + zeta + real(pReal), dimension(thermal_source_maxSizeDotState) :: & + r ! state residuum + real(pReal), dimension(thermal_source_maxSizeDotState,2,maxval(thermal_Nsources)) :: source_dotState + logical :: & + broken, converged_ + + + ph = material_phaseAt(co,el) + me = material_phaseMemberAt(co,ip,el) + + converged_ = .true. + broken = constitutive_thermal_collectDotState(ph,me) + if(broken) return + + do so = 1, thermal_Nsources(ph) + size_so(so) = thermalState(ph)%p(so)%sizeDotState + thermalState(ph)%p(so)%state(1:size_so(so),me) = thermalState(ph)%p(so)%subState0(1:size_so(so),me) & + + thermalState(ph)%p(so)%dotState (1:size_so(so),me) * dt + source_dotState(1:size_so(so),2,so) = 0.0_pReal + enddo + + iteration: do NiterationState = 1, num%nState + + do so = 1, thermal_Nsources(ph) + if(nIterationState > 1) source_dotState(1:size_so(so),2,so) = source_dotState(1:size_so(so),1,so) + source_dotState(1:size_so(so),1,so) = thermalState(ph)%p(so)%dotState(:,me) + enddo + + broken = constitutive_thermal_collectDotState(ph,me) + if(broken) exit iteration + + do so = 1, thermal_Nsources(ph) + zeta = damper(thermalState(ph)%p(so)%dotState(:,me), & + source_dotState(1:size_so(so),1,so),& + source_dotState(1:size_so(so),2,so)) + thermalState(ph)%p(so)%dotState(:,me) = thermalState(ph)%p(so)%dotState(:,me) * zeta & + + source_dotState(1:size_so(so),1,so)* (1.0_pReal - zeta) + r(1:size_so(so)) = thermalState(ph)%p(so)%state (1:size_so(so),me) & + - thermalState(ph)%p(so)%subState0(1:size_so(so),me) & + - thermalState(ph)%p(so)%dotState (1:size_so(so),me) * dt + thermalState(ph)%p(so)%state(1:size_so(so),me) = thermalState(ph)%p(so)%state(1:size_so(so),me) & + - r(1:size_so(so)) + converged_ = converged_ .and. converged(r(1:size_so(so)), & + thermalState(ph)%p(so)%state(1:size_so(so),me), & + thermalState(ph)%p(so)%atol(1:size_so(so))) + enddo + + if(converged_) exit iteration + + enddo iteration + + broken = broken .or. .not. converged_ + + + contains + + !-------------------------------------------------------------------------------------------------- + !> @brief calculate the damping for correction of state and dot state + !-------------------------------------------------------------------------------------------------- + real(pReal) pure function damper(current,previous,previous2) + + real(pReal), dimension(:), intent(in) ::& + current, previous, previous2 + + real(pReal) :: dot_prod12, dot_prod22 + + dot_prod12 = dot_product(current - previous, previous - previous2) + dot_prod22 = dot_product(previous - previous2, previous - previous2) + if ((dot_product(current,previous) < 0.0_pReal .or. dot_prod12 < 0.0_pReal) .and. dot_prod22 > 0.0_pReal) then + damper = 0.75_pReal + 0.25_pReal * tanh(2.0_pReal + 4.0_pReal * dot_prod12 / dot_prod22) + else + damper = 1.0_pReal + endif + + end function damper + +end function integrateThermalState + + module subroutine thermal_initializeRestorationPoints(ph,me) @@ -153,7 +306,7 @@ module subroutine thermal_initializeRestorationPoints(ph,me) integer :: so - do so = 1, size(sourceState(ph)%p) + do so = 1, size(thermalState(ph)%p) thermalState(ph)%p(so)%partitionedState0(:,me) = thermalState(ph)%p(so)%state0(:,me) enddo @@ -168,7 +321,7 @@ module subroutine thermal_windForward(ph,me) integer :: so - do so = 1, size(sourceState(ph)%p) + do so = 1, size(thermalState(ph)%p) thermalState(ph)%p(so)%partitionedState0(:,me) = thermalState(ph)%p(so)%state(:,me) enddo @@ -181,7 +334,7 @@ module subroutine thermal_forward() do ph = 1, size(thermalState) - do so = 1, size(sourceState(ph)%p) + do so = 1, size(thermalState(ph)%p) thermalState(ph)%p(so)%state0 = thermalState(ph)%p(so)%state enddo enddo @@ -200,7 +353,7 @@ module subroutine thermal_restore(ip,el) ph = material_phaseAt(co,el) me = material_phaseMemberAt(co,ip,el) - do so = 1, size(sourceState(ph)%p) + do so = 1, size(thermalState(ph)%p) thermalState(ph)%p(so)%state(:,me) = thermalState(ph)%p(so)%partitionedState0(:,me) enddo @@ -237,4 +390,39 @@ module subroutine constitutive_thermal_setT(T,co,ip,el) end subroutine constitutive_thermal_setT + +!-------------------------------------------------------------------------------------------------- +!> @brief checks if a source mechanism is active or not +!-------------------------------------------------------------------------------------------------- +function thermal_active(source_label,src_length) result(active_source) + + character(len=*), intent(in) :: source_label !< name of source mechanism + integer, intent(in) :: src_length !< max. number of sources in system + logical, dimension(:,:), allocatable :: active_source + + class(tNode), pointer :: & + phases, & + phase, & + sources, thermal, & + src + integer :: p,s + + phases => config_material%get('phase') + allocate(active_source(src_length,phases%length), source = .false. ) + do p = 1, phases%length + phase => phases%get(p) + if (phase%contains('thermal')) then + thermal => phase%get('thermal',defaultVal=emptyList) + sources => thermal%get('source',defaultVal=emptyList) + do s = 1, sources%length + src => sources%get(s) + if(src%get_asString('type') == source_label) active_source(s,p) = .true. + enddo + endif + enddo + + +end function thermal_active + + end submodule constitutive_thermal diff --git a/src/constitutive_thermal_dissipation.f90 b/src/constitutive_thermal_dissipation.f90 index 27653a9ef..44227536c 100644 --- a/src/constitutive_thermal_dissipation.f90 +++ b/src/constitutive_thermal_dissipation.f90 @@ -62,7 +62,7 @@ module function source_thermal_dissipation_init(source_length) result(mySources) src => sources%get(sourceOffset) prm%kappa = src%get_asFloat('kappa') Nconstituents = count(material_phaseAt==p) * discretization_nIPs - call constitutive_allocateState(sourceState(p)%p(sourceOffset),Nconstituents,0,0,0) + call constitutive_allocateState(thermalState(p)%p(sourceOffset),Nconstituents,0,0,0) end associate endif diff --git a/src/constitutive_thermal_externalheat.f90 b/src/constitutive_thermal_externalheat.f90 index 3ef96790e..de1617efa 100644 --- a/src/constitutive_thermal_externalheat.f90 +++ b/src/constitutive_thermal_externalheat.f90 @@ -13,7 +13,7 @@ submodule(constitutive:constitutive_thermal) source_externalheat type :: tParameters !< container type for internal constitutive parameters real(pReal), dimension(:), allocatable :: & - t_n, & + t_n, & f_T integer :: & nIntervals @@ -31,19 +31,20 @@ contains !-------------------------------------------------------------------------------------------------- module function source_thermal_externalheat_init(source_length) result(mySources) - integer, intent(in) :: source_length + integer, intent(in) :: source_length logical, dimension(:,:), allocatable :: mySources class(tNode), pointer :: & phases, & phase, & - sources, & - src + sources, thermal, & + src integer :: Ninstances,sourceOffset,Nconstituents,p print'(/,a)', ' <<<+- source_thermal_externalHeat init -+>>>' - mySources = source_active('thermal_externalheat',source_length) + mySources = thermal_active('externalheat',source_length) + Ninstances = count(mySources) print'(a,i2)', ' # instances: ',Ninstances; flush(IO_STDOUT) if(Ninstances == 0) return @@ -54,15 +55,16 @@ module function source_thermal_externalheat_init(source_length) result(mySources allocate(source_thermal_externalheat_instance(phases%length), source=0) do p = 1, phases%length - phase => phases%get(p) + phase => phases%get(p) if(any(mySources(:,p))) source_thermal_externalheat_instance(p) = count(mySources(:,1:p)) if(count(mySources(:,p)) == 0) cycle - sources => phase%get('source') + thermal => phase%get('thermal') + sources => thermal%get('source') do sourceOffset = 1, sources%length if(mySources(sourceOffset,p)) then source_thermal_externalheat_offset(p) = sourceOffset associate(prm => param(source_thermal_externalheat_instance(p))) - src => sources%get(sourceOffset) + src => sources%get(sourceOffset) prm%t_n = src%get_asFloats('t_n') prm%nIntervals = size(prm%t_n) - 1 @@ -70,7 +72,7 @@ module function source_thermal_externalheat_init(source_length) result(mySources prm%f_T = src%get_asFloats('f_T',requiredSize = size(prm%t_n)) Nconstituents = count(material_phaseAt==p) * discretization_nIPs - call constitutive_allocateState(sourceState(p)%p(sourceOffset),Nconstituents,1,1,0) + call constitutive_allocateState(thermalState(p)%p(sourceOffset),Nconstituents,1,1,0) end associate endif @@ -95,7 +97,7 @@ module subroutine source_thermal_externalheat_dotState(phase, of) sourceOffset = source_thermal_externalheat_offset(phase) - sourceState(phase)%p(sourceOffset)%dotState(1,of) = 1.0_pReal ! state is current time + thermalState(phase)%p(sourceOffset)%dotState(1,of) = 1.0_pReal ! state is current time end subroutine source_thermal_externalheat_dotState @@ -121,7 +123,7 @@ module subroutine source_thermal_externalheat_getRateAndItsTangent(TDot, dTDot_d associate(prm => param(source_thermal_externalheat_instance(phase))) do interval = 1, prm%nIntervals ! scan through all rate segments - frac_time = (sourceState(phase)%p(sourceOffset)%state(1,of) - prm%t_n(interval)) & + frac_time = (thermalState(phase)%p(sourceOffset)%state(1,of) - prm%t_n(interval)) & / (prm%t_n(interval+1) - prm%t_n(interval)) ! fractional time within segment if ( (frac_time < 0.0_pReal .and. interval == 1) & .or. (frac_time >= 1.0_pReal .and. interval == prm%nIntervals) & diff --git a/src/homogenization.f90 b/src/homogenization.f90 index df7369096..8738ba6f1 100644 --- a/src/homogenization.f90 +++ b/src/homogenization.f90 @@ -69,6 +69,13 @@ module homogenization el !< element number end subroutine mech_partition + module subroutine thermal_partition(T,ip,el) + real(pReal), intent(in) :: T + integer, intent(in) :: & + ip, & !< integration point + el !< element number + end subroutine thermal_partition + module subroutine mech_homogenize(dt,ip,el) real(pReal), intent(in) :: dt integer, intent(in) :: & @@ -131,9 +138,10 @@ subroutine homogenization_init call mech_init(num_homog) + call thermal_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(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 diff --git a/src/homogenization_thermal.f90 b/src/homogenization_thermal.f90 index 59e7357b6..f181d97fa 100644 --- a/src/homogenization_thermal.f90 +++ b/src/homogenization_thermal.f90 @@ -11,9 +11,11 @@ contains !-------------------------------------------------------------------------------------------------- module subroutine thermal_init() + print'(/,a)', ' <<<+- homogenization_thermal init -+>>>' - allocate(homogenization_T(discretization_nIPs*discretization_Nelems), source=0.0_pReal) + allocate(homogenization_T(discretization_nIPs*discretization_Nelems)) + end subroutine thermal_init diff --git a/src/thermal_conduction.f90 b/src/thermal_conduction.f90 index f98d36d3b..0cd1678e0 100644 --- a/src/thermal_conduction.f90 +++ b/src/thermal_conduction.f90 @@ -10,6 +10,7 @@ module thermal_conduction use results use constitutive use YAML_types + use discretization implicit none private @@ -38,25 +39,28 @@ contains !> @brief module initialization !> @details reads in material parameters, allocates arrays, and does sanity checks !-------------------------------------------------------------------------------------------------- -subroutine thermal_conduction_init +subroutine thermal_conduction_init(T) - integer :: Ninstances,Nmaterialpoints,h + 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 h = 1, size(material_name_homogenization) - if (thermal_type(h) /= THERMAL_conduction_ID) cycle - homog => material_homogenization%get(h) + 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(h))) + associate(prm => param(thermal_typeInstance(ho))) #if defined (__GFORTRAN__) prm%output = output_asStrings(homogThermal) @@ -64,14 +68,23 @@ subroutine thermal_conduction_init prm%output = homogThermal%get_asStrings('output',defaultVal=emptyStringArray) #endif - Nmaterialpoints=count(material_homogenizationAt==h) + Nmaterialpoints=count(material_homogenizationAt==ho) - allocate (temperature (h)%p(Nmaterialpoints), source=thermal_initialT(h)) - allocate (temperatureRate(h)%p(Nmaterialpoints), source=0.0_pReal) + 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 @@ -89,12 +102,12 @@ subroutine thermal_conduction_getSourceAndItsTangent(Tdot, dTdot_dT, T, ip, el) Tdot, dTdot_dT integer :: & homog - + Tdot = 0.0_pReal dTdot_dT = 0.0_pReal homog = material_homogenizationAt(el) - call constitutive_thermal_getRateAndItsTangents(TDot, dTDot_dT, T, ip, el) + call constitutive_thermal_getRateAndItsTangents(TDot, dTDot_dT, T, ip, el) Tdot = Tdot/real(homogenization_Nconstituents(homog),pReal) dTdot_dT = dTdot_dT/real(homogenization_Nconstituents(homog),pReal) @@ -112,13 +125,13 @@ function thermal_conduction_getConductivity(ip,el) 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))) @@ -168,7 +181,7 @@ function thermal_conduction_getMassDensity(ip,el) el !< element number real(pReal) :: & thermal_conduction_getMassDensity - + integer :: & co diff --git a/src/thermal_isothermal.f90 b/src/thermal_isothermal.f90 index 2a41ada49..09e35931e 100644 --- a/src/thermal_isothermal.f90 +++ b/src/thermal_isothermal.f90 @@ -6,6 +6,7 @@ module thermal_isothermal use prec use config use material + use discretization implicit none public @@ -15,22 +16,33 @@ contains !-------------------------------------------------------------------------------------------------- !> @brief allocates fields, reads information from material configuration file !-------------------------------------------------------------------------------------------------- -subroutine thermal_isothermal_init +subroutine thermal_isothermal_init(T) - integer :: h,Nmaterialpoints + real(pReal), dimension(:), intent(inout) :: T + + integer :: Ninstances,Nmaterialpoints,ho,ip,el,ce print'(/,a)', ' <<<+- thermal_isothermal init -+>>>'; flush(6) - do h = 1, size(material_name_homogenization) - if (thermal_type(h) /= THERMAL_isothermal_ID) cycle + do ho = 1, size(thermal_type) + if (thermal_type(ho) /= THERMAL_isothermal_ID) cycle - Nmaterialpoints = count(material_homogenizationAt == h) + Nmaterialpoints = count(material_homogenizationAt == ho) - allocate(temperature (h)%p(Nmaterialpoints),source=thermal_initialT(h)) - allocate(temperatureRate(h)%p(Nmaterialpoints),source = 0.0_pReal) + 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