From 369ea31a4b9bbf38db32f997eed0526f7aa087aa Mon Sep 17 00:00:00 2001
From: Martin Diehl <mail@martin-diehl.net>
Date: Tue, 24 Mar 2020 15:02:55 +0100
Subject: [PATCH 01/21] name unification for simple copy and paste

---
 src/crystallite.f90 | 56 ++++++++++++++++++++++-----------------------
 1 file changed, 28 insertions(+), 28 deletions(-)

diff --git a/src/crystallite.f90 b/src/crystallite.f90
index a359f09b2..79b67c878 100644
--- a/src/crystallite.f90
+++ b/src/crystallite.f90
@@ -1423,7 +1423,7 @@ subroutine integrateStateRKCK45
      13525.0_pReal/55296.0_pReal, 277.0_pReal/14336.0_pReal,      0.25_pReal]                       !< coefficients in Butcher tableau (used for final integration and error estimate)
 
  real(pReal), dimension(5), parameter :: &
-   C = [0.2_pReal, 0.3_pReal, 0.6_pReal, 1.0_pReal, 0.875_pReal]                                    !< coefficients in Butcher tableau (fractions of original time step in stages 2 to 6)
+   CC = [0.2_pReal, 0.3_pReal, 0.6_pReal, 1.0_pReal, 0.875_pReal]                                    !< coefficients in Butcher tableau (fractions of original time step in stages 2 to 6)
 
  integer :: &
    e, &                                                                                             ! element index in element loop
@@ -1432,7 +1432,7 @@ subroutine integrateStateRKCK45
    stage, &                                                                                         ! stage index in integration stage loop
    n, &
    p, &
-   cc, &
+   c, &
    s, &
    sizeDotState
 
@@ -1455,27 +1455,27 @@ subroutine integrateStateRKCK45
 
    ! --- state update ---
 
-   !$OMP PARALLEL DO PRIVATE(p,cc)
+   !$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)) then
-         p = material_phaseAt(g,e); cc = material_phaseMemberAt(g,i,e)
+         p = material_phaseAt(g,e); c = material_phaseMemberAt(g,i,e)
 
-         plasticState(p)%RKCK45dotState(stage,:,cc) = plasticState(p)%dotState(:,cc)
-         plasticState(p)%dotState(:,cc) = A(1,stage) * plasticState(p)%RKCK45dotState(1,:,cc)
+         plasticState(p)%RKCK45dotState(stage,:,c) = plasticState(p)%dotState(:,c)
+         plasticState(p)%dotState(:,c) = A(1,stage) * plasticState(p)%RKCK45dotState(1,:,c)
 
          do s = 1, phase_Nsources(p)
-           sourceState(p)%p(s)%RKCK45dotState(stage,:,cc) = sourceState(p)%p(s)%dotState(:,cc)
-           sourceState(p)%p(s)%dotState(:,cc) = A(1,stage) * sourceState(p)%p(s)%RKCK45dotState(1,:,cc)
+           sourceState(p)%p(s)%RKCK45dotState(stage,:,c) = sourceState(p)%p(s)%dotState(:,c)
+           sourceState(p)%p(s)%dotState(:,c) = A(1,stage) * sourceState(p)%p(s)%RKCK45dotState(1,:,c)
          enddo
 
          do n = 2, stage
-           plasticState(p)%dotState(:,cc) = plasticState(p)%dotState(:,cc) &
-                                          + A(n,stage) * plasticState(p)%RKCK45dotState(n,:,cc)
+           plasticState(p)%dotState(:,c) = plasticState(p)%dotState(:,c) &
+                                          + A(n,stage) * plasticState(p)%RKCK45dotState(n,:,c)
            do s = 1, phase_Nsources(p)
-             sourceState(p)%p(s)%dotState(:,cc) = sourceState(p)%p(s)%dotState(:,cc) &
-                                                + A(n,stage) * sourceState(p)%p(s)%RKCK45dotState(n,:,cc)
+             sourceState(p)%p(s)%dotState(:,c) = sourceState(p)%p(s)%dotState(:,c) &
+                                                + A(n,stage) * sourceState(p)%p(s)%RKCK45dotState(n,:,c)
            enddo
          enddo
 
@@ -1486,8 +1486,8 @@ subroutine integrateStateRKCK45
     call update_state(1.0_pReal) !MD: 1.0 correct?
     call update_deltaState
     call update_dependentState
-    call update_stress(C(stage))
-    call update_dotState(C(stage))
+    call update_stress(CC(stage))
+    call update_dotState(CC(stage))
 
  enddo
 
@@ -1495,35 +1495,35 @@ subroutine integrateStateRKCK45
 !--------------------------------------------------------------------------------------------------
 ! --- STATE UPDATE WITH ERROR ESTIMATE FOR STATE ---
 
- !$OMP PARALLEL DO PRIVATE(sizeDotState,p,cc)
+ !$OMP PARALLEL DO PRIVATE(sizeDotState,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)) then
-       p = material_phaseAt(g,e); cc = material_phaseMemberAt(g,i,e)
+       p = material_phaseAt(g,e); c = material_phaseMemberAt(g,i,e)
 
        sizeDotState = plasticState(p)%sizeDotState
 
-       plasticState(p)%RKCK45dotState(6,:,cc) = plasticState (p)%dotState(:,cc)
+       plasticState(p)%RKCK45dotState(6,:,c) = plasticState (p)%dotState(:,c)
 
        residuum_plastic(1:sizeDotState,g,i,e) = &
-         matmul(transpose(plasticState(p)%RKCK45dotState(1:6,1:sizeDotState,cc)),DB) &              ! why transpose? Better to transpose constant DB
+         matmul(transpose(plasticState(p)%RKCK45dotState(1:6,1:sizeDotState,c)),DB) &              ! why transpose? Better to transpose constant DB
        * crystallite_subdt(g,i,e)
 
-        plasticState(p)%dotState(:,cc) =  &
-         matmul(transpose(plasticState(p)%RKCK45dotState(1:6,1:sizeDotState,cc)), B)                ! why transpose? Better to transpose constant B
+        plasticState(p)%dotState(:,c) =  &
+         matmul(transpose(plasticState(p)%RKCK45dotState(1:6,1:sizeDotState,c)), B)                ! why transpose? Better to transpose constant B
 
        do s = 1, phase_Nsources(p)
          sizeDotState = sourceState(p)%p(s)%sizeDotState
 
-         sourceState(p)%p(s)%RKCK45dotState(6,:,cc) = sourceState(p)%p(s)%dotState(:,cc)
+         sourceState(p)%p(s)%RKCK45dotState(6,:,c) = sourceState(p)%p(s)%dotState(:,c)
 
          residuum_source(1:sizeDotState,s,g,i,e) = &
-           matmul(transpose(sourceState(p)%p(s)%RKCK45dotState(1:6,1:sizeDotState,cc)),DB) &
+           matmul(transpose(sourceState(p)%p(s)%RKCK45dotState(1:6,1:sizeDotState,c)),DB) &
          * crystallite_subdt(g,i,e)
 
-         sourceState(p)%p(s)%dotState(:,cc)  = &
-           matmul(transpose(sourceState(p)%p(s)%RKCK45dotState(1:6,1:sizeDotState,cc)),B)
+         sourceState(p)%p(s)%dotState(:,c)  = &
+           matmul(transpose(sourceState(p)%p(s)%RKCK45dotState(1:6,1:sizeDotState,c)),B)
        enddo
 
      endif
@@ -1534,17 +1534,17 @@ subroutine integrateStateRKCK45
 
  ! --- relative residui and state convergence ---
 
- !$OMP PARALLEL DO PRIVATE(sizeDotState,p,cc)
+ !$OMP PARALLEL DO PRIVATE(sizeDotState,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)) then
-         p  = material_phaseAt(g,e); cc = material_phaseMemberAt(g,i,e)
+         p  = material_phaseAt(g,e); c = material_phaseMemberAt(g,i,e)
 
          sizeDotState = plasticState(p)%sizeDotState
 
          crystallite_todo(g,i,e) = converged(residuum_plastic(1:sizeDotState,g,i,e), &
-                                             plasticState(p)%state(1:sizeDotState,cc), &
+                                             plasticState(p)%state(1:sizeDotState,c), &
                                              plasticState(p)%atol(1:sizeDotState))
 
          do s = 1, phase_Nsources(p)
@@ -1552,7 +1552,7 @@ subroutine integrateStateRKCK45
 
                   crystallite_todo(g,i,e) = &
                   crystallite_todo(g,i,e) .and. converged(residuum_source(1:sizeDotState,s,g,i,e), &
-                                                          sourceState(p)%p(s)%state(1:sizeDotState,cc), &
+                                                          sourceState(p)%p(s)%state(1:sizeDotState,c), &
                                                           sourceState(p)%p(s)%atol(1:sizeDotState))
          enddo
      endif

From eb6fe8a3a237dd47eb8d40617213a7a50a3c8ae0 Mon Sep 17 00:00:00 2001
From: Martin Diehl <mail@martin-diehl.net>
Date: Tue, 24 Mar 2020 16:14:14 +0100
Subject: [PATCH 02/21] merge into one loop

