nonlocal stiffness calculation:

rather perturb all components at once (and optionally decrease the frequency of the Jacobian update with the iJaco parameter) than perturbing only a single component per cycle

code/crystallite.f90

@@ -925,40 +925,9 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
       
     elseif (any(.not. crystallite_localConstitution)) then                                                                          ! if any nonlocal grain present, we have to do a full loop over all grains after each perturbance
       
-      do e = FEsolving_execElem(1),FEsolving_execElem(2)
+      do k = 1,3
-        myNgrains = homogenization_Ngrains(mesh_element(3,e))
+        do l = 1,3
-        do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e)
+          crystallite_subF(k,l,:,:,:) = crystallite_subF(k,l,:,:,:) + pert_Fg                                                       ! perturb single component
-          do g = 1,myNgrains
-            selectiveDebugger = (e == debug_e .and. i == debug_i .and. g == debug_g)
-            mask = .true.
-            do comp = 1,9
-              kl(:,comp,g,i,e) = maxloc(abs(crystallite_subF(:,:,g,i,e)-crystallite_F0(:,:,g,i,e)), mask)                           ! map from component to array indices for F (sorted in descending order of abs(deltaF))
-              mask(kl(1,comp,g,i,e),kl(2,comp,g,i,e)) = .false.
-            enddo
-            k = kl(1,mod(cycleCounter/iJacoStiffness,9)+1,g,i,e)                                                                                     ! perturb components in the descending order of change in F (-> component with biggest change in F is perturbed in first cycle, component with second biggest change in next cycle, ...)
-            l = kl(2,mod(cycleCounter/iJacoStiffness,9)+1,g,i,e)      
-            if (verboseDebugger .and. selectiveDebugger) then
-              !$OMP CRITICAL (write2out)
-                write (6,*) 'perturb component ',k,l
-              !$OMPEND CRITICAL (write2out)
-            endif
-            crystallite_subF(k,l,g,i,e) = crystallite_subF(k,l,g,i,e) + pert_Fg                                                     ! perturb single component
-      enddo; enddo; enddo
-      
-      !$OMP PARALLEL DO
-        do e = FEsolving_execElem(1),FEsolving_execElem(2) 
-          myNgrains = homogenization_Ngrains(mesh_element(3,e))
-          do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e)
-            do g = 1,myNgrains
-              if (crystallite_todo(g,i,e)) then
-                crystallite_stateConverged(g,i,e) = crystallite_updateState(g,i,e)                                                  ! preguess for state
-                crystallite_temperatureConverged(g,i,e) = crystallite_updateTemperature(g,i,e)
-                if (      .not. crystallite_localConstitution(g,i,e) &
-                    .and. .not. crystallite_todo(g,i,e)) &                                                                          ! if broken non-local... 
-                  crystallite_todo = crystallite_todo .and. crystallite_localConstitution                                           ! ...all non-locals skipped
-              endif
-        enddo; enddo; enddo
-      !$OMPEND PARALLEL DO
                 
           NiterationState = 0_pInt
           crystallite_todo = .true.
@@ -970,6 +939,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
                 myNgrains = homogenization_Ngrains(mesh_element(3,e))
                 do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e)
                   do g = 1,myNgrains
+                    selectiveDebugger = (e == debug_e .and. i == debug_i .and. g == debug_g)
                     if (crystallite_todo(g,i,e)) then
                       crystallite_todo(g,i,e) = crystallite_integrateStress(g,i,e)                                                      ! stress integration
                       if (      .not. crystallite_localConstitution(g,i,e) & 
@@ -994,6 +964,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
                 myNgrains = homogenization_Ngrains(mesh_element(3,e))
                 do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e)
                   do g = 1,myNgrains
+                    selectiveDebugger = (e == debug_e .and. i == debug_i .and. g == debug_g)
                     if (crystallite_todo(g,i,e)) then
                       call constitutive_collectDotState(crystallite_Tstar_v(:,g,i,e), crystallite_subTstar0_v(:,g,i,e), &
                                                         crystallite_Fe, crystallite_Fp, crystallite_Temperature(g,i,e), & 
@@ -1005,8 +976,10 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
                       dot_prod22 = dot_product(delta_dotState2, delta_dotState2)
                       if (      dot_prod22 > 0.0_pReal &
                           .and. (     dot_prod12 < 0.0_pReal &
-                             .or. dot_product(constitutive_dotState(g,i,e)%p, constitutive_previousDotState(g,i,e)%p) < 0.0_pReal))&
+                                 .or. dot_product(constitutive_dotState(g,i,e)%p, &
-                    crystallite_statedamper(g,i,e) = 0.75_pReal + 0.25_pReal * tanh(2.0_pReal + 4.0_pReal * dot_prod12 / dot_prod22)
+                                                  constitutive_previousDotState(g,i,e)%p) < 0.0_pReal) ) &
+                        crystallite_statedamper(g,i,e) = 0.75_pReal &
+                                                       + 0.25_pReal * tanh(2.0_pReal + 4.0_pReal * dot_prod12 / dot_prod22)
                     endif
               enddo; enddo; enddo
             !$OMPEND PARALLEL DO
@@ -1016,6 +989,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
                 myNgrains = homogenization_Ngrains(mesh_element(3,e))
                 do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e)
                   do g = 1,myNgrains
+                    selectiveDebugger = (e == debug_e .and. i == debug_i .and. g == debug_g)
                     if (crystallite_todo(g,i,e)) then
                       crystallite_stateConverged(g,i,e) = crystallite_updateState(g,i,e)                                                ! update state
                       crystallite_temperatureConverged(g,i,e) = crystallite_updateTemperature(g,i,e)                                    ! update temperature
@@ -1040,8 +1014,6 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
             do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e)
               do g = 1,myNgrains
                 if (crystallite_converged(g,i,e)) then                                                                                  ! if stiffness calculation converged...
-              k = kl(1,mod(cycleCounter/iJacoStiffness,9)+1,g,i,e)
-              l = kl(2,mod(cycleCounter/iJacoStiffness,9)+1,g,i,e)      
                   crystallite_dPdF(:,:,k,l,g,i,e) = (crystallite_P(:,:,g,i,e) - storedP(:,:,g,i,e))/pert_Fg                             ! ... use tangent dP_ij/dFg_kl
                 elseif (.not. storedConvergenceFlag(g,i,e)) then                                                                        ! if crystallite didnÕt converge before...
                   crystallite_dPdF(:,:,:,:,g,i,e) = crystallite_fallbackdPdF(:,:,:,:,g,i,e)                                             ! ... use (elastic) fallback
@@ -1072,6 +1044,8 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
     
           crystallite_converged = storedConvergenceFlag
       
+      enddo;enddo                 ! k,l loop
+      
     endif
 
   endif                           ! jacobian calculation