Merge branch 'FPI-one-loop' into 'development'

FPI in one loop

See merge request damask/DAMASK!144
This commit is contained in:
Franz Roters 2020-03-25 22:11:35 +01:00
commit b6103cc7cb
1 changed files with 220 additions and 271 deletions

View File

@ -1003,184 +1003,144 @@ end function integrateStress
!--------------------------------------------------------------------------------------------------
subroutine integrateStateFPI
integer :: &
NiterationState, & !< number of iterations in state loop
e, & !< element index in element loop
i, & !< integration point index in ip loop
g, & !< grain index in grain loop
p, &
c, &
s, &
sizeDotState
real(pReal) :: &
zeta
real(pReal), dimension(constitutive_plasticity_maxSizeDotState) :: &
residuum_plastic ! residuum for plastic state
real(pReal), dimension(constitutive_source_maxSizeDotState) :: &
residuum_source ! residuum for source state
logical :: &
doneWithIntegration
integer :: &
NiterationState, & !< number of iterations in state loop
e, & !< element index in element loop
i, & !< integration point index in ip loop
g, & !< grain index in grain loop
p, &
c, &
s, &
sizeDotState
real(pReal) :: &
zeta
real(pReal), dimension(constitutive_plasticity_maxSizeDotState) :: &
residuum_plastic ! residuum for plastic state
real(pReal), dimension(constitutive_source_maxSizeDotState) :: &
residuum_source ! residuum for source state
logical :: &
nonlocalBroken
! --+>> PREGUESS FOR STATE <<+--
call update_dotState(1.0_pReal)
call update_state(1.0_pReal)
! --+>> PREGUESS FOR STATE <<+--
call update_dotState(1.0_pReal)
call update_state(1.0_pReal)
NiterationState = 0
doneWithIntegration = .false.
crystalliteLooping: do while (.not. doneWithIntegration .and. NiterationState < num%nState)
NiterationState = NiterationState + 1
nonlocalBroken = .false.
!$OMP PARALLEL DO PRIVATE(sizeDotState,residuum_plastic,residuum_source,zeta,p,c)
do e = FEsolving_execElem(1),FEsolving_execElem(2)
do i = FEsolving_execIP(1),FEsolving_execIP(2)
do g = 1,homogenization_Ngrains(material_homogenizationAt(e))
if(crystallite_todo(g,i,e) .and. .not. crystallite_converged(g,i,e) .and. &
(.not. nonlocalBroken .or. crystallite_localPlasticity(g,i,e)) ) then
#ifdef DEBUG
if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0) &
write(6,'(a,i6)') '<< CRYST stateFPI >> state iteration ',NiterationState
#endif
iteration: do NiterationState = 1, num%nState
! store previousDotState and previousDotState2
p = material_phaseAt(g,e); c = material_phaseMemberAt(g,i,e)
!$OMP PARALLEL DO PRIVATE(p,c)
do e = FEsolving_execElem(1),FEsolving_execElem(2)
do i = FEsolving_execIP(1),FEsolving_execIP(2)
do g = 1,homogenization_Ngrains(material_homogenizationAt(e))
if (crystallite_todo(g,i,e) .and. .not. crystallite_converged(g,i,e)) then
p = material_phaseAt(g,e); c = material_phaseMemberAt(g,i,e)
plasticState(p)%previousDotState2(:,c) = merge(plasticState(p)%previousDotState(:,c),&
0.0_pReal,&
NiterationState > 1)
plasticState(p)%previousDotState (:,c) = plasticState(p)%dotState(:,c)
do s = 1, phase_Nsources(p)
sourceState(p)%p(s)%previousDotState2(:,c) = merge(sourceState(p)%p(s)%previousDotState(:,c),&
0.0_pReal, &
NiterationState > 1)
sourceState(p)%p(s)%previousDotState (:,c) = sourceState(p)%p(s)%dotState(:,c)
enddo
plasticState(p)%previousDotState2(:,c) = merge(plasticState(p)%previousDotState(:,c),&
0.0_pReal,&
NiterationState > 1)
plasticState(p)%previousDotState (:,c) = plasticState(p)%dotState(:,c)
do s = 1, phase_Nsources(p)
sourceState(p)%p(s)%previousDotState2(:,c) = merge(sourceState(p)%p(s)%previousDotState(:,c),&
0.0_pReal, &
NiterationState > 1)
sourceState(p)%p(s)%previousDotState (:,c) = sourceState(p)%p(s)%dotState(:,c)
enddo
endif
enddo
enddo
enddo
!$OMP END PARALLEL DO
call constitutive_dependentState(crystallite_partionedF(1:3,1:3,g,i,e), &
crystallite_Fp(1:3,1:3,g,i,e), &
g, i, e)
call update_dependentState
call update_stress(1.0_pReal)
call update_dotState(1.0_pReal)
crystallite_todo(g,i,e) = integrateStress(g,i,e,1.0_pReal)
if(.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) &
nonlocalBroken = .true.
if(.not. crystallite_todo(g,i,e)) exit iteration
!$OMP PARALLEL
!$OMP DO PRIVATE(sizeDotState,residuum_plastic,residuum_source,zeta,p,c)
do e = FEsolving_execElem(1),FEsolving_execElem(2)
do i = FEsolving_execIP(1),FEsolving_execIP(2)
do g = 1,homogenization_Ngrains(material_homogenizationAt(e))
if (crystallite_todo(g,i,e) .and. .not. crystallite_converged(g,i,e)) then
p = material_phaseAt(g,e); c = material_phaseMemberAt(g,i,e)
sizeDotState = plasticState(p)%sizeDotState
call constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), &
crystallite_partionedF0, &
crystallite_Fi(1:3,1:3,g,i,e), &
crystallite_partionedFp0, &
crystallite_subdt(g,i,e), g,i,e)
crystallite_todo(g,i,e) = all(.not. IEEE_is_NaN(plasticState(p)%dotState(:,c)))
do s = 1, phase_Nsources(p)
crystallite_todo(g,i,e) = crystallite_todo(g,i,e) .and. all(.not. IEEE_is_NaN(sourceState(p)%p(s)%dotState(:,c)))
enddo
if(.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) &
nonlocalBroken = .true.
if(.not. crystallite_todo(g,i,e)) exit iteration
zeta = damper(plasticState(p)%dotState (:,c), &
plasticState(p)%previousDotState (:,c), &
plasticState(p)%previousDotState2(:,c))
sizeDotState = plasticState(p)%sizeDotState
zeta = damper(plasticState(p)%dotState (:,c), &
plasticState(p)%previousDotState (:,c), &
plasticState(p)%previousDotState2(:,c))
plasticState(p)%dotState(:,c) = plasticState(p)%dotState(:,c) * zeta &
+ plasticState(p)%previousDotState(:,c) * (1.0_pReal - zeta)
residuum_plastic(1:SizeDotState) = plasticState(p)%state (1:sizeDotState,c) &
- plasticState(p)%subState0(1:sizeDotState,c) &
- plasticState(p)%dotState (1:sizeDotState,c) &
* crystallite_subdt(g,i,e)
plasticState(p)%state(1:sizeDotState,c) = plasticState(p)%state(1:sizeDotState,c) &
- residuum_plastic(1:sizeDotState)
crystallite_converged(g,i,e) = converged(residuum_plastic(1:sizeDotState), &
plasticState(p)%state(1:sizeDotState,c), &
plasticState(p)%atol(1:sizeDotState))
do s = 1, phase_Nsources(p)
sizeDotState = sourceState(p)%p(s)%sizeDotState
zeta = damper(sourceState(p)%p(s)%dotState (:,c), &
sourceState(p)%p(s)%previousDotState (:,c), &
sourceState(p)%p(s)%previousDotState2(:,c))
sourceState(p)%p(s)%dotState(:,c) = sourceState(p)%p(s)%dotState(:,c) * zeta &
+ sourceState(p)%p(s)%previousDotState(:,c)* (1.0_pReal - zeta)
residuum_source(1:sizeDotState) = sourceState(p)%p(s)%state (1:sizeDotState,c) &
- sourceState(p)%p(s)%subState0(1:sizeDotState,c) &
- sourceState(p)%p(s)%dotState (1:sizeDotState,c) &
* crystallite_subdt(g,i,e)
sourceState(p)%p(s)%state(1:sizeDotState,c) = sourceState(p)%p(s)%state(1:sizeDotState,c) &
- residuum_source(1:sizeDotState)
crystallite_converged(g,i,e) = &
crystallite_converged(g,i,e) .and. converged(residuum_source(1:sizeDotState), &
sourceState(p)%p(s)%state(1:sizeDotState,c), &
sourceState(p)%p(s)%atol(1:sizeDotState))
enddo
residuum_plastic(1:SizeDotState) = plasticState(p)%state (1:sizeDotState,c) &
- plasticState(p)%subState0(1:sizeDotState,c) &
- ( plasticState(p)%dotState (:,c) * zeta &
+ plasticState(p)%previousDotState(:,c) * (1.0_pReal-zeta) &
) * crystallite_subdt(g,i,e)
if(crystallite_converged(g,i,e)) then
crystallite_todo(g,i,e) = stateJump(g,i,e)
if(.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) &
nonlocalBroken = .true.
exit iteration
endif
plasticState(p)%state(1:sizeDotState,c) = plasticState(p)%state(1:sizeDotState,c) &
- residuum_plastic(1:sizeDotState)
plasticState(p)%dotState(:,c) = plasticState(p)%dotState(:,c) * zeta &
+ plasticState(p)%previousDotState(:,c) * (1.0_pReal - zeta)
enddo iteration
crystallite_converged(g,i,e) = converged(residuum_plastic(1:sizeDotState), &
plasticState(p)%state(1:sizeDotState,c), &
plasticState(p)%atol(1:sizeDotState))
endif
enddo; enddo; enddo
!$OMP END PARALLEL DO
if(nonlocalBroken) where(.not. crystallite_localPlasticity) crystallite_todo = .false.
if (any(plasticState(:)%nonlocal)) call nonlocalConvergenceCheck
do s = 1, phase_Nsources(p)
sizeDotState = sourceState(p)%p(s)%sizeDotState
contains
zeta = damper(sourceState(p)%p(s)%dotState (:,c), &
sourceState(p)%p(s)%previousDotState (:,c), &
sourceState(p)%p(s)%previousDotState2(:,c))
!--------------------------------------------------------------------------------------------------
!> @brief calculate the damping for correction of state and dot state
!--------------------------------------------------------------------------------------------------
real(pReal) pure function damper(current,previous,previous2)
residuum_source(1:sizeDotState) = sourceState(p)%p(s)%state (1:sizeDotState,c) &
- sourceState(p)%p(s)%subState0(1:sizeDotState,c) &
- ( sourceState(p)%p(s)%dotState (:,c) * zeta &
+ sourceState(p)%p(s)%previousDotState(:,c) * (1.0_pReal - zeta) &
) * crystallite_subdt(g,i,e)
real(pReal), dimension(:), intent(in) ::&
current, previous, previous2
sourceState(p)%p(s)%state(1:sizeDotState,c) = sourceState(p)%p(s)%state(1:sizeDotState,c) &
- residuum_source(1:sizeDotState)
sourceState(p)%p(s)%dotState(:,c) = sourceState(p)%p(s)%dotState(:,c) * zeta &
+ sourceState(p)%p(s)%previousDotState(:,c)* (1.0_pReal - zeta)
real(pReal) :: dot_prod12, dot_prod22
crystallite_converged(g,i,e) = &
crystallite_converged(g,i,e) .and. converged(residuum_source(1:sizeDotState), &
sourceState(p)%p(s)%state(1:sizeDotState,c), &
sourceState(p)%p(s)%atol(1:sizeDotState))
enddo
endif
enddo; enddo; enddo
!$OMP ENDDO
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
!$OMP DO
do e = FEsolving_execElem(1),FEsolving_execElem(2)
do i = FEsolving_execIP(1),FEsolving_execIP(2)
do g = 1,homogenization_Ngrains(material_homogenizationAt(e))
!$OMP FLUSH(crystallite_todo)
if (crystallite_todo(g,i,e) .and. crystallite_converged(g,i,e)) then ! converged and still alive...
crystallite_todo(g,i,e) = stateJump(g,i,e)
!$OMP FLUSH(crystallite_todo)
if (.not. crystallite_todo(g,i,e)) then ! if state jump fails, then convergence is broken
crystallite_converged(g,i,e) = .false.
if (.not. crystallite_localPlasticity(g,i,e)) then ! if broken non-local...
!$OMP CRITICAL (checkTodo)
crystallite_todo = crystallite_todo .and. crystallite_localPlasticity ! ...all non-locals skipped
!$OMP END CRITICAL (checkTodo)
endif
endif
endif
enddo; enddo; enddo
!$OMP ENDDO
!$OMP END PARALLEL
if (any(plasticState(:)%nonlocal)) call nonlocalConvergenceCheck
! --- CHECK IF DONE WITH INTEGRATION ---
doneWithIntegration = .true.
do e = FEsolving_execElem(1),FEsolving_execElem(2)
do i = FEsolving_execIP(1),FEsolving_execIP(2)
do g = 1,homogenization_Ngrains(material_homogenizationAt(e))
if (crystallite_todo(g,i,e) .and. .not. crystallite_converged(g,i,e)) then
doneWithIntegration = .false.
exit
endif
enddo; enddo
enddo
enddo crystalliteLooping
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 damper
end subroutine integrateStateFPI
@ -1578,32 +1538,34 @@ end subroutine setConvergenceFlag
!--------------------------------------------------------------------------------------------------
subroutine update_stress(timeFraction)
real(pReal), intent(in) :: &
timeFraction
integer :: &
e, & !< element index in element loop
i, & !< integration point index in ip loop
g
real(pReal), intent(in) :: &
timeFraction
integer :: &
e, & !< element index in element loop
i, & !< integration point index in ip loop
g
logical :: &
nonlocalBroken
!$OMP PARALLEL DO
do e = FEsolving_execElem(1),FEsolving_execElem(2)
do i = FEsolving_execIP(1),FEsolving_execIP(2)
do g = 1,homogenization_Ngrains(material_homogenizationAt(e))
!$OMP FLUSH(crystallite_todo)
if (crystallite_todo(g,i,e) .and. .not. crystallite_converged(g,i,e)) then
crystallite_todo(g,i,e) = integrateStress(g,i,e,timeFraction)
!$OMP FLUSH(crystallite_todo)
if (.not. crystallite_todo(g,i,e) .and. .not. crystallite_localPlasticity(g,i,e)) then ! if broken non-local...
!$OMP CRITICAL (checkTodo)
crystallite_todo = crystallite_todo .and. crystallite_localPlasticity ! ...all non-locals skipped
!$OMP END CRITICAL (checkTodo)
endif
endif
enddo; enddo; enddo
!$OMP END PARALLEL DO
nonlocalBroken = .false.
!$OMP PARALLEL DO
do e = FEsolving_execElem(1),FEsolving_execElem(2)
do i = FEsolving_execIP(1),FEsolving_execIP(2)
do g = 1,homogenization_Ngrains(material_homogenizationAt(e))
if(crystallite_todo(g,i,e) .and. .not. crystallite_converged(g,i,e) .and. &
(.not. nonlocalBroken .or. crystallite_localPlasticity(g,i,e)) ) then
crystallite_todo(g,i,e) = integrateStress(g,i,e,timeFraction)
if (.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) &
nonlocalBroken = .true.
endif
enddo; enddo; enddo
!$OMP END PARALLEL DO
if(nonlocalBroken) where(.not. crystallite_localPlasticity) crystallite_todo = .false.
end subroutine update_stress
!--------------------------------------------------------------------------------------------------
!> @brief tbd
!--------------------------------------------------------------------------------------------------
@ -1617,7 +1579,7 @@ subroutine update_dependentState
do i = FEsolving_execIP(1),FEsolving_execIP(2)
do g = 1,homogenization_Ngrains(material_homogenizationAt(e))
if (crystallite_todo(g,i,e) .and. .not. crystallite_converged(g,i,e)) &
call constitutive_dependentState(crystallite_Fe(1:3,1:3,g,i,e), &
call constitutive_dependentState(crystallite_partionedF(1:3,1:3,g,i,e), &
crystallite_Fp(1:3,1:3,g,i,e), &
g, i, e)
enddo; enddo; enddo
@ -1666,56 +1628,56 @@ subroutine update_state(timeFraction)
end subroutine update_state
!--------------------------------------------------------------------------------------------------
!> @brief triggers calculation of all new rates
!---------------------------------------------------------------------------------------------------
!> @brief Trigger calculation of all new rates
!> if NaN occurs, crystallite_todo is set to FALSE. Any NaN in a nonlocal propagates to all others
!--------------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------------
subroutine update_dotState(timeFraction)
real(pReal), intent(in) :: &
timeFraction
integer :: &
e, & !< element index in element loop
i, & !< integration point index in ip loop
g, & !< grain index in grain loop
p, &
c, &
s
logical :: &
NaN, &
nonlocalStop
real(pReal), intent(in) :: &
timeFraction
integer :: &
e, & !< element index in element loop
i, & !< integration point index in ip loop
g, & !< grain index in grain loop
p, &
c, &
s
logical :: &
nonlocalBroken
nonlocalStop = .false.
nonlocalBroken = .false.
!$OMP PARALLEL DO PRIVATE (p,c)
do e = FEsolving_execElem(1),FEsolving_execElem(2)
do i = FEsolving_execIP(1),FEsolving_execIP(2)
do g = 1,homogenization_Ngrains(material_homogenizationAt(e))
if(crystallite_todo(g,i,e) .and. .not. crystallite_converged(g,i,e) .and. &
(.not. nonlocalBroken .or. crystallite_localPlasticity(g,i,e)) ) then
call constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), &
crystallite_partionedF0, &
crystallite_Fi(1:3,1:3,g,i,e), &
crystallite_partionedFp0, &
crystallite_subdt(g,i,e)*timeFraction, g,i,e)
p = material_phaseAt(g,e); c = material_phaseMemberAt(g,i,e)
crystallite_todo(g,i,e) = all(.not. IEEE_is_NaN(plasticState(p)%dotState(:,c)))
do s = 1, phase_Nsources(p)
crystallite_todo(g,i,e) = crystallite_todo(g,i,e) .and. all(.not. IEEE_is_NaN(sourceState(p)%p(s)%dotState(:,c)))
enddo
if (.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) &
nonlocalBroken = .true.
endif
enddo; enddo; enddo
!$OMP END PARALLEL DO
!$OMP PARALLEL DO PRIVATE (p,c,NaN)
do e = FEsolving_execElem(1),FEsolving_execElem(2)
do i = FEsolving_execIP(1),FEsolving_execIP(2)
do g = 1,homogenization_Ngrains(material_homogenizationAt(e))
!$OMP FLUSH(nonlocalStop)
if ((crystallite_todo(g,i,e) .and. .not. crystallite_converged(g,i,e)) .and. .not. nonlocalStop) then
call constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), &
crystallite_partionedF0, &
crystallite_Fi(1:3,1:3,g,i,e), &
crystallite_partionedFp0, &
crystallite_subdt(g,i,e)*timeFraction, g,i,e)
p = material_phaseAt(g,e); c = material_phaseMemberAt(g,i,e)
NaN = any(IEEE_is_NaN(plasticState(p)%dotState(:,c)))
do s = 1, phase_Nsources(p)
NaN = NaN .or. any(IEEE_is_NaN(sourceState(p)%p(s)%dotState(:,c)))
enddo
if (NaN) then
crystallite_todo(g,i,e) = .false. ! this one done (and broken)
if (.not. crystallite_localPlasticity(g,i,e)) nonlocalStop = .True.
endif
endif
enddo; enddo; enddo
!$OMP END PARALLEL DO
if (nonlocalStop) crystallite_todo = crystallite_todo .and. crystallite_localPlasticity
if(nonlocalBroken) where(.not. crystallite_localPlasticity) crystallite_todo = .false.
end subroutine update_DotState
!---------------------------------------------------------------------------------------------------
!> @brief Trigger calculation of all new sudden state change
!> if NaN occurs, crystallite_todo is set to FALSE. Any NaN in a nonlocal propagates to all others
!---------------------------------------------------------------------------------------------------
subroutine update_deltaState
integer :: &
@ -1783,62 +1745,49 @@ end subroutine update_deltaState
!--------------------------------------------------------------------------------------------------
logical function stateJump(ipc,ip,el)
integer, intent(in):: &
el, & ! element index
ip, & ! integration point index
ipc ! grain index
integer, intent(in):: &
el, & ! element index
ip, & ! integration point index
ipc ! grain index
integer :: &
c, &
p, &
mySource, &
myOffset, &
mySize
integer :: &
c, &
p, &
mySource, &
myOffset, &
mySize
c = material_phaseMemberAt(ipc,ip,el)
p = material_phaseAt(ipc,el)
c = material_phaseMemberAt(ipc,ip,el)
p = material_phaseAt(ipc,el)
call constitutive_collectDeltaState(crystallite_S(1:3,1:3,ipc,ip,el), &
crystallite_Fe(1:3,1:3,ipc,ip,el), &
crystallite_Fi(1:3,1:3,ipc,ip,el), &
ipc,ip,el)
call constitutive_collectDeltaState(crystallite_S(1:3,1:3,ipc,ip,el), &
crystallite_Fe(1:3,1:3,ipc,ip,el), &
crystallite_Fi(1:3,1:3,ipc,ip,el), &
ipc,ip,el)
myOffset = plasticState(p)%offsetDeltaState
mySize = plasticState(p)%sizeDeltaState
myOffset = plasticState(p)%offsetDeltaState
mySize = plasticState(p)%sizeDeltaState
if( any(IEEE_is_NaN(plasticState(p)%deltaState(1:mySize,c)))) then ! NaN occured in deltaState
stateJump = .false.
return
endif
if( any(IEEE_is_NaN(plasticState(p)%deltaState(1:mySize,c)))) then
stateJump = .false.
return
endif
plasticState(p)%state(myOffset + 1:myOffset + mySize,c) = &
plasticState(p)%state(myOffset + 1:myOffset + mySize,c) + plasticState(p)%deltaState(1:mySize,c)
plasticState(p)%state(myOffset + 1:myOffset + mySize,c) = &
plasticState(p)%state(myOffset + 1:myOffset + mySize,c) + plasticState(p)%deltaState(1:mySize,c)
do mySource = 1, phase_Nsources(p)
myOffset = sourceState(p)%p(mySource)%offsetDeltaState
mySize = sourceState(p)%p(mySource)%sizeDeltaState
if (any(IEEE_is_NaN(sourceState(p)%p(mySource)%deltaState(1:mySize,c)))) then ! NaN occured in deltaState
stateJump = .false.
return
endif
sourceState(p)%p(mySource)%state(myOffset + 1: myOffset + mySize,c) = &
sourceState(p)%p(mySource)%state(myOffset + 1: myOffset + mySize,c) + sourceState(p)%p(mySource)%deltaState(1:mySize,c)
enddo
do mySource = 1, phase_Nsources(p)
myOffset = sourceState(p)%p(mySource)%offsetDeltaState
mySize = sourceState(p)%p(mySource)%sizeDeltaState
if (any(IEEE_is_NaN(sourceState(p)%p(mySource)%deltaState(1:mySize,c)))) then
stateJump = .false.
return
endif
sourceState(p)%p(mySource)%state(myOffset + 1: myOffset + mySize,c) = &
sourceState(p)%p(mySource)%state(myOffset + 1: myOffset + mySize,c) + sourceState(p)%p(mySource)%deltaState(1:mySize,c)
enddo
#ifdef DEBUG
if (any(dNeq0(plasticState(p)%deltaState(1:mySize,c))) &
.and. iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0 &
.and. ((el == debug_e .and. ip == debug_i .and. ipc == debug_g) &
.or. .not. iand(debug_level(debug_crystallite), debug_levelSelective) /= 0)) then
write(6,'(a,i8,1x,i2,1x,i3, /)') '<< CRYST >> update state at el ip ipc ',el,ip,ipc
write(6,'(a,/,(12x,12(e12.5,1x)),/)') '<< CRYST >> deltaState', plasticState(p)%deltaState(1:mySize,c)
write(6,'(a,/,(12x,12(e12.5,1x)),/)') '<< CRYST >> new state', &
plasticState(p)%state(myOffset + 1 : &
myOffset + mySize,c)
endif
#endif
stateJump = .true.
stateJump = .true.
end function stateJump