---
 src/crystallite.f90 | 34 +++++++++++++++++++++++++---------
 1 file changed, 25 insertions(+), 9 deletions(-)

diff --git a/src/crystallite.f90 b/src/crystallite.f90
index 79b67c878..12994589b 100644
--- a/src/crystallite.f90
+++ b/src/crystallite.f90
@@ -1435,31 +1435,32 @@ subroutine integrateStateRKCK45
    c, &
    s, &
    sizeDotState
-
-   ! ToDo: MD: once all constitutives use allocate state, attach residuum arrays to the state in case of RKCK45
+ logical :: &
+   nonlocalBroken
 
  real(pReal), dimension(constitutive_plasticity_maxSizeDotState,            &
                         homogenization_maxNgrains,discretization_nIP,discretization_nElem) :: &
-   residuum_plastic                                                                         ! relative residuum from evolution in microstructure
+   residuum_plastic
  real(pReal), dimension(constitutive_source_maxSizeDotState, &
                         maxval(phase_Nsources), &
                         homogenization_maxNgrains,discretization_nIP,discretization_nElem) :: &
-   residuum_source                                                                   ! relative residuum from evolution in microstructure
+   residuum_source
 
 
  call update_dotState(1.0_pReal)
 
  ! --- SECOND TO SIXTH RUNGE KUTTA STEP ---
-
+ nonlocalBroken = .false.
  do stage = 1,5
 
    ! --- state update ---
 
-   !$OMP PARALLEL DO PRIVATE(p,c)
+   !$OMP PARALLEL DO PRIVATE(sizeDotState,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)) then
+
          p = material_phaseAt(g,e); c = material_phaseMemberAt(g,i,e)
 
          plasticState(p)%RKCK45dotState(stage,:,c) = plasticState(p)%dotState(:,c)
@@ -1479,12 +1480,27 @@ subroutine integrateStateRKCK45
            enddo
          enddo
 
+         sizeDotState = plasticState(p)%sizeDotState
+         plasticState(p)%state(1:sizeDotState,c) = plasticState(p)%subState0(1:sizeDotState,c) &
+                                                 + plasticState(p)%dotState (1:sizeDotState,c) &
+                                                 * crystallite_subdt(g,i,e)
+         do s = 1, phase_Nsources(p)
+           sizeDotState = sourceState(p)%p(s)%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)
+         enddo
+
+         crystallite_todo(g,i,e) = stateJump(g,i,e)
+         if(.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) &
+           nonlocalBroken = .true.
+         if(.not. crystallite_todo(g,i,e)) cycle
+
        endif
-     enddo; enddo; enddo
+   enddo; enddo; enddo
    !$OMP END PARALLEL DO
 
-    call update_state(1.0_pReal) !MD: 1.0 correct?
-    call update_deltaState
+  if(nonlocalBroken) where(.not. crystallite_localPlasticity) crystallite_todo = .false.
     call update_dependentState
     call update_stress(CC(stage))
     call update_dotState(CC(stage))

From 44e24a9c4f3eec300b1cb6724e54874b8cfd54de Mon Sep 17 00:00:00 2001
From: Martin Diehl <mail@martin-diehl.net>
Date: Tue, 24 Mar 2020 20:25:29 +0100
Subject: [PATCH 03/21] merging into one loop

---
 src/crystallite.f90 | 11 +++++++++--
 1 file changed, 9 insertions(+), 2 deletions(-)

diff --git a/src/crystallite.f90 b/src/crystallite.f90
index 12994589b..09c782deb 100644
--- a/src/crystallite.f90
+++ b/src/crystallite.f90
@@ -1496,13 +1496,20 @@ subroutine integrateStateRKCK45
            nonlocalBroken = .true.
          if(.not. crystallite_todo(g,i,e)) cycle
 
+         call constitutive_dependentState(crystallite_Fe(1:3,1:3,g,i,e), &
+                                          crystallite_Fp(1:3,1:3,g,i,e), &
+                                          g, i, e)
+
+         crystallite_todo(g,i,e) = integrateStress(g,i,e,CC(stage))
+         if(.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) &
+           nonlocalBroken = .true.
+         if(.not. crystallite_todo(g,i,e)) cycle
+
        endif
    enddo; enddo; enddo
    !$OMP END PARALLEL DO
 
   if(nonlocalBroken) where(.not. crystallite_localPlasticity) crystallite_todo = .false.
-    call update_dependentState
-    call update_stress(CC(stage))
     call update_dotState(CC(stage))
 
  enddo

From 0740c9f339cc445b2c08f33cd8a8bc9017078df0 Mon Sep 17 00:00:00 2001
From: Martin Diehl <mail@martin-diehl.net>
Date: Tue, 24 Mar 2020 20:33:26 +0100
Subject: [PATCH 04/21] dot state in loop

---
 src/crystallite.f90 | 17 ++++++++++++++---
 1 file changed, 14 insertions(+), 3 deletions(-)

diff --git a/src/crystallite.f90 b/src/crystallite.f90
index 09c782deb..2b60224b8 100644
--- a/src/crystallite.f90
+++ b/src/crystallite.f90
@@ -1505,15 +1505,26 @@ subroutine integrateStateRKCK45
            nonlocalBroken = .true.
          if(.not. crystallite_todo(g,i,e)) cycle
 
+         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)*CC(stage), 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)) cycle
+
        endif
    enddo; enddo; enddo
    !$OMP END PARALLEL DO
 
-  if(nonlocalBroken) where(.not. crystallite_localPlasticity) crystallite_todo = .false.
-    call update_dotState(CC(stage))
-
  enddo
 
+  if(nonlocalBroken) where(.not. crystallite_localPlasticity) crystallite_todo = .false.
 
 !--------------------------------------------------------------------------------------------------
 ! --- STATE UPDATE WITH ERROR ESTIMATE FOR STATE ---

From 66aa20ad3974045dc3e74804db79ccc493e3253d Mon Sep 17 00:00:00 2001
From: Martin Diehl <mail@martin-diehl.net>
Date: Tue, 24 Mar 2020 20:35:27 +0100
Subject: [PATCH 05/21] extra check for nonlocal needed

---
 src/crystallite.f90 | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/crystallite.f90 b/src/crystallite.f90
index 2b60224b8..882eb3d01 100644
--- a/src/crystallite.f90
+++ b/src/crystallite.f90
@@ -1459,7 +1459,7 @@ subroutine integrateStateRKCK45
    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)) then
+        if(crystallite_todo(g,i,e) .and. (.not. nonlocalBroken .or. crystallite_localPlasticity(g,i,e)) ) then
 
          p = material_phaseAt(g,e); c = material_phaseMemberAt(g,i,e)
 

From e7c585a02eefd1f2fb9f6ea6af8d963ef5515765 Mon Sep 17 00:00:00 2001
From: Martin Diehl <mail@martin-diehl.net>
Date: Tue, 24 Mar 2020 20:47:41 +0100
Subject: [PATCH 06/21] loop order that allows more memory efficient code

---
 src/crystallite.f90 | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/crystallite.f90 b/src/crystallite.f90
index 882eb3d01..77fd00fe7 100644
--- a/src/crystallite.f90
+++ b/src/crystallite.f90
@@ -1451,7 +1451,6 @@ subroutine integrateStateRKCK45
 
  ! --- SECOND TO SIXTH RUNGE KUTTA STEP ---
  nonlocalBroken = .false.
- do stage = 1,5
 
    ! --- state update ---
 
@@ -1461,6 +1460,7 @@ subroutine integrateStateRKCK45
        do g = 1,homogenization_Ngrains(material_homogenizationAt(e))
         if(crystallite_todo(g,i,e) .and. (.not. nonlocalBroken .or. crystallite_localPlasticity(g,i,e)) ) then
 
+ do stage = 1,5
          p = material_phaseAt(g,e); c = material_phaseMemberAt(g,i,e)
 
          plasticState(p)%RKCK45dotState(stage,:,c) = plasticState(p)%dotState(:,c)
@@ -1518,11 +1518,11 @@ subroutine integrateStateRKCK45
            nonlocalBroken = .true.
          if(.not. crystallite_todo(g,i,e)) cycle
 
+ enddo
        endif
    enddo; enddo; enddo
    !$OMP END PARALLEL DO
 
- enddo
 
   if(nonlocalBroken) where(.not. crystallite_localPlasticity) crystallite_todo = .false.
 

From 8f6dc054d0895ec603e4840d70f22592d79beb37 Mon Sep 17 00:00:00 2001
From: Martin Diehl <mail@martin-diehl.net>
Date: Tue, 24 Mar 2020 22:25:40 +0100
Subject: [PATCH 07/21] move up the nonlocal skip

---
 src/crystallite.f90 | 7 +++++--
 1 file changed, 5 insertions(+), 2 deletions(-)

diff --git a/src/crystallite.f90 b/src/crystallite.f90
index 77fd00fe7..01ef8bd27 100644
--- a/src/crystallite.f90
+++ b/src/crystallite.f90
@@ -1524,7 +1524,6 @@ subroutine integrateStateRKCK45
    !$OMP END PARALLEL DO
 
 
