not needed as global variable

This commit is contained in:
Martin Diehl 2020-12-27 16:13:31 +01:00
parent 609d69a3e7
commit 2eed6fdfdb
2 changed files with 41 additions and 60 deletions

View File

@ -64,10 +64,6 @@ module constitutive
real(pReal), dimension(:,:,:,:,:), allocatable, public :: & real(pReal), dimension(:,:,:,:,:), allocatable, public :: &
crystallite_partitionedF !< def grad to be reached at end of homog inc crystallite_partitionedF !< def grad to be reached at end of homog inc
logical, dimension(:,:,:), allocatable :: &
crystallite_converged !< convergence flag
type :: tTensorContainer type :: tTensorContainer
real(pReal), dimension(:,:,:), allocatable :: data real(pReal), dimension(:,:,:), allocatable :: data
end type end type
@ -185,10 +181,10 @@ module constitutive
! == cleaned:end =================================================================================== ! == cleaned:end ===================================================================================
module function crystallite_stress(dt,co,ip,el) module function crystallite_stress(dt,co,ip,el) result(converged_)
real(pReal), intent(in) :: dt real(pReal), intent(in) :: dt
integer, intent(in) :: co, ip, el integer, intent(in) :: co, ip, el
logical :: crystallite_stress logical :: converged_
end function crystallite_stress end function crystallite_stress
module function constitutive_homogenizedC(co,ip,el) result(C) module function constitutive_homogenizedC(co,ip,el) result(C)
@ -872,10 +868,8 @@ subroutine crystallite_init
source = crystallite_partitionedF) source = crystallite_partitionedF)
allocate(crystallite_subdt(cMax,iMax,eMax),source=0.0_pReal) allocate(crystallite_subdt(cMax,iMax,eMax),source=0.0_pReal)
allocate(crystallite_orientation(cMax,iMax,eMax)) allocate(crystallite_orientation(cMax,iMax,eMax))
allocate(crystallite_converged(cMax,iMax,eMax), source=.true.)
num_crystallite => config_numerics%get('crystallite',defaultVal=emptyDict) num_crystallite => config_numerics%get('crystallite',defaultVal=emptyDict)
@ -1253,7 +1247,7 @@ end function crystallite_push33ToRef
!> @brief integrate stress, state with adaptive 1st order explicit Euler method !> @brief integrate stress, state with adaptive 1st order explicit Euler method
!> using Fixed Point Iteration to adapt the stepsize !> using Fixed Point Iteration to adapt the stepsize
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine integrateSourceState(co,ip,el) function integrateSourceState(co,ip,el) result(broken)
integer, intent(in) :: & integer, intent(in) :: &
el, & !< element index in element loop el, & !< element index in element loop
@ -1273,12 +1267,13 @@ subroutine integrateSourceState(co,ip,el)
r ! state residuum r ! state residuum
real(pReal), dimension(constitutive_source_maxSizeDotState,2,maxval(phase_Nsources)) :: source_dotState real(pReal), dimension(constitutive_source_maxSizeDotState,2,maxval(phase_Nsources)) :: source_dotState
logical :: & logical :: &
broken broken, converged_
ph = material_phaseAt(co,el) ph = material_phaseAt(co,el)
me = material_phaseMemberAt(co,ip,el) me = material_phaseMemberAt(co,ip,el)
converged_ = .true.
broken = constitutive_thermal_collectDotState(ph,me) broken = constitutive_thermal_collectDotState(ph,me)
broken = broken .or. constitutive_damage_collectDotState(crystallite_S(1:3,1:3,co,ip,el), co,ip,el,ph,me) broken = broken .or. constitutive_damage_collectDotState(crystallite_S(1:3,1:3,co,ip,el), co,ip,el,ph,me)
if(broken) return if(broken) return
@ -1313,19 +1308,20 @@ subroutine integrateSourceState(co,ip,el)
- sourceState(ph)%p(so)%dotState (1:size_so(so),me) * crystallite_subdt(co,ip,el) - sourceState(ph)%p(so)%dotState (1:size_so(so),me) * crystallite_subdt(co,ip,el)
sourceState(ph)%p(so)%state(1:size_so(so),me) = sourceState(ph)%p(so)%state(1:size_so(so),me) & sourceState(ph)%p(so)%state(1:size_so(so),me) = sourceState(ph)%p(so)%state(1:size_so(so),me) &
- r(1:size_so(so)) - r(1:size_so(so))
crystallite_converged(co,ip,el) = & converged_ = converged_ .and. converged(r(1:size_so(so)), &
crystallite_converged(co,ip,el) .and. converged(r(1:size_so(so)), & sourceState(ph)%p(so)%state(1:size_so(so),me), &
sourceState(ph)%p(so)%state(1:size_so(so),me), & sourceState(ph)%p(so)%atol(1:size_so(so)))
sourceState(ph)%p(so)%atol(1:size_so(so)))
enddo enddo
if(crystallite_converged(co,ip,el)) then if(converged_) then
broken = constitutive_damage_deltaState(crystallite_Fe(1:3,1:3,co,ip,el),co,ip,el,ph,me) broken = constitutive_damage_deltaState(crystallite_Fe(1:3,1:3,co,ip,el),co,ip,el,ph,me)
exit iteration exit iteration
endif endif
enddo iteration enddo iteration
broken = broken .or. .not. converged_
contains contains
@ -1349,7 +1345,7 @@ subroutine integrateSourceState(co,ip,el)
end function damper end function damper
end subroutine integrateSourceState end function integrateSourceState
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------

