no need for separate loop

src/constitutive.f90

@@ -398,6 +398,7 @@ module constitutive
     converged, &
     crystallite_init, &
     crystallite_stress, &
+    crystallite_stress2, &
     crystallite_stressTangent, &
     crystallite_orientations, &
     crystallite_push33ToRef, &
@@ -1152,6 +1153,138 @@ function crystallite_stress()
 end function crystallite_stress
 
 
+!--------------------------------------------------------------------------------------------------
+!> @brief calculate stress (P)
+!--------------------------------------------------------------------------------------------------
+function crystallite_stress2(co,ip,el)
+
+  integer, intent(in) :: &
+    co, &
+    ip, &
+    el
+
+  logical :: crystallite_stress2
+
+  real(pReal) :: &
+    formerSubStep
+  integer :: &
+    NiterationCrystallite, &                                                                        ! number of iterations in crystallite loop
+    s, ph, me
+  logical :: todo
+  real(pReal) :: subFrac     !ToDo: need to set some values to false for different Ngrains
+  real(pReal), dimension(3,3) :: &
+    subLp0, &                                                                                       !< plastic velocity grad at start of crystallite inc
+    subLi0                                                                                          !< intermediate velocity grad at start of crystallite inc
+
+
+
+
+
+!--------------------------------------------------------------------------------------------------
+! initialize to starting condition
+  crystallite_subStep(co,ip,el) = 0.0_pReal
+
+      ph = material_phaseAt(co,el)
+      me = material_phaseMemberAt(co,ip,el)
+      subLi0 = constitutive_mech_partionedLi0(ph)%data(1:3,1:3,me)
+      subLp0 = crystallite_partitionedLp0(1:3,1:3,co,ip,el)
+      homogenizationRequestsCalculation: if (crystallite_requested(co,ip,el)) then
+        plasticState    (material_phaseAt(co,el))%subState0(      :,material_phaseMemberAt(co,ip,el)) = &
+        plasticState    (material_phaseAt(co,el))%partitionedState0(:,material_phaseMemberAt(co,ip,el))
+
+        do s = 1, phase_Nsources(material_phaseAt(co,el))
+          sourceState(material_phaseAt(co,el))%p(s)%subState0(      :,material_phaseMemberAt(co,ip,el)) = &
+          sourceState(material_phaseAt(co,el))%p(s)%partitionedState0(:,material_phaseMemberAt(co,ip,el))
+        enddo
+        crystallite_subFp0(1:3,1:3,co,ip,el) = constitutive_mech_partionedFp0(ph)%data(1:3,1:3,me)
+        crystallite_subFi0(1:3,1:3,co,ip,el) = constitutive_mech_partionedFi0(ph)%data(1:3,1:3,me)
+        crystallite_subF0(1:3,1:3,co,ip,el)  = crystallite_partitionedF0(1:3,1:3,co,ip,el)
+        subFrac = 0.0_pReal
+        crystallite_subStep(co,ip,el) = 1.0_pReal/num%subStepSizeCryst
+        todo = .true.
+        crystallite_converged(co,ip,el) = .false.                                                      ! pretend failed step of 1/subStepSizeCryst
+      endif homogenizationRequestsCalculation
+
+  todo = .true.
+  NiterationCrystallite = 0
+  cutbackLooping: do while (todo)
+    NiterationCrystallite = NiterationCrystallite + 1
+
+!--------------------------------------------------------------------------------------------------
+!  wind forward
+          if (crystallite_converged(co,ip,el)) then
+            formerSubStep = crystallite_subStep(co,ip,el)
+            subFrac = subFrac + crystallite_subStep(co,ip,el)
+            crystallite_subStep(co,ip,el) = min(1.0_pReal - subFrac, &
+                                             num%stepIncreaseCryst * crystallite_subStep(co,ip,el))
+
+            todo = crystallite_subStep(co,ip,el) > 0.0_pReal                        ! still time left to integrate on?
+
+            if (todo) then
+              crystallite_subF0 (1:3,1:3,co,ip,el) = crystallite_subF(1:3,1:3,co,ip,el)
+              subLp0 = crystallite_Lp  (1:3,1:3,co,ip,el)
+              subLi0 = constitutive_mech_Li(ph)%data(1:3,1:3,me)
+              crystallite_subFp0(1:3,1:3,co,ip,el) = constitutive_mech_Fp(ph)%data(1:3,1:3,me)
+              crystallite_subFi0(1:3,1:3,co,ip,el) = constitutive_mech_Fi(ph)%data(1:3,1:3,me)
+              plasticState(    material_phaseAt(co,el))%subState0(:,material_phaseMemberAt(co,ip,el)) &
+                = plasticState(material_phaseAt(co,el))%state(    :,material_phaseMemberAt(co,ip,el))
+              do s = 1, phase_Nsources(material_phaseAt(co,el))
+                sourceState(    material_phaseAt(co,el))%p(s)%subState0(:,material_phaseMemberAt(co,ip,el)) &
+                  = sourceState(material_phaseAt(co,el))%p(s)%state(    :,material_phaseMemberAt(co,ip,el))
+              enddo
+            endif
+
+!--------------------------------------------------------------------------------------------------
+!  cut back (reduced time and restore)
+          else
+            crystallite_subStep(co,ip,el)       = num%subStepSizeCryst * crystallite_subStep(co,ip,el)
+            constitutive_mech_Fp(ph)%data(1:3,1:3,me) =            crystallite_subFp0(1:3,1:3,co,ip,el)
+            constitutive_mech_Fi(ph)%data(1:3,1:3,me) =            crystallite_subFi0(1:3,1:3,co,ip,el)
+            crystallite_S    (1:3,1:3,co,ip,el) =            crystallite_S0    (1:3,1:3,co,ip,el)
+            if (crystallite_subStep(co,ip,el) < 1.0_pReal) then                                        ! actual (not initial) cutback
+              crystallite_Lp (1:3,1:3,co,ip,el) =            subLp0
+              constitutive_mech_Li(ph)%data(1:3,1:3,me) =            subLi0
+            endif
+            plasticState    (material_phaseAt(co,el))%state(    :,material_phaseMemberAt(co,ip,el)) &
+              = plasticState(material_phaseAt(co,el))%subState0(:,material_phaseMemberAt(co,ip,el))
+            do s = 1, phase_Nsources(material_phaseAt(co,el))
+              sourceState(    material_phaseAt(co,el))%p(s)%state(    :,material_phaseMemberAt(co,ip,el)) &
+                = sourceState(material_phaseAt(co,el))%p(s)%subState0(:,material_phaseMemberAt(co,ip,el))
+            enddo
+
+                                                                                                    ! cant restore dotState here, since not yet calculated in first cutback after initialization
+            todo = crystallite_subStep(co,ip,el) > num%subStepMinCryst                          ! still on track or already done (beyond repair)
+          endif
+
+!--------------------------------------------------------------------------------------------------
+!  prepare for integration
+          if (todo) then
+            crystallite_subF(1:3,1:3,co,ip,el) = crystallite_subF0(1:3,1:3,co,ip,el) &
+                                            + crystallite_subStep(co,ip,el) *( crystallite_partitionedF (1:3,1:3,co,ip,el) &
+                                                                           -crystallite_partitionedF0(1:3,1:3,co,ip,el))
+            crystallite_Fe(1:3,1:3,co,ip,el) = matmul(crystallite_subF(1:3,1:3,co,ip,el), &
+                                                   math_inv33(matmul(constitutive_mech_Fi(ph)%data(1:3,1:3,me), &
+                                                                     constitutive_mech_Fp(ph)%data(1:3,1:3,me))))
+            crystallite_subdt(co,ip,el) = crystallite_subStep(co,ip,el) * crystallite_dt(co,ip,el)
+            crystallite_converged(co,ip,el) = .false.
+            call integrateState(co,ip,el)
+            call integrateSourceState(co,ip,el)
+          endif
+
+
+
+!--------------------------------------------------------------------------------------------------
+!  integrate --- requires fully defined state array (basic + dependent state)
+    if (.not. crystallite_converged(co,ip,el) .and. crystallite_subStep(co,ip,el) > num%subStepMinCryst) &            ! do not try non-converged but fully cutbacked any further
+      todo = .true.
+  enddo cutbackLooping
+
+! return whether converged or not
+  crystallite_stress2 = crystallite_converged(co,ip,el)
+
+end function crystallite_stress2
+
+
 !--------------------------------------------------------------------------------------------------
 !> @brief Backup data for homog cutback.
 !--------------------------------------------------------------------------------------------------