-  if(nonlocalBroken) where(.not. crystallite_localPlasticity) crystallite_todo = .false.
 
 !--------------------------------------------------------------------------------------------------
 ! --- STATE UPDATE WITH ERROR ESTIMATE FOR STATE ---
@@ -1533,7 +1532,9 @@ subroutine integrateStateRKCK45
    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)) then
+        
+        if(crystallite_todo(g,i,e) .and. (.not. nonlocalBroken .or. crystallite_localPlasticity(g,i,e)) ) then
+       
        p = material_phaseAt(g,e); c = material_phaseMemberAt(g,i,e)
 
        sizeDotState = plasticState(p)%sizeDotState
@@ -1564,6 +1565,8 @@ subroutine integrateStateRKCK45
    enddo; enddo; enddo
  !$OMP END PARALLEL DO
 
+  if(nonlocalBroken) where(.not. crystallite_localPlasticity) crystallite_todo = .false.
+ 
  call update_state(1.0_pReal)
 
  ! --- relative residui and state convergence ---

From 2dd3ccdad1f5915fe406289517208866cbd464dc Mon Sep 17 00:00:00 2001
From: Martin Diehl <mail@martin-diehl.net>
Date: Wed, 25 Mar 2020 10:22:21 +0100
Subject: [PATCH 08/21] no need to transpose

---
 src/crystallite.f90 | 14 +++++++-------
 1 file changed, 7 insertions(+), 7 deletions(-)

diff --git a/src/crystallite.f90 b/src/crystallite.f90
index 01ef8bd27..70910725d 100644
--- a/src/crystallite.f90
+++ b/src/crystallite.f90
@@ -1532,9 +1532,9 @@ subroutine integrateStateRKCK45
    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. nonlocalBroken .or. crystallite_localPlasticity(g,i,e)) ) then
-       
+
        p = material_phaseAt(g,e); c = material_phaseMemberAt(g,i,e)
 
        sizeDotState = plasticState(p)%sizeDotState
@@ -1542,11 +1542,11 @@ subroutine integrateStateRKCK45
        plasticState(p)%RKCK45dotState(6,:,c) = plasticState (p)%dotState(:,c)
 
        residuum_plastic(1:sizeDotState,g,i,e) = &
-         matmul(transpose(plasticState(p)%RKCK45dotState(1:6,1:sizeDotState,c)),DB) &              ! why transpose? Better to transpose constant DB
+         matmul(DB,plasticState(p)%RKCK45dotState(1:6,1:sizeDotState,c)) &
        * crystallite_subdt(g,i,e)
 
         plasticState(p)%dotState(:,c) =  &
-         matmul(transpose(plasticState(p)%RKCK45dotState(1:6,1:sizeDotState,c)), B)                ! why transpose? Better to transpose constant B
+         matmul(B,plasticState(p)%RKCK45dotState(1:6,1:sizeDotState,c))
 
        do s = 1, phase_Nsources(p)
          sizeDotState = sourceState(p)%p(s)%sizeDotState
@@ -1554,11 +1554,11 @@ subroutine integrateStateRKCK45
          sourceState(p)%p(s)%RKCK45dotState(6,:,c) = sourceState(p)%p(s)%dotState(:,c)
 
          residuum_source(1:sizeDotState,s,g,i,e) = &
-           matmul(transpose(sourceState(p)%p(s)%RKCK45dotState(1:6,1:sizeDotState,c)),DB) &
+           matmul(DB,sourceState(p)%p(s)%RKCK45dotState(1:6,1:sizeDotState,c)) &
          * crystallite_subdt(g,i,e)
 
          sourceState(p)%p(s)%dotState(:,c)  = &
-           matmul(transpose(sourceState(p)%p(s)%RKCK45dotState(1:6,1:sizeDotState,c)),B)
+           matmul(B,sourceState(p)%p(s)%RKCK45dotState(1:6,1:sizeDotState,c))
        enddo
 
      endif
@@ -1566,7 +1566,7 @@ subroutine integrateStateRKCK45
  !$OMP END PARALLEL DO
 
   if(nonlocalBroken) where(.not. crystallite_localPlasticity) crystallite_todo = .false.
- 
+
  call update_state(1.0_pReal)
 
  ! --- relative residui and state convergence ---

From 52ca7cc43ccb050eb48c540e32e6ba74c38cdf52 Mon Sep 17 00:00:00 2001
From: Martin Diehl <mail@martin-diehl.net>
Date: Wed, 25 Mar 2020 10:27:55 +0100
Subject: [PATCH 09/21] only one stateJump per integration

---
 src/crystallite.f90 | 8 --------
 1 file changed, 8 deletions(-)

diff --git a/src/crystallite.f90 b/src/crystallite.f90
index 70910725d..5fc5ca91d 100644
--- a/src/crystallite.f90
+++ b/src/crystallite.f90
@@ -1449,11 +1449,8 @@ subroutine integrateStateRKCK45
 
  call update_dotState(1.0_pReal)
 
- ! --- SECOND TO SIXTH RUNGE KUTTA STEP ---
  nonlocalBroken = .false.
 
-   ! --- state update ---
-
    !$OMP PARALLEL DO PRIVATE(sizeDotState,p,c)
    do e = FEsolving_execElem(1),FEsolving_execElem(2)
      do i = FEsolving_execIP(1),FEsolving_execIP(2)
@@ -1491,11 +1488,6 @@ subroutine integrateStateRKCK45
                                                        * crystallite_subdt(g,i,e)
          enddo
 
-         crystallite_todo(g,i,e) = stateJump(g,i,e)
-         if(.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) &
-           nonlocalBroken = .true.
-         if(.not. crystallite_todo(g,i,e)) cycle
-
          call constitutive_dependentState(crystallite_Fe(1:3,1:3,g,i,e), &
                                           crystallite_Fp(1:3,1:3,g,i,e), &
                                           g, i, e)

From 02774a89d96c1d0d7f8b6ff7e1a48803581c6270 Mon Sep 17 00:00:00 2001
From: Martin Diehl <mail@martin-diehl.net>
Date: Wed, 25 Mar 2020 10:33:41 +0100
Subject: [PATCH 10/21] 2 space indentation

---
 src/crystallite.f90 | 214 ++++++++++++++++++++++----------------------
 1 file changed, 108 insertions(+), 106 deletions(-)

diff --git a/src/crystallite.f90 b/src/crystallite.f90
index 5fc5ca91d..c0e1da027 100644
--- a/src/crystallite.f90
+++ b/src/crystallite.f90
@@ -1404,116 +1404,118 @@ end subroutine integrateStateRK4
 !> adaptive step size  (use 5th order solution to advance = "local extrapolation")
 !--------------------------------------------------------------------------------------------------
 subroutine integrateStateRKCK45
