Heavy bug fixing: Butcher tableau multiplication in RK45 integration scheme, state updating of adaptive Euler, wrong assignment to evolution rate instead of state in fixed-point iterator.

Exchanged possibly non thread-safe matmul in RKCK45, simplified RK45 integration step 1--3 vs 4 handling.
(Thanks to Pratheek, Luv, and Chen for their help!)
This commit is contained in:
Philip Eisenlohr 2014-09-02 19:46:52 +00:00
parent 5f1e49c053
commit 630d9efffd
1 changed files with 359 additions and 338 deletions

View File

@ -53,9 +53,9 @@ module crystallite
crystallite_Lp, & !< current plastic velocitiy grad (end of converged time step)
crystallite_Lp0, & !< plastic velocitiy grad at start of FE inc
crystallite_partionedLp0,& !< plastic velocity grad at start of homog inc
crystallite_Fe, & !< current "elastic" def grad (end of converged time step)
crystallite_P !< 1st Piola-Kirchhoff stress per grain
real(pReal), dimension(:,:,:,:,:), allocatable, private :: &
crystallite_Fe, & !< current "elastic" def grad (end of converged time step)
crystallite_subFe0,& !< "elastic" def grad at start of crystallite inc
crystallite_invFp, & !< inverse of current plastic def grad (end of converged time step)
crystallite_subFp0,& !< plastic def grad at start of crystallite inc
@ -173,7 +173,7 @@ subroutine crystallite_init(temperature)
use lattice, only: &
lattice_structure
use constitutive, only: &
constitutive_microstructure
constitutive_microstructure ! derived (shortcut) quantities of given state
use constitutive_damage, only: &
constitutive_damage_microstructure
use constitutive_thermal, only: &
@ -424,7 +424,8 @@ subroutine crystallite_init(temperature)
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e)
do g = 1_pInt,myNgrains
call constitutive_microstructure(temperature, &
crystallite_Fe(1:3,1:3,g,i,e), crystallite_Fp(1:3,1:3,g,i,e),g,i,e) ! update dependent state variables to be consistent with basic states
crystallite_Fe(1:3,1:3,g,i,e), &
crystallite_Fp(1:3,1:3,g,i,e),g,i,e) ! update dependent state variables to be consistent with basic states
call constitutive_damage_microstructure(crystallite_Tstar_v(1:6,g,i,e), &
crystallite_Fe(1:3,1:3,g,i,e), &
g,i,e) ! update dependent state variables to be consistent with basic states
@ -907,11 +908,10 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
if (crystallite_subStep(g,i,e) > 0.0_pReal) then
crystallite_subF0(1:3,1:3,g,i,e) = crystallite_subF(1:3,1:3,g,i,e) ! ...def grad
!$OMP FLUSH(crystallite_subF0)
crystallite_subFp0(1:3,1:3,g,i,e) = crystallite_Fp(1:3,1:3,g,i,e) ! ...plastic def grad
crystallite_subFe0(1:3,1:3,g,i,e) = &
math_mul33x33(crystallite_subF(1:3,1:3,g,i,e), crystallite_invFp(1:3,1:3,g,i,e)) ! only needed later on for stiffness calculation
crystallite_subLp0(1:3,1:3,g,i,e) = crystallite_Lp(1:3,1:3,g,i,e) ! ...plastic velocity gradient
crystallite_subFp0(1:3,1:3,g,i,e) = crystallite_Fp(1:3,1:3,g,i,e) ! ...plastic def grad
crystallite_subFe0(1:3,1:3,g,i,e) = math_mul33x33(crystallite_subF(1:3,1:3,g,i,e), &
crystallite_invFp(1:3,1:3,g,i,e)) ! only needed later on for stiffness calculation
!if abbrevation, make c and p private in omp
plasticState(mappingConstitutive(2,g,i,e))%subState0(:,mappingConstitutive(1,g,i,e)) = &
plasticState(mappingConstitutive(2,g,i,e))%state( :,mappingConstitutive(1,g,i,e))
@ -1091,6 +1091,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
enddo
enddo elementLooping5
! --+>> STIFFNESS CALCULATION <<+--
computeJacobian: if(updateJaco) then
@ -1162,16 +1163,18 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
elementLooping7: do e = FEsolving_execElem(1),FEsolving_execElem(2)
myNgrains = homogenization_Ngrains(mesh_element(3,e))
forall (i = FEsolving_execIP(1,e):FEsolving_execIP(2,e), g = 1:myNgrains)
plasticState(mappingConstitutive(2,g,i,e))%state_backup(:,mappingConstitutive(1,g,i,e)) = &
plasticState(mappingConstitutive(2,g,i,e))%state( :,mappingConstitutive(1,g,i,e))
plasticState(mappingConstitutive(2,g,i,e))%dotState_backup(:,mappingConstitutive(1,g,i,e)) = &
plasticState(mappingConstitutive(2,g,i,e))%dotState(:,mappingConstitutive(1,g,i,e))
damageState( mappingConstitutive(2,g,i,e))%state_backup(:,mappingConstitutive(1,g,i,e)) = &
damageState( mappingConstitutive(2,g,i,e))%state( :,mappingConstitutive(1,g,i,e))
damageState(mappingConstitutive(2,g,i,e))%dotState_backup(:,mappingConstitutive(1,g,i,e)) = &
damageState(mappingConstitutive(2,g,i,e))%dotState(:,mappingConstitutive(1,g,i,e))
thermalState(mappingConstitutive(2,g,i,e))%state_backup(:,mappingConstitutive(1,g,i,e)) = &
thermalState(mappingConstitutive(2,g,i,e))%state( :,mappingConstitutive(1,g,i,e))
plasticState(mappingConstitutive(2,g,i,e))%dotState_backup(:,mappingConstitutive(1,g,i,e)) = &
plasticState(mappingConstitutive(2,g,i,e))%dotState( :,mappingConstitutive(1,g,i,e))
damageState( mappingConstitutive(2,g,i,e))%dotState_backup(:,mappingConstitutive(1,g,i,e)) = &
damageState( mappingConstitutive(2,g,i,e))%dotState( :,mappingConstitutive(1,g,i,e))
thermalState(mappingConstitutive(2,g,i,e))%dotState_backup(:,mappingConstitutive(1,g,i,e)) = &
thermalState(mappingConstitutive(2,g,i,e))%dotState( :,mappingConstitutive(1,g,i,e))
@ -1206,18 +1209,21 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
do e = FEsolving_execElem(1),FEsolving_execElem(2)
myNgrains = homogenization_Ngrains(mesh_element(3,e))
forall (i = FEsolving_execIP(1,e):FEsolving_execIP(2,e), g = 1:myNgrains)
plasticState(mappingConstitutive(2,g,i,e))%state( :,mappingConstitutive(1,g,i,e)) = &
plasticState(mappingConstitutive(2,g,i,e))%state_backup(:,mappingConstitutive(1,g,i,e))
plasticState(mappingConstitutive(2,g,i,e))%dotState(:,mappingConstitutive(1,g,i,e)) = &
plasticState(mappingConstitutive(2,g,i,e))%dotState_backup(:,mappingConstitutive(1,g,i,e))
damageState( mappingConstitutive(2,g,i,e))%state( :,mappingConstitutive(1,g,i,e)) = &
damageState( mappingConstitutive(2,g,i,e))%state_backup(:,mappingConstitutive(1,g,i,e))
damageState(mappingConstitutive(2,g,i,e))%dotState(:,mappingConstitutive(1,g,i,e)) = &
damageState(mappingConstitutive(2,g,i,e))%dotState_backup(:,mappingConstitutive(1,g,i,e))
thermalState(mappingConstitutive(2,g,i,e))%state( :,mappingConstitutive(1,g,i,e)) = &
thermalState(mappingConstitutive(2,g,i,e))%state_backup(:,mappingConstitutive(1,g,i,e))
plasticState(mappingConstitutive(2,g,i,e))%dotState( :,mappingConstitutive(1,g,i,e)) = &
plasticState(mappingConstitutive(2,g,i,e))%dotState_backup(:,mappingConstitutive(1,g,i,e))
damageState( mappingConstitutive(2,g,i,e))%dotState( :,mappingConstitutive(1,g,i,e)) = &
damageState( mappingConstitutive(2,g,i,e))%dotState_backup(:,mappingConstitutive(1,g,i,e))
thermalState(mappingConstitutive(2,g,i,e))%dotState( :,mappingConstitutive(1,g,i,e)) = &
thermalState(mappingConstitutive(2,g,i,e))%dotState_backup(:,mappingConstitutive(1,g,i,e))
crystallite_Fp(1:3,1:3,g,i,e) = Fp_backup(1:3,1:3,g,i,e)
crystallite_invFp(1:3,1:3,g,i,e) = InvFp_backup(1:3,1:3,g,i,e)
crystallite_Fe(1:3,1:3,g,i,e) = Fe_backup(1:3,1:3,g,i,e)
@ -1231,18 +1237,21 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
do e = FEsolving_execElem(1),FEsolving_execElem(2)
myNgrains = homogenization_Ngrains(mesh_element(3,e))
forall (i = FEsolving_execIP(1,e):FEsolving_execIP(2,e), g = 1:myNgrains)
plasticState(mappingConstitutive(2,g,i,e))%state( :,mappingConstitutive(1,g,i,e)) = &
plasticState(mappingConstitutive(2,g,i,e))%subState0(:,mappingConstitutive(1,g,i,e))
plasticState(mappingConstitutive(2,g,i,e))%dotState(:,mappingConstitutive(1,g,i,e)) = &
plasticState(mappingConstitutive(2,g,i,e))%dotState_backup(:,mappingConstitutive(1,g,i,e))
damageState( mappingConstitutive(2,g,i,e))%state( :,mappingConstitutive(1,g,i,e)) = &
damageState( mappingConstitutive(2,g,i,e))%subState0(:,mappingConstitutive(1,g,i,e))
damageState(mappingConstitutive(2,g,i,e))%dotState(:,mappingConstitutive(1,g,i,e)) = &
damageState(mappingConstitutive(2,g,i,e))%dotState_backup(:,mappingConstitutive(1,g,i,e))
thermalState(mappingConstitutive(2,g,i,e))%state( :,mappingConstitutive(1,g,i,e)) = &
thermalState(mappingConstitutive(2,g,i,e))%subState0(:,mappingConstitutive(1,g,i,e))
plasticState(mappingConstitutive(2,g,i,e))%dotState( :,mappingConstitutive(1,g,i,e)) = &
plasticState(mappingConstitutive(2,g,i,e))%dotState_backup(:,mappingConstitutive(1,g,i,e))
damageState( mappingConstitutive(2,g,i,e))%dotState( :,mappingConstitutive(1,g,i,e)) = &
damageState( mappingConstitutive(2,g,i,e))%dotState_backup(:,mappingConstitutive(1,g,i,e))
thermalState(mappingConstitutive(2,g,i,e))%dotState( :,mappingConstitutive(1,g,i,e)) = &
thermalState(mappingConstitutive(2,g,i,e))%dotState_backup(:,mappingConstitutive(1,g,i,e))
crystallite_Fp(1:3,1:3,g,i,e) = crystallite_subFp0(1:3,1:3,g,i,e)
crystallite_Fe(1:3,1:3,g,i,e) = crystallite_subFe0(1:3,1:3,g,i,e)
crystallite_Lp(1:3,1:3,g,i,e) = crystallite_subLp0(1:3,1:3,g,i,e)
@ -1328,18 +1337,21 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
elementLooping10: do e = FEsolving_execElem(1),FEsolving_execElem(2)
myNgrains = homogenization_Ngrains(mesh_element(3,e))
forall (i = FEsolving_execIP(1,e):FEsolving_execIP(2,e), g = 1:myNgrains)
plasticState(mappingConstitutive(2,g,i,e))%state( :,mappingConstitutive(1,g,i,e)) = &
plasticState(mappingConstitutive(2,g,i,e))%state_backup(:,mappingConstitutive(1,g,i,e))
plasticState(mappingConstitutive(2,g,i,e))%dotState(:,mappingConstitutive(1,g,i,e)) = &
plasticState(mappingConstitutive(2,g,i,e))%dotState_backup(:,mappingConstitutive(1,g,i,e))
damageState( mappingConstitutive(2,g,i,e))%state( :,mappingConstitutive(1,g,i,e)) = &
damageState( mappingConstitutive(2,g,i,e))%state_backup(:,mappingConstitutive(1,g,i,e))
damageState(mappingConstitutive(2,g,i,e))%dotState(:,mappingConstitutive(1,g,i,e)) = &
damageState(mappingConstitutive(2,g,i,e))%dotState_backup(:,mappingConstitutive(1,g,i,e))
thermalState(mappingConstitutive(2,g,i,e))%state( :,mappingConstitutive(1,g,i,e)) = &
thermalState(mappingConstitutive(2,g,i,e))%state_backup(:,mappingConstitutive(1,g,i,e))
plasticState(mappingConstitutive(2,g,i,e))%dotState( :,mappingConstitutive(1,g,i,e)) = &
plasticState(mappingConstitutive(2,g,i,e))%dotState_backup(:,mappingConstitutive(1,g,i,e))
damageState( mappingConstitutive(2,g,i,e))%dotState( :,mappingConstitutive(1,g,i,e)) = &
damageState( mappingConstitutive(2,g,i,e))%dotState_backup(:,mappingConstitutive(1,g,i,e))
thermalState(mappingConstitutive(2,g,i,e))%dotState( :,mappingConstitutive(1,g,i,e)) = &
thermalState(mappingConstitutive(2,g,i,e))%dotState_backup(:,mappingConstitutive(1,g,i,e))
crystallite_subF(1:3,1:3,g,i,e) = F_backup(1:3,1:3,g,i,e)
crystallite_Fp(1:3,1:3,g,i,e) = Fp_backup(1:3,1:3,g,i,e)
crystallite_invFp(1:3,1:3,g,i,e) = InvFp_backup(1:3,1:3,g,i,e)
@ -1420,7 +1432,7 @@ subroutine crystallite_integrateStateRK4()
real(pReal), dimension(4), parameter :: &
TIMESTEPFRACTION = [0.5_pReal, 0.5_pReal, 1.0_pReal, 1.0_pReal] ! factor giving the fraction of the original timestep used for Runge Kutta Integration
real(pReal), dimension(4), parameter :: &
WEIGHT = [1.0_pReal, 2.0_pReal, 2.0_pReal, 1.0_pReal] ! weight of slope used for Runge Kutta integration
WEIGHT = [1.0_pReal, 2.0_pReal, 2.0_pReal, 1.0_pReal/6.0_pReal] ! weight of slope used for Runge Kutta integration (final weight divided by 6)
integer(pInt) :: e, & ! element index in element loop
i, & ! integration point index in ip loop
@ -1428,7 +1440,6 @@ subroutine crystallite_integrateStateRK4()
p, & ! phase loop
c, &
n, &
mySizeDotState, &
mySizePlasticDotState, &
mySizeDamageDotState, &
mySizeThermalDotState
@ -1514,7 +1525,6 @@ subroutine crystallite_integrateStateRK4()
if (crystallite_todo(g,i,e)) then
p = mappingConstitutive(2,g,i,e)
c = mappingConstitutive(1,g,i,e)
first3steps: if (n < 4) then
plasticState(p)%RK4dotState(:,c) = plasticState(p)%RK4dotState(:,c) &
+ weight(n)*plasticState(p)%dotState(:,c)
@ -1522,21 +1532,11 @@ subroutine crystallite_integrateStateRK4()
+ weight(n)*damageState(p)%dotState(:,c)
thermalState(p)%RK4dotState(:,c) = thermalState(p)%RK4dotState(:,c) &
+ weight(n)*thermalState(p)%dotState(:,c)
else first3steps
plasticState(p)%RK4dotState(:,c) = plasticState(p)%RK4dotState(:,c) &
+ weight(n)*plasticState(p)%dotState(:,c) / 6.0_pReal
damageState(p)%RK4dotState(:,c) = damageState(p)%RK4dotState(:,c) &
+ weight(n)*damageState(p)%dotState(:,c) / 6.0_pReal
thermalState(p)%RK4dotState(:,c) = thermalState(p)%RK4dotState(:,c) &
+ weight(n)*thermalState(p)%dotState(:,c) / 6.0_pReal
endif first3steps
endif
enddo; enddo; enddo
!$OMP ENDDO
!$OMP DO PRIVATE(mySizeDotState,mySizePlasticDotState,mySizeDamageDotState,mySizeThermalDotState,p,c)
!$OMP DO PRIVATE(mySizePlasticDotState,mySizeDamageDotState,mySizeThermalDotState,p,c)
do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
if (crystallite_todo(g,i,e)) then
@ -1554,19 +1554,19 @@ subroutine crystallite_integrateStateRK4()
thermalState(p)%State(1:mySizeThermalDotState,c) = thermalState(p)%subState0(1:mySizeThermalDotState,c) &
+ thermalState(p)%dotState (1:mySizeThermalDotState,c) &
* crystallite_subdt(g,i,e) * timeStepFraction(n)
if (n == 4) then ! final integration step
#ifndef _OPENMP
if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt &
if (n == 4 &
.and. iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt &
.and. ((e == debug_e .and. i == debug_i .and. g == debug_g) &
.or. .not. iand(debug_level(debug_crystallite), debug_levelSelective) /= 0_pInt)) then
.or. .not. iand(debug_level(debug_crystallite), debug_levelSelective) /= 0_pInt)) then ! final integration step
write(6,'(a,i8,1x,i2,1x,i3,/)') '<< CRYST >> updateState at el ip g ',e,i,g
write(6,'(a,/,(12x,12(e12.5,1x)),/)') '<< CRYST >> dotState', plasticState(p)%dotState(:,c)
write(6,'(a,/,(12x,12(e12.5,1x)),/)') '<< CRYST >> new state', plasticState(p)%State(:,c)
write(6,'(a,/,(12x,12(e12.5,1x)),/)') '<< CRYST >> dotState', plasticState(p)%dotState(1:mySizePlasticDotState,c)
write(6,'(a,/,(12x,12(e12.5,1x)),/)') '<< CRYST >> new state', plasticState(p)%state(1:mySizePlasticDotState,c)
endif
#endif
endif
endif
enddo; enddo; enddo
!$OMP ENDDO
@ -1594,7 +1594,8 @@ subroutine crystallite_integrateStateRK4()
!$OMP DO
do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
if (crystallite_todo(g,i,e)) then
call constitutive_microstructure(crystallite_temperature(i,e), crystallite_Fe(1:3,1:3,g,i,e), &
call constitutive_microstructure(crystallite_temperature(i,e), &
crystallite_Fe(1:3,1:3,g,i,e), &
crystallite_Fp(1:3,1:3,g,i,e), g, i, e) ! update dependent state variables to be consistent with basic states
call constitutive_damage_microstructure(crystallite_Tstar_v(1:6,g,i,e), &
crystallite_Fe(1:3,1:3,g,i,e), &
@ -1627,7 +1628,7 @@ subroutine crystallite_integrateStateRK4()
! --- dot state and RK dot state---
first3steps2: if (n < 4) then
first3steps: if (n < 4) then
!$OMP DO
do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
if (crystallite_todo(g,i,e)) then
@ -1667,7 +1668,7 @@ subroutine crystallite_integrateStateRK4()
endif
enddo; enddo; enddo
!$OMP ENDDO
endif first3steps2
endif first3steps
!$OMP END PARALLEL
enddo
@ -1774,8 +1775,7 @@ subroutine crystallite_integrateStateRKCK45()
s, & ! state index
p, &
cc, &
mySizeDotState, & ! size of dot State
mySizePlasticDotState, &
mySizePlasticDotState, & ! size of dot States
mySizeDamageDotState, &
mySizeThermalDotState
integer(pInt), dimension(2) :: &
@ -1785,13 +1785,13 @@ subroutine crystallite_integrateStateRKCK45()
gIter ! bounds for grain iteration
real(pReal), dimension(constitutive_maxSizeDotState,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: &
stateResiduum, & ! residuum from evolution in micrstructure
stateResiduum, & ! residuum from evolution in microstructure
relStateResiduum ! relative residuum from evolution in microstructure
real(pReal), dimension(constitutive_damage_maxSizeDotState,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: &
damageStateResiduum, & ! residuum from evolution in micrstructure
damageStateResiduum, & ! residuum from evolution in microstructure
relDamageStateResiduum ! relative residuum from evolution in microstructure
real(pReal), dimension(constitutive_thermal_maxSizeDotState,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: &
thermalStateResiduum, & ! residuum from evolution in micrstructure
thermalStateResiduum, & ! residuum from evolution in microstructure
relThermalStateResiduum ! relative residuum from evolution in microstructure
logical :: &
singleRun ! flag indicating computation for single (g,i,e) triple
@ -1885,9 +1885,9 @@ subroutine crystallite_integrateStateRKCK45()
plasticState(p)%dotState(:,cc) = A(1,2) * plasticState(p)%RKCK45dotState(1,:,cc) &
+ A(2,2) * plasticState(p)%RKCK45dotState(2,:,cc)
damageState(p)%dotState(:,cc) = A(1,1) * damageState(p)%RKCK45dotState(1,:,cc) &
damageState(p)%dotState(:,cc) = A(1,2) * damageState( p)%RKCK45dotState(1,:,cc) &
+ A(2,2) * damageState( p)%RKCK45dotState(2,:,cc)
thermalState(p)%dotState(:,cc) = A(1,1) * thermalState(p)%RKCK45dotState(1,:,cc) &
thermalState(p)%dotState(:,cc) = A(1,2) * thermalState(p)%RKCK45dotState(1,:,cc) &
+ A(2,2) * thermalState(p)%RKCK45dotState(2,:,cc)
elseif (n == 3) then
@ -1895,10 +1895,10 @@ subroutine crystallite_integrateStateRKCK45()
plasticState(p)%dotState(:,cc) = A(1,3) * plasticState(p)%RKCK45dotState(1,:,cc) &
+ A(2,3) * plasticState(p)%RKCK45dotState(2,:,cc) &
+ A(3,3) * plasticState(p)%RKCK45dotState(3,:,cc)
damageState(p)%dotState(:,cc) = A(1,1) * damageState(p)%RKCK45dotState(1,:,cc) &
damageState(p)%dotState(:,cc) = A(1,3) * damageState( p)%RKCK45dotState(1,:,cc) &
+ A(2,3) * damageState( p)%RKCK45dotState(2,:,cc) &
+ A(3,3) * damageState( p)%RKCK45dotState(3,:,cc)
thermalState(p)%dotState(:,cc) = A(1,1) * thermalState(p)%RKCK45dotState(1,:,cc) &
thermalState(p)%dotState(:,cc) = A(1,3) * thermalState(p)%RKCK45dotState(1,:,cc) &
+ A(2,3) * thermalState(p)%RKCK45dotState(2,:,cc) &
+ A(3,3) * thermalState(p)%RKCK45dotState(3,:,cc)
@ -1908,11 +1908,11 @@ subroutine crystallite_integrateStateRKCK45()
+ A(2,4) * plasticState(p)%RKCK45dotState(2,:,cc) &
+ A(3,4) * plasticState(p)%RKCK45dotState(3,:,cc) &
+ A(4,4) * plasticState(p)%RKCK45dotState(4,:,cc)
damageState(p)%dotState(:,cc) = A(1,1) * damageState(p)%RKCK45dotState(1,:,cc) &
damageState(p)%dotState(:,cc) = A(1,4) * damageState( p)%RKCK45dotState(1,:,cc) &
+ A(2,4) * damageState( p)%RKCK45dotState(2,:,cc) &
+ A(3,4) * damageState( p)%RKCK45dotState(3,:,cc) &
+ A(4,4) * damageState( p)%RKCK45dotState(4,:,cc)
thermalState(p)%dotState(:,cc) = A(1,1) * thermalState(p)%RKCK45dotState(1,:,cc) &
thermalState(p)%dotState(:,cc) = A(1,4) * thermalState(p)%RKCK45dotState(1,:,cc) &
+ A(2,4) * thermalState(p)%RKCK45dotState(2,:,cc) &
+ A(3,4) * thermalState(p)%RKCK45dotState(3,:,cc) &
+ A(4,4) * thermalState(p)%RKCK45dotState(4,:,cc)
@ -1923,12 +1923,12 @@ subroutine crystallite_integrateStateRKCK45()
+ A(3,5) * plasticState(p)%RKCK45dotState(3,:,cc) &
+ A(4,5) * plasticState(p)%RKCK45dotState(4,:,cc) &
+ A(5,5) * plasticState(p)%RKCK45dotState(5,:,cc)
damageState(p)%dotState(:,cc) = A(1,1) * damageState(p)%RKCK45dotState(1,:,cc) &
damageState(p)%dotState(:,cc) = A(1,5) * damageState( p)%RKCK45dotState(1,:,cc) &
+ A(2,5) * damageState( p)%RKCK45dotState(2,:,cc) &
+ A(3,5) * damageState( p)%RKCK45dotState(3,:,cc) &
+ A(4,5) * damageState( p)%RKCK45dotState(4,:,cc) &
+ A(5,5) * damageState( p)%RKCK45dotState(5,:,cc)
thermalState(p)%dotState(:,cc) = A(1,1) * thermalState(p)%RKCK45dotState(1,:,cc) &
thermalState(p)%dotState(:,cc) = A(1,5) * thermalState(p)%RKCK45dotState(1,:,cc) &
+ A(2,5) * thermalState(p)%RKCK45dotState(2,:,cc) &
+ A(3,5) * thermalState(p)%RKCK45dotState(3,:,cc) &
+ A(4,5) * thermalState(p)%RKCK45dotState(4,:,cc) &
@ -1938,7 +1938,7 @@ subroutine crystallite_integrateStateRKCK45()
endif
enddo; enddo; enddo
!$OMP ENDDO
!$OMP DO PRIVATE(mySizeDotState,mySizePlasticDotState,mySizeDamageDotState,mySizeThermalDotState,p,cc)
!$OMP DO PRIVATE(mySizePlasticDotState,mySizeDamageDotState,mySizeThermalDotState,p,cc)
do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
if (crystallite_todo(g,i,e)) then
p = mappingConstitutive(2,g,i,e)
@ -2082,7 +2082,7 @@ subroutine crystallite_integrateStateRKCK45()
enddo; enddo; enddo
!$OMP ENDDO
!$OMP DO PRIVATE(mySizeDotState,mySizePlasticDotState,mySizeDamageDotState,mySizeThermalDotState,p,cc)
!$OMP DO PRIVATE(mySizePlasticDotState,mySizeDamageDotState,mySizeThermalDotState,p,cc)
do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
if (crystallite_todo(g,i,e)) then
p = mappingConstitutive(2,g,i,e)
@ -2095,12 +2095,30 @@ subroutine crystallite_integrateStateRKCK45()
! NEED TO DO THE ADDITION IN THIS LENGTHY WAY BECAUSE OF PARALLELIZATION
! CAN'T USE A REDUCTION CLAUSE ON A POINTER OR USER DEFINED TYPE
stateResiduum(1:mySizePlasticDotState,g,i,e) = matmul(DB, &
plasticState(p)%RKCK45dotState(1:6,1:mySizePlasticDotState,cc))*crystallite_subdt(g,i,e)
damageStateResiduum(1:mySizeDamageDotState,g,i,e) = matmul(DB, &
damageState(p)%RKCK45dotState(1:6,1:mySizeDamageDotState,cc))*crystallite_subdt(g,i,e)
thermalStateResiduum(1:mySizeThermalDotState,g,i,e) = matmul(DB, &
thermalState(p)%RKCK45dotState(1:6,1:mySizeThermalDotState,cc))*crystallite_subdt(g,i,e)
stateResiduum(1:mySizePlasticDotState,g,i,e) = &
( DB(1) * plasticState(p)%RKCK45dotState(1,1:mySizePlasticDotState,cc) &
+ DB(2) * plasticState(p)%RKCK45dotState(2,1:mySizePlasticDotState,cc) &
+ DB(3) * plasticState(p)%RKCK45dotState(3,1:mySizePlasticDotState,cc) &
+ DB(4) * plasticState(p)%RKCK45dotState(4,1:mySizePlasticDotState,cc) &
+ DB(5) * plasticState(p)%RKCK45dotState(5,1:mySizePlasticDotState,cc) &
+ DB(6) * plasticState(p)%RKCK45dotState(6,1:mySizePlasticDotState,cc) &
) * crystallite_subdt(g,i,e)
damageStateResiduum(1:mySizeDamageDotState,g,i,e) = &
( DB(1) * damageState(p)%RKCK45dotState(1,1:mySizeDamageDotState,cc) &
+ DB(2) * damageState(p)%RKCK45dotState(2,1:mySizeDamageDotState,cc) &
+ DB(3) * damageState(p)%RKCK45dotState(3,1:mySizeDamageDotState,cc) &
+ DB(4) * damageState(p)%RKCK45dotState(4,1:mySizeDamageDotState,cc) &
+ DB(5) * damageState(p)%RKCK45dotState(5,1:mySizeDamageDotState,cc) &
+ DB(6) * damageState(p)%RKCK45dotState(6,1:mySizeDamageDotState,cc) &
) * crystallite_subdt(g,i,e)
thermalStateResiduum(1:mySizethermalDotState,g,i,e) = &
( DB(1) * thermalState(p)%RKCK45dotState(1,1:mySizeThermalDotState,cc) &
+ DB(2) * thermalState(p)%RKCK45dotState(2,1:mySizeThermalDotState,cc) &
+ DB(3) * thermalState(p)%RKCK45dotState(3,1:mySizeThermalDotState,cc) &
+ DB(4) * thermalState(p)%RKCK45dotState(4,1:mySizeThermalDotState,cc) &
+ DB(5) * thermalState(p)%RKCK45dotState(5,1:mySizeThermalDotState,cc) &
+ DB(6) * thermalState(p)%RKCK45dotState(6,1:mySizeThermalDotState,cc) &
) * crystallite_subdt(g,i,e)
! --- dot state ---
@ -2128,7 +2146,7 @@ subroutine crystallite_integrateStateRKCK45()
! --- state and update ---
!$OMP DO PRIVATE(mySizeDotState,mySizePlasticDotState,mySizeDamageDotState,mySizeThermalDotState,p,cc)
!$OMP DO PRIVATE(mySizePlasticDotState,mySizeDamageDotState,mySizeThermalDotState,p,cc)
do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
if (crystallite_todo(g,i,e)) then
@ -2152,7 +2170,7 @@ subroutine crystallite_integrateStateRKCK45()
! --- relative residui and state convergence ---
!$OMP DO PRIVATE(mySizeDotState,mySizePlasticDotState,mySizeDamageDotState,mySizeThermalDotState,p,cc)
!$OMP DO PRIVATE(mySizePlasticDotState,mySizeDamageDotState,mySizeThermalDotState,p,cc,s)
do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
if (crystallite_todo(g,i,e)) then
p = mappingConstitutive(2,g,i,e)
@ -2167,6 +2185,9 @@ subroutine crystallite_integrateStateRKCK45()
forall (s = 1_pInt:mySizeThermalDotState, abs(thermalState(p)%state(s,cc)) > 0.0_pReal) &
relThermalStateResiduum(s,g,i,e) = thermalStateResiduum(s,g,i,e) / thermalState(p)%state(s,cc)
!$OMP FLUSH(relStateResiduum)
!$OMP FLUSH(relDamageStateResiduum)
!$OMP FLUSH(relThermalStateResiduum)
! @Martin: do we need flushing? why..?
crystallite_todo(g,i,e) = &
( all(abs(relStateResiduum(1:mySizePlasticDotState,g,i,e)) < &
rTol_crystalliteState .or. &
@ -2187,9 +2208,9 @@ subroutine crystallite_integrateStateRKCK45()
.or. .not. iand(debug_level(debug_crystallite), debug_levelSelective) /= 0_pInt)) then
write(6,'(a,i8,1x,i3,1x,i3,/)') '<< CRYST >> updateState at el ip g ',e,i,g
write(6,'(a,/,(12x,12(f12.1,1x)),/)') '<< CRYST >> absolute residuum tolerance', &
stateResiduum(1:mySizeDotState,g,i,e) / plasticState(p)%aTolState(1:mySizePlasticDotState)
stateResiduum(1:mySizePlasticDotState,g,i,e) / plasticState(p)%aTolState(1:mySizePlasticDotState)
write(6,'(a,/,(12x,12(f12.1,1x)),/)') '<< CRYST >> relative residuum tolerance', &
relStateResiduum(1:mySizeDotState,g,i,e) / rTol_crystalliteState
relStateResiduum(1:mySizePlasticDotState,g,i,e) / rTol_crystalliteState
write(6,'(a,/,(12x,12(e12.5,1x)),/)') '<< CRYST >> dotState', &
plasticState(p)%dotState(1:mySizePlasticDotState,cc)
write(6,'(a,/,(12x,12(e12.5,1x)),/)') '<< CRYST >> new state', &
@ -2337,8 +2358,7 @@ subroutine crystallite_integrateStateAdaptiveEuler()
s, & ! state index
p, &
c, &
mySizeDotState, & ! size of dot State
mySizePlasticDotState, &
mySizePlasticDotState, & ! size of dot states
mySizeDamageDotState, &
mySizeThermalDotState
integer(pInt), dimension(2) :: &
@ -2423,7 +2443,7 @@ subroutine crystallite_integrateStateAdaptiveEuler()
! --- STATE UPDATE (EULER INTEGRATION) ---
!$OMP DO PRIVATE(mySizeDotState,mySizePlasticDotState,mySizeDamageDotState,mySizeThermalDotState,p,c)
!$OMP DO PRIVATE(mySizePlasticDotState,mySizeDamageDotState,mySizeThermalDotState,p,c)
do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
if (crystallite_todo(g,i,e)) then
p = mappingConstitutive(2,g,i,e)
@ -2559,7 +2579,7 @@ subroutine crystallite_integrateStateAdaptiveEuler()
relThermalStateResiduum = 0.0_pReal
!$OMP END SINGLE
!$OMP DO PRIVATE(mySizeDotState,mySizePlasticDotState,mySizeDamageDotState,mySizeThermalDotState,p,c)
!$OMP DO PRIVATE(mySizePlasticDotState,mySizeDamageDotState,mySizeThermalDotState,p,c,s)
do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
if (crystallite_todo(g,i,e)) then
p = mappingConstitutive(2,g,i,e)
@ -2580,6 +2600,8 @@ subroutine crystallite_integrateStateAdaptiveEuler()
* crystallite_subdt(g,i,e) ! contribution to absolute residuum in state
!$OMP FLUSH(stateResiduum)
!$OMP FLUSH(damageStateResiduum)
!$OMP FLUSH(thermalStateResiduum)
! --- relative residui ---
forall (s = 1_pInt:mySizePlasticDotState, abs(plasticState(p)%dotState(s,c)) > 0.0_pReal) &
@ -2589,6 +2611,8 @@ subroutine crystallite_integrateStateAdaptiveEuler()
forall (s = 1_pInt:mySizeThermalDotState, abs(thermalState(p)%dotState(s,c)) > 0.0_pReal) &
relThermalStateResiduum(s,g,i,e) = thermalStateResiduum(s,g,i,e) / thermalState(p)%dotState(s,c)
!$OMP FLUSH(relStateResiduum)
!$OMP FLUSH(relDamageStateResiduum)
!$OMP FLUSH(relthermalStateResiduum)
#ifndef _OPENMP
@ -2710,7 +2734,6 @@ subroutine crystallite_integrateStateEuler()
g, & ! grain index in grain loop
p, & ! phase loop
c, &
mySizeDotState, &
mySizePlasticDotState, &
mySizeDamageDotState, &
mySizeThermalDotState
@ -2775,7 +2798,7 @@ eIter = FEsolving_execElem(1:2)
! --- UPDATE STATE ---
!$OMP DO PRIVATE(mySizeDotState,mySizePlasticDotState,mySizeDamageDotState,mySizeThermalDotState,p,c)
!$OMP DO PRIVATE(mySizePlasticDotState,mySizeDamageDotState,mySizeThermalDotState,p,c)
do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
if (crystallite_todo(g,i,e) .and. .not. crystallite_converged(g,i,e)) then
p = mappingConstitutive(2,g,i,e)
@ -2783,13 +2806,13 @@ eIter = FEsolving_execElem(1:2)
mySizePlasticDotState = plasticState(p)%sizeDotState
mySizeDamageDotState = damageState( p)%sizeDotState
mySizeThermalDotState = thermalState(p)%sizeDotState
plasticState(p)%State(1:mySizePlasticDotState,c) = plasticState(p)%subState0(1:mySizePlasticDotState,c) &
plasticState(p)%state(1:mySizePlasticDotState,c) = plasticState(p)%state( 1:mySizePlasticDotState,c) &
+ plasticState(p)%dotState(1:mySizePlasticDotState,c) &
* crystallite_subdt(g,i,e)
damageState(p)%State(1:mySizeDamageDotState,c) = damageState(p)%subState0(1:mySizeDamageDotState,c) &
damageState( p)%state(1:mySizeDamageDotState,c) = damageState( p)%state( 1:mySizeDamageDotState,c) &
+ damageState( p)%dotState(1:mySizeDamageDotState,c) &
* crystallite_subdt(g,i,e)
thermalState(p)%State(1:mySizeThermalDotState,c) = thermalState(p)%subState0(1:mySizeThermalDotState,c) &
thermalState(p)%state(1:mySizeThermalDotState,c) = thermalState(p)%state( 1:mySizeThermalDotState,c) &
+ thermalState(p)%dotState(1:mySizeThermalDotState,c) &
* crystallite_subdt(g,i,e)
@ -2950,8 +2973,7 @@ subroutine crystallite_integrateStateFPI()
g, & !< grain index in grain loop
p, &
c, &
mySizeDotState, & ! size of dot State
mySizePlasticDotState, &
mySizePlasticDotState, & ! size of dot states
mySizeDamageDotState, &
mySizeThermalDotState, &
ss
@ -3008,7 +3030,6 @@ subroutine crystallite_integrateStateFPI()
do g = gIter(1,e), gIter(2,e)
p = mappingConstitutive(2,g,i,e)
c = mappingConstitutive(1,g,i,e)
plasticState(p)%previousDotState (:,c) = 0.0_pReal
plasticState(p)%previousDotState2(:,c) = 0.0_pReal
damageState( p)%previousDotState (:,c) = 0.0_pReal
@ -3064,7 +3085,7 @@ subroutine crystallite_integrateStateFPI()
! --- UPDATE STATE ---
!$OMP DO PRIVATE(mySizeDotState,mySizePlasticDotState,mySizeDamageDotState,mySizeThermalDotState,p,c)
!$OMP DO PRIVATE(mySizePlasticDotState,mySizeDamageDotState,mySizeThermalDotState,p,c)
do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
if (crystallite_todo(g,i,e)) then
p = mappingConstitutive(2,g,i,e)
@ -3072,10 +3093,10 @@ subroutine crystallite_integrateStateFPI()
mySizePlasticDotState = plasticState(p)%sizeDotState
mySizeDamageDotState = damageState( p)%sizeDotState
mySizeThermalDotState = thermalState(p)%sizeDotState
plasticState(p)%State(1:mySizePlasticDotState,c) = plasticState(p)%subState0(1:mySizePlasticDotState,c) &
plasticState(p)%state(1:mySizePlasticDotState,c) = plasticState(p)%subState0(1:mySizePlasticDotState,c) &
+ plasticState(p)%dotState (1:mySizePlasticDotState,c) &
* crystallite_subdt(g,i,e)
damageState(p)%State(1:mySizeDamageDotState,c) = damageState(p)%subState0(1:mySizeDamageDotState,c) &
damageState( p)%state(1:mySizeDamageDotState,c) = damageState( p)%subState0(1:mySizeDamageDotState,c) &
+ damageState( p)%dotState (1:mySizeDamageDotState,c) &
* crystallite_subdt(g,i,e)
thermalState(p)%State(1:mySizeThermalDotState,c) = thermalState(p)%subState0(1:mySizeThermalDotState,c) &
@ -3100,7 +3121,7 @@ subroutine crystallite_integrateStateFPI()
!$OMP DO PRIVATE(p,c)
do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
if (crystallite_todo(g,i,e) .and. .not. crystallite_converged(g,i,e)) then !ToDo: should this be independent of todo/converged?
if (crystallite_todo(g,i,e) .and. .not. crystallite_converged(g,i,e)) then
call constitutive_microstructure(crystallite_temperature(i,e), crystallite_Fe(1:3,1:3,g,i,e), &
crystallite_Fp(1:3,1:3,g,i,e), g, i, e) ! update dependent state variables to be consistent with basic states
call constitutive_damage_microstructure(crystallite_Tstar_v(1:6,g,i,e), &
@ -3113,10 +3134,10 @@ endif
p = mappingConstitutive(2,g,i,e)
c = mappingConstitutive(1,g,i,e)
plasticState(p)%previousDotState2(:,c) = plasticState(p)%previousDotState(:,c)
plasticState(p)%previousDotState(:,c) = plasticState(p)%dotState(:,c)
damageState( p)%previousDotState2(:,c) = damageState( p)%previousDotState(:,c)
damageState(p)%previousDotState(:,c) = damageState(p)%dotState(:,c)
thermalState(p)%previousDotState2(:,c) = thermalState(p)%previousDotState(:,c)
plasticState(p)%previousDotState (:,c) = plasticState(p)%dotState(:,c)
damageState( p)%previousDotState (:,c) = damageState( p)%dotState(:,c)
thermalState(p)%previousDotState (:,c) = thermalState(p)%dotState(:,c)
enddo; enddo; enddo
!$OMP ENDDO
@ -3150,7 +3171,7 @@ endif
!$OMP DO
do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
if (crystallite_todo(g,i,e) .and. .not. crystallite_converged(g,i,e)) then !ToDo: should this be independent of todo/converged?
if (crystallite_todo(g,i,e) .and. .not. crystallite_converged(g,i,e)) then
call constitutive_collectDotState(crystallite_Tstar_v(1:6,g,i,e), crystallite_Fe, &
crystallite_Fp, crystallite_temperature(i,e), &
crystallite_subdt(g,i,e), crystallite_subFrac, g,i,e)
@ -3193,7 +3214,7 @@ endif
!$OMP& mySizePlasticDotState,mySizeDamageDotState,mySizeThermalDotState, &
!$OMP& damageStateResiduum,thermalStateResiduum,damageStateDamper,thermalStateDamper, &
!$OMP& tempDamageState,tempThermalState,p,c, &
!$OMP& statedamper,mySizeDotState,stateResiduum,tempState)
!$OMP& statedamper,stateResiduum,tempState)
do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
if (crystallite_todo(g,i,e) .and. .not. crystallite_converged(g,i,e)) then
@ -3202,6 +3223,7 @@ endif
mySizePlasticDotState = plasticState(p)%sizeDotState
mySizeDamageDotState = damageState( p)%sizeDotState
mySizeThermalDotState = thermalState(p)%sizeDotState
dot_prod12 = dot_product(plasticState(p)%dotState(:,c) - plasticState(p)%previousDotState(:,c), &
plasticState(p)%previousDotState(:,c) - plasticState(p)%previousDotState2(:,c))
dot_prod22 = dot_product(plasticState(p)%previousDotState(:,c) - plasticState(p)%previousDotState2(:,c), &
@ -3214,6 +3236,7 @@ endif
else
stateDamper = 1.0_pReal
endif
dot_prod12 = dot_product(damageState(p)%dotState(:,c) - damageState(p)%previousDotState(:,c), &
damageState(p)%previousDotState(:,c) - damageState(p)%previousDotState2(:,c))
dot_prod22 = dot_product(damageState(p)%previousDotState(:,c) - damageState(p)%previousDotState2(:,c), &
@ -3226,6 +3249,7 @@ endif
else
damageStateDamper = 1.0_pReal
endif
dot_prod12 = dot_product(thermalState(p)%dotState(:,c) - thermalState(p)%previousDotState(:,c), &
thermalState(p)%previousDotState(:,c) - thermalState(p)%previousDotState2(:,c))
dot_prod22 = dot_product(thermalState(p)%previousDotState(:,c) - thermalState(p)%previousDotState2(:,c), &
@ -3304,9 +3328,9 @@ endif
!$OMP END CRITICAL (distributionState)
endif
endif
plasticState(p)%dotState(1:mySizePlasticDotState,c) = tempState(1:mySizePlasticDotState) ! copy local backup to global pointer
damageState(p)%dotState (1:mySizeDamageDotState, c) = tempDamageState(1:mySizeDamageDotState) ! copy local backup to global pointer
thermalState(p)%dotState(1:mySizeThermalDotState,c) = tempThermalState(1:mySizeThermalDotState) ! copy local backup to global pointer
plasticState(p)%state(1:mySizePlasticDotState,c) = tempState(1:mySizePlasticDotState)
damageState( p)%state(1:mySizeDamageDotState, c) = tempDamageState(1:mySizeDamageDotState)
thermalState(p)%state(1:mySizeThermalDotState,c) = tempThermalState(1:mySizeThermalDotState)
endif
enddo; enddo; enddo
!$OMP ENDDO
@ -3402,10 +3426,7 @@ logical function crystallite_stateJump(g,i,e)
integer(pInt) :: &
c, &
p, &
mySizeDotState, &
mySizePlasticDotState, &
mySizeDamageDotState, &
mySizeThermalDotState
mySizePlasticDotState
c = mappingConstitutive(1,g,i,e)
p = mappingConstitutive(2,g,i,e)
@ -3802,7 +3823,7 @@ subroutine crystallite_orientations
use material, only: &
material_phase, &
homogenization_Ngrains, &
phase_localPlasticity
plasticState
use mesh, only: &
mesh_element, &
mesh_ipNeighborhood, &
@ -3865,8 +3886,8 @@ subroutine crystallite_orientations
!$OMP PARALLEL DO PRIVATE(myPhase,neighboring_e,neighboring_i,neighboringPhase)
do e = FEsolving_execElem(1),FEsolving_execElem(2)
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e)
myPhase = material_phase(1,i,e) ! get my phase
if (.not. phase_localPlasticity(myPhase)) then ! if nonlocal model
myPhase = material_phase(1,i,e) ! get my phase (non-local models make no sense with more than one grain per material point)
if (plasticState(myPhase)%nonLocal) then ! if nonlocal model
! --- calculate disorientation between me and my neighbor ---
do n = 1_pInt,FE_NipNeighbors(FE_celltype(FE_geomtype(mesh_element(2,e)))) ! loop through my neighbors
@ -3874,7 +3895,7 @@ subroutine crystallite_orientations
neighboring_i = mesh_ipNeighborhood(2,n,i,e)
if (neighboring_e > 0 .and. neighboring_i > 0) then ! if neighbor exists
neighboringPhase = material_phase(1,neighboring_i,neighboring_e) ! get my neighbor's phase
if (.not. phase_localPlasticity(neighboringPhase)) then ! neighbor got also nonlocal plasticity
if (plasticState(neighboringPhase)%nonLocal) then ! neighbor got also nonlocal plasticity
if (lattice_structure(myPhase) == lattice_structure(neighboringPhase)) then ! if my neighbor has same crystal structure like me
crystallite_disorientation(:,n,1,i,e) = &
lattice_qDisorientation( crystallite_orientation(1:4,1,i,e), &
@ -3937,8 +3958,8 @@ function crystallite_postResults(ipc, ip, el)
material_texture, &
homogenization_Ngrains
use constitutive, only: &
constitutive_postResults, &
constitutive_homogenizedC
constitutive_homogenizedC, &
constitutive_postResults
use constitutive_damage, only: &
constitutive_damage_postResults
use constitutive_thermal, only: &