src/homogenization.f90

@@ -180,7 +180,7 @@ subroutine materialpoint_stressAndItsTangent(dt)
     NiterationMPstate, &
     i, &                                                                                            !< integration point number
     e, &                                                                                            !< element number
-    myNgrains
+    myNgrains, co
   real(pReal), dimension(discretization_nIPs,discretization_Nelems) :: &
     subFrac, &
     subStep
@@ -285,7 +285,7 @@ subroutine materialpoint_stressAndItsTangent(dt)
 
 !--------------------------------------------------------------------------------------------------
 ! deformation partitioning
-      !$OMP PARALLEL DO PRIVATE(myNgrains,m)
+      !$OMP PARALLEL DO PRIVATE(myNgrains,m,co)
       elementLooping2: do e = FEsolving_execElem(1),FEsolving_execElem(2)
         myNgrains = homogenization_Nconstituents(material_homogenizationAt(e))
         IpLooping2: do i = FEsolving_execIP(1),FEsolving_execIP(2)
@@ -300,19 +300,12 @@ subroutine materialpoint_stressAndItsTangent(dt)
           else
             crystallite_requested(1:myNgrains,i,e) = .false.                                        ! calculation for constituents not required anymore
           endif
-        enddo IpLooping2
+          converged(i,e) = .true.
-      enddo elementLooping2
+          do co = 1, myNgrains
-      !$OMP END PARALLEL DO
+            converged(i,e) = converged(i,e) .and. crystallite_stress2(co,i,e)
+          enddo
 