-
- real(pReal), dimension(5,5), parameter :: &
-   A = reshape([&
-     .2_pReal, .075_pReal,   .3_pReal, -11.0_pReal/54.0_pReal,  1631.0_pReal/55296.0_pReal, &
-     .0_pReal, .225_pReal,  -.9_pReal,   2.5_pReal,              175.0_pReal/512.0_pReal, &
-     .0_pReal,   .0_pReal,  1.2_pReal, -70.0_pReal/27.0_pReal,   575.0_pReal/13824.0_pReal, &
-     .0_pReal,   .0_pReal,   .0_pReal,  35.0_pReal/27.0_pReal, 44275.0_pReal/110592.0_pReal, &
-     .0_pReal,   .0_pReal,   .0_pReal,    .0_pReal,              253.0_pReal/4096.0_pReal], &
-     [5,5], order=[2,1])                                                                            !< coefficients in Butcher tableau (used for preliminary integration in stages 2 to 6)
-
- real(pReal), dimension(6), parameter :: &
-   B = &
-     [37.0_pReal/378.0_pReal, .0_pReal, 250.0_pReal/621.0_pReal, &
-     125.0_pReal/594.0_pReal, .0_pReal, 512.0_pReal/1771.0_pReal], &                                !< coefficients in Butcher tableau (used for final integration and error estimate)
-   DB = B - &
-     [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,      0.25_pReal]                       !< coefficients in Butcher tableau (used for final integration and error estimate)
-
- real(pReal), dimension(5), parameter :: &
-   CC = [0.2_pReal, 0.3_pReal, 0.6_pReal, 1.0_pReal, 0.875_pReal]                                    !< coefficients in Butcher tableau (fractions of original time step in stages 2 to 6)
-
- 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
-   n, &
-   p, &
-   c, &
-   s, &
-   sizeDotState
- logical :: &
-   nonlocalBroken
-
- real(pReal), dimension(constitutive_plasticity_maxSizeDotState,            &
-                        homogenization_maxNgrains,discretization_nIP,discretization_nElem) :: &
-   residuum_plastic
- real(pReal), dimension(constitutive_source_maxSizeDotState, &
-                        maxval(phase_Nsources), &
-                        homogenization_maxNgrains,discretization_nIP,discretization_nElem) :: &
-   residuum_source
-
-
- call update_dotState(1.0_pReal)
-
- nonlocalBroken = .false.
-
-   !$OMP PARALLEL DO PRIVATE(sizeDotState,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. nonlocalBroken .or. crystallite_localPlasticity(g,i,e)) ) then
-
- do stage = 1,5
+ 
+  real(pReal), dimension(5,5), parameter :: &
+    A = reshape([&
+      .2_pReal, .075_pReal,   .3_pReal, -11.0_pReal/54.0_pReal,  1631.0_pReal/55296.0_pReal, &
+      .0_pReal, .225_pReal,  -.9_pReal,   2.5_pReal,              175.0_pReal/512.0_pReal, &
+      .0_pReal,   .0_pReal,  1.2_pReal, -70.0_pReal/27.0_pReal,   575.0_pReal/13824.0_pReal, &
+      .0_pReal,   .0_pReal,   .0_pReal,  35.0_pReal/27.0_pReal, 44275.0_pReal/110592.0_pReal, &
+      .0_pReal,   .0_pReal,   .0_pReal,    .0_pReal,              253.0_pReal/4096.0_pReal], &
+      [5,5], order=[2,1])                                                                           !< coefficients in Butcher tableau (used for preliminary integration in stages 2 to 6)
+ 
+  real(pReal), dimension(6), parameter :: &
+    B = &
+      [37.0_pReal/378.0_pReal, .0_pReal, 250.0_pReal/621.0_pReal, &
+      125.0_pReal/594.0_pReal, .0_pReal, 512.0_pReal/1771.0_pReal], &                               !< coefficients in Butcher tableau (used for final integration and error estimate)
+    DB = B - &
+      [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,      0.25_pReal]                      !< coefficients in Butcher tableau (used for final integration and error estimate)
+ 
+  real(pReal), dimension(5), parameter :: &
+    CC = [0.2_pReal, 0.3_pReal, 0.6_pReal, 1.0_pReal, 0.875_pReal]                                  !< coefficients in Butcher tableau (fractions of original time step in stages 2 to 6)
+ 
+  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
+    n, &
+    p, &
+    c, &
+    s, &
+    sizeDotState
+  logical :: &
+    nonlocalBroken
+ 
+  real(pReal), dimension(constitutive_plasticity_maxSizeDotState,            &
+                         homogenization_maxNgrains,discretization_nIP,discretization_nElem) :: &
+    residuum_plastic
+  real(pReal), dimension(constitutive_source_maxSizeDotState, &
+                         maxval(phase_Nsources), &
+                         homogenization_maxNgrains,discretization_nIP,discretization_nElem) :: &
+    residuum_source
+ 
+ 
+  call update_dotState(1.0_pReal)
+ 
+  nonlocalBroken = .false.
+ 
+  !$OMP PARALLEL DO PRIVATE(sizeDotState,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. nonlocalBroken .or. crystallite_localPlasticity(g,i,e)) ) then
+ 
          p = material_phaseAt(g,e); c = material_phaseMemberAt(g,i,e)
-
-         plasticState(p)%RKCK45dotState(stage,:,c) = plasticState(p)%dotState(:,c)
-         plasticState(p)%dotState(:,c) = A(1,stage) * plasticState(p)%RKCK45dotState(1,:,c)
-
-         do s = 1, phase_Nsources(p)
-           sourceState(p)%p(s)%RKCK45dotState(stage,:,c) = sourceState(p)%p(s)%dotState(:,c)
-           sourceState(p)%p(s)%dotState(:,c) = A(1,stage) * sourceState(p)%p(s)%RKCK45dotState(1,:,c)
-         enddo
-
-         do n = 2, stage
-           plasticState(p)%dotState(:,c) = plasticState(p)%dotState(:,c) &
-                                          + A(n,stage) * plasticState(p)%RKCK45dotState(n,:,c)
+         
+         do stage = 1,5
+ 
+           plasticState(p)%RKCK45dotState(stage,:,c) = plasticState(p)%dotState(:,c)
+           plasticState(p)%dotState(:,c) = A(1,stage) * plasticState(p)%RKCK45dotState(1,:,c)
+ 
            do s = 1, phase_Nsources(p)
-             sourceState(p)%p(s)%dotState(:,c) = sourceState(p)%p(s)%dotState(:,c) &
-                                                + A(n,stage) * sourceState(p)%p(s)%RKCK45dotState(n,:,c)
+             sourceState(p)%p(s)%RKCK45dotState(stage,:,c) = sourceState(p)%p(s)%dotState(:,c)
+             sourceState(p)%p(s)%dotState(:,c) = A(1,stage) * sourceState(p)%p(s)%RKCK45dotState(1,:,c)
            enddo
+ 
+           do n = 2, stage
+             plasticState(p)%dotState(:,c) = plasticState(p)%dotState(:,c) &
+                                            + A(n,stage) * plasticState(p)%RKCK45dotState(n,:,c)
+             do s = 1, phase_Nsources(p)
+               sourceState(p)%p(s)%dotState(:,c) = sourceState(p)%p(s)%dotState(:,c) &
+                                                  + A(n,stage) * sourceState(p)%p(s)%RKCK45dotState(n,:,c)
+             enddo
+           enddo
+ 
+           sizeDotState = plasticState(p)%sizeDotState
+           plasticState(p)%state(1:sizeDotState,c) = plasticState(p)%subState0(1:sizeDotState,c) &
+                                                   + plasticState(p)%dotState (1:sizeDotState,c) &
+                                                   * crystallite_subdt(g,i,e)
+           do s = 1, phase_Nsources(p)
+             sizeDotState = sourceState(p)%p(s)%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)
+           enddo
+ 
+           call constitutive_dependentState(crystallite_Fe(1:3,1:3,g,i,e), &
+                                            crystallite_Fp(1:3,1:3,g,i,e), &
+                                            g, i, e)
+ 
+           crystallite_todo(g,i,e) = integrateStress(g,i,e,CC(stage))
+           if(.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) &
+             nonlocalBroken = .true.
+           if(.not. crystallite_todo(g,i,e)) cycle
+ 
+           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)*CC(stage), 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)) cycle
+ 
          enddo
-
-         sizeDotState = plasticState(p)%sizeDotState
-         plasticState(p)%state(1:sizeDotState,c) = plasticState(p)%subState0(1:sizeDotState,c) &
-                                                 + plasticState(p)%dotState (1:sizeDotState,c) &
-                                                 * crystallite_subdt(g,i,e)
-         do s = 1, phase_Nsources(p)
-           sizeDotState = sourceState(p)%p(s)%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)
-         enddo
-
-         call constitutive_dependentState(crystallite_Fe(1:3,1:3,g,i,e), &
-                                          crystallite_Fp(1:3,1:3,g,i,e), &
-                                          g, i, e)
-
-         crystallite_todo(g,i,e) = integrateStress(g,i,e,CC(stage))
-         if(.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) &
-           nonlocalBroken = .true.
-         if(.not. crystallite_todo(g,i,e)) cycle
-
-         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)*CC(stage), 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)) cycle
-
- enddo
-       endif
-   enddo; enddo; enddo
-   !$OMP END PARALLEL DO
+      
+      endif
+  enddo; enddo; enddo
+  !$OMP END PARALLEL DO
 
 
 

From 9a188784e2498743772906afa18354f4cfd8d372 Mon Sep 17 00:00:00 2001
From: Martin Diehl <mail@martin-diehl.net>
Date: Wed, 25 Mar 2020 10:40:53 +0100
Subject: [PATCH 11/21] no need for an extra loop

---
 src/crystallite.f90 | 219 ++++++++++++++++++++------------------------
 1 file changed, 101 insertions(+), 118 deletions(-)

diff --git a/src/crystallite.f90 b/src/crystallite.f90
index c0e1da027..f230aebeb 100644
--- a/src/crystallite.f90
+++ b/src/crystallite.f90
@@ -1455,146 +1455,129 @@ subroutine integrateStateRKCK45
   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. nonlocalBroken .or. crystallite_localPlasticity(g,i,e)) ) then
+        if(crystallite_todo(g,i,e) .and. (.not. nonlocalBroken .or. crystallite_localPlasticity(g,i,e)) ) then
  
-         p = material_phaseAt(g,e); c = material_phaseMemberAt(g,i,e)
-         
-         do stage = 1,5
+          p = material_phaseAt(g,e); c = material_phaseMemberAt(g,i,e)
+          
+          do stage = 1,5
  
-           plasticState(p)%RKCK45dotState(stage,:,c) = plasticState(p)%dotState(:,c)
-           plasticState(p)%dotState(:,c) = A(1,stage) * plasticState(p)%RKCK45dotState(1,:,c)
+            plasticState(p)%RKCK45dotState(stage,:,c) = plasticState(p)%dotState(:,c)
+            plasticState(p)%dotState(:,c) = A(1,stage) * plasticState(p)%RKCK45dotState(1,:,c)
  
