223031012
/
DAMASK_EICMD
mirror of https://github.com/achalhp/DAMASK_EICMD.git
1
0
Fork 0
Code Issues Projects Releases Wiki Activity

Merge branch 'simplify-crystallite' into 'development' 

Simplify crystallite

See merge request damask/DAMASK!242
This commit is contained in:
Franz Roters 2020-10-02 14:14:40 +02:00
parent 58229b8851 bd5d557fbb
commit 365cf9a222
1 changed files with 277 additions and 343 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

204
src/crystallite.f90
View File

@ -54,12 +54,10 @@ module crystallite
! !
crystallite_Lp0, & !< plastic velocitiy grad at start of FE inc crystallite_Lp0, & !< plastic velocitiy grad at start of FE inc
crystallite_partionedLp0, & !< plastic velocity grad at start of homog inc crystallite_partionedLp0, & !< plastic velocity grad at start of homog inc
crystallite_subLp0,& !< plastic velocity grad at start of crystallite inc
! !
crystallite_Li, & !< current intermediate velocitiy grad (end of converged time step) crystallite_Li, & !< current intermediate velocitiy grad (end of converged time step)
crystallite_Li0, & !< intermediate velocitiy grad at start of FE inc crystallite_Li0, & !< intermediate velocitiy grad at start of FE inc
crystallite_partionedLi0, & !< intermediate velocity grad at start of homog inc crystallite_partionedLi0, & !< intermediate velocity grad at start of homog inc
crystallite_subLi0, & !< intermediate velocity grad at start of crystallite inc
! !
crystallite_S0, & !< 2nd Piola-Kirchhoff stress vector at start of FE inc crystallite_S0, & !< 2nd Piola-Kirchhoff stress vector at start of FE inc
crystallite_partionedS0 !< 2nd Piola-Kirchhoff stress vector at start of homog inc crystallite_partionedS0 !< 2nd Piola-Kirchhoff stress vector at start of homog inc
@ -183,7 +181,6 @@ subroutine crystallite_init
crystallite_Li,crystallite_Lp, & crystallite_Li,crystallite_Lp, &
crystallite_subF,crystallite_subF0, & crystallite_subF,crystallite_subF0, &
crystallite_subFp0,crystallite_subFi0, & crystallite_subFp0,crystallite_subFi0, &
crystallite_subLi0,crystallite_subLp0, &
source = crystallite_partionedF) source = crystallite_partionedF)
allocate(crystallite_dPdF(3,3,3,3,cMax,iMax,eMax),source=0.0_pReal) allocate(crystallite_dPdF(3,3,3,3,cMax,iMax,eMax),source=0.0_pReal)
@ -326,34 +323,17 @@ function crystallite_stress()
startIP, endIP, & startIP, endIP, &
s s
logical, dimension(homogenization_maxNgrains,discretization_nIP,discretization_nElem) :: todo !ToDo: need to set some values to false for different Ngrains logical, dimension(homogenization_maxNgrains,discretization_nIP,discretization_nElem) :: todo !ToDo: need to set some values to false for different Ngrains
real(pReal), dimension(:,:,:,:,:), allocatable :: &
subLp0,& !< plastic velocity grad at start of crystallite inc
subLi0 !< intermediate velocity grad at start of crystallite inc
todo = .false. todo = .false.
#ifdef DEBUG subLp0 = crystallite_partionedLp0
if (debugCrystallite%selective & subLi0 = crystallite_partionedLi0
.and. FEsolving_execElem(1) <= debugCrystallite%element &
.and. debugCrystallite%element <= FEsolving_execElem(2)) then
print'(/,a,i8,1x,i2,1x,i3)', '<< CRYST stress >> boundary and initial values at el ip ipc ', &
debugCrystallite%element,debugCrystallite%ip, debugCrystallite%grain
print'(a,/,3(12x,3(f14.9,1x)/))', '<< CRYST stress >> F ', &
transpose(crystallite_partionedF(1:3,1:3,debugCrystallite%grain, &
debugCrystallite%ip,debugCrystallite%element))
print'(a,/,3(12x,3(f14.9,1x)/))', '<< CRYST stress >> F0 ', &
transpose(crystallite_partionedF0(1:3,1:3,debugCrystallite%grain, &
debugCrystallite%ip,debugCrystallite%element))
print'(a,/,3(12x,3(f14.9,1x)/))', '<< CRYST stress >> Fp0', &
transpose(crystallite_partionedFp0(1:3,1:3,debugCrystallite%grain, &
debugCrystallite%ip,debugCrystallite%element))
print'(a,/,3(12x,3(f14.9,1x)/))', '<< CRYST stress >> Fi0', &
transpose(crystallite_partionedFi0(1:3,1:3,debugCrystallite%grain, &
debugCrystallite%ip,debugCrystallite%element))
print'(a,/,3(12x,3(f14.9,1x)/))', '<< CRYST stress >> Lp0', &
transpose(crystallite_partionedLp0(1:3,1:3,debugCrystallite%grain, &
debugCrystallite%ip,debugCrystallite%element))
print'(a,/,3(12x,3(f14.9,1x)/))', '<< CRYST stress >> Li0', &
transpose(crystallite_partionedLi0(1:3,1:3,debugCrystallite%grain, &
debugCrystallite%ip,debugCrystallite%element))
endif
#endif
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! initialize to starting condition ! initialize to starting condition
@ -370,9 +350,7 @@ function crystallite_stress()
sourceState(material_phaseAt(c,e))%p(s)%partionedState0(:,material_phaseMemberAt(c,i,e)) sourceState(material_phaseAt(c,e))%p(s)%partionedState0(:,material_phaseMemberAt(c,i,e))
enddo enddo
crystallite_subFp0(1:3,1:3,c,i,e) = crystallite_partionedFp0(1:3,1:3,c,i,e) crystallite_subFp0(1:3,1:3,c,i,e) = crystallite_partionedFp0(1:3,1:3,c,i,e)
crystallite_subLp0(1:3,1:3,c,i,e) = crystallite_partionedLp0(1:3,1:3,c,i,e)
crystallite_subFi0(1:3,1:3,c,i,e) = crystallite_partionedFi0(1:3,1:3,c,i,e) crystallite_subFi0(1:3,1:3,c,i,e) = crystallite_partionedFi0(1:3,1:3,c,i,e)
crystallite_subLi0(1:3,1:3,c,i,e) = crystallite_partionedLi0(1:3,1:3,c,i,e)
crystallite_subF0(1:3,1:3,c,i,e) = crystallite_partionedF0(1:3,1:3,c,i,e) crystallite_subF0(1:3,1:3,c,i,e) = crystallite_partionedF0(1:3,1:3,c,i,e)
crystallite_subFrac(c,i,e) = 0.0_pReal crystallite_subFrac(c,i,e) = 0.0_pReal
crystallite_subStep(c,i,e) = 1.0_pReal/num%subStepSizeCryst crystallite_subStep(c,i,e) = 1.0_pReal/num%subStepSizeCryst
@ -415,8 +393,8 @@ function crystallite_stress()
todo(c,i,e) = crystallite_subStep(c,i,e) > 0.0_pReal ! still time left to integrate on? todo(c,i,e) = crystallite_subStep(c,i,e) > 0.0_pReal ! still time left to integrate on?
if (todo(c,i,e)) then if (todo(c,i,e)) then
crystallite_subF0 (1:3,1:3,c,i,e) = crystallite_subF(1:3,1:3,c,i,e) crystallite_subF0 (1:3,1:3,c,i,e) = crystallite_subF(1:3,1:3,c,i,e)
crystallite_subLp0(1:3,1:3,c,i,e) = crystallite_Lp (1:3,1:3,c,i,e) subLp0(1:3,1:3,c,i,e) = crystallite_Lp (1:3,1:3,c,i,e)
crystallite_subLi0(1:3,1:3,c,i,e) = crystallite_Li (1:3,1:3,c,i,e) subLi0(1:3,1:3,c,i,e) = crystallite_Li (1:3,1:3,c,i,e)
crystallite_subFp0(1:3,1:3,c,i,e) = crystallite_Fp (1:3,1:3,c,i,e) crystallite_subFp0(1:3,1:3,c,i,e) = crystallite_Fp (1:3,1:3,c,i,e)
crystallite_subFi0(1:3,1:3,c,i,e) = crystallite_Fi (1:3,1:3,c,i,e) crystallite_subFi0(1:3,1:3,c,i,e) = crystallite_Fi (1:3,1:3,c,i,e)
plasticState( material_phaseAt(c,e))%subState0(:,material_phaseMemberAt(c,i,e)) & plasticState( material_phaseAt(c,e))%subState0(:,material_phaseMemberAt(c,i,e)) &
@ -435,8 +413,8 @@ function crystallite_stress()
crystallite_Fi (1:3,1:3,c,i,e) = crystallite_subFi0(1:3,1:3,c,i,e) crystallite_Fi (1:3,1:3,c,i,e) = crystallite_subFi0(1:3,1:3,c,i,e)
crystallite_S (1:3,1:3,c,i,e) = crystallite_S0 (1:3,1:3,c,i,e) crystallite_S (1:3,1:3,c,i,e) = crystallite_S0 (1:3,1:3,c,i,e)
if (crystallite_subStep(c,i,e) < 1.0_pReal) then ! actual (not initial) cutback if (crystallite_subStep(c,i,e) < 1.0_pReal) then ! actual (not initial) cutback
crystallite_Lp (1:3,1:3,c,i,e) = crystallite_subLp0(1:3,1:3,c,i,e) crystallite_Lp (1:3,1:3,c,i,e) = subLp0(1:3,1:3,c,i,e)
crystallite_Li (1:3,1:3,c,i,e) = crystallite_subLi0(1:3,1:3,c,i,e) crystallite_Li (1:3,1:3,c,i,e) = subLi0(1:3,1:3,c,i,e)
endif endif
plasticState (material_phaseAt(c,e))%state( :,material_phaseMemberAt(c,i,e)) & plasticState (material_phaseAt(c,e))%state( :,material_phaseMemberAt(c,i,e)) &
= plasticState(material_phaseAt(c,e))%subState0(:,material_phaseMemberAt(c,i,e)) = plasticState(material_phaseAt(c,e))%subState0(:,material_phaseMemberAt(c,i,e))
@ -460,6 +438,7 @@ function crystallite_stress()
math_inv33(crystallite_Fi(1:3,1:3,c,i,e))) math_inv33(crystallite_Fi(1:3,1:3,c,i,e)))
crystallite_subdt(c,i,e) = crystallite_subStep(c,i,e) * crystallite_dt(c,i,e) crystallite_subdt(c,i,e) = crystallite_subStep(c,i,e) * crystallite_dt(c,i,e)
crystallite_converged(c,i,e) = .false. crystallite_converged(c,i,e) = .false.
call integrateState(c,i,e)
endif endif
enddo enddo
@ -467,9 +446,10 @@ function crystallite_stress()
enddo elementLooping3 enddo elementLooping3
!$OMP END PARALLEL DO !$OMP END PARALLEL DO
call nonlocalConvergenceCheck
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! integrate --- requires fully defined state array (basic + dependent state) ! integrate --- requires fully defined state array (basic + dependent state)
if (any(todo)) call integrateState(todo) ! TODO: unroll into proper elementloop to avoid N^2 for single point evaluation
where(.not. crystallite_converged .and. crystallite_subStep > num%subStepMinCryst) & ! do not try non-converged but fully cutbacked any further where(.not. crystallite_converged .and. crystallite_subStep > num%subStepMinCryst) & ! do not try non-converged but fully cutbacked any further
todo = .true. ! TODO: again unroll this into proper elementloop to avoid N^2 for single point evaluation todo = .true. ! TODO: again unroll this into proper elementloop to avoid N^2 for single point evaluation
@ -1125,14 +1105,14 @@ 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(todo) subroutine integrateStateFPI(g,i,e)
logical, dimension(:,:,:), intent(in) :: todo integer, intent(in) :: &
integer :: &
NiterationState, & !< number of iterations in state loop
e, & !< element index in element loop e, & !< element index in element loop
i, & !< integration point index in ip loop i, & !< integration point index in ip loop
g, & !< grain index in grain loop g !< grain index in grain loop
integer :: &
NiterationState, & !< number of iterations in state loop
p, & p, &
c, & c, &
s, & s, &
@ -1147,16 +1127,9 @@ subroutine integrateStateFPI(todo)
plastic_dotState plastic_dotState
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 :: &
nonlocalBroken, broken broken
nonlocalBroken = .false.
!$OMP PARALLEL DO PRIVATE(size_pl,size_so,r,zeta,p,c,plastic_dotState,source_dotState,broken)
do e = FEsolving_execElem(1),FEsolving_execElem(2)
do i = FEsolving_execIP(1),FEsolving_execIP(2)
do g = 1,homogenization_Ngrains(material_homogenizationAt(e))
p = material_phaseAt(g,e) p = material_phaseAt(g,e)
if(todo(g,i,e) .and. .not. (nonlocalBroken .and. plasticState(p)%nonlocal)) then
c = material_phaseMemberAt(g,i,e) c = material_phaseMemberAt(g,i,e)
broken = constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & broken = constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), &
@ -1164,8 +1137,7 @@ subroutine integrateStateFPI(todo)
crystallite_Fi(1:3,1:3,g,i,e), & crystallite_Fi(1:3,1:3,g,i,e), &
crystallite_partionedFp0, & crystallite_partionedFp0, &
crystallite_subdt(g,i,e), g,i,e,p,c) crystallite_subdt(g,i,e), g,i,e,p,c)
if(broken .and. plasticState(p)%nonlocal) nonlocalBroken = .true. if(broken) return
if(broken) cycle
size_pl = plasticState(p)%sizeDotState size_pl = plasticState(p)%sizeDotState
plasticState(p)%state(1:size_pl,c) = plasticState(p)%subState0(1:size_pl,c) & plasticState(p)%state(1:size_pl,c) = plasticState(p)%subState0(1:size_pl,c) &
@ -1236,12 +1208,7 @@ subroutine integrateStateFPI(todo)
endif endif
enddo iteration enddo iteration
if(broken .and. plasticState(p)%nonlocal) nonlocalBroken = .true.
endif
enddo; enddo; enddo
!$OMP END PARALLEL DO
if (nonlocalBroken) call nonlocalConvergenceCheck
contains contains
@ -1271,29 +1238,21 @@ end subroutine integrateStateFPI
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief integrate state with 1st order explicit Euler method !> @brief integrate state with 1st order explicit Euler method
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine integrateStateEuler(todo) subroutine integrateStateEuler(g,i,e)
logical, dimension(:,:,:), intent(in) :: todo integer, intent(in) :: &
integer :: &
e, & !< element index in element loop e, & !< element index in element loop
i, & !< integration point index in ip loop i, & !< integration point index in ip loop
g, & !< grain index in grain loop g !< grain index in grain loop
integer :: &
p, & p, &
c, & c, &
s, & s, &
sizeDotState sizeDotState
logical :: & logical :: &
nonlocalBroken, broken broken
nonlocalBroken = .false.
!$OMP PARALLEL DO PRIVATE (sizeDotState,p,c,broken)
do e = FEsolving_execElem(1),FEsolving_execElem(2)
do i = FEsolving_execIP(1),FEsolving_execIP(2)
do g = 1,homogenization_Ngrains(material_homogenizationAt(e))
p = material_phaseAt(g,e) p = material_phaseAt(g,e)
if(todo(g,i,e) .and. .not. (nonlocalBroken .and. plasticState(p)%nonlocal)) then
c = material_phaseMemberAt(g,i,e) c = material_phaseMemberAt(g,i,e)
broken = constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & broken = constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), &
@ -1301,8 +1260,7 @@ subroutine integrateStateEuler(todo)
crystallite_Fi(1:3,1:3,g,i,e), & crystallite_Fi(1:3,1:3,g,i,e), &
crystallite_partionedFp0, & crystallite_partionedFp0, &
crystallite_subdt(g,i,e), g,i,e,p,c) crystallite_subdt(g,i,e), g,i,e,p,c)
if(broken .and. plasticState(p)%nonlocal) nonlocalBroken = .true. if(broken) return
if(broken) cycle
sizeDotState = plasticState(p)%sizeDotState sizeDotState = plasticState(p)%sizeDotState
plasticState(p)%state(1:sizeDotState,c) = plasticState(p)%subState0(1:sizeDotState,c) & plasticState(p)%state(1:sizeDotState,c) = plasticState(p)%subState0(1:sizeDotState,c) &
@ -1318,17 +1276,10 @@ subroutine integrateStateEuler(todo)
broken = constitutive_deltaState(crystallite_S(1:3,1:3,g,i,e), & broken = constitutive_deltaState(crystallite_S(1:3,1:3,g,i,e), &
crystallite_Fe(1:3,1:3,g,i,e), & crystallite_Fe(1:3,1:3,g,i,e), &
crystallite_Fi(1:3,1:3,g,i,e),g,i,e,p,c) crystallite_Fi(1:3,1:3,g,i,e),g,i,e,p,c)
if(broken .and. plasticState(p)%nonlocal) nonlocalBroken = .true. if(broken) return
if(broken) cycle
broken = integrateStress(g,i,e) broken = integrateStress(g,i,e)
if(broken .and. plasticState(p)%nonlocal) nonlocalBroken = .true.
crystallite_converged(g,i,e) = .not. broken crystallite_converged(g,i,e) = .not. broken
endif
enddo; enddo; enddo
!$OMP END PARALLEL DO
if (nonlocalBroken) call nonlocalConvergenceCheck
end subroutine integrateStateEuler end subroutine integrateStateEuler
@ -1336,33 +1287,25 @@ end subroutine 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(todo) subroutine integrateStateAdaptiveEuler(g,i,e)
logical, dimension(:,:,:), intent(in) :: todo
integer, intent(in) :: &
e, & !< element index in element loop
i, & !< integration point index in ip loop
g !< grain index in grain loop
integer :: & integer :: &
e, & ! element index in element loop
i, & ! integration point index in ip loop
g, & ! grain index in grain loop
p, & p, &
c, & c, &
s, & s, &
sizeDotState sizeDotState
logical :: & logical :: &
nonlocalBroken, broken broken
real(pReal), dimension(constitutive_plasticity_maxSizeDotState) :: residuum_plastic real(pReal), dimension(constitutive_plasticity_maxSizeDotState) :: residuum_plastic
real(pReal), dimension(constitutive_source_maxSizeDotState,maxval(phase_Nsources)) :: residuum_source real(pReal), dimension(constitutive_source_maxSizeDotState,maxval(phase_Nsources)) :: residuum_source
nonlocalBroken = .false.
!$OMP PARALLEL DO PRIVATE(sizeDotState,p,c,residuum_plastic,residuum_source,broken)
do e = FEsolving_execElem(1),FEsolving_execElem(2)
do i = FEsolving_execIP(1),FEsolving_execIP(2)
do g = 1,homogenization_Ngrains(material_homogenizationAt(e))
broken = .false.
p = material_phaseAt(g,e)
if(todo(g,i,e) .and. .not. (nonlocalBroken .and. plasticState(p)%nonlocal)) then
p = material_phaseAt(g,e)
c = material_phaseMemberAt(g,i,e) c = material_phaseMemberAt(g,i,e)
broken = constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & broken = constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), &
@ -1370,7 +1313,7 @@ subroutine integrateStateAdaptiveEuler(todo)
crystallite_Fi(1:3,1:3,g,i,e), & crystallite_Fi(1:3,1:3,g,i,e), &
crystallite_partionedFp0, & crystallite_partionedFp0, &
crystallite_subdt(g,i,e), g,i,e,p,c) crystallite_subdt(g,i,e), g,i,e,p,c)
if(broken) cycle if(broken) return
sizeDotState = plasticState(p)%sizeDotState sizeDotState = plasticState(p)%sizeDotState
@ -1389,17 +1332,17 @@ subroutine integrateStateAdaptiveEuler(todo)
broken = constitutive_deltaState(crystallite_S(1:3,1:3,g,i,e), & broken = constitutive_deltaState(crystallite_S(1:3,1:3,g,i,e), &
crystallite_Fe(1:3,1:3,g,i,e), & crystallite_Fe(1:3,1:3,g,i,e), &
crystallite_Fi(1:3,1:3,g,i,e),g,i,e,p,c) crystallite_Fi(1:3,1:3,g,i,e),g,i,e,p,c)
if(broken) cycle if(broken) return
broken = integrateStress(g,i,e) broken = integrateStress(g,i,e)
if(broken) cycle if(broken) return
broken = constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & broken = constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), &
crystallite_partionedF0, & crystallite_partionedF0, &
crystallite_Fi(1:3,1:3,g,i,e), & crystallite_Fi(1:3,1:3,g,i,e), &
crystallite_partionedFp0, & crystallite_partionedFp0, &
crystallite_subdt(g,i,e), g,i,e,p,c) crystallite_subdt(g,i,e), g,i,e,p,c)
if(broken) cycle if(broken) return
sizeDotState = plasticState(p)%sizeDotState sizeDotState = plasticState(p)%sizeDotState
@ -1417,22 +1360,15 @@ subroutine integrateStateAdaptiveEuler(todo)
sourceState(p)%p(s)%atol(1:sizeDotState)) sourceState(p)%p(s)%atol(1:sizeDotState))
enddo enddo
endif
if(broken .and. plasticState(p)%nonlocal) nonlocalBroken = .true.
enddo; enddo; enddo
!$OMP END PARALLEL DO
if (nonlocalBroken) call nonlocalConvergenceCheck
end subroutine integrateStateAdaptiveEuler 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(todo) subroutine integrateStateRK4(g,i,e)
logical, dimension(:,:,:), intent(in) :: todo integer, intent(in) :: g,i,e
real(pReal), dimension(3,3), parameter :: & real(pReal), dimension(3,3), parameter :: &
A = reshape([& A = reshape([&
@ -1445,7 +1381,7 @@ subroutine integrateStateRK4(todo)
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(todo,A,B,C) call integrateStateRK(g,i,e,A,B,C)
end subroutine integrateStateRK4 end subroutine integrateStateRK4
@ -1453,9 +1389,9 @@ end subroutine 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(todo) subroutine integrateStateRKCK45(g,i,e)
logical, dimension(:,:,:), intent(in) :: todo integer, intent(in) :: g,i,e
real(pReal), dimension(5,5), parameter :: & real(pReal), dimension(5,5), parameter :: &
A = reshape([& A = reshape([&
@ -1475,7 +1411,7 @@ subroutine integrateStateRKCK45(todo)
[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(todo,A,B,C,DB) call integrateStateRK(g,i,e,A,B,C,DB)
end subroutine integrateStateRKCK45 end subroutine integrateStateRKCK45
@ -1484,18 +1420,18 @@ end subroutine 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(todo,A,B,CC,DB) subroutine integrateStateRK(g,i,e,A,B,CC,DB)
logical, dimension(:,:,:), intent(in) :: todo
real(pReal), dimension(:,:), intent(in) :: A real(pReal), dimension(:,:), intent(in) :: A
real(pReal), dimension(:), intent(in) :: B, CC real(pReal), dimension(:), intent(in) :: B, CC
real(pReal), dimension(:), intent(in), optional :: DB real(pReal), dimension(:), intent(in), optional :: DB
integer, intent(in) :: &
e, & !< element index in element loop
i, & !< integration point index in ip loop
g !< grain index in grain loop
integer :: & integer :: &
e, & ! element index in element loop
i, & ! integration point index in ip loop
g, & ! grain index in grain loop
stage, & ! stage index in integration stage loop stage, & ! stage index in integration stage loop
n, & n, &
p, & p, &
@ -1503,19 +1439,11 @@ subroutine integrateStateRK(todo,A,B,CC,DB)
s, & s, &
sizeDotState sizeDotState
logical :: & logical :: &
nonlocalBroken, broken broken
real(pReal), dimension(constitutive_source_maxSizeDotState,size(B),maxval(phase_Nsources)) :: source_RKdotState real(pReal), dimension(constitutive_source_maxSizeDotState,size(B),maxval(phase_Nsources)) :: source_RKdotState
real(pReal), dimension(constitutive_plasticity_maxSizeDotState,size(B)) :: plastic_RKdotState real(pReal), dimension(constitutive_plasticity_maxSizeDotState,size(B)) :: plastic_RKdotState
nonlocalBroken = .false.
!$OMP PARALLEL DO PRIVATE(sizeDotState,p,c,plastic_RKdotState,source_RKdotState,broken)
do e = FEsolving_execElem(1),FEsolving_execElem(2)
do i = FEsolving_execIP(1),FEsolving_execIP(2)
do g = 1,homogenization_Ngrains(material_homogenizationAt(e))
broken = .false.
p = material_phaseAt(g,e) p = material_phaseAt(g,e)
if(todo(g,i,e) .and. .not. (nonlocalBroken .and. plasticState(p)%nonlocal)) then
c = material_phaseMemberAt(g,i,e) c = material_phaseMemberAt(g,i,e)
broken = constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & broken = constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), &
@ -1523,7 +1451,7 @@ subroutine integrateStateRK(todo,A,B,CC,DB)
crystallite_Fi(1:3,1:3,g,i,e), & crystallite_Fi(1:3,1:3,g,i,e), &
crystallite_partionedFp0, & crystallite_partionedFp0, &
crystallite_subdt(g,i,e), g,i,e,p,c) crystallite_subdt(g,i,e), g,i,e,p,c)
if(broken) cycle if(broken) return
do stage = 1,size(A,1) do stage = 1,size(A,1)
sizeDotState = plasticState(p)%sizeDotState sizeDotState = plasticState(p)%sizeDotState
@ -1568,7 +1496,7 @@ subroutine integrateStateRK(todo,A,B,CC,DB)
if(broken) exit if(broken) exit
enddo enddo
if(broken) cycle if(broken) return
sizeDotState = plasticState(p)%sizeDotState sizeDotState = plasticState(p)%sizeDotState
@ -1597,22 +1525,16 @@ subroutine integrateStateRK(todo,A,B,CC,DB)
sourceState(p)%p(s)%state(1:sizeDotState,c), & sourceState(p)%p(s)%state(1:sizeDotState,c), &
sourceState(p)%p(s)%atol(1:sizeDotState)) sourceState(p)%p(s)%atol(1:sizeDotState))
enddo enddo
if(broken) cycle if(broken) return
broken = constitutive_deltaState(crystallite_S(1:3,1:3,g,i,e), & broken = constitutive_deltaState(crystallite_S(1:3,1:3,g,i,e), &
crystallite_Fe(1:3,1:3,g,i,e), & crystallite_Fe(1:3,1:3,g,i,e), &
crystallite_Fi(1:3,1:3,g,i,e),g,i,e,p,c) crystallite_Fi(1:3,1:3,g,i,e),g,i,e,p,c)
if(broken) cycle if(broken) return
broken = integrateStress(g,i,e) broken = integrateStress(g,i,e)
crystallite_converged(g,i,e) = .not. broken crystallite_converged(g,i,e) = .not. broken
endif
if(broken .and. plasticState(p)%nonlocal) nonlocalBroken = .true.
enddo; enddo; enddo
!$OMP END PARALLEL DO
if(nonlocalBroken) call nonlocalConvergenceCheck
end subroutine integrateStateRK end subroutine integrateStateRK
@ -1624,7 +1546,19 @@ end subroutine integrateStateRK
subroutine nonlocalConvergenceCheck subroutine nonlocalConvergenceCheck
integer :: e,i,p integer :: e,i,p
logical :: nonlocal_broken
nonlocal_broken = .false.
!$OMP PARALLEL DO PRIVATE(p)
do e = FEsolving_execElem(1),FEsolving_execElem(2)
p = material_phaseAt(1,e)
do i = FEsolving_execIP(1),FEsolving_execIP(2)
if(plasticState(p)%nonlocal .and. .not. crystallite_converged(1,i,e)) nonlocal_broken = .true.
enddo
enddo
!$OMP END PARALLEL DO
if(.not. nonlocal_broken) return
!$OMP PARALLEL DO PRIVATE(p) !$OMP PARALLEL DO PRIVATE(p)
do e = FEsolving_execElem(1),FEsolving_execElem(2) do e = FEsolving_execElem(1),FEsolving_execElem(2)
p = material_phaseAt(1,e) p = material_phaseAt(1,e)