-!--------------------------------------------------------------------------------------------------
-! crystallite integration
-      converged = crystallite_stress() !ToDo: MD not sure if that is the best logic
 
-!--------------------------------------------------------------------------------------------------
-! state update
-     !$OMP PARALLEL DO PRIVATE(m)
-      elementLooping3: do e = FEsolving_execElem(1),FEsolving_execElem(2)
-        IpLooping3: do i = FEsolving_execIP(1),FEsolving_execIP(2)
           if (requested(i,e) .and. .not. doneAndHappy(1,i,e)) then
             if (.not. converged(i,e)) then
               doneAndHappy(1:2,i,e) = [.true.,.false.]
@@ -326,8 +319,8 @@ subroutine materialpoint_stressAndItsTangent(dt)
               converged(i,e) = all(doneAndHappy(1:2,i,e))                                           ! converged if done and happy
             endif
           endif
-        enddo IpLooping3
+        enddo IpLooping2
-      enddo elementLooping3
+      enddo elementLooping2
       !$OMP END PARALLEL DO
 
     enddo convergenceLooping
@@ -339,11 +332,11 @@ subroutine materialpoint_stressAndItsTangent(dt)
   if (.not. terminallyIll ) then
     call crystallite_orientations()                                                                 ! calculate crystal orientations
     !$OMP PARALLEL DO
-    elementLooping4: do e = FEsolving_execElem(1),FEsolving_execElem(2)
+    elementLooping3: do e = FEsolving_execElem(1),FEsolving_execElem(2)
-      IpLooping4: do i = FEsolving_execIP(1),FEsolving_execIP(2)
+      IpLooping3: do i = FEsolving_execIP(1),FEsolving_execIP(2)
         call mech_homogenize(i,e)
-      enddo IpLooping4
+      enddo IpLooping3
-    enddo elementLooping4
+    enddo elementLooping3
     !$OMP END PARALLEL DO
   else
     print'(/,a,/)', ' << HOMOG >> Material Point terminally ill'