-           do s = 1, phase_Nsources(p)
-             sourceState(p)%p(s)%RKCK45dotState(stage,:,c) = sourceState(p)%p(s)%dotState(:,c)
-             sourceState(p)%p(s)%dotState(:,c) = A(1,stage) * sourceState(p)%p(s)%RKCK45dotState(1,:,c)
-           enddo
+            do s = 1, phase_Nsources(p)
+              sourceState(p)%p(s)%RKCK45dotState(stage,:,c) = sourceState(p)%p(s)%dotState(:,c)
+              sourceState(p)%p(s)%dotState(:,c) = A(1,stage) * sourceState(p)%p(s)%RKCK45dotState(1,:,c)
+            enddo
  
-           do n = 2, stage
-             plasticState(p)%dotState(:,c) = plasticState(p)%dotState(:,c) &
-                                            + A(n,stage) * plasticState(p)%RKCK45dotState(n,:,c)
-             do s = 1, phase_Nsources(p)
-               sourceState(p)%p(s)%dotState(:,c) = sourceState(p)%p(s)%dotState(:,c) &
-                                                  + A(n,stage) * sourceState(p)%p(s)%RKCK45dotState(n,:,c)
-             enddo
-           enddo
+            do n = 2, stage
+              plasticState(p)%dotState(:,c) = plasticState(p)%dotState(:,c) &
+                                             + A(n,stage) * plasticState(p)%RKCK45dotState(n,:,c)
+              do s = 1, phase_Nsources(p)
+                sourceState(p)%p(s)%dotState(:,c) = sourceState(p)%p(s)%dotState(:,c) &
+                                                   + A(n,stage) * sourceState(p)%p(s)%RKCK45dotState(n,:,c)
+              enddo
+            enddo
  
-           sizeDotState = plasticState(p)%sizeDotState
-           plasticState(p)%state(1:sizeDotState,c) = plasticState(p)%subState0(1:sizeDotState,c) &
-                                                   + plasticState(p)%dotState (1:sizeDotState,c) &
-                                                   * crystallite_subdt(g,i,e)
-           do s = 1, phase_Nsources(p)
-             sizeDotState = sourceState(p)%p(s)%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)
-           enddo
+            sizeDotState = plasticState(p)%sizeDotState
+            plasticState(p)%state(1:sizeDotState,c) = plasticState(p)%subState0(1:sizeDotState,c) &
+                                                    + plasticState(p)%dotState (1:sizeDotState,c) &
+                                                    * crystallite_subdt(g,i,e)
+            do s = 1, phase_Nsources(p)
+              sizeDotState = sourceState(p)%p(s)%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)
+            enddo
  
-           call constitutive_dependentState(crystallite_Fe(1:3,1:3,g,i,e), &
-                                            crystallite_Fp(1:3,1:3,g,i,e), &
-                                            g, i, e)
+            call constitutive_dependentState(crystallite_Fe(1:3,1:3,g,i,e), &
+                                             crystallite_Fp(1:3,1:3,g,i,e), &
+                                             g, i, e)
  
-           crystallite_todo(g,i,e) = integrateStress(g,i,e,CC(stage))
-           if(.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) &
-             nonlocalBroken = .true.
-           if(.not. crystallite_todo(g,i,e)) cycle
+            crystallite_todo(g,i,e) = integrateStress(g,i,e,CC(stage))
+            if(.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) &
+              nonlocalBroken = .true.
+            if(.not. crystallite_todo(g,i,e)) cycle
  
-           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)*CC(stage), 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)) cycle
+            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)*CC(stage), 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)) cycle
  
-         enddo
+          enddo
       
-      endif
+          if(.not. crystallite_todo(g,i,e)) exit
+
+          sizeDotState = plasticState(p)%sizeDotState
+          plasticState(p)%RKCK45dotState(6,:,c) = plasticState (p)%dotState(:,c)
+
+          residuum_plastic(1:sizeDotState,g,i,e) = &
+            matmul(DB,plasticState(p)%RKCK45dotState(1:6,1:sizeDotState,c)) &
+          * crystallite_subdt(g,i,e)
+
+           plasticState(p)%dotState(:,c) =  &
+            matmul(B,plasticState(p)%RKCK45dotState(1:6,1:sizeDotState,c))
+
+          do s = 1, phase_Nsources(p)
+            sizeDotState = sourceState(p)%p(s)%sizeDotState
+            sourceState(p)%p(s)%RKCK45dotState(6,:,c) = sourceState(p)%p(s)%dotState(:,c)
+
+            residuum_source(1:sizeDotState,s,g,i,e) = &
+              matmul(DB,sourceState(p)%p(s)%RKCK45dotState(1:6,1:sizeDotState,c)) &
+            * crystallite_subdt(g,i,e)
+
+            sourceState(p)%p(s)%dotState(:,c)  = &
+              matmul(B,sourceState(p)%p(s)%RKCK45dotState(1:6,1:sizeDotState,c))
+          enddo
+
+        endif
   enddo; enddo; enddo
   !$OMP END PARALLEL DO
 
-
-
-!--------------------------------------------------------------------------------------------------
-! --- STATE UPDATE WITH ERROR ESTIMATE FOR STATE ---
-
- !$OMP PARALLEL DO PRIVATE(sizeDotState,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. nonlocalBroken .or. crystallite_localPlasticity(g,i,e)) ) then
-
-       p = material_phaseAt(g,e); c = material_phaseMemberAt(g,i,e)
-
-       sizeDotState = plasticState(p)%sizeDotState
-
-       plasticState(p)%RKCK45dotState(6,:,c) = plasticState (p)%dotState(:,c)
-
-       residuum_plastic(1:sizeDotState,g,i,e) = &
-         matmul(DB,plasticState(p)%RKCK45dotState(1:6,1:sizeDotState,c)) &
-       * crystallite_subdt(g,i,e)
-
-        plasticState(p)%dotState(:,c) =  &
-         matmul(B,plasticState(p)%RKCK45dotState(1:6,1:sizeDotState,c))
-
-       do s = 1, phase_Nsources(p)
-         sizeDotState = sourceState(p)%p(s)%sizeDotState
-
-         sourceState(p)%p(s)%RKCK45dotState(6,:,c) = sourceState(p)%p(s)%dotState(:,c)
-
-         residuum_source(1:sizeDotState,s,g,i,e) = &
-           matmul(DB,sourceState(p)%p(s)%RKCK45dotState(1:6,1:sizeDotState,c)) &
-         * crystallite_subdt(g,i,e)
-
-         sourceState(p)%p(s)%dotState(:,c)  = &
-           matmul(B,sourceState(p)%p(s)%RKCK45dotState(1:6,1:sizeDotState,c))
-       enddo
-
-     endif
-   enddo; enddo; enddo
- !$OMP END PARALLEL DO
-
   if(nonlocalBroken) where(.not. crystallite_localPlasticity) crystallite_todo = .false.
 
- call update_state(1.0_pReal)
+  call update_state(1.0_pReal)
 
- ! --- relative residui and state convergence ---
+  ! --- relative residui and state convergence ---
 
- !$OMP PARALLEL DO PRIVATE(sizeDotState,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)) then
-         p  = material_phaseAt(g,e); c = material_phaseMemberAt(g,i,e)
+  !$OMP PARALLEL DO PRIVATE(sizeDotState,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)) then
+          p  = material_phaseAt(g,e); c = material_phaseMemberAt(g,i,e)
 
-         sizeDotState = plasticState(p)%sizeDotState
+          sizeDotState = plasticState(p)%sizeDotState
 
-         crystallite_todo(g,i,e) = converged(residuum_plastic(1:sizeDotState,g,i,e), &
-                                             plasticState(p)%state(1:sizeDotState,c), &
-                                             plasticState(p)%atol(1:sizeDotState))
+          crystallite_todo(g,i,e) = converged(residuum_plastic(1:sizeDotState,g,i,e), &
+                                              plasticState(p)%state(1:sizeDotState,c), &
+                                              plasticState(p)%atol(1:sizeDotState))
 
-         do s = 1, phase_Nsources(p)
-           sizeDotState = sourceState(p)%p(s)%sizeDotState
+          do s = 1, phase_Nsources(p)
+            sizeDotState = sourceState(p)%p(s)%sizeDotState
 
-                  crystallite_todo(g,i,e) = &
-                  crystallite_todo(g,i,e) .and. converged(residuum_source(1:sizeDotState,s,g,i,e), &
-                                                          sourceState(p)%p(s)%state(1:sizeDotState,c), &
-                                                          sourceState(p)%p(s)%atol(1:sizeDotState))
-         enddo
-     endif
-   enddo; enddo; enddo
- !$OMP END PARALLEL DO
+                   crystallite_todo(g,i,e) = &
+                   crystallite_todo(g,i,e) .and. converged(residuum_source(1:sizeDotState,s,g,i,e), &
+                                                           sourceState(p)%p(s)%state(1:sizeDotState,c), &
+                                                           sourceState(p)%p(s)%atol(1:sizeDotState))
+          enddo
+      endif
+    enddo; enddo; enddo
+  !$OMP END PARALLEL DO
 