View File

@ -951,7 +951,7 @@ end function integrateStress
!> @brief integrate stress, state with adaptive 1st order explicit Euler method !> @brief integrate stress, state with adaptive 1st order explicit Euler method
!> using Fixed Point Iteration to adapt the stepsize !> using Fixed Point Iteration to adapt the stepsize
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine integrateStateFPI(F_0,F,Delta_t,co,ip,el) function integrateStateFPI(F_0,F,Delta_t,co,ip,el) result(broken)
real(pReal), intent(in),dimension(3,3) :: F_0,F real(pReal), intent(in),dimension(3,3) :: F_0,F
real(pReal), intent(in) :: Delta_t real(pReal), intent(in) :: Delta_t
@ -1004,11 +1004,7 @@ subroutine integrateStateFPI(F_0,F,Delta_t,co,ip,el)
- plasticState(ph)%dotState (1:size_pl,me) * Delta_t - plasticState(ph)%dotState (1:size_pl,me) * Delta_t
plasticState(ph)%state(1:size_pl,me) = plasticState(ph)%state(1:size_pl,me) & plasticState(ph)%state(1:size_pl,me) = plasticState(ph)%state(1:size_pl,me) &
- r(1:size_pl) - r(1:size_pl)
crystallite_converged(co,ip,el) = converged(r(1:size_pl), & if (converged(r(1:size_pl),plasticState(ph)%state(1:size_pl,me),plasticState(ph)%atol(1:size_pl))) then
plasticState(ph)%state(1:size_pl,me), &
plasticState(ph)%atol(1:size_pl))
if(crystallite_converged(co,ip,el)) then
broken = constitutive_deltaState(crystallite_S(1:3,1:3,co,ip,el), & broken = constitutive_deltaState(crystallite_S(1:3,1:3,co,ip,el), &
constitutive_mech_Fi(ph)%data(1:3,1:3,me),co,ip,el,ph,me) constitutive_mech_Fi(ph)%data(1:3,1:3,me),co,ip,el,ph,me)
exit iteration exit iteration
@ -1016,7 +1012,6 @@ subroutine integrateStateFPI(F_0,F,Delta_t,co,ip,el)
enddo iteration enddo iteration
contains contains
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -1039,13 +1034,13 @@ subroutine integrateStateFPI(F_0,F,Delta_t,co,ip,el)
end function damper end function damper
end subroutine integrateStateFPI end function integrateStateFPI
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief integrate state with 1st order explicit Euler method !> @brief integrate state with 1st order explicit Euler method
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine integrateStateEuler(F_0,F,Delta_t,co,ip,el) function integrateStateEuler(F_0,F,Delta_t,co,ip,el) result(broken)
real(pReal), intent(in),dimension(3,3) :: F_0,F real(pReal), intent(in),dimension(3,3) :: F_0,F
real(pReal), intent(in) :: Delta_t real(pReal), intent(in) :: Delta_t
@ -1075,15 +1070,14 @@ subroutine integrateStateEuler(F_0,F,Delta_t,co,ip,el)
if(broken) return if(broken) return
broken = integrateStress(F,Delta_t,co,ip,el) broken = integrateStress(F,Delta_t,co,ip,el)
crystallite_converged(co,ip,el) = .not. broken
end subroutine integrateStateEuler end function integrateStateEuler
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief integrate stress, state with 1st order Euler method with adaptive step size !> @brief integrate stress, state with 1st order Euler method with adaptive step size
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine integrateStateAdaptiveEuler(F_0,F,Delta_t,co,ip,el) function integrateStateAdaptiveEuler(F_0,F,Delta_t,co,ip,el) result(broken)
real(pReal), intent(in),dimension(3,3) :: F_0,F real(pReal), intent(in),dimension(3,3) :: F_0,F
real(pReal), intent(in) :: Delta_t real(pReal), intent(in) :: Delta_t
@ -1123,24 +1117,22 @@ subroutine integrateStateAdaptiveEuler(F_0,F,Delta_t,co,ip,el)
broken = mech_collectDotState(Delta_t, co,ip,el,ph,me) broken = mech_collectDotState(Delta_t, co,ip,el,ph,me)
if(broken) return if(broken) return
broken = .not. converged(residuum_plastic(1:sizeDotState) + 0.5_pReal * plasticState(ph)%dotState(:,me) * Delta_t, &
plasticState(ph)%state(1:sizeDotState,me), &
plasticState(ph)%atol(1:sizeDotState))
sizeDotState = plasticState(ph)%sizeDotState end function integrateStateAdaptiveEuler
crystallite_converged(co,ip,el) = converged(residuum_plastic(1:sizeDotState) &
+ 0.5_pReal * plasticState(ph)%dotState(:,me) * Delta_t, &
plasticState(ph)%state(1:sizeDotState,me), &
plasticState(ph)%atol(1:sizeDotState))
end subroutine integrateStateAdaptiveEuler
!--------------------------------------------------------------------------------------------------- !---------------------------------------------------------------------------------------------------
!> @brief Integrate state (including stress integration) with the classic Runge Kutta method !> @brief Integrate state (including stress integration) with the classic Runge Kutta method
!--------------------------------------------------------------------------------------------------- !---------------------------------------------------------------------------------------------------
subroutine integrateStateRK4(F_0,F,Delta_t,co,ip,el) function integrateStateRK4(F_0,F,Delta_t,co,ip,el) result(broken)
real(pReal), intent(in),dimension(3,3) :: F_0,F real(pReal), intent(in),dimension(3,3) :: F_0,F
real(pReal), intent(in) :: Delta_t real(pReal), intent(in) :: Delta_t
integer, intent(in) :: co,ip,el integer, intent(in) :: co,ip,el
logical :: broken
real(pReal), dimension(3,3), parameter :: & real(pReal), dimension(3,3), parameter :: &
A = reshape([& A = reshape([&
@ -1153,19 +1145,20 @@ subroutine integrateStateRK4(F_0,F,Delta_t,co,ip,el)
real(pReal), dimension(4), parameter :: & real(pReal), dimension(4), parameter :: &
B = [1.0_pReal/6.0_pReal, 1.0_pReal/3.0_pReal, 1.0_pReal/3.0_pReal, 1.0_pReal/6.0_pReal] B = [1.0_pReal/6.0_pReal, 1.0_pReal/3.0_pReal, 1.0_pReal/3.0_pReal, 1.0_pReal/6.0_pReal]
call integrateStateRK(F_0,F,Delta_t,co,ip,el,A,B,C) broken = integrateStateRK(F_0,F,Delta_t,co,ip,el,A,B,C)
end subroutine integrateStateRK4 end function integrateStateRK4
!--------------------------------------------------------------------------------------------------- !---------------------------------------------------------------------------------------------------
!> @brief Integrate state (including stress integration) with the Cash-Carp method !> @brief Integrate state (including stress integration) with the Cash-Carp method
!--------------------------------------------------------------------------------------------------- !---------------------------------------------------------------------------------------------------
subroutine integrateStateRKCK45(F_0,F,Delta_t,co,ip,el) function integrateStateRKCK45(F_0,F,Delta_t,co,ip,el) result(broken)
real(pReal), intent(in),dimension(3,3) :: F_0,F real(pReal), intent(in),dimension(3,3) :: F_0,F
real(pReal), intent(in) :: Delta_t real(pReal), intent(in) :: Delta_t
integer, intent(in) :: co,ip,el integer, intent(in) :: co,ip,el
logical :: broken
real(pReal), dimension(5,5), parameter :: & real(pReal), dimension(5,5), parameter :: &
A = reshape([& A = reshape([&
@ -1185,16 +1178,16 @@ subroutine integrateStateRKCK45(F_0,F,Delta_t,co,ip,el)
[2825.0_pReal/27648.0_pReal, .0_pReal, 18575.0_pReal/48384.0_pReal,& [2825.0_pReal/27648.0_pReal, .0_pReal, 18575.0_pReal/48384.0_pReal,&
13525.0_pReal/55296.0_pReal, 277.0_pReal/14336.0_pReal, 1._pReal/4._pReal] 13525.0_pReal/55296.0_pReal, 277.0_pReal/14336.0_pReal, 1._pReal/4._pReal]
call integrateStateRK(F_0,F,Delta_t,co,ip,el,A,B,C,DB) broken = integrateStateRK(F_0,F,Delta_t,co,ip,el,A,B,C,DB)
end subroutine integrateStateRKCK45 end function integrateStateRKCK45
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Integrate state (including stress integration) with an explicit Runge-Kutta method or an !> @brief Integrate state (including stress integration) with an explicit Runge-Kutta method or an
!! embedded explicit Runge-Kutta method !! embedded explicit Runge-Kutta method
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine integrateStateRK(F_0,F,Delta_t,co,ip,el,A,B,C,DB) function integrateStateRK(F_0,F,Delta_t,co,ip,el,A,B,C,DB) result(broken)
real(pReal), intent(in),dimension(3,3) :: F_0,F real(pReal), intent(in),dimension(3,3) :: F_0,F
real(pReal), intent(in) :: Delta_t real(pReal), intent(in) :: Delta_t
@ -1205,15 +1198,14 @@ subroutine integrateStateRK(F_0,F,Delta_t,co,ip,el,A,B,C,DB)
el, & !< element index in element loop el, & !< element index in element loop
ip, & !< integration point index in ip loop ip, & !< integration point index in ip loop
co !< grain index in grain loop co !< grain index in grain loop
logical :: broken
integer :: & integer :: &
stage, & ! stage index in integration stage loop stage, & ! stage index in integration stage loop
n, & n, &
ph, & ph, &
me, & me, &
sizeDotState sizeDotState
logical :: &
broken
real(pReal), dimension(constitutive_plasticity_maxSizeDotState,size(B)) :: plastic_RKdotState real(pReal), dimension(constitutive_plasticity_maxSizeDotState,size(B)) :: plastic_RKdotState
@ -1266,10 +1258,8 @@ subroutine integrateStateRK(F_0,F,Delta_t,co,ip,el,A,B,C,DB)
if(broken) return if(broken) return
broken = integrateStress(F,Delta_t,co,ip,el) broken = integrateStress(F,Delta_t,co,ip,el)
crystallite_converged(co,ip,el) = .not. broken
end function integrateStateRK
end subroutine integrateStateRK
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -1479,15 +1469,14 @@ end function constitutive_homogenizedC
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief calculate stress (P) !> @brief calculate stress (P)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module function crystallite_stress(dt,co,ip,el) module function crystallite_stress(dt,co,ip,el) result(converged_)
real(pReal), intent(in) :: dt real(pReal), intent(in) :: dt
integer, intent(in) :: & integer, intent(in) :: &
co, & co, &
ip, & ip, &
el el
logical :: converged_
logical :: crystallite_stress
real(pReal) :: & real(pReal) :: &
formerSubStep formerSubStep
@ -1519,7 +1508,7 @@ module function crystallite_stress(dt,co,ip,el)
subFrac = 0.0_pReal subFrac = 0.0_pReal
subStep = 1.0_pReal/num%subStepSizeCryst subStep = 1.0_pReal/num%subStepSizeCryst
todo = .true. todo = .true.
crystallite_converged(co,ip,el) = .false. ! pretend failed step of 1/subStepSizeCryst converged_ = .false. ! pretend failed step of 1/subStepSizeCryst
todo = .true. todo = .true.
NiterationCrystallite = 0 NiterationCrystallite = 0
@ -1528,7 +1517,7 @@ module function crystallite_stress(dt,co,ip,el)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! wind forward ! wind forward
if (crystallite_converged(co,ip,el)) then if (converged_) then
formerSubStep = subStep formerSubStep = subStep
subFrac = subFrac + subStep subFrac = subFrac + subStep
subStep = min(1.0_pReal - subFrac, num%stepIncreaseCryst * subStep) subStep = min(1.0_pReal - subFrac, num%stepIncreaseCryst * subStep)
@ -1579,17 +1568,13 @@ module function crystallite_stress(dt,co,ip,el)
math_inv33(matmul(constitutive_mech_Fi(ph)%data(1:3,1:3,me), & math_inv33(matmul(constitutive_mech_Fi(ph)%data(1:3,1:3,me), &
constitutive_mech_Fp(ph)%data(1:3,1:3,me)))) constitutive_mech_Fp(ph)%data(1:3,1:3,me))))
crystallite_subdt(co,ip,el) = subStep * dt crystallite_subdt(co,ip,el) = subStep * dt
crystallite_converged(co,ip,el) = .false. converged_ = .not. integrateState(subF0,crystallite_subF(1:3,1:3,co,ip,el),&
call integrateState(subF0,crystallite_subF(1:3,1:3,co,ip,el),& crystallite_subdt(co,ip,el),co,ip,el)
crystallite_subdt(co,ip,el),co,ip,el) converged_ = converged_ .and. .not. integrateSourceState(co,ip,el)
call integrateSourceState(co,ip,el)
endif endif
enddo cutbackLooping enddo cutbackLooping
! return whether converged or not
crystallite_stress = crystallite_converged(co,ip,el)
end function crystallite_stress end function crystallite_stress
end submodule constitutive_mech end submodule constitutive_mech