Merge branch 'development' into misc-improvements

CMakeLists.txt

@@ -148,7 +148,7 @@ else ()
     set(OPENMP "${OPENMP}")
 endif ()
 
-# syntax check only (mainly for pre-receive hook, works only with gfortran)
+# syntax check only (mainly for pre-receive hook)
 if (CMAKE_BUILD_TYPE STREQUAL "SYNTAXONLY")
     set (BUILDCMD_POST "${BUILDCMD_POST} -fsyntax-only")
 endif ()

VERSION

@@ -1 +1 @@
-v2.0.3-2204-gdb1f2151
+v2.0.3-2243-gbb03483b

processing/pre/geom_fromVoronoiTessellation.py

@@ -21,36 +21,42 @@ def findClosestSeed(seeds, weights, point):
     return np.argmin(np.sum((np.broadcast_to(point,(len(seeds),3))-seeds)**2,axis=1) - weights)
 
 
-def Laguerre_tessellation(grid, seeds, grains, size, periodic, weights, cpus):
+def Laguerre_tessellation(grid, size, seeds, weights, origin = np.zeros(3), periodic = True, cpus = 2):
 
     if periodic:
         weights_p = np.tile(weights,27).flatten(order='F')                                          # Laguerre weights (1,2,3,1,2,3,...,1,2,3)
         seeds_p = np.vstack((seeds  -np.array([size[0],0.,0.]),seeds,  seeds  +np.array([size[0],0.,0.])))
         seeds_p = np.vstack((seeds_p-np.array([0.,size[1],0.]),seeds_p,seeds_p+np.array([0.,size[1],0.])))
         seeds_p = np.vstack((seeds_p-np.array([0.,0.,size[2]]),seeds_p,seeds_p+np.array([0.,0.,size[2]])))
+        coords = damask.grid_filters.cell_coord0(grid*3,size*3,-origin-size).reshape(-1,3,order='F')
     else:
         weights_p = weights.flatten()
         seeds_p   = seeds
+        coords = damask.grid_filters.cell_coord0(grid,size,-origin).reshape(-1,3,order='F')
 
     if cpus > 1:
-        default_args = partial(findClosestSeed,seeds_p,weights_p)
+        pool = multiprocessing.Pool(processes = cpus)
-        pool = multiprocessing.Pool(processes = cpus)                                               # initialize workers
+        result = pool.map_async(partial(findClosestSeed,seeds_p,weights_p), [coord for coord in coords])
-        result = pool.map_async(default_args, [point for point in grid])                            # evaluate function in parallel
         pool.close()
         pool.join()
-        closestSeeds = np.array(result.get()).flatten()
+        closest_seed = np.array(result.get())
     else:
-        closestSeeds= np.array([findClosestSeed(seeds_p,weights_p,point) for point in grid])
+        closest_seed= np.array([findClosestSeed(seeds_p,weights_p,coord) for coord in coords])
 
-    return grains[closestSeeds%seeds.shape[0]]
+    if periodic:
+        closest_seed = closest_seed.reshape(grid*3)
+        return closest_seed[grid[0]:grid[0]*2,grid[1]:grid[1]*2,grid[2]:grid[2]*2]%seeds.shape[0]
+    else:
+        return closest_seed
 
 
-def Voronoi_tessellation(grid, seeds, grains, size, periodic = True):
+def Voronoi_tessellation(grid, size, seeds, origin = np.zeros(3), periodic = True):
 
+    coords = damask.grid_filters.cell_coord0(grid,size,-origin).reshape(-1,3,order='F')
     KDTree = spatial.cKDTree(seeds,boxsize=size) if periodic else spatial.cKDTree(seeds)
-    devNull,closestSeeds = KDTree.query(grid)
+    devNull,closest_seed = KDTree.query(coords)
 
-    return grains[closestSeeds]
+    return closest_seed
 
 
 # --------------------------------------------------------------------
@@ -186,18 +192,15 @@ for name in filenames:
 
     grains = table.get(options.microstructure) if options.microstructure in table.labels else np.arange(len(seeds))+1
     grainIDs  = np.unique(grains).astype('i')
-    NgrainIDs = len(grainIDs)
 
     if options.eulers in table.labels:
         eulers = table.get(options.eulers)
 
-    coords = damask.grid_filters.cell_coord0(grid,size,-origin).reshape(-1,3,order='F')
-
     if options.laguerre:
-        indices = Laguerre_tessellation(coords,seeds,grains,size,options.periodic,
+        indices = grains[Laguerre_tessellation(grid,size,seeds,table.get(options.weight),origin,
-                                        table.get(options.weight),options.cpus)
+                                               options.periodic,options.cpus)]
     else:
-        indices = Voronoi_tessellation (coords,seeds,grains,size,options.periodic)
+        indices = grains[Voronoi_tessellation (grid,size,seeds,origin,options.periodic)]
 
     config_header = []
     if options.config:

processing/pre/geom_translate.py

@@ -52,7 +52,7 @@ for name in filenames:
   geom = damask.Geom.from_file(StringIO(''.join(sys.stdin.read())) if name is None else name)
 
   substituted = geom.get_microstructure()
-  for old,new in zip(sub[0::2],sub[1::2]): substituted[substituted==old] = new                      # substitute microstructure indices
+  for old,new in zip(sub[0::2],sub[1::2]): substituted[geom.microstructure==old] = new              # substitute microstructure indices
   substituted += options.microstructure                                                             # constant shift
 
   damask.util.croak(geom.update(substituted,origin=geom.get_origin()+options.origin))

python/.coveragerc