- call update_deltaState
- call update_dependentState
- call update_stress(1.0_pReal)
- call setConvergenceFlag
- if (any(plasticState(:)%nonlocal)) call nonlocalConvergenceCheck
+  call update_deltaState
+  call update_dependentState
+  call update_stress(1.0_pReal)
+  
+  call setConvergenceFlag
+  if (any(plasticState(:)%nonlocal)) call nonlocalConvergenceCheck
 
 end subroutine integrateStateRKCK45
 

From 32b9a5ab1587e1ff86a1cd82d34ea106d1a76933 Mon Sep 17 00:00:00 2001
From: Martin Diehl <mail@martin-diehl.net>
Date: Wed, 25 Mar 2020 11:07:47 +0100
Subject: [PATCH 12/21] all in one loop

---
 src/crystallite.f90 | 55 ++++++++++++++++++++++++++-------------------
 1 file changed, 32 insertions(+), 23 deletions(-)

diff --git a/src/crystallite.f90 b/src/crystallite.f90
index f230aebeb..3d8e4b841 100644
--- a/src/crystallite.f90
+++ b/src/crystallite.f90
@@ -1404,7 +1404,7 @@ end subroutine integrateStateRK4
 !> adaptive step size  (use 5th order solution to advance = "local extrapolation")
 !--------------------------------------------------------------------------------------------------
 subroutine integrateStateRKCK45
- 
+
   real(pReal), dimension(5,5), parameter :: &
     A = reshape([&
       .2_pReal, .075_pReal,   .3_pReal, -11.0_pReal/54.0_pReal,  1631.0_pReal/55296.0_pReal, &
@@ -1413,7 +1413,7 @@ subroutine integrateStateRKCK45
       .0_pReal,   .0_pReal,   .0_pReal,  35.0_pReal/27.0_pReal, 44275.0_pReal/110592.0_pReal, &
       .0_pReal,   .0_pReal,   .0_pReal,    .0_pReal,              253.0_pReal/4096.0_pReal], &
       [5,5], order=[2,1])                                                                           !< coefficients in Butcher tableau (used for preliminary integration in stages 2 to 6)
- 
+
   real(pReal), dimension(6), parameter :: &
     B = &
       [37.0_pReal/378.0_pReal, .0_pReal, 250.0_pReal/621.0_pReal, &
@@ -1421,10 +1421,10 @@ subroutine integrateStateRKCK45
     DB = B - &
       [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,      0.25_pReal]                      !< coefficients in Butcher tableau (used for final integration and error estimate)
- 
+
   real(pReal), dimension(5), parameter :: &
     CC = [0.2_pReal, 0.3_pReal, 0.6_pReal, 1.0_pReal, 0.875_pReal]                                  !< coefficients in Butcher tableau (fractions of original time step in stages 2 to 6)
- 
+
   integer :: &
     e, &                                                                                            ! element index in element loop
     i, &                                                                                            ! integration point index in ip loop
@@ -1437,7 +1437,7 @@ subroutine integrateStateRKCK45
     sizeDotState
   logical :: &
     nonlocalBroken
- 
+
   real(pReal), dimension(constitutive_plasticity_maxSizeDotState,            &
                          homogenization_maxNgrains,discretization_nIP,discretization_nElem) :: &
     residuum_plastic
@@ -1445,30 +1445,30 @@ subroutine integrateStateRKCK45
                          maxval(phase_Nsources), &
                          homogenization_maxNgrains,discretization_nIP,discretization_nElem) :: &
     residuum_source
- 
- 
+
+
   call update_dotState(1.0_pReal)
- 
+
   nonlocalBroken = .false.
- 
+
   !$OMP PARALLEL DO PRIVATE(sizeDotState,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. nonlocalBroken .or. crystallite_localPlasticity(g,i,e)) ) then
- 
+
           p = material_phaseAt(g,e); c = material_phaseMemberAt(g,i,e)
-          
+
           do stage = 1,5
- 
+
             plasticState(p)%RKCK45dotState(stage,:,c) = plasticState(p)%dotState(:,c)
             plasticState(p)%dotState(:,c) = A(1,stage) * plasticState(p)%RKCK45dotState(1,:,c)
- 
+
             do s = 1, phase_Nsources(p)
               sourceState(p)%p(s)%RKCK45dotState(stage,:,c) = sourceState(p)%p(s)%dotState(:,c)
               sourceState(p)%p(s)%dotState(:,c) = A(1,stage) * sourceState(p)%p(s)%RKCK45dotState(1,:,c)
             enddo
- 
+
             do n = 2, stage
               plasticState(p)%dotState(:,c) = plasticState(p)%dotState(:,c) &
                                              + A(n,stage) * plasticState(p)%RKCK45dotState(n,:,c)
@@ -1477,7 +1477,7 @@ subroutine integrateStateRKCK45
                                                    + A(n,stage) * sourceState(p)%p(s)%RKCK45dotState(n,:,c)
               enddo
             enddo
- 
+
             sizeDotState = plasticState(p)%sizeDotState
             plasticState(p)%state(1:sizeDotState,c) = plasticState(p)%subState0(1:sizeDotState,c) &
                                                     + plasticState(p)%dotState (1:sizeDotState,c) &
@@ -1488,16 +1488,16 @@ subroutine integrateStateRKCK45
                                                           + sourceState(p)%p(s)%dotState (1:sizeDotState,c) &
                                                           * crystallite_subdt(g,i,e)
             enddo
- 
+
             call constitutive_dependentState(crystallite_Fe(1:3,1:3,g,i,e), &
                                              crystallite_Fp(1:3,1:3,g,i,e), &
                                              g, i, e)
- 
+
             crystallite_todo(g,i,e) = integrateStress(g,i,e,CC(stage))
             if(.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) &
               nonlocalBroken = .true.
             if(.not. crystallite_todo(g,i,e)) cycle
- 
+
             call constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), &
                                               crystallite_partionedF0, &
                                               crystallite_Fi(1:3,1:3,g,i,e), &
@@ -1510,9 +1510,9 @@ subroutine integrateStateRKCK45
             if(.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) &
               nonlocalBroken = .true.
             if(.not. crystallite_todo(g,i,e)) cycle
- 
+
           enddo
-      
+
           if(.not. crystallite_todo(g,i,e)) exit
 
           sizeDotState = plasticState(p)%sizeDotState
@@ -1537,14 +1537,23 @@ subroutine integrateStateRKCK45
               matmul(B,sourceState(p)%p(s)%RKCK45dotState(1:6,1:sizeDotState,c))
           enddo
 
+          sizeDotState = plasticState(p)%sizeDotState
+          plasticState(p)%state(1:sizeDotState,c) = plasticState(p)%subState0(1:sizeDotState,c) &
+                                                  + plasticState(p)%dotState (1:sizeDotState,c) &
+                                                    * crystallite_subdt(g,i,e)
+          do s = 1, phase_Nsources(p)
+            sizeDotState = sourceState(p)%p(s)%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)
+          enddo
+
         endif
   enddo; enddo; enddo
   !$OMP END PARALLEL DO
 
   if(nonlocalBroken) where(.not. crystallite_localPlasticity) crystallite_todo = .false.
 
-  call update_state(1.0_pReal)
-
   ! --- relative residui and state convergence ---
 
   !$OMP PARALLEL DO PRIVATE(sizeDotState,p,c)
@@ -1575,7 +1584,7 @@ subroutine integrateStateRKCK45
   call update_deltaState
   call update_dependentState
   call update_stress(1.0_pReal)
-  
+
   call setConvergenceFlag
   if (any(plasticState(:)%nonlocal)) call nonlocalConvergenceCheck
 

From 67f64cf7e1c708f7c222f0d9b6fe89784b1c962d Mon Sep 17 00:00:00 2001
From: Martin Diehl <mail@martin-diehl.net>
Date: Wed, 25 Mar 2020 11:13:46 +0100
Subject: [PATCH 13/21] correct iteration skipping

---
 src/crystallite.f90 | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/src/crystallite.f90 b/src/crystallite.f90
index 3d8e4b841..f32bd3c36 100644
--- a/src/crystallite.f90
+++ b/src/crystallite.f90
@@ -1496,7 +1496,7 @@ subroutine integrateStateRKCK45
             crystallite_todo(g,i,e) = integrateStress(g,i,e,CC(stage))
             if(.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) &
               nonlocalBroken = .true.
-            if(.not. crystallite_todo(g,i,e)) cycle
+            if(.not. crystallite_todo(g,i,e)) exit
 
             call constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), &
                                               crystallite_partionedF0, &
@@ -1509,11 +1509,11 @@ subroutine integrateStateRKCK45
             enddo
             if(.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) &
               nonlocalBroken = .true.
-            if(.not. crystallite_todo(g,i,e)) cycle
+            if(.not. crystallite_todo(g,i,e)) exit
 
           enddo
 
-          if(.not. crystallite_todo(g,i,e)) exit
+          if(.not. crystallite_todo(g,i,e)) cycle
 
           sizeDotState = plasticState(p)%sizeDotState
           plasticState(p)%RKCK45dotState(6,:,c) = plasticState (p)%dotState(:,c)

From 652846cdc9c48b2341b40928bd393d7431eb20cd Mon Sep 17 00:00:00 2001
From: Martin Diehl <mail@martin-diehl.net>
Date: Wed, 25 Mar 2020 11:16:26 +0100
Subject: [PATCH 14/21] no need for extra loop

---
 src/crystallite.f90 | 25 ++++++++-----------------
 1 file changed, 8 insertions(+), 17 deletions(-)

diff --git a/src/crystallite.f90 b/src/crystallite.f90
index f32bd3c36..6d7d6e80b 100644
--- a/src/crystallite.f90
+++ b/src/crystallite.f90
@@ -1548,21 +1548,6 @@ subroutine integrateStateRKCK45
                                                           * crystallite_subdt(g,i,e)
           enddo
 
-        endif
-  enddo; enddo; enddo
-  !$OMP END PARALLEL DO
-
-  if(nonlocalBroken) where(.not. crystallite_localPlasticity) crystallite_todo = .false.
-
-  ! --- relative residui and state convergence ---
-
-  !$OMP PARALLEL DO PRIVATE(sizeDotState,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)) then
-          p  = material_phaseAt(g,e); c = material_phaseMemberAt(g,i,e)
-
           sizeDotState = plasticState(p)%sizeDotState
 
           crystallite_todo(g,i,e) = converged(residuum_plastic(1:sizeDotState,g,i,e), &
@@ -1577,10 +1562,16 @@ subroutine integrateStateRKCK45
                                                            sourceState(p)%p(s)%state(1:sizeDotState,c), &
                                                            sourceState(p)%p(s)%atol(1:sizeDotState))
           enddo
-      endif
-    enddo; enddo; enddo
+          if(.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) &
+            nonlocalBroken = .true.
+          if(.not. crystallite_todo(g,i,e)) cycle
+
+        endif
+  enddo; enddo; enddo
   !$OMP END PARALLEL DO
 
+  if(nonlocalBroken) where(.not. crystallite_localPlasticity) crystallite_todo = .false.
+
   call update_deltaState
   call update_dependentState
   call update_stress(1.0_pReal)

From b54a109d990db79ae8eabe51d4648f253a2fda31 Mon Sep 17 00:00:00 2001
From: Martin Diehl <mail@martin-diehl.net>
Date: Wed, 25 Mar 2020 11:20:39 +0100
Subject: [PATCH 15/21] do plasticState and sourceState at once

---
 src/crystallite.f90 | 15 ++++-----------
 1 file changed, 4 insertions(+), 11 deletions(-)

diff --git a/src/crystallite.f90 b/src/crystallite.f90
index 6d7d6e80b..64998794e 100644
--- a/src/crystallite.f90
+++ b/src/crystallite.f90
@@ -1541,27 +1541,20 @@ subroutine integrateStateRKCK45
           plasticState(p)%state(1:sizeDotState,c) = plasticState(p)%subState0(1:sizeDotState,c) &
                                                   + plasticState(p)%dotState (1:sizeDotState,c) &
                                                     * crystallite_subdt(g,i,e)
+          crystallite_todo(g,i,e) = converged(residuum_plastic(1:sizeDotState,g,i,e), &
+                                              plasticState(p)%state(1:sizeDotState,c), &
+                                              plasticState(p)%atol(1:sizeDotState))
           do s = 1, phase_Nsources(p)
             sizeDotState = sourceState(p)%p(s)%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)
-          enddo
-
-          sizeDotState = plasticState(p)%sizeDotState
-
-          crystallite_todo(g,i,e) = converged(residuum_plastic(1:sizeDotState,g,i,e), &
-                                              plasticState(p)%state(1:sizeDotState,c), &
-                                              plasticState(p)%atol(1:sizeDotState))
-
-          do s = 1, phase_Nsources(p)
-            sizeDotState = sourceState(p)%p(s)%sizeDotState
-
                    crystallite_todo(g,i,e) = &
                    crystallite_todo(g,i,e) .and. converged(residuum_source(1:sizeDotState,s,g,i,e), &
                                                            sourceState(p)%p(s)%state(1:sizeDotState,c), &
                                                            sourceState(p)%p(s)%atol(1:sizeDotState))
           enddo
+
           if(.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) &
             nonlocalBroken = .true.
           if(.not. crystallite_todo(g,i,e)) cycle

From 86abba477a285c79b54f70793167fb99844ff7ed Mon Sep 17 00:00:00 2001
From: Martin Diehl <mail@martin-diehl.net>
Date: Wed, 25 Mar 2020 11:27:03 +0100
Subject: [PATCH 16/21] use same loop

---
 src/crystallite.f90 | 16 +++++++++++++---
 1 file changed, 13 insertions(+), 3 deletions(-)

diff --git a/src/crystallite.f90 b/src/crystallite.f90
index 64998794e..26607a381 100644
--- a/src/crystallite.f90
+++ b/src/crystallite.f90
@@ -1447,10 +1447,7 @@ subroutine integrateStateRKCK45
     residuum_source
 
 
-  call update_dotState(1.0_pReal)
-
   nonlocalBroken = .false.
-
   !$OMP PARALLEL DO PRIVATE(sizeDotState,p,c)
   do e = FEsolving_execElem(1),FEsolving_execElem(2)
     do i = FEsolving_execIP(1),FEsolving_execIP(2)
@@ -1459,6 +1456,19 @@ subroutine integrateStateRKCK45
 
           p = material_phaseAt(g,e); c = material_phaseMemberAt(g,i,e)
 
+          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)) cycle
+
           do stage = 1,5
 
             plasticState(p)%RKCK45dotState(stage,:,c) = plasticState(p)%dotState(:,c)

From 939d2af1d9309462b069b0e4fcd2e3915ee5d3d6 Mon Sep 17 00:00:00 2001
From: Martin Diehl <mail@martin-diehl.net>
Date: Wed, 25 Mar 2020 18:58:58 +0100
Subject: [PATCH 17/21] bugfix: explicit nonlocal requires F, not Fe

---
 src/crystallite.f90 | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/crystallite.f90 b/src/crystallite.f90
index 2386d0aee..f5f5a7ca4 100644
--- a/src/crystallite.f90
+++ b/src/crystallite.f90
@@ -1499,7 +1499,7 @@ subroutine integrateStateRKCK45
                                                           * crystallite_subdt(g,i,e)
             enddo
 
-            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)
 

From 599de26dad312ac57e80485a4c5b2ea8dea4b3db Mon Sep 17 00:00:00 2001
From: Martin Diehl <mail@martin-diehl.net>
Date: Wed, 25 Mar 2020 19:44:51 +0100
Subject: [PATCH 18/21] further integration

---
 src/crystallite.f90 | 16 +++++++++++++---
 1 file changed, 13 insertions(+), 3 deletions(-)

diff --git a/src/crystallite.f90 b/src/crystallite.f90
index f5f5a7ca4..16666f490 100644
--- a/src/crystallite.f90
+++ b/src/crystallite.f90
@@ -1564,19 +1564,29 @@ subroutine integrateStateRKCK45
                                                            sourceState(p)%p(s)%state(1:sizeDotState,c), &
                                                            sourceState(p)%p(s)%atol(1:sizeDotState))
           enddo
-
           if(.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) &
             nonlocalBroken = .true.
           if(.not. crystallite_todo(g,i,e)) cycle
 
+          crystallite_todo(g,i,e) = stateJump(g,i,e)
+          if(.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) &
+            nonlocalBroken = .true.
+          if(.not. crystallite_todo(g,i,e)) cycle
+
+          call constitutive_dependentState(crystallite_partionedF(1:3,1:3,g,i,e), &
+                                           crystallite_Fp(1:3,1:3,g,i,e), &
+                                           g, i, e)
+
+          if(.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) &
+            nonlocalBroken = .true.
+          if(.not. crystallite_todo(g,i,e)) cycle
+        
         endif
   enddo; enddo; enddo
   !$OMP END PARALLEL DO
 
   if(nonlocalBroken) where(.not. crystallite_localPlasticity) crystallite_todo = .false.
 
-  call update_deltaState
-  call update_dependentState
   call update_stress(1.0_pReal)
 
   call setConvergenceFlag

From aa2d4401079bb3baaa559733c8f0f663b512a09d Mon Sep 17 00:00:00 2001
From: Martin Diehl <mail@martin-diehl.net>
Date: Wed, 25 Mar 2020 20:37:37 +0100
Subject: [PATCH 19/21] finally one loop left

---
 src/crystallite.f90 | 11 ++++++-----
 1 file changed, 6 insertions(+), 5 deletions(-)

diff --git a/src/crystallite.f90 b/src/crystallite.f90
index 16666f490..dfabf51d6 100644
--- a/src/crystallite.f90
+++ b/src/crystallite.f90
@@ -1580,16 +1580,17 @@ subroutine integrateStateRKCK45
           if(.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) &
             nonlocalBroken = .true.
           if(.not. crystallite_todo(g,i,e)) cycle