@@ -0,0 +1,2 @@
+[run]
+omit = tests/*

python/.gitignore

@@ -1,3 +1,4 @@
 *.pyc
 dist
 damask.egg-info
+.coverage

python/damask/_geom.py

@@ -1,11 +1,15 @@
 import sys
 from io import StringIO
+import multiprocessing
+from functools import partial
 
 import numpy as np
-from scipy import ndimage
+from scipy import ndimage,spatial
 
 from . import VTK
 from . import util
+from . import Environment
+from . import grid_filters
 
 
 class Geom:
@@ -50,7 +54,7 @@ class Geom:
 
     def update(self,microstructure=None,size=None,origin=None,rescale=False):
         """
-        Updates microstructure and size.
+        Update microstructure and size.
 
         Parameters
         ----------
@@ -113,7 +117,7 @@ class Geom:
 
     def set_comments(self,comments):
         """
-        Replaces all existing comments.
+        Replace all existing comments.
 
         Parameters
         ----------
@@ -127,7 +131,7 @@ class Geom:
 
     def add_comments(self,comments):
         """
-        Appends comments to existing comments.
+        Append comments to existing comments.
 
         Parameters
         ----------
@@ -140,7 +144,7 @@ class Geom:
 
     def set_microstructure(self,microstructure):
         """
-        Replaces the existing microstructure representation.
+        Replace the existing microstructure representation.
 
         Parameters
         ----------
@@ -159,7 +163,7 @@ class Geom:
 
     def set_size(self,size):
         """
-        Replaces the existing size information.
+        Replace the existing size information.
 
         Parameters
         ----------
@@ -179,7 +183,7 @@ class Geom:
 
     def set_origin(self,origin):
         """
-        Replaces the existing origin information.
+        Replace the existing origin information.
 
         Parameters
         ----------
@@ -196,7 +200,7 @@ class Geom:
 
     def set_homogenization(self,homogenization):
         """
-        Replaces the existing homogenization index.
+        Replace the existing homogenization index.
 
         Parameters
         ----------
@@ -264,7 +268,7 @@ class Geom:
     @staticmethod
     def from_file(fname):
         """
-        Reads a geom file.
+        Read a geom file.
 
         Parameters
         ----------
@@ -325,6 +329,81 @@ class Geom:
         return Geom(microstructure.reshape(grid),size,origin,homogenization,comments)
 
 
+    @staticmethod
+    def _find_closest_seed(seeds, weights, point):
+        return np.argmin(np.sum((np.broadcast_to(point,(len(seeds),3))-seeds)**2,axis=1) - weights)
+    @staticmethod
+    def from_Laguerre_tessellation(grid,size,seeds,weights,periodic=True):
+        """
+        Generate geometry from Laguerre tessellation.
+
+        Parameters
+        ----------
+        grid : numpy.ndarray of shape (3)
+            number of grid points in x,y,z direction.
+        size : list or numpy.ndarray of shape (3)
+            physical size of the microstructure in meter.
+        seeds : numpy.ndarray of shape (:,3)
+            position of the seed points in meter. All points need to lay within the box.
+        weights : numpy.ndarray of shape (seeds.shape[0])
+            weights of the seeds. Setting all weights to 1.0 gives a standard Voronoi tessellation.
+        periodic : Boolean, optional
+            perform a periodic tessellation. Defaults to True.
+
+        """
+        if periodic:
+            weights_p = np.tile(weights,27).flatten(order='F')                                      # Laguerre weights (1,2,3,1,2,3,...,1,2,3)
+            seeds_p = np.vstack((seeds  -np.array([size[0],0.,0.]),seeds,  seeds  +np.array([size[0],0.,0.])))
+            seeds_p = np.vstack((seeds_p-np.array([0.,size[1],0.]),seeds_p,seeds_p+np.array([0.,size[1],0.])))
+            seeds_p = np.vstack((seeds_p-np.array([0.,0.,size[2]]),seeds_p,seeds_p+np.array([0.,0.,size[2]])))
+            coords  = grid_filters.cell_coord0(grid*3,size*3,-size).reshape(-1,3,order='F')
+
+        else:
+            weights_p = weights.flatten()
+            seeds_p   = seeds
+            coords    = grid_filters.cell_coord0(grid,size).reshape(-1,3,order='F')
+
+        pool = multiprocessing.Pool(processes = int(Environment().options['DAMASK_NUM_THREADS']))
+        result = pool.map_async(partial(Geom._find_closest_seed,seeds_p,weights_p), [coord for coord in coords])
+        pool.close()
+        pool.join()
+        microstructure = np.array(result.get())
+
+        if periodic:
+            microstructure = microstructure.reshape(grid*3)
+            microstructure = microstructure[grid[0]:grid[0]*2,grid[1]:grid[1]*2,grid[2]:grid[2]*2]%seeds.shape[0]
+        else:
+            microstructure = microstructure.reshape(grid)
+
+        #comments = 'geom.py:from_Laguerre_tessellation v{}'.format(version)
+        return Geom(microstructure+1,size,homogenization=1)
+
+
+    @staticmethod
+    def from_Voronoi_tessellation(grid,size,seeds,periodic=True):
+        """
+        Generate geometry from Voronoi tessellation.
+
+        Parameters
+        ----------
+        grid : numpy.ndarray of shape (3)
+            number of grid points in x,y,z direction.
+        size : list or numpy.ndarray of shape (3)
+            physical size of the microstructure in meter.
+        seeds : numpy.ndarray of shape (:,3)
+            position of the seed points in meter. All points need to lay within the box.
+        periodic : Boolean, optional
+            perform a periodic tessellation. Defaults to True.
+
+        """
+        coords = grid_filters.cell_coord0(grid,size).reshape(-1,3,order='F')
+        KDTree = spatial.cKDTree(seeds,boxsize=size) if periodic else spatial.cKDTree(seeds)
+        devNull,microstructure = KDTree.query(coords)
+
+        #comments = 'geom.py:from_Voronoi_tessellation v{}'.format(version)
+        return Geom(microstructure.reshape(grid)+1,size,homogenization=1)
+
+
     def to_file(self,fname,pack=None):
         """
         Writes a geom file.

python/damask/_vtk.py

@@ -72,7 +72,7 @@ class VTK:
         connectivity : numpy.ndarray of np.dtype = int
             Cell connectivity (0-based), first dimension determines #Cells, second dimension determines #Nodes/Cell.
         cell_type : str
-            Name of the vtk.vtkCell subclass. Tested for TRIANGLE, QUAD, and HEXAHEDRON.
+            Name of the vtk.vtkCell subclass. Tested for TRIANGLE, QUAD, TETRA, and HEXAHEDRON.
 
         """
         vtk_nodes = vtk.vtkPoints()

python/tests/test_Geom.py

@@ -72,7 +72,7 @@ class TestGeom:
         if update: modified.to_file(reference)
         assert geom_equal(modified,Geom.from_file(reference))
 
-    @pytest.mark.parametrize('stencil',[(1),(2),(3),(4)])
+    @pytest.mark.parametrize('stencil',[1,2,3,4])
     def test_clean(self,default,update,reference_dir,stencil):
         modified = copy.deepcopy(default)
         modified.clean(stencil)
@@ -82,12 +82,12 @@ class TestGeom:
         assert geom_equal(modified,Geom.from_file(reference))
 
     @pytest.mark.parametrize('grid',[
-                                     ((10,11,10)),
+                                     (10,11,10),
-                                     ([10,13,10]),
+                                     [10,13,10],
-                                     (np.array((10,10,10))),
+                                     np.array((10,10,10)),
-                                     (np.array((8, 10,12))),
+                                     np.array((8, 10,12)),
-                                     (np.array((5, 4, 20))),
+                                     np.array((5, 4, 20)),
-                                     (np.array((10,20,2)) )
+                                     np.array((10,20,2))
                                     ]
                             )
     def test_scale(self,default,update,reference_dir,grid):
@@ -97,3 +97,37 @@ class TestGeom:
         reference = os.path.join(reference_dir,'scale_{}.geom'.format(tag))
         if update: modified.to_file(reference)
         assert geom_equal(modified,Geom.from_file(reference))
+
+    @pytest.mark.parametrize('periodic',[True,False])
+    def test_tessellation_approaches(self,periodic):
+        grid   = np.random.randint(10,20,3)
+        size   = np.random.random(3) + 1.0
+        N_seeds= np.random.randint(10,30)
+        seeds  = np.random.rand(N_seeds,3) * np.broadcast_to(size,(N_seeds,3))
+        Voronoi  = Geom.from_Voronoi_tessellation( grid,size,seeds,                 periodic)
+        Laguerre = Geom.from_Laguerre_tessellation(grid,size,seeds,np.ones(N_seeds),periodic)
+        assert geom_equal(Laguerre,Voronoi)
+
+    def test_Laguerre_weights(self):
+        grid   = np.random.randint(10,20,3)
+        size   = np.random.random(3) + 1.0
+        N_seeds= np.random.randint(10,30)
+        seeds  = np.random.rand(N_seeds,3) * np.broadcast_to(size,(N_seeds,3))
+        weights= np.full((N_seeds),-np.inf)
+        ms     = np.random.randint(1, N_seeds+1)
+        weights[ms-1] = np.random.random()
+        Laguerre = Geom.from_Laguerre_tessellation(grid,size,seeds,weights,np.random.random()>0.5)
+        assert np.all(Laguerre.microstructure == ms)
+
+    @pytest.mark.parametrize('approach',['Laguerre','Voronoi'])
+    def test_tessellate_bicrystal(self,approach):
+        grid  = np.random.randint(5,10,3)*2
+        size  = grid.astype(np.float)
+        seeds = np.vstack((size*np.array([0.5,0.25,0.5]),size*np.array([0.5,0.75,0.5])))
+        microstructure = np.ones(grid)
+        microstructure[:,grid[1]//2:,:] = 2
+        if   approach == 'Laguerre':
+            geom = Geom.from_Laguerre_tessellation(grid,size,seeds,np.ones(2),np.random.random()>0.5)
+        elif approach == 'Voronoi':
+            geom = Geom.from_Voronoi_tessellation(grid,size,seeds,            np.random.random()>0.5)
+        assert np.all(geom.microstructure == microstructure)

python/tests/test_VTK.py

@@ -0,0 +1,47 @@
+import os
+
+import pytest
+import numpy as np
+
+from damask import VTK
+
+@pytest.fixture
+def reference_dir(reference_dir_base):
+    """Directory containing reference results."""
+    return os.path.join(reference_dir_base,'Result')
+
+class TestVTK:
+
+    def test_rectilinearGrid(self,tmp_path):
+        grid   = np.random.randint(5,10,3)*2
+        size   = np.random.random(3) + 1.0
+        origin = np.random.random(3)
+        v = VTK.from_rectilinearGrid(grid,size,origin)
+        s = v.__repr__()
+        v.write(os.path.join(tmp_path,'rectilinearGrid'))
+        v = VTK.from_file(os.path.join(tmp_path,'rectilinearGrid.vtr'))
+        assert(s == v.__repr__())
+
+    def test_polyData(self,tmp_path):
+        points = np.random.rand(3,100)
+        v = VTK.from_polyData(points)
+        s = v.__repr__()
+        v.write(os.path.join(tmp_path,'polyData'))
+        v = VTK.from_file(os.path.join(tmp_path,'polyData.vtp'))
+        assert(s == v.__repr__())
+
+    @pytest.mark.parametrize('cell_type,n',[
+                                            ('VTK_hexahedron',8),
+                                            ('TETRA',4),
+                                            ('quad',4),
+                                            ('VTK_TRIANGLE',3)
+                                            ]
+                            )
+    def test_unstructuredGrid(self,tmp_path,cell_type,n):
+        nodes = np.random.rand(n,3)
+        connectivity = np.random.choice(np.arange(n),n,False).reshape(-1,n)
+        v = VTK.from_unstructuredGrid(nodes,connectivity,cell_type)
+        s = v.__repr__()
+        v.write(os.path.join(tmp_path,'unstructuredGrid'))
+        v = VTK.from_file(os.path.join(tmp_path,'unstructuredGrid.vtu'))
+        assert(s == v.__repr__())

src/DAMASK_interface.f90

@@ -25,7 +25,7 @@ module DAMASK_interface
   
   implicit none
   private
-  logical,                       public, protected :: &
+  logical,          volatile,    public, protected :: &
     SIGTERM, &                                                                                      !< termination signal 
     SIGUSR1, &                                                                                      !< 1. user-defined signal
     SIGUSR2                                                                                         !< 2. user-defined signal

src/DAMASK_marc.f90

@@ -43,6 +43,9 @@ module DAMASK_interface
   logical,          protected, public :: symmetricSolver
   character(len=*), parameter, public :: INPUTFILEEXTENSION = '.dat'
 
+  logical, dimension(:,:), public, allocatable :: &
+    calcMode                                                                                        !< calculate or collect (ping pong scheme)
+
   public :: &
     DAMASK_interface_init, &
     getSolverJobName
@@ -102,7 +105,6 @@ function getSolverJobName()
   character(len=*), parameter   :: pathSep = achar(47)//achar(92)                                   ! forward and backward slash
   integer :: extPos
 
-  inputName=''
   inquire(5, name=inputName)                                                                        ! determine inputfile
   extPos = len_trim(inputName)-4
   getSolverJobName=inputName(scan(inputName,pathSep,back=.true.)+1:extPos)
@@ -179,6 +181,7 @@ subroutine hypela2(d,g,e,de,s,t,dt,ngens,m,nn,kcus,matus,ndi,nshear,disp, &
  use FEsolving
  use debug
  use discretization_marc
+ use homogenization
  use CPFEM
 
  implicit none

src/FEsolving.f90

@@ -8,16 +8,8 @@ module FEsolving
   implicit none
   public
 
-  logical :: &
-    terminallyIll = .false.                                                                         !< at least one material point is terminally ill
-
   integer, dimension(2) :: &
     FEsolving_execElem, &                                                                           !< for ping-pong scheme always whole range, otherwise one specific element
     FEsolving_execIP                                                                                !< for ping-pong scheme always range to max IP, otherwise one specific IP
 
-#if defined(Marc4DAMASK)
-  logical, dimension(:,:), allocatable :: &
-    calcMode                                                                                        !< do calculation or simply collect when using ping pong scheme
-#endif
-
 end module FEsolving

src/constitutive_plastic_dislotwin.f90

@@ -568,7 +568,22 @@ module subroutine plastic_dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,instance,of)
                        + ddot_gamma_dtau_slip(i) * prm%P_sl(k,l,i) * prm%P_sl(m,n,i)
   enddo slipContribution
 
-  !ToDo: Why do this before shear banding?
+  call kinetics_twin(Mp,T,dot_gamma_sl,instance,of,dot_gamma_twin,ddot_gamma_dtau_twin)
+  twinContibution: do i = 1, prm%sum_N_tw
+    Lp = Lp + dot_gamma_twin(i)*prm%P_tw(1:3,1:3,i)
+    forall (k=1:3,l=1:3,m=1:3,n=1:3) &
+      dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) &
+                       + ddot_gamma_dtau_twin(i)* prm%P_tw(k,l,i)*prm%P_tw(m,n,i)
+  enddo twinContibution
+
+  call kinetics_trans(Mp,T,dot_gamma_sl,instance,of,dot_gamma_tr,ddot_gamma_dtau_trans)
+  transContibution: do i = 1, prm%sum_N_tr
+    Lp = Lp + dot_gamma_tr(i)*prm%P_tr(1:3,1:3,i)
+    forall (k=1:3,l=1:3,m=1:3,n=1:3) &
+      dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) &
+                       + ddot_gamma_dtau_trans(i)* prm%P_tr(k,l,i)*prm%P_tr(m,n,i)
+  enddo transContibution
+
   Lp      = Lp      * f_unrotated
   dLp_dMp = dLp_dMp * f_unrotated
 
@@ -598,23 +613,6 @@ module subroutine plastic_dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,instance,of)
 
   endif shearBandingContribution
 
-  call kinetics_twin(Mp,T,dot_gamma_sl,instance,of,dot_gamma_twin,ddot_gamma_dtau_twin)
-  twinContibution: do i = 1, prm%sum_N_tw
-    Lp = Lp + dot_gamma_twin(i)*prm%P_tw(1:3,1:3,i) * f_unrotated
-    forall (k=1:3,l=1:3,m=1:3,n=1:3) &
-      dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) &
-                       + ddot_gamma_dtau_twin(i)* prm%P_tw(k,l,i)*prm%P_tw(m,n,i) * f_unrotated
-  enddo twinContibution
-
-  call kinetics_trans(Mp,T,dot_gamma_sl,instance,of,dot_gamma_tr,ddot_gamma_dtau_trans)
-  transContibution: do i = 1, prm%sum_N_tr
-    Lp = Lp + dot_gamma_tr(i)*prm%P_tr(1:3,1:3,i) * f_unrotated
-    forall (k=1:3,l=1:3,m=1:3,n=1:3) &
-      dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) &
-                       + ddot_gamma_dtau_trans(i)* prm%P_tr(k,l,i)*prm%P_tr(m,n,i) * f_unrotated
-  enddo transContibution
-
-
   end associate
 
 end subroutine plastic_dislotwin_LpAndItsTangent

src/constitutive_plastic_nonlocal.f90

@@ -16,15 +16,15 @@ submodule(constitutive) plastic_nonlocal
     kB = 1.38e-23_pReal                                                                             !< Boltzmann constant in J/Kelvin
 
   ! storage order of dislocation types
-  integer, dimension(8), parameter :: &
+  integer, dimension(*), parameter :: &
     sgl = [1,2,3,4,5,6,7,8]                                                                         !< signed (single)
-  integer, dimension(5), parameter :: &
+  integer, dimension(*), parameter :: &
     edg = [1,2,5,6,9], &                                                                            !< edge
     scr = [3,4,7,8,10]                                                                              !< screw
-  integer, dimension(4), parameter :: &
+  integer, dimension(*), parameter :: &
     mob = [1,2,3,4], &                                                                              !< mobile
     imm = [5,6,7,8]                                                                                 !< immobile (blocked)
-  integer, dimension(2), parameter :: &
+  integer, dimension(*), parameter :: &
     dip = [9,10], &                                                                                 !< dipole
     imm_edg = imm(1:2), &                                                                           !< immobile edge
     imm_scr = imm(3:4)                                                                              !< immobile screw
@@ -1611,7 +1611,7 @@ end subroutine stateInit
 !--------------------------------------------------------------------------------------------------
 !> @brief calculates kinetics
 !--------------------------------------------------------------------------------------------------
-subroutine kinetics(v, dv_dtau, dv_dtauNS, tau, tauNS, tauThreshold, c, Temperature, instance)
+pure subroutine kinetics(v, dv_dtau, dv_dtauNS, tau, tauNS, tauThreshold, c, Temperature, instance)
 
   integer, intent(in) :: &
     c, &                                                                                            !< dislocation character (1:edge, 2:screw)
@@ -1726,7 +1726,7 @@ end subroutine kinetics
 !> @brief returns copy of current dislocation densities from state
 !> @details raw values is rectified
 !--------------------------------------------------------------------------------------------------
-function getRho(instance,of,ip,el)
+pure function getRho(instance,of,ip,el)
 
   integer, intent(in) :: instance, of,ip,el
   real(pReal), dimension(param(instance)%sum_N_sl,10) :: getRho
@@ -1751,7 +1751,7 @@ end function getRho
 !> @brief returns copy of current dislocation densities from state
 !> @details raw values is rectified
 !--------------------------------------------------------------------------------------------------
-function getRho0(instance,of,ip,el)
+pure function getRho0(instance,of,ip,el)
 
   integer, intent(in) :: instance, of,ip,el
   real(pReal), dimension(param(instance)%sum_N_sl,10) :: getRho0

src/crystallite.f90

@@ -804,7 +804,6 @@ logical function integrateStress(ipc,ip,el,timeFraction)
                                       residuumLi_old, &                                             ! last residuum of intermediate velocity gradient
                                       deltaLi, &                                                    ! direction of next guess
                                       Fe, &                                                         ! elastic deformation gradient
-                                      Fe_new, &
                                       S, &                                                          ! 2nd Piola-Kirchhoff Stress in plastic (lattice) configuration
                                       A, &
                                       B, &
@@ -911,21 +910,19 @@ logical function integrateStress(ipc,ip,el,timeFraction)
         cycle LpLoop
       endif
 
-     !* calculate Jacobian for correction term
+      calculateJacobiLi: if (mod(jacoCounterLp, num%iJacoLpresiduum) == 0) then
-      if (mod(jacoCounterLp, num%iJacoLpresiduum) == 0) then
         jacoCounterLp = jacoCounterLp + 1
 
         do o=1,3; do p=1,3
-          dFe_dLp(o,1:3,p,1:3) = A(o,p)*transpose(invFi_new)                                        ! dFe_dLp(i,j,k,l) = -dt * A(i,k) invFi(l,j)
+          dFe_dLp(o,1:3,p,1:3) = - dt * A(o,p)*transpose(invFi_new)                                 ! dFe_dLp(i,j,k,l) = -dt * A(i,k) invFi(l,j)
         enddo; enddo
-        dFe_dLp  = - dt * dFe_dLp
         dRLp_dLp = math_identity2nd(9) &
                  - math_3333to99(math_mul3333xx3333(math_mul3333xx3333(dLp_dS,dS_dFe),dFe_dLp))
         temp_9 = math_33to9(residuumLp)
         call dgesv(9,1,dRLp_dLp,9,devNull_9,temp_9,9,ierr)                                          ! solve dRLp/dLp * delta Lp = -res for delta Lp
         if (ierr /= 0) return ! error
         deltaLp = - math_9to33(temp_9)
-      endif
+      endif calculateJacobiLi
 
       Lpguess = Lpguess &
               + deltaLp * steplengthLp
@@ -953,8 +950,7 @@ logical function integrateStress(ipc,ip,el,timeFraction)
       cycle LiLoop
     endif
 
-    !* calculate Jacobian for correction term
+    calculateJacobiLp: if (mod(jacoCounterLi, num%iJacoLpresiduum) == 0) then
-    if (mod(jacoCounterLi, num%iJacoLpresiduum) == 0) then
       jacoCounterLi = jacoCounterLi + 1
 
       temp_33 = matmul(matmul(A,B),invFi_current)
@@ -973,7 +969,7 @@ logical function integrateStress(ipc,ip,el,timeFraction)
       call dgesv(9,1,dRLi_dLi,9,devNull_9,temp_9,9,ierr)                                            ! solve dRLi/dLp * delta Li = -res for delta Li
       if (ierr /= 0) return ! error
       deltaLi = - math_9to33(temp_9)
-    endif
+    endif calculateJacobiLp
 
     Liguess = Liguess &
             + deltaLi * steplengthLi
@@ -982,17 +978,15 @@ logical function integrateStress(ipc,ip,el,timeFraction)
   invFp_new = matmul(invFp_current,B)
   call math_invert33(Fp_new,devNull,error,invFp_new)
   if (error) return ! error
-  Fp_new = Fp_new / math_det33(Fp_new)**(1.0_pReal/3.0_pReal)                                       ! regularize
-  Fe_new = matmul(matmul(F,invFp_new),invFi_new)
 
   integrateStress = .true.
   crystallite_P    (1:3,1:3,ipc,ip,el) = matmul(matmul(F,invFp_new),matmul(S,transpose(invFp_new)))
   crystallite_S    (1:3,1:3,ipc,ip,el) = S
   crystallite_Lp   (1:3,1:3,ipc,ip,el) = Lpguess
   crystallite_Li   (1:3,1:3,ipc,ip,el) = Liguess
-  crystallite_Fp   (1:3,1:3,ipc,ip,el) = Fp_new
+  crystallite_Fp   (1:3,1:3,ipc,ip,el) = Fp_new / math_det33(Fp_new)**(1.0_pReal/3.0_pReal)         ! regularize
   crystallite_Fi   (1:3,1:3,ipc,ip,el) = Fi_new
-  crystallite_Fe   (1:3,1:3,ipc,ip,el) = Fe_new
+  crystallite_Fe   (1:3,1:3,ipc,ip,el) = matmul(matmul(F,invFp_new),invFi_new)
 
 end function integrateStress
 
@@ -1014,15 +1008,15 @@ subroutine integrateStateFPI
     sizeDotState
   real(pReal) :: &
     zeta
-  real(pReal), dimension(constitutive_plasticity_maxSizeDotState) :: &
+  real(pReal), dimension(max(constitutive_plasticity_maxSizeDotState,constitutive_source_maxSizeDotState)) :: &
-    residuum_plastic                                                                                ! residuum for plastic state
+    r                                                                                               ! state residuum
-  real(pReal), dimension(constitutive_source_maxSizeDotState) :: &
+  real(pReal), dimension(:), allocatable :: plastic_dotState_p1, plastic_dotState_p2
-    residuum_source                                                                                 ! residuum for source state
+  real(pReal), dimension(constitutive_source_maxSizeDotState,2,maxval(phase_Nsources)) :: source_dotState
   logical :: &
     nonlocalBroken
 
   nonlocalBroken = .false.
-  !$OMP PARALLEL DO PRIVATE(sizeDotState,residuum_plastic,residuum_source,zeta,p,c)
+  !$OMP PARALLEL DO PRIVATE(sizeDotState,r,zeta,p,c,plastic_dotState_p1, plastic_dotState_p2,source_dotState)
   do e = FEsolving_execElem(1),FEsolving_execElem(2)
     do i = FEsolving_execIP(1),FEsolving_execIP(2)
       do g = 1,homogenization_Ngrains(material_homogenizationAt(e))
@@ -1047,24 +1041,23 @@ subroutine integrateStateFPI
           plasticState(p)%state(1:sizeDotState,c) = plasticState(p)%subState0(1:sizeDotState,c) &
                                                   + plasticState(p)%dotState (1:sizeDotState,c) &
                                                     * crystallite_subdt(g,i,e)
+          plastic_dotState_p2 = 0.0_pReal * plasticState(p)%dotState (1:sizeDotState,c)             ! ToDo can be done smarter/clearer
           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)
+            source_dotState(1:sizeDotState,2,s) = 0.0_pReal
           enddo
 
           iteration: do NiterationState = 1, num%nState
 
-            plasticState(p)%previousDotState2(:,c) = merge(plasticState(p)%previousDotState(:,c),&
+            if(nIterationState > 1) plastic_dotState_p2 = plastic_dotState_p1
-                                                           0.0_pReal,&
+            plastic_dotState_p1 = plasticState(p)%dotState(:,c)
-                                                           NiterationState > 1)
-            plasticState(p)%previousDotState (:,c) = plasticState(p)%dotState(:,c)
             do s = 1, phase_Nsources(p)
-              sourceState(p)%p(s)%previousDotState2(:,c) = merge(sourceState(p)%p(s)%previousDotState(:,c),&
+              sizeDotState = sourceState(p)%p(s)%sizeDotState
-                                                                 0.0_pReal, &
+              if(nIterationState > 1) source_dotState(1:sizeDotState,2,s) = source_dotState(1:sizeDotState,1,s)
-                                                                 NiterationState > 1)
+              source_dotState(1:sizeDotState,1,s) = sourceState(p)%p(s)%dotState(:,c)
-              sourceState(p)%p(s)%previousDotState (:,c) = sourceState(p)%p(s)%dotState(:,c)
             enddo
 
             call constitutive_dependentState(crystallite_partionedF(1:3,1:3,g,i,e), &
@@ -1072,8 +1065,6 @@ subroutine integrateStateFPI
                                              g, i, e)
 
             crystallite_todo(g,i,e) = integrateStress(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)) exit iteration
 
             call constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), &
@@ -1085,59 +1076,52 @@ subroutine integrateStateFPI
             do s = 1, phase_Nsources(p)
               crystallite_todo(g,i,e) = crystallite_todo(g,i,e) .and. all(.not. IEEE_is_NaN(sourceState(p)%p(s)%dotState(:,c)))
             enddo
-            if(.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) &
-              nonlocalBroken = .true.
             if(.not. crystallite_todo(g,i,e)) exit iteration
 
             sizeDotState = plasticState(p)%sizeDotState
-            zeta = damper(plasticState(p)%dotState         (:,c), &
+            zeta = damper(plasticState(p)%dotState(:,c),plastic_dotState_p1,plastic_dotState_p2)
-                          plasticState(p)%previousDotState (:,c), &
-                          plasticState(p)%previousDotState2(:,c))
             plasticState(p)%dotState(:,c) = plasticState(p)%dotState(:,c) * zeta &
-                                          + plasticState(p)%previousDotState(:,c) * (1.0_pReal - zeta)
+                                          + plastic_dotState_p1 * (1.0_pReal - zeta)
-            residuum_plastic(1:SizeDotState) = plasticState(p)%state    (1:sizeDotState,c) &
+            r(1:SizeDotState) = plasticState(p)%state    (1:sizeDotState,c) &
                               - plasticState(p)%subState0(1:sizeDotState,c)  &
-                                             - plasticState(p)%dotState (1:sizeDotState,c)  &
+                              - plasticState(p)%dotState (1:sizeDotState,c) * crystallite_subdt(g,i,e)
-                                               * crystallite_subdt(g,i,e)
             plasticState(p)%state(1:sizeDotState,c) = plasticState(p)%state(1:sizeDotState,c) &
-                                                    - residuum_plastic(1:sizeDotState)
+                                                    - r(1:sizeDotState)
-            crystallite_converged(g,i,e) = converged(residuum_plastic(1:sizeDotState), &
+            crystallite_converged(g,i,e) = converged(r(1:sizeDotState), &
                                                      plasticState(p)%state(1:sizeDotState,c), &
                                                      plasticState(p)%atol(1:sizeDotState))
             do s = 1, phase_Nsources(p)
               sizeDotState  = sourceState(p)%p(s)%sizeDotState
               zeta = damper(sourceState(p)%p(s)%dotState(:,c), &
-                            sourceState(p)%p(s)%previousDotState (:,c), &
+                            source_dotState(1:sizeDotState,1,s),&
-                            sourceState(p)%p(s)%previousDotState2(:,c))
+                            source_dotState(1:sizeDotState,2,s))
               sourceState(p)%p(s)%dotState(:,c) = sourceState(p)%p(s)%dotState(:,c) * zeta &
-                                                + sourceState(p)%p(s)%previousDotState(:,c)* (1.0_pReal - zeta)
+                                                + source_dotState(1:sizeDotState,1,s)* (1.0_pReal - zeta)
-              residuum_source(1:sizeDotState) = sourceState(p)%p(s)%state    (1:sizeDotState,c)  &
+              r(1:sizeDotState) = sourceState(p)%p(s)%state    (1:sizeDotState,c)  &
                                 - sourceState(p)%p(s)%subState0(1:sizeDotState,c)  &
-                                              - sourceState(p)%p(s)%dotState (1:sizeDotState,c)  &
+                                - sourceState(p)%p(s)%dotState (1:sizeDotState,c) * crystallite_subdt(g,i,e)
-                                                * crystallite_subdt(g,i,e)
               sourceState(p)%p(s)%state(1:sizeDotState,c) = sourceState(p)%p(s)%state(1:sizeDotState,c) &
-                                                          - residuum_source(1:sizeDotState)
+                                                          - r(1:sizeDotState)
               crystallite_converged(g,i,e) = &
-              crystallite_converged(g,i,e) .and. converged(residuum_source(1:sizeDotState), &
+              crystallite_converged(g,i,e) .and. converged(r(1:sizeDotState), &
                                                            sourceState(p)%p(s)%state(1:sizeDotState,c), &
                                                            sourceState(p)%p(s)%atol(1:sizeDotState))
             enddo
 
             if(crystallite_converged(g,i,e)) then
               crystallite_todo(g,i,e) = stateJump(g,i,e)
-              if(.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) &
-                nonlocalBroken = .true.
               exit iteration
             endif
 
           enddo iteration
+          if(.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) &
+            nonlocalBroken = .true.
 
         endif
   enddo; enddo; enddo
   !$OMP END PARALLEL DO
 
-  if(nonlocalBroken) where(.not. crystallite_localPlasticity) crystallite_todo = .false.
+  if (nonlocalBroken) call nonlocalConvergenceCheck
-  if (any(plasticState(:)%nonlocal)) call nonlocalConvergenceCheck
 
   contains
 
@@ -1232,8 +1216,7 @@ subroutine integrateStateEuler
   enddo; enddo; enddo
   !$OMP END PARALLEL DO
 
-  if(nonlocalBroken) where(.not. crystallite_localPlasticity) crystallite_todo = .false.
+  if (nonlocalBroken) call nonlocalConvergenceCheck
-  if (any(plasticState(:)%nonlocal)) call nonlocalConvergenceCheck
 
 end subroutine integrateStateEuler
 
@@ -1254,17 +1237,11 @@ subroutine integrateStateAdaptiveEuler
   logical :: &
     nonlocalBroken
 
-  real(pReal), dimension(constitutive_plasticity_maxSizeDotState,            &
+  real(pReal), dimension(constitutive_plasticity_maxSizeDotState) :: residuum_plastic
-                         homogenization_maxNgrains,discretization_nIP,discretization_nElem) :: &
+  real(pReal), dimension(constitutive_source_maxSizeDotState,maxval(phase_Nsources)) :: residuum_source
-    residuum_plastic
-  real(pReal), dimension(constitutive_source_maxSizeDotState,&
-                         maxval(phase_Nsources), &
-                         homogenization_maxNgrains,discretization_nIP,discretization_nElem) :: &
-    residuum_source
-
 
   nonlocalBroken = .false.
-  !$OMP PARALLEL DO PRIVATE(sizeDotState,p,c)
+  !$OMP PARALLEL DO PRIVATE(sizeDotState,p,c,residuum_plastic,residuum_source)
   do e = FEsolving_execElem(1),FEsolving_execElem(2)
     do i = FEsolving_execIP(1),FEsolving_execIP(2)
       do g = 1,homogenization_Ngrains(material_homogenizationAt(e))
@@ -1287,15 +1264,14 @@ subroutine integrateStateAdaptiveEuler
 
           sizeDotState = plasticState(p)%sizeDotState
 
-          residuum_plastic(1:sizeDotState,g,i,e) = plasticState(p)%dotstate(1:sizeDotState,c) &
+          residuum_plastic(1:sizeDotState) = - plasticState(p)%dotstate(1:sizeDotState,c) * 0.5_pReal * crystallite_subdt(g,i,e)
-                                                 * (- 0.5_pReal * crystallite_subdt(g,i,e))
           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
 
-            residuum_source(1:sizeDotState,s,g,i,e) = sourceState(p)%p(s)%dotstate(1:sizeDotState,c) &
+            residuum_source(1:sizeDotState,s)  = - sourceState(p)%p(s)%dotstate(1:sizeDotState,c) &
-                                                    * (- 0.5_pReal * crystallite_subdt(g,i,e))
+                                               * 0.5_pReal * crystallite_subdt(g,i,e)
             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
@@ -1330,21 +1306,17 @@ subroutine integrateStateAdaptiveEuler
 
           sizeDotState = plasticState(p)%sizeDotState
 
-          residuum_plastic(1:sizeDotState,g,i,e) = residuum_plastic(1:sizeDotState,g,i,e) &
+          crystallite_converged(g,i,e) = converged(residuum_plastic(1:sizeDotState) &
-                                                 + 0.5_pReal * plasticState(p)%dotState(:,c) * crystallite_subdt(g,i,e)
+                                                   + 0.5_pReal * plasticState(p)%dotState(:,c) * crystallite_subdt(g,i,e), &
-
-          crystallite_converged(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
 
-            residuum_source(1:sizeDotState,s,g,i,e) = &
-            residuum_source(1:sizeDotState,s,g,i,e) + 0.5_pReal * sourceState(p)%p(s)%dotState(:,c) * crystallite_subdt(g,i,e)
-
             crystallite_converged(g,i,e) = &
-            crystallite_converged(g,i,e) .and. converged(residuum_source(1:sizeDotState,s,g,i,e), &
+            crystallite_converged(g,i,e) .and. converged(residuum_source(1:sizeDotState,s) &
+                                                         + 0.5_pReal*sourceState(p)%p(s)%dotState(:,c)*crystallite_subdt(g,i,e), &
                                                          sourceState(p)%p(s)%state(1:sizeDotState,c), &
                                                          sourceState(p)%p(s)%atol(1:sizeDotState))
            enddo
@@ -1353,7 +1325,7 @@ subroutine integrateStateAdaptiveEuler
   enddo; enddo; enddo
   !$OMP END PARALLEL DO
 
-  if (any(plasticState(:)%nonlocal)) call nonlocalConvergenceCheck
+  if (nonlocalBroken) call nonlocalConvergenceCheck
 
 end subroutine integrateStateAdaptiveEuler
 
@@ -1387,8 +1359,10 @@ subroutine integrateStateRK4
   logical :: &
     nonlocalBroken
 
+  real(pReal), dimension(constitutive_plasticity_maxSizeDotState,4) :: plastic_RK4dotState
+  real(pReal), dimension(constitutive_source_maxSizeDotState,4,maxval(phase_Nsources)) :: source_RK4dotState
   nonlocalBroken = .false.
-  !$OMP PARALLEL DO PRIVATE(sizeDotState,p,c)
+  !$OMP PARALLEL DO PRIVATE(sizeDotState,p,c,plastic_RK4dotState,source_RK4dotState)
   do e = FEsolving_execElem(1),FEsolving_execElem(2)
     do i = FEsolving_execIP(1),FEsolving_execIP(2)
       do g = 1,homogenization_Ngrains(material_homogenizationAt(e))
@@ -1410,20 +1384,23 @@ subroutine integrateStateRK4
           if(.not. crystallite_todo(g,i,e)) cycle
 
           do stage = 1,3
-
+            sizeDotState = plasticState(p)%sizeDotState
-            plasticState(p)%RK4dotState(stage,:,c) = plasticState(p)%dotState(:,c)
+            plastic_RK4dotState(1:sizeDotState,stage) = plasticState(p)%dotState(:,c)
-            plasticState(p)%dotState(:,c) = A(1,stage) * plasticState(p)%RK4dotState(1,:,c)
+            plasticState(p)%dotState(:,c) = A(1,stage) * plastic_RK4dotState(1:sizeDotState,1)
             do s = 1, phase_Nsources(p)
-              sourceState(p)%p(s)%RK4dotState(stage,:,c) = sourceState(p)%p(s)%dotState(:,c)
+              sizeDotState = sourceState(p)%p(s)%sizeDotState
-              sourceState(p)%p(s)%dotState(:,c) = A(1,stage) * sourceState(p)%p(s)%RK4dotState(1,:,c)
+              source_RK4dotState(1:sizeDotState,stage,s) = sourceState(p)%p(s)%dotState(:,c)
+              sourceState(p)%p(s)%dotState(:,c) = A(1,stage) * source_RK4dotState(1:sizeDotState,1,s)
             enddo
 
             do n = 2, stage
+              sizeDotState = plasticState(p)%sizeDotState
               plasticState(p)%dotState(:,c) = plasticState(p)%dotState(:,c) &
-                                            + A(n,stage) * plasticState(p)%RK4dotState(n,:,c)
+                                            + A(n,stage) * plastic_RK4dotState(1:sizeDotState,n)
               do s = 1, phase_Nsources(p)
+                sizeDotState = sourceState(p)%p(s)%sizeDotState
                 sourceState(p)%p(s)%dotState(:,c) = sourceState(p)%p(s)%dotState(:,c) &
-                                                  + A(n,stage) * sourceState(p)%p(s)%RK4dotState(n,:,c)
+                                                  + A(n,stage) * source_RK4dotState(1:sizeDotState,n,s)
               enddo
             enddo
 
@@ -1466,9 +1443,9 @@ subroutine integrateStateRK4
 
           sizeDotState = plasticState(p)%sizeDotState
 
-          plasticState(p)%RK4dotState(4,:,c) = plasticState (p)%dotState(:,c)
+          plastic_RK4dotState(1:sizeDotState,4) = plasticState (p)%dotState(:,c)
 
-          plasticState(p)%dotState(:,c) =  matmul(B,plasticState(p)%RK4dotState(1:4,1:sizeDotState,c))
+          plasticState(p)%dotState(:,c) =  matmul(plastic_RK4dotState(1:sizeDotState,1:4),B)
           plasticState(p)%state(1:sizeDotState,c) = plasticState(p)%subState0(1:sizeDotState,c) &
                                                   + plasticState(p)%dotState (1:sizeDotState,c) &
                                                     * crystallite_subdt(g,i,e)
@@ -1476,9 +1453,9 @@ subroutine integrateStateRK4
           do s = 1, phase_Nsources(p)
             sizeDotState = sourceState(p)%p(s)%sizeDotState
 
-            sourceState(p)%p(s)%RK4dotState(4,:,c) = sourceState(p)%p(s)%dotState(:,c)
+            source_RK4dotState(1:sizeDotState,4,s) = sourceState(p)%p(s)%dotState(:,c)
 
-            sourceState(p)%p(s)%dotState(:,c)  = matmul(B,sourceState(p)%p(s)%RK4dotState(1:4,1:sizeDotState,c))
+            sourceState(p)%p(s)%dotState(:,c)  = matmul(source_RK4dotState(1:sizeDotState,1:4,s),B)
             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)
@@ -1506,8 +1483,7 @@ subroutine integrateStateRK4
   enddo; enddo; enddo
   !$OMP END PARALLEL DO
 
-  if(nonlocalBroken) where(.not. crystallite_localPlasticity) crystallite_todo = .false.
+  if (nonlocalBroken) call nonlocalConvergenceCheck
-  if (any(plasticState(:)%nonlocal)) call nonlocalConvergenceCheck
 
 end subroutine integrateStateRK4
 
@@ -1550,18 +1526,11 @@ subroutine integrateStateRKCK45
     sizeDotState
   logical :: &
     nonlocalBroken
-
+  real(pReal), dimension(constitutive_plasticity_maxSizeDotState,6) :: plastic_RKdotState
-  real(pReal), dimension(constitutive_plasticity_maxSizeDotState,            &
+  real(pReal), dimension(constitutive_source_maxSizeDotState,6,maxval(phase_Nsources)) :: source_RKdotState
-                         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
-
 
   nonlocalBroken = .false.
-  !$OMP PARALLEL DO PRIVATE(sizeDotState,p,c)
+  !$OMP PARALLEL DO PRIVATE(sizeDotState,p,c,plastic_RKdotState,source_RKdotState)
   do e = FEsolving_execElem(1),FEsolving_execElem(2)
     do i = FEsolving_execIP(1),FEsolving_execIP(2)
       do g = 1,homogenization_Ngrains(material_homogenizationAt(e))
@@ -1583,20 +1552,23 @@ subroutine integrateStateRKCK45
           if(.not. crystallite_todo(g,i,e)) cycle
 
           do stage = 1,5
-
+            sizeDotState = plasticState(p)%sizeDotState
-            plasticState(p)%RKCK45dotState(stage,:,c) = plasticState(p)%dotState(:,c)
+            plastic_RKdotState(1:sizeDotState,stage) = plasticState(p)%dotState(:,c)
-            plasticState(p)%dotState(:,c) = A(1,stage) * plasticState(p)%RKCK45dotState(1,:,c)
+            plasticState(p)%dotState(:,c) = A(1,stage) * plastic_RKdotState(1:sizeDotState,1)
             do s = 1, phase_Nsources(p)
-              sourceState(p)%p(s)%RKCK45dotState(stage,:,c) = sourceState(p)%p(s)%dotState(:,c)
+              sizeDotState = sourceState(p)%p(s)%sizeDotState
-              sourceState(p)%p(s)%dotState(:,c) = A(1,stage) * sourceState(p)%p(s)%RKCK45dotState(1,:,c)
+              source_RKdotState(1:sizeDotState,stage,s) = sourceState(p)%p(s)%dotState(:,c)
+              sourceState(p)%p(s)%dotState(:,c) = A(1,stage) * source_RKdotState(1:sizeDotState,1,s)
             enddo
 
             do n = 2, stage
+              sizeDotState = plasticState(p)%sizeDotState
               plasticState(p)%dotState(:,c) = plasticState(p)%dotState(:,c) &
-                                            + A(n,stage) * plasticState(p)%RKCK45dotState(n,:,c)
+                                            + A(n,stage) * plastic_RKdotState(1:sizeDotState,n)
               do s = 1, phase_Nsources(p)
+                sizeDotState = sourceState(p)%p(s)%sizeDotState
                 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) * source_RKdotState(1:sizeDotState,n,s)
               enddo
             enddo
 
@@ -1639,29 +1611,27 @@ subroutine integrateStateRKCK45
 
           sizeDotState = plasticState(p)%sizeDotState
 
-          plasticState(p)%RKCK45dotState(6,:,c) = plasticState (p)%dotState(:,c)
+          plastic_RKdotState(1:sizeDotState,6) = plasticState (p)%dotState(:,c)
-          residuum_plastic(1:sizeDotState,g,i,e) = matmul(DB,plasticState(p)%RKCK45dotState(1:6,1:sizeDotState,c)) &
+          plasticState(p)%dotState(:,c) = matmul(plastic_RKdotState(1:sizeDotState,1:6),B)
-                                                   * 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)
-          crystallite_todo(g,i,e) = converged(residuum_plastic(1:sizeDotState,g,i,e), &
+          crystallite_todo(g,i,e) = converged( matmul(plastic_RKdotState(1:sizeDotState,1:6),DB) &
+                                                   * crystallite_subdt(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)
+            source_RKdotState(1:sizeDotState,6,s) = 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)) &
+            sourceState(p)%p(s)%dotState(:,c)  = matmul(source_RKdotState(1:sizeDotState,1:6,s),B)
-                                                      * 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), &
+                                      converged(matmul(source_RKdotState(1:sizeDotState,1:6,s),DB) &
+                                                      * crystallite_subdt(g,i,e), &
                                                 sourceState(p)%p(s)%state(1:sizeDotState,c), &
                                                 sourceState(p)%p(s)%atol(1:sizeDotState))
           enddo
@@ -1687,8 +1657,7 @@ subroutine integrateStateRKCK45
   enddo; enddo; enddo
   !$OMP END PARALLEL DO
 
-  if(nonlocalBroken) where(.not. crystallite_localPlasticity) crystallite_todo = .false.
+  if (nonlocalBroken) call nonlocalConvergenceCheck
-  if (any(plasticState(:)%nonlocal)) call nonlocalConvergenceCheck
 
 end subroutine integrateStateRKCK45
 
@@ -1699,7 +1668,6 @@ end subroutine integrateStateRKCK45
 !--------------------------------------------------------------------------------------------------
 subroutine nonlocalConvergenceCheck
 
- if (any(.not. crystallite_converged .and. .not. crystallite_localPlasticity)) &                    ! any non-local not yet converged (or broken)...
   where( .not. crystallite_localPlasticity) crystallite_converged = .false.
 
 end subroutine nonlocalConvergenceCheck

src/grid/spectral_utilities.f90

@@ -1056,14 +1056,15 @@ subroutine utilities_updateCoords(F)
   call MPI_Irecv(IPfluct_padded(:,:,:,grid3+2),c,MPI_DOUBLE,rank_t,0,PETSC_COMM_WORLD,r,ierr)
   if(ierr /=0) call IO_error(894, ext_msg='update_IPcoords/MPI_Irecv')
   call MPI_Wait(r,s,ierr)
+  if(ierr /=0) call IO_error(894, ext_msg='update_IPcoords/MPI_Wait')
 
   ! send top layer to process above
-  if(ierr /=0) call IO_error(894, ext_msg='update_IPcoords/MPI_Wait')
   call MPI_Isend(IPfluct_padded(:,:,:,grid3+1),c,MPI_DOUBLE,rank_t,0,PETSC_COMM_WORLD,r,ierr)
   if(ierr /=0) call IO_error(894, ext_msg='update_IPcoords/MPI_Isend')
   call MPI_Irecv(IPfluct_padded(:,:,:,1),      c,MPI_DOUBLE,rank_b,0,PETSC_COMM_WORLD,r,ierr)
   if(ierr /=0) call IO_error(894, ext_msg='update_IPcoords/MPI_Irecv')
   call MPI_Wait(r,s,ierr)
+  if(ierr /=0) call IO_error(894, ext_msg='update_IPcoords/MPI_Wait')
 
  !--------------------------------------------------------------------------------------------------
  ! calculate nodal displacements

src/homogenization.f90

@@ -29,6 +29,8 @@ module homogenization
 
 !--------------------------------------------------------------------------------------------------
 ! General variables for the homogenization at a  material point
+  logical, public :: &
+    terminallyIll = .false.                                                                         !< at least one material point is terminally ill
   real(pReal),   dimension(:,:,:,:),   allocatable, public :: &
     materialpoint_F0, &                                                                             !< def grad of IP at start of FE increment
     materialpoint_F, &                                                                              !< def grad of IP to be reached at end of FE increment
@@ -171,22 +173,6 @@ subroutine homogenization_init
 
   write(6,'(/,a)') ' <<<+-  homogenization init  -+>>>'; flush(6)
 
-  if (iand(debug_level(debug_homogenization), debug_levelBasic) /= 0) then
-    write(6,'(a32,1x,7(i8,1x))')   'materialpoint_dPdF:             ', shape(materialpoint_dPdF)
-    write(6,'(a32,1x,7(i8,1x))')   'materialpoint_F0:               ', shape(materialpoint_F0)
-    write(6,'(a32,1x,7(i8,1x))')   'materialpoint_F:                ', shape(materialpoint_F)
-    write(6,'(a32,1x,7(i8,1x))')   'materialpoint_subF0:            ', shape(materialpoint_subF0)
-    write(6,'(a32,1x,7(i8,1x))')   'materialpoint_subF:             ', shape(materialpoint_subF)
-    write(6,'(a32,1x,7(i8,1x))')   'materialpoint_P:                ', shape(materialpoint_P)
-    write(6,'(a32,1x,7(i8,1x))')   'materialpoint_subFrac:          ', shape(materialpoint_subFrac)
-    write(6,'(a32,1x,7(i8,1x))')   'materialpoint_subStep:          ', shape(materialpoint_subStep)
-    write(6,'(a32,1x,7(i8,1x))')   'materialpoint_subdt:            ', shape(materialpoint_subdt)
-    write(6,'(a32,1x,7(i8,1x))')   'materialpoint_requested:        ', shape(materialpoint_requested)
-    write(6,'(a32,1x,7(i8,1x))')   'materialpoint_converged:        ', shape(materialpoint_converged)
-    write(6,'(a32,1x,7(i8,1x),/)') 'materialpoint_doneAndHappy:     ', shape(materialpoint_doneAndHappy)
-  endif
-  flush(6)
-
   if (debug_g < 1 .or. debug_g > homogenization_Ngrains(material_homogenizationAt(debug_e))) &
     call IO_error(602,ext_msg='constituent', el=debug_e, g=debug_g)
 

src/lattice.f90

@@ -18,16 +18,16 @@ module lattice
 
 !--------------------------------------------------------------------------------------------------
 ! face centered cubic
-  integer, dimension(2), parameter :: &
+  integer, dimension(*), parameter :: &
     FCC_NSLIPSYSTEM = [12, 6]                                                                       !< # of slip systems per family for fcc
 
-  integer, dimension(1), parameter :: &
+  integer, dimension(*), parameter :: &
     FCC_NTWINSYSTEM = [12]                                                                          !< # of twin systems per family for fcc
 
-  integer, dimension(1), parameter :: &
+  integer, dimension(*), parameter :: &
     FCC_NTRANSSYSTEM = [12]                                                                         !< # of transformation systems per family for fcc
 
-  integer, dimension(1), parameter :: &
+  integer, dimension(*), parameter :: &
     FCC_NCLEAVAGESYSTEM = [3]                                                                       !< # of cleavage systems per family for fcc
 
   integer, parameter  :: &
@@ -109,13 +109,13 @@ module lattice
 
 !--------------------------------------------------------------------------------------------------
 ! body centered cubic
-  integer, dimension(2), parameter :: &
+  integer, dimension(*), parameter :: &
     BCC_NSLIPSYSTEM = [12, 12]                                                                      !< # of slip systems per family for bcc
 
-  integer, dimension(1), parameter :: &
+  integer, dimension(*), parameter :: &
     BCC_NTWINSYSTEM = [12]                                                                          !< # of twin systems per family for bcc
 
-  integer, dimension(1), parameter :: &
+  integer, dimension(*), parameter :: &
     BCC_NCLEAVAGESYSTEM = [3]                                                                       !< # of cleavage systems per family for bcc
 
   integer, parameter  :: &
@@ -187,10 +187,10 @@ module lattice
 
 !--------------------------------------------------------------------------------------------------
 ! hexagonal
-  integer, dimension(6), parameter :: &
+  integer, dimension(*), parameter :: &
     HEX_NSLIPSYSTEM = [3, 3, 3, 6, 12, 6]                                                           !< # of slip systems per family for hex
 
-  integer, dimension(4), parameter :: &
+  integer, dimension(*), parameter :: &
     HEX_NTWINSYSTEM = [6, 6, 6, 6]                                                                  !< # of slip systems per family for hex
 
   integer, parameter  :: &
@@ -280,7 +280,7 @@ module lattice
 
 !--------------------------------------------------------------------------------------------------
 ! body centered tetragonal
-  integer, dimension(13), parameter :: &
+  integer, dimension(*), parameter :: &
     BCT_NSLIPSYSTEM = [2, 2, 2, 4, 2, 4, 2, 2, 4, 8, 4, 8, 8 ]                                      !< # of slip systems per family for bct (Sn) Bieler J. Electr Mater 2009
 
   integer, parameter :: &
@@ -362,7 +362,7 @@ module lattice
 
 !--------------------------------------------------------------------------------------------------
 ! orthorhombic
-  integer, dimension(3), parameter :: &
+  integer, dimension(*), parameter :: &
     ORT_NCLEAVAGESYSTEM = [1, 1, 1]                                                                 !< # of cleavage systems per family for ortho
 
   integer, parameter  :: &

src/material.f90

@@ -724,14 +724,6 @@ subroutine material_allocatePlasticState(phase,NipcMyPhase,&
   allocate(plasticState(phase)%state               (sizeState,NipcMyPhase),   source=0.0_pReal)
 
   allocate(plasticState(phase)%dotState            (sizeDotState,NipcMyPhase),source=0.0_pReal)
-  if (numerics_integrator == 1) then
-    allocate(plasticState(phase)%previousDotState  (sizeDotState,NipcMyPhase),source=0.0_pReal)
-    allocate(plasticState(phase)%previousDotState2 (sizeDotState,NipcMyPhase),source=0.0_pReal)
-  endif
-  if (numerics_integrator == 4) &
-    allocate(plasticState(phase)%RK4dotState       (4,sizeDotState,NipcMyPhase),source=0.0_pReal)
-  if (numerics_integrator == 5) &
-    allocate(plasticState(phase)%RKCK45dotState  (6,sizeDotState,NipcMyPhase),source=0.0_pReal)
 
   allocate(plasticState(phase)%deltaState        (sizeDeltaState,NipcMyPhase),source=0.0_pReal)
 
@@ -762,14 +754,6 @@ subroutine material_allocateSourceState(phase,of,NipcMyPhase,&
   allocate(sourceState(phase)%p(of)%state               (sizeState,NipcMyPhase),   source=0.0_pReal)
 
   allocate(sourceState(phase)%p(of)%dotState            (sizeDotState,NipcMyPhase),source=0.0_pReal)
-  if (numerics_integrator == 1) then
-    allocate(sourceState(phase)%p(of)%previousDotState  (sizeDotState,NipcMyPhase),source=0.0_pReal)
-    allocate(sourceState(phase)%p(of)%previousDotState2 (sizeDotState,NipcMyPhase),source=0.0_pReal)
-  endif
-  if (numerics_integrator == 4) &
-    allocate(sourceState(phase)%p(of)%RK4dotState     (4,sizeDotState,NipcMyPhase),source=0.0_pReal)
-  if (numerics_integrator == 5) &
-    allocate(sourceState(phase)%p(of)%RKCK45dotState  (6,sizeDotState,NipcMyPhase),source=0.0_pReal)
 
   allocate(sourceState(phase)%p(of)%deltaState        (sizeDeltaState,NipcMyPhase),source=0.0_pReal)
 

src/math.f90

@@ -30,14 +30,12 @@ module math
       1.0_pReal,0.0_pReal,0.0_pReal, &
       0.0_pReal,1.0_pReal,0.0_pReal, &
       0.0_pReal,0.0_pReal,1.0_pReal  &
-      ],[3,3])                                                                                      !< 3x3 Identity
+      ],shape(math_I3))                                                                             !< 3x3 Identity
 
-  real(pReal), dimension(6), parameter, private :: &
+  real(pReal), dimension(*), parameter, private :: &
-    NRMMANDEL = [&
+    NRMMANDEL = [1.0_pReal, 1.0_pReal,1.0_pReal, sqrt(2.0_pReal), sqrt(2.0_pReal), sqrt(2.0_pReal)] !< forward weighting for Mandel notation
-      1.0_pReal,       1.0_pReal,       1.0_pReal, &
-      sqrt(2.0_pReal), sqrt(2.0_pReal), sqrt(2.0_pReal) ]                                           !< forward weighting for Mandel notation
 
-  real(pReal), dimension(6), parameter, private :: &
+  real(pReal), dimension(*), parameter, private :: &
     INVNRMMANDEL = 1.0_pReal/NRMMANDEL                                                              !< backward weighting for Mandel notation
 
   integer, dimension (2,6), parameter, private :: &
@@ -48,7 +46,7 @@ module math
       1,2, &
       2,3, &
       1,3  &
-      ],[2,6])                                                                                      !< arrangement in Nye notation.
+      ],shape(MAPNYE))                                                                              !< arrangement in Nye notation.
 
   integer, dimension (2,6), parameter, private :: &
     MAPVOIGT = reshape([&
@@ -58,7 +56,7 @@ module math
       2,3, &
       1,3, &
       1,2  &
-      ],[2,6])                                                                                      !< arrangement in Voigt notation
+      ],shape(MAPVOIGT))                                                                            !< arrangement in Voigt notation
 
   integer, dimension (2,9), parameter, private :: &
     MAPPLAIN = reshape([&
@@ -71,7 +69,7 @@ module math
       3,1, &
       3,2, &
       3,3  &
-      ],[2,9])                                                                                      !< arrangement in Plain notation
+      ],shape(MAPPLAIN))                                                                            !< arrangement in Plain notation
 
   interface math_eye
     module procedure math_identity2nd
@@ -487,21 +485,19 @@ function math_invSym3333(A)
 
   real(pReal),dimension(3,3,3,3),intent(in) :: A
 
-  integer :: ierr
   integer,     dimension(6)        :: ipiv6
   real(pReal), dimension(6,6)      :: temp66
   real(pReal), dimension(6*(64+2)) :: work
-  logical                          :: error
+  integer                          :: ierr_i, ierr_f
   external :: &
     dgetrf, &
     dgetri
 
   temp66 = math_sym3333to66(A)
-  call dgetrf(6,6,temp66,6,ipiv6,ierr)
+  call dgetrf(6,6,temp66,6,ipiv6,ierr_i)
-  error = (ierr /= 0)
+  call dgetri(6,temp66,6,ipiv6,work,size(work,1),ierr_f)
-  call dgetri(6,temp66,6,ipiv6,work,size(work,1),ierr)
+
-  error = error .or. (ierr /= 0)
+  if (ierr_i /= 0 .or. ierr_f /= 0) then
-  if (error) then
     call IO_error(400, ext_msg = 'math_invSym3333')
   else
     math_invSym3333 = math_66toSym3333(temp66)

src/prec.f90

@@ -30,10 +30,6 @@ module prec
     real(pReal), dimension(:), pointer :: p
   end type group_float
 
-  type :: group_int
-    integer, dimension(:), pointer :: p
-  end type group_int
-
   ! http://stackoverflow.com/questions/3948210/can-i-have-a-pointer-to-an-item-in-an-allocatable-array
   type :: tState
     integer :: &
@@ -50,12 +46,7 @@ module prec
       deltaState                                                                                    !< increment of state change
     real(pReal), allocatable, dimension(:,:) :: &
       partionedState0, &
-      subState0, &
+      subState0
-      previousDotState, &
-      previousDotState2
-    real(pReal), allocatable, dimension(:,:,:) :: &
-      RK4dotState, &
-      RKCK45dotState
   end type
 
   type, extends(tState) :: tPlasticState