-        
+
+          crystallite_todo(g,i,e) = integrateStress(g,i,e)
+          if(.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) &
+            nonlocalBroken = .true.
+          crystallite_converged(g,i,e) = crystallite_todo(g,i,e)                                    ! consider converged if not broken
+
         endif
   enddo; enddo; enddo
   !$OMP END PARALLEL DO
 
   if(nonlocalBroken) where(.not. crystallite_localPlasticity) crystallite_todo = .false.
-
-  call update_stress(1.0_pReal)
-
-  call setConvergenceFlag
   if (any(plasticState(:)%nonlocal)) call nonlocalConvergenceCheck
 
 end subroutine integrateStateRKCK45

From 3ce8902245b88ed0196cbae30469458fbff95b03 Mon Sep 17 00:00:00 2001
From: Martin Diehl <mail@martin-diehl.net>
Date: Wed, 25 Mar 2020 20:47:10 +0100
Subject: [PATCH 20/21] no need to split up

---
 src/crystallite.f90 | 18 +++++++-----------
 1 file changed, 7 insertions(+), 11 deletions(-)

diff --git a/src/crystallite.f90 b/src/crystallite.f90
index dfabf51d6..821122c24 100644
--- a/src/crystallite.f90
+++ b/src/crystallite.f90
@@ -1535,6 +1535,13 @@ subroutine integrateStateRKCK45
            plasticState(p)%dotState(:,c) =  &
             matmul(B,plasticState(p)%RKCK45dotState(1:6,1:sizeDotState,c))
 
+          plasticState(p)%state(1:sizeDotState,c) = plasticState(p)%subState0(1:sizeDotState,c) &
+                                                  + plasticState(p)%dotState (1:sizeDotState,c) &
+                                                    * crystallite_subdt(g,i,e)
+          crystallite_todo(g,i,e) = converged(residuum_plastic(1:sizeDotState,g,i,e), &
+                                              plasticState(p)%state(1:sizeDotState,c), &
+                                              plasticState(p)%atol(1:sizeDotState))
+
           do s = 1, phase_Nsources(p)
             sizeDotState = sourceState(p)%p(s)%sizeDotState
             sourceState(p)%p(s)%RKCK45dotState(6,:,c) = sourceState(p)%p(s)%dotState(:,c)
@@ -1545,17 +1552,6 @@ subroutine integrateStateRKCK45
 
             sourceState(p)%p(s)%dotState(:,c)  = &
               matmul(B,sourceState(p)%p(s)%RKCK45dotState(1:6,1:sizeDotState,c))
-          enddo
-
-          sizeDotState = plasticState(p)%sizeDotState
-          plasticState(p)%state(1:sizeDotState,c) = plasticState(p)%subState0(1:sizeDotState,c) &
-                                                  + plasticState(p)%dotState (1:sizeDotState,c) &
-                                                    * crystallite_subdt(g,i,e)
-          crystallite_todo(g,i,e) = converged(residuum_plastic(1:sizeDotState,g,i,e), &
-                                              plasticState(p)%state(1:sizeDotState,c), &
-                                              plasticState(p)%atol(1:sizeDotState))
-          do s = 1, phase_Nsources(p)
-            sizeDotState = sourceState(p)%p(s)%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)

From a4674a6461df3e4f45960817ef4108f766b7fd82 Mon Sep 17 00:00:00 2001
From: Martin Diehl <mail@martin-diehl.net>
Date: Wed, 25 Mar 2020 20:55:37 +0100
Subject: [PATCH 21/21] whitespace adjustments

---
 src/crystallite.f90 | 40 ++++++++++++++++------------------------
 1 file changed, 16 insertions(+), 24 deletions(-)

diff --git a/src/crystallite.f90 b/src/crystallite.f90
index 821122c24..47b675f33 100644
--- a/src/crystallite.f90
+++ b/src/crystallite.f90
@@ -1473,7 +1473,6 @@ subroutine integrateStateRKCK45
 
             plasticState(p)%RKCK45dotState(stage,:,c) = plasticState(p)%dotState(:,c)
             plasticState(p)%dotState(:,c) = A(1,stage) * plasticState(p)%RKCK45dotState(1,:,c)
-
             do s = 1, phase_Nsources(p)
               sourceState(p)%p(s)%RKCK45dotState(stage,:,c) = sourceState(p)%p(s)%dotState(:,c)
               sourceState(p)%p(s)%dotState(:,c) = A(1,stage) * sourceState(p)%p(s)%RKCK45dotState(1,:,c)
@@ -1481,22 +1480,22 @@ subroutine integrateStateRKCK45
 
             do n = 2, stage
               plasticState(p)%dotState(:,c) = plasticState(p)%dotState(:,c) &
-                                             + A(n,stage) * plasticState(p)%RKCK45dotState(n,:,c)
+                                            + A(n,stage) * plasticState(p)%RKCK45dotState(n,:,c)
               do s = 1, phase_Nsources(p)
                 sourceState(p)%p(s)%dotState(:,c) = sourceState(p)%p(s)%dotState(:,c) &
-                                                   + A(n,stage) * sourceState(p)%p(s)%RKCK45dotState(n,:,c)
+                                                  + A(n,stage) * sourceState(p)%p(s)%RKCK45dotState(n,:,c)
               enddo
             enddo
 
             sizeDotState = plasticState(p)%sizeDotState
             plasticState(p)%state(1:sizeDotState,c) = plasticState(p)%subState0(1:sizeDotState,c) &
                                                     + plasticState(p)%dotState (1:sizeDotState,c) &
-                                                    * crystallite_subdt(g,i,e)
+                                                      * crystallite_subdt(g,i,e)
             do s = 1, phase_Nsources(p)
               sizeDotState = sourceState(p)%p(s)%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)
+                                                            * crystallite_subdt(g,i,e)
             enddo
 
             call constitutive_dependentState(crystallite_partionedF(1:3,1:3,g,i,e), &
@@ -1526,15 +1525,11 @@ subroutine integrateStateRKCK45
           if(.not. crystallite_todo(g,i,e)) cycle
 
           sizeDotState = plasticState(p)%sizeDotState
+
           plasticState(p)%RKCK45dotState(6,:,c) = plasticState (p)%dotState(:,c)
-
-          residuum_plastic(1:sizeDotState,g,i,e) = &
-            matmul(DB,plasticState(p)%RKCK45dotState(1:6,1:sizeDotState,c)) &
-          * crystallite_subdt(g,i,e)
-
-           plasticState(p)%dotState(:,c) =  &
-            matmul(B,plasticState(p)%RKCK45dotState(1:6,1:sizeDotState,c))
-
+          residuum_plastic(1:sizeDotState,g,i,e) = matmul(DB,plasticState(p)%RKCK45dotState(1:6,1:sizeDotState,c)) &
+                                                   * crystallite_subdt(g,i,e)
+          plasticState(p)%dotState(:,c) =  matmul(B,plasticState(p)%RKCK45dotState(1:6,1:sizeDotState,c))
           plasticState(p)%state(1:sizeDotState,c) = plasticState(p)%subState0(1:sizeDotState,c) &
                                                   + plasticState(p)%dotState (1:sizeDotState,c) &
                                                     * crystallite_subdt(g,i,e)
@@ -1544,21 +1539,18 @@ subroutine integrateStateRKCK45
 
           do s = 1, phase_Nsources(p)
             sizeDotState = sourceState(p)%p(s)%sizeDotState
+
             sourceState(p)%p(s)%RKCK45dotState(6,:,c) = sourceState(p)%p(s)%dotState(:,c)
-
-            residuum_source(1:sizeDotState,s,g,i,e) = &
-              matmul(DB,sourceState(p)%p(s)%RKCK45dotState(1:6,1:sizeDotState,c)) &
-            * crystallite_subdt(g,i,e)
-
-            sourceState(p)%p(s)%dotState(:,c)  = &
-              matmul(B,sourceState(p)%p(s)%RKCK45dotState(1:6,1:sizeDotState,c))
+            residuum_source(1:sizeDotState,s,g,i,e) = matmul(DB,sourceState(p)%p(s)%RKCK45dotState(1:6,1:sizeDotState,c)) &
+                                                      * crystallite_subdt(g,i,e)
+            sourceState(p)%p(s)%dotState(:,c)  = matmul(B,sourceState(p)%p(s)%RKCK45dotState(1:6,1:sizeDotState,c))
             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)
-                   crystallite_todo(g,i,e) = &
-                   crystallite_todo(g,i,e) .and. converged(residuum_source(1:sizeDotState,s,g,i,e), &
-                                                           sourceState(p)%p(s)%state(1:sizeDotState,c), &
-                                                           sourceState(p)%p(s)%atol(1:sizeDotState))
+            crystallite_todo(g,i,e) = crystallite_todo(g,i,e) .and. &
+                                      converged(residuum_source(1:sizeDotState,s,g,i,e), &
+                                                sourceState(p)%p(s)%state(1:sizeDotState,c), &
+                                                sourceState(p)%p(s)%atol(1:sizeDotState))
           enddo
           if(.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) &
             nonlocalBroken = .true.