Merge branch 'development' into empty-table-init

CMakeLists.txt

@@ -108,7 +108,7 @@ file(STRINGS "$ENV{PETSC_DIR}/$ENV{PETSC_ARCH}/lib/petsc/conf/petscvariables" PE
 string(REPLACE "PETSC_FC_INCLUDES = " "" PETSC_INCLUDES "${PETSC_INCLUDES}")
 message("PETSC_INCLUDES:\n${PETSC_INCLUDES}\n")
 
-set(CMAKE_Fortran_FLAGS_${CMAKE_BUILD_TYPE} "${BUILDCMD_PRE} ${OPENMP_FLAGS} ${STANDARD_CHECK} ${OPTIMIZATION_FLAGS} ${COMPILE_FLAGS} ${PRECISION_FLAGS}")
+set(CMAKE_Fortran_FLAGS_${CMAKE_BUILD_TYPE} "${BUILDCMD_PRE} ${OPENMP_FLAGS} ${STANDARD_CHECK} ${OPTIMIZATION_FLAGS} ${COMPILE_FLAGS}")
 set(CMAKE_Fortran_LINK_EXECUTABLE           "${BUILDCMD_PRE} ${CMAKE_Fortran_COMPILER} ${OPENMP_FLAGS} ${OPTIMIZATION_FLAGS} ${LINKER_FLAGS}")
 
 if(CMAKE_BUILD_TYPE STREQUAL "DEBUG")

LICENSE

@@ -1,9 +1,9 @@
 Copyright 2011-2022 Max-Planck-Institut für Eisenforschung GmbH
 
 DAMASK is free software: you can redistribute it and/or modify
-it under the terms of the GNU Affero General Public License as published by
-the Free Software Foundation, either version 3 of the License, or
-(at your option) any later version.
+it under the terms of the GNU Affero General Public License as
+published by the Free Software Foundation, either version 3 of the
+License, or (at your option) any later version.
 
 This program is distributed in the hope that it will be useful,
 but WITHOUT ANY WARRANTY; without even the implied warranty of

PRIVATE

@@ -1 +1 @@
-Subproject commit 0e82975b23a1bd4310c523388f1cadf1b8e03dd0
+Subproject commit cdca8ab0c14b637c18279e0ea236caa148d15e5e

VERSION

@@ -1 +1 @@
-3.0.0-alpha6-395-g9696c1967
+3.0.0-alpha6-523-g66f129273

cmake/Compiler-GNU.cmake

@@ -135,10 +135,3 @@ set (DEBUG_FLAGS "${DEBUG_FLAGS} -fsanitize=undefined")
 # detect undefined behavior
 # Additional options
 # -fsanitize=address,leak,thread
-
-#------------------------------------------------------------------------------------------------
-#  precision settings
-set (PRECISION_FLAGS "${PRECISION_FLAGS} -fdefault-real-8")
-# set precision to 8 bytes for standard real (=8 for pReal). Will set size of double to 16 bytes as long as -fdefault-double-8 is not set
-set (PRECISION_FLAGS "${PRECISION_FLAGS} -fdefault-double-8")
-# set precision to 8 bytes for double real, would be 16 bytes if -fdefault-real-8 is used

cmake/Compiler-Intel.cmake

@@ -118,8 +118,3 @@ set (DEBUG_FLAGS "${DEBUG_FLAGS} -debug all")
 # -check:                  Checks at runtime, where
 #    arg_temp_created:       will cause a lot of warnings because we create a bunch of temporary arrays (performance?)
 #    stack:
-
-#------------------------------------------------------------------------------------------------
-#  precision settings
-set (PRECISION_FLAGS "${PRECISION_FLAGS} -real-size 64")
-# set precision for standard real to 32 | 64 | 128 (= 4 | 8 | 16 bytes, type pReal is always 8 bytes)

cmake/Compiler-IntelLLVM.cmake

@@ -117,8 +117,3 @@ set (DEBUG_FLAGS "${DEBUG_FLAGS} -debug all")
 # -check:                  Checks at runtime, where
 #    arg_temp_created:       will cause a lot of warnings because we create a bunch of temporary arrays (performance?)
 #    stack:
-
-#------------------------------------------------------------------------------------------------
-#  precision settings
-set (PRECISION_FLAGS "${PRECISION_FLAGS} -real-size 64")
-# set precision for standard real to 32 | 64 | 128 (= 4 | 8 | 16 bytes, type pReal is always 8 bytes)

examples/config/phase/damage/anisobrittle_cubic.yaml

@@ -9,5 +9,5 @@ s_crit: [0.006666]
 dot_o: 1.e-3
 q: 20
 
-D_11: 1.0
+l_c: 1.0
 mu: 0.001

examples/config/phase/damage/isobrittle_generic.yaml

@@ -2,6 +2,6 @@ type: isobrittle
 
 output: [f_phi]
 
-W_crit: 1400000.0
-D_11: 1.0
+G_crit: 1400000.0
+l_c: 1.0
 mu: 0.001

examples/config/phase/mechanical/plastic/dislotwin_alpha-Brass-shearbanding.yaml

@@ -0,0 +1,15 @@
+type: dislotwin
+
+references:
+  - N. Jia et al.,
+    Acta Materialia 60(3):1099-1115, 2012,
+    https://doi.org/10.1016/j.actamat.2011.10.047
+  - N. Jia et al.,
+    Acta Materialia 60:3415-3434, 2012,
+    https://doi.org/10.1016/j.actamat.2012.03.005
+
+gamma_0_sb: 0.0001
+tau_sb: 180.0e6 # tau_hat_sb
+Q_sb: 4.0e-19 # Q_0
+p_sb: 1.15
+q_sb: 1.0

install/MarcMentat/2021.3.1/Marc_tools/include_linux64.patch

@@ -1,18 +1,22 @@
 ---
 +++
-@@ -119,6 +119,11 @@ if test "$MSCCOSIM_VERSION" = ""; then
+@@ -119,6 +119,15 @@ if test "$MSCCOSIM_VERSION" = ""; then
    MSCCOSIM_VERSION="2020"
  fi
 
 +# DAMASK uses the HDF5 compiler wrapper around the Intel compiler
-+H5FC="$(h5fc -shlib -show)"
-+HDF5_LIB=${H5FC//ifort/}
++H5FC=$(h5fc -shlib -show)
++if [[ "$H5FC" == *"$dir is"* ]]; then
++  H5FC=$(echo $(echo "$H5FC" | tail -n1) | sed -e "s/\-shlib/-fPIC -integer-size 64 -real-size 64 -qopenmp/g")
++  H5FC=${H5FC%-lmpifort*}
++fi
++HDF5_LIB=${H5FC//*"ifort"/}
 +FCOMP="$H5FC"
 +
  # AEM
  if test "$MARCDLLOUTDIR" = ""; then
   DLLOUTDIR="$MARC_LIB"
-@@ -439,7 +444,7 @@ if test "$MARC_INTEGER_SIZE" = "i4" ; then
+@@ -439,7 +448,7 @@ if test "$MARC_INTEGER_SIZE" = "i4" ; then
    I8DEFINES=
    I8CDEFINES=
  else
@@ -22,7 +26,7 @@
    I8CDEFINES="-U_DOUBLE -D_SINGLE"
  fi
 
-@@ -556,7 +561,7 @@ then
+@@ -556,7 +565,7 @@ then
  PROFILE=" $PROFILE -pg"
  fi
 
@@ -31,7 +35,7 @@
  if test "$MTHREAD" = "OPENMP"
  then
    FORT_OPT=" $FORT_OPT -qopenmp"
-@@ -569,7 +574,7 @@ else
+@@ -569,7 +578,7 @@ else
    FORT_OPT=" $FORT_OPT -save -zero"
  fi
  if test "$MARCHDF_HDF" = "HDF"; then
@@ -40,7 +44,7 @@
  fi
 
  FORTLOW="$FCOMP $FORT_OPT $PROFILE -O0 $I8FFLAGS -I$MARC_SOURCE/common \
-@@ -583,6 +588,30 @@ FORTNA="$FCOMP $FORT_OPT -fno-alias -O3 $I8FFLAGS -I$MARC_SOURCE/common \
+@@ -583,6 +592,30 @@ FORTNA="$FCOMP $FORT_OPT -fno-alias -O3 $I8FFLAGS -I$MARC_SOURCE/common \
  #  for compiling free form f90 files. high opt, integer(4)
  FORTF90="$FCOMP -c -O3"
 
@@ -71,7 +75,7 @@
  if test "$MARCDEBUG" = "ON"
  then
    FORTLOW="$FCOMP $FORT_OPT $PROFILE $I8FFLAGS -I$MARC_SOURCE/common \
-@@ -739,7 +768,7 @@ SECLIBS="-L$MARC_LIB -llapi"
+@@ -739,7 +772,7 @@ SECLIBS="-L$MARC_LIB -llapi"
 
  SOLVERLIBS="${BCSSOLVERLIBS} ${VKISOLVERLIBS} ${CASISOLVERLIBS} ${MF2SOLVERLIBS} \
    -L$MARC_MKL \
@@ -80,7 +84,7 @@
 
  SOLVERLIBS_DLL=${SOLVERLIBS}
  if test "$AEM_DLL" -eq 1
-@@ -762,7 +791,7 @@ then
+@@ -762,7 +795,7 @@ then
   OPENSSL=NONE
  fi
 
@@ -89,7 +93,7 @@
 
  # Uncomment the following lines to turn on the trace and comment out the next 4 lines
  # if test $MPITYPE = intelmpi
-@@ -772,7 +801,7 @@ SYSLIBS=" $OPENMP -lpthread -shared-intel -cxxlib"
+@@ -772,7 +805,7 @@ SYSLIBS=" $OPENMP -lpthread -shared-intel -cxxlib"
  # fi
  if test $MPITYPE = intelmpi
  then

install/MarcMentat/2022.1/Marc_tools/include_linux64.patch

@@ -5,14 +5,18 @@
  fi
  
 +# DAMASK uses the HDF5 compiler wrapper around the Intel compiler
-+H5FC="$(h5fc -shlib -show)"
-+HDF5_LIB=${H5FC//ifort/}
++H5FC=$(h5fc -shlib -show)
++if [[ "$H5FC" == *"$dir is"* ]]; then
++  H5FC=$(echo $(echo "$H5FC" | tail -n1) | sed -e "s/\-shlib/-fPIC -integer-size 64 -real-size 64 -qopenmp/g")
++  H5FC=${H5FC%-lmpifort*}
++fi
++HDF5_LIB=${H5FC//*"ifort"/}
 +FCOMP="$H5FC"
 +
  # AEM
  if test "$MARCDLLOUTDIR" = ""; then
   DLLOUTDIR="$MARC_LIB"
-@@ -477,8 +482,8 @@ if test "$MARC_INTEGER_SIZE" = "i4" ; then
+@@ -477,8 +486,8 @@ if test "$MARC_INTEGER_SIZE" = "i4" ; then
    I8DEFINES=
    I8CDEFINES=
  else
@@ -22,7 +26,7 @@
    I8CDEFINES="-U_DOUBLE -D_SINGLE"
  fi
  
-@@ -594,7 +599,7 @@ then
+@@ -594,7 +605,7 @@ then
  PROFILE=" $PROFILE -pg"
  fi
  
@@ -31,7 +35,7 @@
  if test "$MTHREAD" = "OPENMP"
  then
    FORT_OPT=" $FORT_OPT -qopenmp"
-@@ -607,7 +612,7 @@ else
+@@ -607,7 +616,7 @@ else
    FORT_OPT=" $FORT_OPT -save -zero"
  fi
  if test "$MARCHDF_HDF" = "HDF"; then
@@ -40,7 +44,7 @@
  fi
  
  FORTLOW="$FCOMP $FORT_OPT $PROFILE -O0 $I8FFLAGS -I$MARC_SOURCE/common \
-@@ -621,6 +626,29 @@ FORTNA="$FCOMP $FORT_OPT -fno-alias -O3 $I8FFLAGS -I$MARC_SOURCE/common \
+@@ -621,6 +630,29 @@ FORTNA="$FCOMP $FORT_OPT -fno-alias -O3 $I8FFLAGS -I$MARC_SOURCE/common \
  #  for compiling free form f90 files. high opt, integer(4)
  FORTF90="$FCOMP -c -O3"
 

install/MarcMentat/apply_DAMASK_modifications.py

@@ -6,15 +6,18 @@ import glob
 import argparse
 import shutil
 from pathlib import Path
+import subprocess
+import shlex
 
-import damask
+sys.path.append(str(Path(__file__).parents[2]/'python/damask'))
+import solver
 
 def copy_and_patch(patch,orig,editor):
     try:
         shutil.copyfile(orig,orig.parent/patch.stem)
     except shutil.SameFileError:
         pass
-    damask.util.run(f'patch {orig.parent/patch.stem} {patch} --backup --forward')
+    subprocess.run(shlex.split(f'patch {orig.parent/patch.stem} {patch} --backup --forward'))
     with open(orig.parent/patch.stem) as f_in:
         content = f_in.read()
     with open(orig.parent/patch.stem,'w') as f_out:
@@ -28,15 +31,16 @@ parser = argparse.ArgumentParser(
 parser.add_argument('--editor', dest='editor', metavar='string', default='vi',
                     help='Name of the editor (executable) used by Marc Mentat')
 parser.add_argument('--marc-root', dest='marc_root', metavar='string',
-                    default=damask.solver._marc._marc_root,
+                    default=solver._marc._marc_root,
                     help='Marc root directory')
 parser.add_argument('--marc-version', dest='marc_version', metavar='string',
-                    default=damask.solver._marc._marc_version,
+                    default=solver._marc._marc_version,
                     help='Marc version')
 parser.add_argument('--damask-root', dest='damask_root', metavar = 'string',
-                    default=damask.solver._marc._damask_root,
+                    default=solver._marc._damask_root,
                     help='DAMASK root directory')
 
+
 args = parser.parse_args()
 marc_root = Path(args.marc_root).expanduser()
 damask_root = Path(args.damask_root).expanduser()
@@ -52,7 +56,7 @@ matches = {'Marc_tools':  [['comp_user','comp_damask_*mp'],
 
 for cmd in ['patch','xvfb-run']:
     try:
-        damask.util.run(f'{cmd} --help')
+        subprocess.run(shlex.split(f'{cmd} --help'))
     except FileNotFoundError:
         print(f'"{cmd}" not found, please install')
         sys.exit()
@@ -71,7 +75,7 @@ print('compiling Mentat menu binaries...')
 
 executable = marc_root/f'mentat{marc_version}/bin/mentat'
 menu_file  = marc_root/f'mentat{marc_version}/menus/linux64/main.msb'
-damask.util.run(f'xvfb-run -a {executable} -compile {menu_file}')
+subprocess.run(shlex.split(f'xvfb-run -a {executable} -compile {menu_file}'))
 
 print('setting file access rights...')
 

python/damask/_colormap.py

@@ -48,7 +48,12 @@ class Colormap(mpl.colors.ListedColormap):
 
     def __eq__(self,
                other: object) -> bool:
-        """Test equality of colormaps."""
+        """
+        Return self==other.
+
+        Test equality of other.
+
+        """
         if not isinstance(other, Colormap):
             return NotImplemented
         return         len(self.colors) == len(other.colors) \
@@ -56,31 +61,61 @@ class Colormap(mpl.colors.ListedColormap):
 
     def __add__(self,
                 other: 'Colormap') -> 'Colormap':
-        """Concatenate."""
+        """
+        Return self+other.
+
+        Concatenate.
+
+        """
         return Colormap(np.vstack((self.colors,other.colors)),
                         f'{self.name}+{other.name}')
 
     def __iadd__(self,
                  other: 'Colormap') -> 'Colormap':
-        """Concatenate (in-place)."""
+        """
+        Return self+=other.
+
+        Concatenate (in-place).
+
+        """
         return self.__add__(other)
 
     def __mul__(self,
                 factor: int) -> 'Colormap':
-        """Repeat."""
+        """
+        Return self*other.
+
+        Repeat.
+
+        """
         return Colormap(np.vstack([self.colors]*factor),f'{self.name}*{factor}')
 
     def __imul__(self,
                  factor: int) -> 'Colormap':
-        """Repeat (in-place)."""
+        """
+        Return self*=other.
+
+        Repeat (in-place).
+
+        """
         return self.__mul__(factor)
 
     def __invert__(self) -> 'Colormap':
-        """Reverse."""
+        """
+        Return ~self.
+
+        Reverse.
+
+        """
         return self.reversed()
 
     def __repr__(self) -> str:
-        """Show as matplotlib figure."""
+        """
+        Return repr(self).
+
+        Show as matplotlib figure.
+
+        """
         fig = plt.figure(self.name,figsize=(5,.5))
         ax1 = fig.add_axes([0, 0, 1, 1])
         ax1.set_axis_off()
@@ -385,7 +420,7 @@ class Colormap(mpl.colors.ListedColormap):
         GOM_str = '1 1 {name} 9 {name} '.format(name=self.name.replace(" ","_")) \
                 +  '0 1 0 3 0 0 -1 9 \\ 0 0 0 255 255 255 0 0 255 ' \
                 + f'30 NO_UNIT 1 1 64 64 64 255 1 0 0 0 0 0 0 3 0 {self.N}' \
-                + ' '.join([f' 0 {c[0]} {c[1]} {c[2]} 255 1' for c in reversed((self.colors*255).astype(int))]) \
+                + ' '.join([f' 0 {c[0]} {c[1]} {c[2]} 255 1' for c in reversed((self.colors*255).astype(np.int64))]) \
                 + '\n'
 
         self._get_file_handle(fname,'.legend').write(GOM_str)

python/damask/_config.py

@@ -64,7 +64,12 @@ class Config(dict):
         super().__init__(**kwargs)
 
     def __repr__(self) -> str:
-        """Show as in file."""
+        """
+        Return repr(self).
+
+        Show as in file.
+
+        """
         output = StringIO()
         self.save(output)
         output.seek(0)
@@ -72,7 +77,12 @@ class Config(dict):
 
 
     def __copy__(self: MyType) -> MyType:
-        """Create deep copy."""
+        """
+        Return deepcopy(self).
+
+        Create deep copy.
+
+        """
         return copy.deepcopy(self)
 
     copy = __copy__
@@ -81,6 +91,8 @@ class Config(dict):
     def __or__(self: MyType,
                other) -> MyType:
         """
+        Return self|other.
+
         Update configuration with contents of other.
 
         Parameters
@@ -105,7 +117,12 @@ class Config(dict):
 
     def __ior__(self: MyType,
                 other) -> MyType:
-        """Update configuration with contents of other."""
+        """
+        Return self|=other.
+
+        Update configuration with contents of other.
+
+        """
         return self.__or__(other)
 
 

python/damask/_crystal.py

@@ -403,7 +403,12 @@ class Crystal():
 
 
     def __repr__(self):
-        """Give short human-readable summary."""
+        """
+        Return repr(self).
+
+        Give short human-readable summary.
+
+        """
         family = f'Crystal family: {self.family}'
         return family if self.lattice is None else \
                util.srepr([family,
@@ -415,7 +420,9 @@ class Crystal():
     def __eq__(self,
                other: object) -> bool:
         """
-        Equal to other.
+        Return self==other.
+
+        Test equality of other.
 
         Parameters
         ----------

python/damask/_grid.py

@@ -62,7 +62,12 @@ class Grid:
         self.comments = [] if comments_ is None else [str(c) for c in comments_]
 
     def __repr__(self) -> str:
-        """Give short human-readable summary."""
+        """
+        Return repr(self).
+
+        Give short human-readable summary.
+
+        """
         mat_min = np.nanmin(self.material)
         mat_max = np.nanmax(self.material)
         mat_N   = self.N_materials
@@ -76,7 +81,12 @@ class Grid:
 
 
     def __copy__(self) -> 'Grid':
-        """Create deep copy."""
+        """
+        Return deepcopy(self).
+
+        Create deep copy.
+
+        """
         return copy.deepcopy(self)
 
     copy = __copy__
@@ -85,6 +95,8 @@ class Grid:
     def __eq__(self,
                other: object) -> bool:
         """
+        Return self==other.
+
         Test equality of other.
 
         Parameters
@@ -117,8 +129,8 @@ class Grid:
         self._material = np.copy(material)
 
         if self.material.dtype in np.sctypes['float'] and \
-           np.all(self.material == self.material.astype(int).astype(float)):
-            self._material = self.material.astype(int)
+           np.all(self.material == self.material.astype(np.int64).astype(float)):
+            self._material = self.material.astype(np.int64)
 
 
     @property
@@ -285,7 +297,7 @@ class Grid:
             raise TypeError(f'mismatch between {cells.prod()} expected entries and {i} found')
 
         if not np.any(np.mod(material,1) != 0.0):                                                   # no float present
-            material = material.astype(int) - (1 if material.min() > 0 else 0)
+            material = material.astype(np.int64) - (1 if material.min() > 0 else 0)
 
         return Grid(material = material.reshape(cells,order='F'),
                     size     = size,
@@ -916,7 +928,7 @@ class Grid:
                                            cval=np.nanmax(self.material) + 1 if fill is None else fill)
             # avoid scipy interpolation errors for rotations close to multiples of 90°
             material = material_temp if np.prod(material_temp.shape) != np.prod(material.shape) else \
-                       np.rot90(material,k=np.rint(angle/90.).astype(int),axes=axes)
+                       np.rot90(material,k=np.rint(angle/90.).astype(np.int64),axes=axes)
 
         origin = self.origin-(np.asarray(material.shape)-self.cells)*.5 * self.size/self.cells
 
@@ -1094,7 +1106,7 @@ class Grid:
 
         rng = np.random.default_rng(rng_seed)
 
-        d = np.floor(distance).astype(int)
+        d = np.floor(distance).astype(np.int64)
         ext = np.linspace(-d,d,1+2*d,dtype=float),
         xx,yy,zz = np.meshgrid(ext,ext,ext)
         footprint = xx**2+yy**2+zz**2 <= distance**2+distance*1e-8
@@ -1197,7 +1209,7 @@ class Grid:
             mask = np.sum(np.power(coords_rot/r,2.0**np.array(exponent)),axis=-1) > 1.0
 
         if periodic:                                                                                # translate back to center
-            mask = np.roll(mask,((c/self.size-0.5)*self.cells).round().astype(int),(0,1,2))
+            mask = np.roll(mask,((c/self.size-0.5)*self.cells).round().astype(np.int64),(0,1,2))
 
         return Grid(material = np.where(np.logical_not(mask) if inverse else mask,
                                         self.material,
@@ -1249,7 +1261,7 @@ class Grid:
             return np.any(stencil != me if selection is None else
                           np.in1d(stencil,np.array(list(selection - {me}))))
 
-        d = np.floor(distance).astype(int)
+        d = np.floor(distance).astype(np.int64)
         ext = np.linspace(-d,d,1+2*d,dtype=float),
         xx,yy,zz = np.meshgrid(ext,ext,ext)
         footprint = xx**2+yy**2+zz**2 <= distance**2+distance*1e-8

python/damask/_orientation.py

@@ -120,14 +120,24 @@ class Orientation(Rotation,Crystal):
 
 
     def __repr__(self) -> str:
-        """Give short human-readable summary."""
+        """
+        Return repr(self).
+
+        Give short human-readable summary.
+
+        """
         return util.srepr([Crystal.__repr__(self),
                            Rotation.__repr__(self)])
 
 
     def __copy__(self: MyType,
                  rotation: Union[FloatSequence, Rotation] = None) -> MyType:
-        """Create deep copy."""
+        """
+        Return deepcopy(self).
+
+        Create deep copy.
+
+        """
         dup = copy.deepcopy(self)
         if rotation is not None:
             dup.quaternion = Rotation(rotation).quaternion
@@ -140,7 +150,9 @@ class Orientation(Rotation,Crystal):
     def __eq__(self,
                other: object) -> bool:
         """
-        Equal to other.
+        Return self==other.
+
+        Test equality of other.
 
         Parameters
         ----------
@@ -158,7 +170,9 @@ class Orientation(Rotation,Crystal):
     def __ne__(self,
                other: object) -> bool:
         """
-        Not equal to other.
+        Return self!=other.
+
+        Test inequality of other.
 
         Parameters
         ----------
@@ -448,9 +462,12 @@ class Orientation(Rotation,Crystal):
             elif self.family == 'orthorhombic':
                 return (np.prod(1. >= rho_abs,axis=-1)).astype(bool)
             elif self.family == 'monoclinic':
-                return (1. >= rho_abs[...,1]).astype(bool)
+                return np.logical_or(   1. >= rho_abs[...,1],
+                                     np.isnan(rho_abs[...,1]))
+            elif self.family == 'triclinic':
+                return np.ones(rho_abs.shape[:-1]).astype(bool)
             else:
-                return np.all(np.isfinite(rho_abs),axis=-1)
+                raise TypeError(f'unknown symmetry "{self.family}"')
 
 
     @property

python/damask/_result.py

@@ -83,7 +83,7 @@ class Result:
 
     >>> import damask
     >>> r = damask.Result('my_file.hdf5')
-    >>> r.add_Cauchy()
+    >>> r.add_stress_Cauchy()
     >>> r.add_equivalent_Mises('sigma')
     >>> r.export_VTK()
     >>> r_last = r.view(increments=-1)
@@ -152,14 +152,24 @@ class Result:
 
 
     def __copy__(self) -> "Result":
-        """Create deep copy."""
+        """
+        Return deepcopy(self).
+
+        Create deep copy.
+
+        """
         return copy.deepcopy(self)
 
     copy = __copy__
 
 
     def __repr__(self) -> str:
-        """Give short human-readable summary."""
+        """
+        Return repr(self).
+
+        Give short human-readable summary.
+
+        """
         with h5py.File(self.fname,'r') as f:
             header = [f'Created by {f.attrs["creator"]}',
                       f'        on {f.attrs["created"]}',
@@ -334,7 +344,7 @@ class Result:
 
         >>> import damask
         >>> r = damask.Result('my_file.hdf5')
-        >>> r_first = r.view(increment=0)
+        >>> r_first = r.view(increments=0)
 
         Get a view that shows all results between simulation times of 10 to 40:
 

python/damask/_rotation.py

@@ -88,14 +88,24 @@ class Rotation:
 
 
     def __repr__(self) -> str:
-        """Give short human-readable summary."""
+        """
+        Return repr(self).
+
+        Give short human-readable summary.
+
+        """
         return f'Quaternion{" " if self.quaternion.shape == (4,) else "s of shape "+str(self.quaternion.shape[:-1])+chr(10)}'\
                + str(self.quaternion)
 
 
     def __copy__(self: MyType,
                  rotation: Union[FloatSequence, 'Rotation'] = None) -> MyType:
-        """Create deep copy."""
+        """
+        Return deepcopy(self).
+
+        Create deep copy.
+
+        """
         dup = copy.deepcopy(self)
         if rotation is not None:
             dup.quaternion = Rotation(rotation).quaternion
@@ -106,7 +116,12 @@ class Rotation:
 
     def __getitem__(self,
                     item: Union[Tuple[int], int, bool, np.bool_, np.ndarray]):
-        """Return slice according to item."""
+        """
+        Return self[item].
+
+        Return slice according to item.
+
+        """
         return self.copy() if self.shape == () else \
                self.copy(self.quaternion[item+(slice(None),)] if isinstance(item,tuple) else self.quaternion[item])
 
@@ -114,7 +129,9 @@ class Rotation:
     def __eq__(self,
                other: object) -> bool:
         """
-        Equal to other.
+        Return self==other.
+
+        Test equality of other.
 
         Parameters
         ----------
@@ -130,7 +147,9 @@ class Rotation:
     def __ne__(self,
                other: object) -> bool:
         """
-        Not equal to other.
+        Return self!=other.
+
+        Test inequality of other.
 
         Parameters
         ----------
@@ -214,12 +233,22 @@ class Rotation:
 
 
     def __len__(self) -> int:
-        """Length of leading/leftmost dimension of array."""
+        """
+        Return len(self).
+
+        Length of leading/leftmost dimension of array.
+
+        """
         return 0 if self.shape == () else self.shape[0]
 
 
     def __invert__(self: MyType) -> MyType:
-        """Inverse rotation (backward rotation)."""
+        """
+        Return ~self.
+
+        Inverse rotation (backward rotation).
+
+        """
         dup = self.copy()
         dup.quaternion[...,1:] *= -1
         return dup
@@ -228,6 +257,8 @@ class Rotation:
     def __pow__(self: MyType,
                 exp: Union[float, int]) -> MyType:
         """
+        Return self**exp.
+
         Perform the rotation 'exp' times.
 
         Parameters
@@ -243,6 +274,8 @@ class Rotation:
     def __ipow__(self: MyType,
                  exp: Union[float, int]) -> MyType:
         """
+        Return self**=exp.
+
         Perform the rotation 'exp' times (in-place).
 
         Parameters
@@ -257,6 +290,8 @@ class Rotation:
     def __mul__(self: MyType,
                 other: MyType) -> MyType:
         """
+        Return self*other.
+
         Compose with other.
 
         Parameters
@@ -284,6 +319,8 @@ class Rotation:
     def __imul__(self: MyType,
                  other: MyType) -> MyType:
         """
+        Return self*=other.
+
         Compose with other (in-place).
 
         Parameters
@@ -298,6 +335,8 @@ class Rotation:
     def __truediv__(self: MyType,
                     other: MyType) -> MyType:
         """
+        Return self/other.
+
         Compose with inverse of other.
 
         Parameters
@@ -319,6 +358,8 @@ class Rotation:
     def __itruediv__(self: MyType,
                      other: MyType) -> MyType:
         """
+        Return self/=other.
+
         Compose with inverse of other (in-place).
 
         Parameters
@@ -333,7 +374,9 @@ class Rotation:
     def __matmul__(self,
                    other: np.ndarray) -> np.ndarray:
         """
-        Rotate vector, second order tensor, or fourth order tensor.
+        Return self@other.
+
+        Rotate vector, second-order tensor, or fourth-order tensor.
 
         Parameters
         ----------
@@ -365,7 +408,7 @@ class Rotation:
                 R = self.as_matrix()
                 return np.einsum('...im,...jn,...ko,...lp,...mnop',R,R,R,R,other)
             else:
-                raise ValueError('can only rotate vectors, 2nd order tensors, and 4th order tensors')
+                raise ValueError('can only rotate vectors, second-order tensors, and fourth-order tensors')
         elif isinstance(other, Rotation):
             raise TypeError('use "R1*R2", i.e. multiplication, to compose rotations "R1" and "R2"')
         else:
@@ -1372,7 +1415,7 @@ class Rotation:
         w[np.isclose(w[...,0],1.0+0.0j),1:] = 0.
         w[np.isclose(w[...,1],1.0+0.0j),2:] = 0.
         vr = np.swapaxes(vr,-1,-2)
-        ax = np.where(np.abs(diag_delta)<1e-12,
+        ax = np.where(np.abs(diag_delta)<1e-13,
                              np.real(vr[np.isclose(w,1.0+0.0j)]).reshape(om.shape[:-2]+(3,)),
                       np.abs(np.real(vr[np.isclose(w,1.0+0.0j)]).reshape(om.shape[:-2]+(3,))) \
                       *np.sign(diag_delta))
@@ -1581,14 +1624,13 @@ class Rotation:
     @staticmethod
     def _ho2ax(ho: np.ndarray) -> np.ndarray:
         """Homochoric vector to axis–angle pair."""
-        tfit = np.array([+1.0000000000018852,      -0.5000000002194847,
-                         -0.024999992127593126,    -0.003928701544781374,
-                         -0.0008152701535450438,   -0.0002009500426119712,
-                         -0.00002397986776071756,  -0.00008202868926605841,
-                         +0.00012448715042090092,  -0.0001749114214822577,
-                         +0.0001703481934140054,   -0.00012062065004116828,
-                         +0.000059719705868660826, -0.00001980756723965647,
-                         +0.000003953714684212874, -0.00000036555001439719544])
+        tfit = np.array([+0.9999999999999968,     -0.49999999999986866,     -0.025000000000632055,
+                         -0.003928571496460683,   -0.0008164666077062752,   -0.00019411896443261646,
+                         -0.00004985822229871769, -0.000014164962366386031, -1.9000248160936107e-6,
+                         -5.72184549898506e-6,    +7.772149920658778e-6,    -0.00001053483452909705,
+                         +9.528014229335313e-6,   -5.660288876265125e-6,    +1.2844901692764126e-6,
+                         +1.1255185726258763e-6,  -1.3834391419956455e-6,   +7.513691751164847e-7,
+                         -2.401996891720091e-7,   +4.386887017466388e-8,    -3.5917775353564864e-9])
         hmag_squared = np.sum(ho**2.,axis=-1,keepdims=True)
         s = np.sum(tfit*hmag_squared**np.arange(len(tfit)),axis=-1,keepdims=True)
         with np.errstate(invalid='ignore'):
@@ -1679,7 +1721,7 @@ class Rotation:
 
         """
         with np.errstate(invalid='ignore',divide='ignore'):
-            # get pyramide and scale by grid parameter ratio
+            # get pyramid and scale by grid parameter ratio
             XYZ = np.take_along_axis(cu,Rotation._get_pyramid_order(cu,'forward'),-1) * _sc
             order = np.abs(XYZ[...,1:2]) <= np.abs(XYZ[...,0:1])
             q = np.pi/12.0 * np.where(order,XYZ[...,1:2],XYZ[...,0:1]) \

python/damask/_table.py

@@ -37,7 +37,12 @@ class Table:
 
 
     def __repr__(self) -> str:
-        """Give short human-readable summary."""
+        """
+        Return repr(self).
+
+        Give short human-readable summary.
+
+        """
         self._relabel('shapes')
         data_repr = self.data.__repr__()
         self._relabel('uniform')
@@ -46,7 +51,12 @@ class Table:
 
     def __eq__(self,
                other: object) -> bool:
-        """Compare to other Table."""
+        """
+        Return self==other.
+
+        Test equality of other.
+
+        """
         return NotImplemented if not isinstance(other,Table) else \
                self.shapes == other.shapes and self.data.equals(other.data)
 
@@ -54,7 +64,9 @@ class Table:
     def __getitem__(self,
                     item: Union[slice, Tuple[slice, ...]]) -> 'Table':
         """
-        Slice the Table according to item.
+        Return self[item].
+
+        Return slice according to item.
 
         Parameters
         ----------
@@ -102,12 +114,22 @@ class Table:
 
 
     def __len__(self) -> int:
-        """Number of rows."""
+        """
+        Return len(self).
+
+        Number of rows.
+
+        """
         return len(self.data)
 
 
     def __copy__(self) -> 'Table':
-        """Create deep copy."""
+        """
+        Return deepcopy(self).
+
+        Create deep copy.
+
+        """
         return copy.deepcopy(self)
 
     copy = __copy__
@@ -134,7 +156,7 @@ class Table:
         labels = []
         for label in what:
             shape = self.shapes[label]
-            size = np.prod(shape,dtype=int)
+            size = np.prod(shape,dtype=np.int64)
             if   how == 'uniform':
                 labels += [label] * size
             elif how == 'shapes':
@@ -168,7 +190,7 @@ class Table:
                      shape: Tuple[int, ...],
                      info: str = None):
         if info is not None:
-            specific = f'{label}{" "+str(shape) if np.prod(shape,dtype=int) > 1 else ""}: {info}'
+            specific = f'{label}{" "+str(shape) if np.prod(shape,dtype=np.int64) > 1 else ""}: {info}'
             general  = util.execution_stamp('Table')
             self.comments.append(f'{specific} / {general}')
 
@@ -401,7 +423,7 @@ class Table:
         else:
 
             dup.shapes[label] = data.shape[1:] if len(data.shape) > 1 else (1,)
-            size = np.prod(data.shape[1:],dtype=int)
+            size = np.prod(data.shape[1:],dtype=np.int64)
             new = pd.DataFrame(data=data.reshape(-1,size),
                                columns=[label]*size,
                               )

python/damask/_vtk.py

@@ -39,7 +39,12 @@ class VTK:
 
 
     def __repr__(self) -> str:
-        """Give short human-readable summary."""
+        """
+        Return repr(self).
+
+        Give short human-readable summary.
+
+        """
         info = [self.vtk_data.__vtkname__]
 
         for data in ['Cell Data', 'Point Data']:
@@ -54,7 +59,9 @@ class VTK:
     def __eq__(self,
                other: object) -> bool:
         """
-        Equal to other.
+        Return self==other.
+
+        Test equality of other.
 
         Parameters
         ----------
@@ -187,7 +194,7 @@ class VTK:
         ----------
         nodes : numpy.ndarray, shape (:,3)
             Spatial position of the nodes.
-        connectivity : numpy.ndarray of np.dtype = int
+        connectivity : numpy.ndarray of np.dtype = np.int64
             Cell connectivity (0-based), first dimension determines #Cells,
             second dimension determines #Nodes/Cell.
         cell_type : str

python/damask/seeds.py

@@ -45,7 +45,7 @@ def from_random(size: _FloatSequence,
     else:
         grid_coords = _grid_filters.coordinates0_point(cells,size).reshape(-1,3,order='F')
         coords = grid_coords[rng.choice(_np.prod(cells),N_seeds, replace=False)] \
-               + _np.broadcast_to(size_/_np.array(cells,int),(N_seeds,3))*(rng.random((N_seeds,3))*.5-.25) # wobble w/o leaving grid
+               + _np.broadcast_to(size_/_np.array(cells,_np.int64),(N_seeds,3))*(rng.random((N_seeds,3))*.5-.25) # wobble w/o leaving grid
 
     return coords
 

python/damask/solver/_marc.py

@@ -10,17 +10,25 @@ _damask_root = str(Path(__file__).parents[3])
 class Marc:
     """Wrapper to run DAMASK with MSC Marc."""
 
-    def __init__(self,marc_version=_marc_version,marc_root=_marc_root,damask_root=_damask_root):
+    def __init__(self,
+                 version: str = _marc_version,
+                 marc_root: str = _marc_root,
+                 damask_root: str = _damask_root):
         """
         Create a Marc solver object.
 
         Parameters
         ----------
-        version : string
-            Marc version
+        version : str, optional
+            Marc version. Defaults to latest supported Marc version.
+        marc_root : str, optional
+            Marc root location. Defaults to /opt/msc.
+        damask_root : str, optional
+            DAMASK root location.
+            Default is autodected based on location of the Python library.
 
         """
-        self.marc_version = marc_version
+        self.marc_version = version
         self.marc_root = Path(marc_root)
         self.damask_root = Path(damask_root)
 
@@ -44,9 +52,10 @@ class Marc:
         return path_tools
 
 
-    def submit_job(self, model, job,
-                   compile      = False,
-                   optimization = ''):
+    def submit_job(self, model: str, job: str,
+                   compile: bool = False,
+                   optimization: str = '',
+                   env = None):
         """
         Assemble command line arguments and call Marc executable.
 
@@ -62,6 +71,8 @@ class Marc:
         optimization : str, optional
             Optimization level '' (-O0), 'l' (-O1), or 'h' (-O3).
             Defaults to ''.
+        env : dict, optional
+            Environment for execution.
 
         """
         usersub = (self.damask_root/'src/Marc/DAMASK_Marc').with_suffix('.f90' if compile else '.marc')
@@ -73,18 +84,16 @@ class Marc:
 
         cmd = f'{self.tools_path/script} -jid {model}_{job} -nprocd 1 -autorst 0 -ci n -cr n -dcoup 0 -b no -v no ' \
             + (f'-u {usersub} -save y' if compile else f'-prog {usersub.with_suffix("")}')
-
         print(cmd)
 
-        ret = subprocess.run(shlex.split(cmd),capture_output=True)
+        ret = subprocess.run(shlex.split(cmd),capture_output=True,env=env)
 
-        try:
-            v = int(re.search('Exit number ([0-9]+)',ret.stderr.decode()).group(1))
-            if 3004 != v:
+        if (m := re.search('Exit number ([0-9]+)',ret.stderr.decode())) is not None:
+            if 3004 != (v := int(m.group(1))):
                 print(ret.stderr.decode())
                 print(ret.stdout.decode())
                 raise RuntimeError(f'Marc simulation failed ({v})')
-        except (AttributeError,ValueError):
+        else:
             print(ret.stderr.decode())
             print(ret.stdout.decode())
             raise RuntimeError('Marc simulation failed (unknown return value)')

python/damask/util.py

@@ -1,43 +1,25 @@
 """Miscellaneous helper functionality."""
 
-import sys
-import datetime
-import os
-import subprocess
-import shlex
-import re
-import signal
-import fractions
-from collections import abc
-from functools import reduce, partial
-from typing import Callable, Union, Iterable, Sequence, Dict, List, Tuple, Literal, Any, Collection, TextIO
-from pathlib import Path
+import sys as _sys
+import datetime as _datetime
+import os as _os
+import subprocess as _subprocess
+import shlex as _shlex
+import re as _re
+import signal as _signal
+import fractions as _fractions
+from collections import abc as _abc
+from functools import reduce as _reduce, partial as _partial
+from typing import Callable as _Callable, Union as _Union, Iterable as _Iterable, Sequence as _Sequence, Dict as _Dict,  \
+                   List as _List, Tuple as _Tuple, Literal as _Literal, Any as _Any, Collection as _Collection, TextIO as _TextIO
+from pathlib import Path as _Path
 
-import numpy as np
-import h5py
+import numpy as _np
+import h5py as _h5py
 
-from . import version
-from ._typehints import FloatSequence, NumpyRngSeed, IntCollection, FileHandle
-
-# limit visibility
-__all__=[
-         'srepr',
-         'emph', 'deemph', 'warn', 'strikeout',
-         'run',
-         'open_text',
-         'natural_sort',
-         'show_progress',
-         'scale_to_coprime',
-         'project_equal_angle', 'project_equal_area',
-         'hybrid_IA',
-         'execution_stamp',
-         'shapeshifter', 'shapeblender',
-         'extend_docstring', 'extended_docstring',
-         'Bravais_to_Miller', 'Miller_to_Bravais',
-         'DREAM3D_base_group', 'DREAM3D_cell_data_group',
-         'dict_prune', 'dict_flatten',
-         'tail_repack',
-        ]
+from . import version as _version
+from ._typehints import FloatSequence as _FloatSequence, NumpyRngSeed as _NumpyRngSeed, IntCollection as _IntCollection, \
+                        FileHandle as _FileHandle
 
 # https://svn.blender.org/svnroot/bf-blender/trunk/blender/build_files/scons/tools/bcolors.py
 # https://stackoverflow.com/questions/287871
@@ -154,8 +136,8 @@ def strikeout(msg) -> str:
 
 def run(cmd: str,
         wd: str = './',
-        env: Dict[str, str] = None,
-        timeout: int = None) -> Tuple[str, str]:
+        env: _Dict[str, str] = None,
+        timeout: int = None) -> _Tuple[str, str]:
     """
     Run a command.
 
@@ -178,26 +160,26 @@ def run(cmd: str,
     """
     def pass_signal(sig,_,proc,default):
         proc.send_signal(sig)
-        signal.signal(sig,default)
-        signal.raise_signal(sig)
+        _signal.signal(sig,default)
+        _signal.raise_signal(sig)
 
-    signals = [signal.SIGINT,signal.SIGTERM]
+    signals = [_signal.SIGINT,_signal.SIGTERM]
 
     print(f"running '{cmd}' in '{wd}'")
-    process = subprocess.Popen(shlex.split(cmd),
-                             stdout = subprocess.PIPE,
-                             stderr = subprocess.PIPE,
-                             env = os.environ if env is None else env,
+    process = _subprocess.Popen(_shlex.split(cmd),
+                                stdout = _subprocess.PIPE,
+                                stderr = _subprocess.PIPE,
+                                env = _os.environ if env is None else env,
                                 cwd = wd,
                                 encoding = 'utf-8')
     # ensure that process is terminated (https://stackoverflow.com/questions/22916783)
-    sig_states = [signal.signal(sig,partial(pass_signal,proc=process,default=signal.getsignal(sig))) for sig in signals]
+    sig_states = [_signal.signal(sig,_partial(pass_signal,proc=process,default=_signal.getsignal(sig))) for sig in signals]
 
     try:
         stdout,stderr = process.communicate(timeout=timeout)
     finally:
         for sig,state in zip(signals,sig_states):
-            signal.signal(sig,state)
+            _signal.signal(sig,state)
 
     if process.returncode != 0:
         print(stdout)
@@ -207,8 +189,8 @@ def run(cmd: str,
     return stdout, stderr
 
 
-def open_text(fname: FileHandle,
-              mode: Literal['r','w'] = 'r') -> TextIO:
+def open_text(fname: _FileHandle,
+              mode: _Literal['r','w'] = 'r') -> _TextIO:                                            # noqa
     """
     Open a text file.
 
@@ -224,11 +206,11 @@ def open_text(fname: FileHandle,
     f : file handle
 
     """
-    return fname if not isinstance(fname, (str,Path)) else \
-           open(Path(fname).expanduser(),mode,newline=('\n' if mode == 'w' else None))
+    return fname if not isinstance(fname, (str,_Path)) else \
+           open(_Path(fname).expanduser(),mode,newline=('\n' if mode == 'w' else None))
 
 
-def natural_sort(key: str) -> List[Union[int, str]]:
+def natural_sort(key: str) -> _List[_Union[int, str]]:
     """
     Natural sort.
 
@@ -240,13 +222,13 @@ def natural_sort(key: str) -> List[Union[int, str]]:
 
     """
     convert = lambda text: int(text) if text.isdigit() else text
-    return [ convert(c) for c in re.split('([0-9]+)', key) ]
+    return [ convert(c) for c in _re.split('([0-9]+)', key) ]
 
 
-def show_progress(iterable: Iterable,
+def show_progress(iterable: _Iterable,
                   N_iter: int = None,
                   prefix: str = '',
-                  bar_length: int = 50) -> Any:
+                  bar_length: int = 50) -> _Any:
     """
     Decorate a loop with a progress bar.
 
@@ -264,7 +246,7 @@ def show_progress(iterable: Iterable,
         Length of progress bar in characters. Defaults to 50.
 
     """
-    if isinstance(iterable,abc.Sequence):
+    if isinstance(iterable,_abc.Sequence):
         if N_iter is None:
             N = len(iterable)
         else:
@@ -285,7 +267,7 @@ def show_progress(iterable: Iterable,
             status.update(i)
 
 
-def scale_to_coprime(v: FloatSequence) -> np.ndarray:
+def scale_to_coprime(v: _FloatSequence) -> _np.ndarray:
     """
     Scale vector to co-prime (relatively prime) integers.
 
@@ -304,30 +286,30 @@ def scale_to_coprime(v: FloatSequence) -> np.ndarray:
 
     def get_square_denominator(x):
         """Denominator of the square of a number."""
-        return fractions.Fraction(x ** 2).limit_denominator(MAX_DENOMINATOR).denominator
+        return _fractions.Fraction(x ** 2).limit_denominator(MAX_DENOMINATOR).denominator
 
     def lcm(a,b):
         """Least common multiple."""
         try:
-            return np.lcm(a,b)                                                                      # numpy > 1.18
+            return _np.lcm(a,b)                                                                      # numpy > 1.18
         except AttributeError:
-            return a * b // np.gcd(a, b)
+            return a * b // _np.gcd(a, b)
 
-    v_ = np.array(v)
-    m = (v_ * reduce(lcm, map(lambda x: int(get_square_denominator(x)),v_))**0.5).astype(int)
-    m = m//reduce(np.gcd,m)
+    v_ = _np.array(v)
+    m = (v_ * _reduce(lcm, map(lambda x: int(get_square_denominator(x)),v_))**0.5).astype(_np.int64)
+    m = m//_reduce(_np.gcd,m)
 
-    with np.errstate(invalid='ignore'):
-        if not np.allclose(np.ma.masked_invalid(v_/m),v_[np.argmax(abs(v_))]/m[np.argmax(abs(v_))]):
+    with _np.errstate(invalid='ignore'):
+        if not _np.allclose(_np.ma.masked_invalid(v_/m),v_[_np.argmax(abs(v_))]/m[_np.argmax(abs(v_))]):
             raise ValueError(f'invalid result "{m}" for input "{v_}"')
 
     return m
 
 
-def project_equal_angle(vector: np.ndarray,
-                        direction: Literal['x', 'y', 'z'] = 'z',
+def project_equal_angle(vector: _np.ndarray,
+                        direction: _Literal['x', 'y', 'z'] = 'z',                                   # noqa
                         normalize: bool = True,
-                        keepdims: bool = False) -> np.ndarray:
+                        keepdims: bool = False) -> _np.ndarray:
     """
     Apply equal-angle projection to vector.
 
@@ -367,15 +349,15 @@ def project_equal_angle(vector: np.ndarray,
 
     """
     shift = 'zyx'.index(direction)
-    v = np.roll(vector/np.linalg.norm(vector,axis=-1,keepdims=True) if normalize else vector,
+    v = _np.roll(vector/_np.linalg.norm(vector,axis=-1,keepdims=True) if normalize else vector,
                  shift,axis=-1)
-    return np.roll(np.block([v[...,:2]/(1.0+np.abs(v[...,2:3])),np.zeros_like(v[...,2:3])]),
+    return _np.roll(_np.block([v[...,:2]/(1.0+_np.abs(v[...,2:3])),_np.zeros_like(v[...,2:3])]),
                     -shift if keepdims else 0,axis=-1)[...,:3 if keepdims else 2]
 
-def project_equal_area(vector: np.ndarray,
-                       direction: Literal['x', 'y', 'z'] = 'z',
+def project_equal_area(vector: _np.ndarray,
+                       direction: _Literal['x', 'y', 'z'] = 'z',                                    # noqa
                        normalize: bool = True,
-                       keepdims: bool = False) -> np.ndarray:
+                       keepdims: bool = False) -> _np.ndarray:
     """
     Apply equal-area projection to vector.
 
@@ -416,22 +398,22 @@ def project_equal_area(vector: np.ndarray,
 
     """
     shift = 'zyx'.index(direction)
-    v = np.roll(vector/np.linalg.norm(vector,axis=-1,keepdims=True) if normalize else vector,
+    v = _np.roll(vector/_np.linalg.norm(vector,axis=-1,keepdims=True) if normalize else vector,
                  shift,axis=-1)
-    return np.roll(np.block([v[...,:2]/np.sqrt(1.0+np.abs(v[...,2:3])),np.zeros_like(v[...,2:3])]),
+    return _np.roll(_np.block([v[...,:2]/_np.sqrt(1.0+_np.abs(v[...,2:3])),_np.zeros_like(v[...,2:3])]),
                     -shift if keepdims else 0,axis=-1)[...,:3 if keepdims else 2]
 
 def execution_stamp(class_name: str,
                     function_name: str = None) -> str:
     """Timestamp the execution of a (function within a) class."""
-    now = datetime.datetime.now().astimezone().strftime('%Y-%m-%d %H:%M:%S%z')
+    now = _datetime.datetime.now().astimezone().strftime('%Y-%m-%d %H:%M:%S%z')
     _function_name = '' if function_name is None else f'.{function_name}'
-    return f'damask.{class_name}{_function_name} v{version} ({now})'
+    return f'damask.{class_name}{_function_name} v{_version} ({now})'
 
 
-def hybrid_IA(dist: np.ndarray,
+def hybrid_IA(dist: _np.ndarray,
               N: int,
-              rng_seed: NumpyRngSeed = None) -> np.ndarray:
+              rng_seed: _NumpyRngSeed = None) -> _np.ndarray:
     """
     Hybrid integer approximation.
 
@@ -446,23 +428,23 @@ def hybrid_IA(dist: np.ndarray,
         If None, then fresh, unpredictable entropy will be pulled from the OS.
 
     """
-    N_opt_samples,N_inv_samples = (max(np.count_nonzero(dist),N),0)                                 # random subsampling if too little samples requested
+    N_opt_samples,N_inv_samples = (max(_np.count_nonzero(dist),N),0)                                 # random subsampling if too little samples requested
 
     scale_,scale,inc_factor = (0.0,float(N_opt_samples),1.0)
-    while (not np.isclose(scale, scale_)) and (N_inv_samples != N_opt_samples):
-        repeats = np.rint(scale*dist).astype(np.int64)
-        N_inv_samples = np.sum(repeats)
+    while (not _np.isclose(scale, scale_)) and (N_inv_samples != N_opt_samples):
+        repeats = _np.rint(scale*dist).astype(_np.int64)
+        N_inv_samples = _np.sum(repeats)
         scale_,scale,inc_factor = (scale,scale+inc_factor*0.5*(scale - scale_), inc_factor*2.0) \
                                    if N_inv_samples < N_opt_samples else \
                                   (scale_,0.5*(scale_ + scale), 1.0)
 
-    return np.repeat(np.arange(len(dist)),repeats)[np.random.default_rng(rng_seed).permutation(N_inv_samples)[:N]]
+    return _np.repeat(_np.arange(len(dist)),repeats)[_np.random.default_rng(rng_seed).permutation(N_inv_samples)[:N]]
 
 
-def shapeshifter(fro: Tuple[int, ...],
-                 to: Tuple[int, ...],
-                 mode: Literal['left','right'] = 'left',
-                 keep_ones: bool = False) -> Tuple[int, ...]:
+def shapeshifter(fro: _Tuple[int, ...],
+                 to: _Tuple[int, ...],
+                 mode: _Literal['left','right'] = 'left',                                           # noqa
+                 keep_ones: bool = False) -> _Tuple[int, ...]:
     """
     Return dimensions that reshape 'fro' to become broadcastable to 'to'.
 
@@ -486,8 +468,8 @@ def shapeshifter(fro: Tuple[int, ...],
     new_dims : tuple
         Dimensions for reshape.
 
-    Example
-    -------
+    Examples
+    --------
     >>> import numpy as np
     >>> from damask import util
     >>> a = np.ones((3,4,2))
@@ -496,36 +478,29 @@ def shapeshifter(fro: Tuple[int, ...],
     >>> (a * np.broadcast_to(b_extended,a.shape)).shape
     (3,4,2)
 
-
     """
-    if len(fro) == 0 and len(to) == 0: return ()
+    if len(fro) == 0 and len(to) == 0: return tuple()
+    _fro = [1] if len(fro) == 0 else list(fro)[::-1 if mode=='left' else 1]
+    _to  = [1] if len(to)  == 0 else list(to) [::-1 if mode=='left' else 1]
 
-    beg = dict(left ='(^.*\\b)',
-               right='(^.*?\\b)')
-    sep = dict(left ='(.*\\b)',
-               right='(.*?\\b)')
-    end = dict(left ='(.*?$)',
-               right='(.*$)')
-    fro = (1,) if len(fro) == 0 else fro
-    to  = (1,) if len(to) == 0 else to
-    try:
-        match = re.match(beg[mode]
-                      +f',{sep[mode]}'.join(map(lambda x: f'{x}'
-                                                          if x>1 or (keep_ones and len(fro)>1) else
-                                                          '\\d+',fro))
-                      +f',{end[mode]}',','.join(map(str,to))+',')
-        assert match
-        grp = match.groups()
-    except AssertionError:
-        raise ValueError(f'shapes cannot be shifted {fro} --> {to}')
-    fill: Any = ()
-    for g,d in zip(grp,fro+(None,)):
-        fill += (1,)*g.count(',')+(d,)
-    return fill[:-1]
+    final_shape: _List[int] = []
+    index = 0
+    for i,item in enumerate(_to):
+        if item==_fro[index]:
+            final_shape.append(item)
+            index+=1
+        else:
+            final_shape.append(1)
+            if _fro[index]==1 and not keep_ones:
+                index+=1
+        if index==len(_fro):
+            final_shape = final_shape+[1]*(len(_to)-i-1)
+            break
+    if index!=len(_fro): raise ValueError(f'shapes cannot be shifted {fro} --> {to}')
+    return tuple(final_shape[::-1] if mode=='left' else final_shape)
 
-
-def shapeblender(a: Tuple[int, ...],
-                 b: Tuple[int, ...]) -> Tuple[int, ...]:
+def shapeblender(a: _Tuple[int, ...],
+                 b: _Tuple[int, ...]) -> _Tuple[int, ...]:
     """
     Return a shape that overlaps the rightmost entries of 'a' with the leftmost of 'b'.
 
@@ -553,7 +528,7 @@ def shapeblender(a: Tuple[int, ...],
     return a + b[i:]
 
 
-def extend_docstring(extra_docstring: str) -> Callable:
+def extend_docstring(extra_docstring: str) -> _Callable:
     """
     Decorator: Append to function's docstring.
 
@@ -569,8 +544,8 @@ def extend_docstring(extra_docstring: str) -> Callable:
     return _decorator
 
 
-def extended_docstring(f: Callable,
-                       extra_docstring: str) -> Callable:
+def extended_docstring(f: _Callable,
+                       extra_docstring: str) -> _Callable:
     """
     Decorator: Combine another function's docstring with a given docstring.
 
@@ -588,7 +563,7 @@ def extended_docstring(f: Callable,
     return _decorator
 
 
-def DREAM3D_base_group(fname: Union[str, Path]) -> str:
+def DREAM3D_base_group(fname: _Union[str, _Path]) -> str:
     """
     Determine the base group of a DREAM.3D file.
 
@@ -606,7 +581,7 @@ def DREAM3D_base_group(fname: Union[str, Path]) -> str:
         Path to the base group.
 
     """
-    with h5py.File(Path(fname).expanduser(),'r') as f:
+    with _h5py.File(_Path(fname).expanduser(),'r') as f:
         base_group = f.visit(lambda path: path.rsplit('/',2)[0] if '_SIMPL_GEOMETRY/SPACING' in path else None)
 
     if base_group is None:
@@ -614,7 +589,7 @@ def DREAM3D_base_group(fname: Union[str, Path]) -> str:
 
     return base_group
 
-def DREAM3D_cell_data_group(fname: Union[str, Path]) -> str:
+def DREAM3D_cell_data_group(fname: _Union[str, _Path]) -> str:
     """
     Determine the cell data group of a DREAM.3D file.
 
@@ -634,10 +609,10 @@ def DREAM3D_cell_data_group(fname: Union[str, Path]) -> str:
 
     """
     base_group = DREAM3D_base_group(fname)
-    with h5py.File(Path(fname).expanduser(),'r') as f:
+    with _h5py.File(_Path(fname).expanduser(),'r') as f:
         cells = tuple(f['/'.join([base_group,'_SIMPL_GEOMETRY','DIMENSIONS'])][()][::-1])
         cell_data_group = f[base_group].visititems(lambda path,obj: path.split('/')[0] \
-                                                   if isinstance(obj,h5py._hl.dataset.Dataset) and np.shape(obj)[:-1] == cells \
+                                                   if isinstance(obj,_h5py._hl.dataset.Dataset) and _np.shape(obj)[:-1] == cells \
                                                    else None)
 
     if cell_data_group is None:
@@ -647,8 +622,8 @@ def DREAM3D_cell_data_group(fname: Union[str, Path]) -> str:
 
 
 def Bravais_to_Miller(*,
-                      uvtw: np.ndarray = None,
-                      hkil: np.ndarray = None) -> np.ndarray:
+                      uvtw: _np.ndarray = None,
+                      hkil: _np.ndarray = None) -> _np.ndarray:
     """
     Transform 4 Miller–Bravais indices to 3 Miller indices of crystal direction [uvw] or plane normal (hkl).
 
@@ -665,19 +640,19 @@ def Bravais_to_Miller(*,
     """
     if (uvtw is not None) ^ (hkil is None):
         raise KeyError('specify either "uvtw" or "hkil"')
-    axis,basis  = (np.array(uvtw),np.array([[1,0,-1,0],
+    axis,basis  = (_np.array(uvtw),_np.array([[1,0,-1,0],
                                               [0,1,-1,0],
                                               [0,0, 0,1]])) \
                   if hkil is None else \
-                  (np.array(hkil),np.array([[1,0,0,0],
+                  (_np.array(hkil),_np.array([[1,0,0,0],
                                               [0,1,0,0],
                                               [0,0,0,1]]))
-    return np.einsum('il,...l',basis,axis)
+    return _np.einsum('il,...l',basis,axis)
 
 
 def Miller_to_Bravais(*,
-                      uvw: np.ndarray = None,
-                      hkl: np.ndarray = None) -> np.ndarray:
+                      uvw: _np.ndarray = None,
+                      hkl: _np.ndarray = None) -> _np.ndarray:
     """
     Transform 3 Miller indices to 4 Miller–Bravais indices of crystal direction [uvtw] or plane normal (hkil).
 
@@ -694,19 +669,19 @@ def Miller_to_Bravais(*,
     """
     if (uvw is not None) ^ (hkl is None):
         raise KeyError('specify either "uvw" or "hkl"')
-    axis,basis  = (np.array(uvw),np.array([[ 2,-1, 0],
+    axis,basis  = (_np.array(uvw),_np.array([[ 2,-1, 0],
                                              [-1, 2, 0],
                                              [-1,-1, 0],
                                              [ 0, 0, 3]])/3) \
                   if hkl is None else \
-                  (np.array(hkl),np.array([[ 1, 0, 0],
+                  (_np.array(hkl),_np.array([[ 1, 0, 0],
                                              [ 0, 1, 0],
                                              [-1,-1, 0],
                                              [ 0, 0, 1]]))
-    return np.einsum('il,...l',basis,axis)
+    return _np.einsum('il,...l',basis,axis)
 
 
-def dict_prune(d: Dict) -> Dict:
+def dict_prune(d: _Dict) -> _Dict:
     """
     Recursively remove empty dictionaries.
 
@@ -732,7 +707,7 @@ def dict_prune(d: Dict) -> Dict:
     return new
 
 
-def dict_flatten(d: Dict) -> Dict:
+def dict_flatten(d: _Dict) -> _Dict:
     """
     Recursively remove keys of single-entry dictionaries.
 
@@ -756,8 +731,8 @@ def dict_flatten(d: Dict) -> Dict:
     return new
 
 
-def tail_repack(extended: Union[str, Sequence[str]],
-                existing: List[str] = []) -> List[str]:
+def tail_repack(extended: _Union[str, _Sequence[str]],
+                existing: _List[str] = []) -> _List[str]:
     """
     Repack tailing characters into single string if all are new.
 
@@ -782,11 +757,11 @@ def tail_repack(extended: Union[str, Sequence[str]],
 
     """
     return [extended] if isinstance(extended,str) else existing + \
-         ([''.join(extended[len(existing):])] if np.prod([len(i) for i in extended[len(existing):]]) == 1 else
+         ([''.join(extended[len(existing):])] if _np.prod([len(i) for i in extended[len(existing):]]) == 1 else
           list(extended[len(existing):]))
 
 
-def aslist(arg: Union[IntCollection,int,None]) -> List:
+def aslist(arg: _Union[_IntCollection, int, None]) -> _List:
     """
     Transform argument to list.
 
@@ -801,7 +776,7 @@ def aslist(arg: Union[IntCollection,int,None]) -> List:
         Entity transformed into list.
 
     """
-    return [] if arg is None else list(arg) if isinstance(arg,(np.ndarray,Collection)) else [arg]
+    return [] if arg is None else list(arg) if isinstance(arg,(_np.ndarray,_Collection)) else [arg]
 
 
 ####################################################################################################
@@ -834,11 +809,11 @@ class ProgressBar:
         self.total = total
         self.prefix = prefix
         self.bar_length = bar_length
-        self.time_start = self.time_last_update = datetime.datetime.now()
+        self.time_start = self.time_last_update = _datetime.datetime.now()
         self.fraction_last = 0.0
 
-        sys.stderr.write(f"{self.prefix} {'░'*self.bar_length}   0% ETA n/a")
-        sys.stderr.flush()
+        _sys.stderr.write(f"{self.prefix} {'░'*self.bar_length}   0% ETA n/a")
+        _sys.stderr.flush()
 
     def update(self,
                iteration: int) -> None:
@@ -846,17 +821,17 @@ class ProgressBar:
         fraction = (iteration+1) / self.total
 
         if (filled_length := int(self.bar_length * fraction)) > int(self.bar_length * self.fraction_last) or \
-            datetime.datetime.now() - self.time_last_update > datetime.timedelta(seconds=10):
-            self.time_last_update = datetime.datetime.now()
+            _datetime.datetime.now() - self.time_last_update > _datetime.timedelta(seconds=10):
+            self.time_last_update = _datetime.datetime.now()
             bar = '█' * filled_length + '░' * (self.bar_length - filled_length)
-            remaining_time = (datetime.datetime.now() - self.time_start) \
+            remaining_time = (_datetime.datetime.now() - self.time_start) \
                            * (self.total - (iteration+1)) / (iteration+1)
-            remaining_time -= datetime.timedelta(microseconds=remaining_time.microseconds)          # remove μs
-            sys.stderr.write(f'\r{self.prefix} {bar} {fraction:>4.0%} ETA {remaining_time}')
-            sys.stderr.flush()
+            remaining_time -= _datetime.timedelta(microseconds=remaining_time.microseconds)          # remove μs
+            _sys.stderr.write(f'\r{self.prefix} {bar} {fraction:>4.0%} ETA {remaining_time}')
+            _sys.stderr.flush()
 
         self.fraction_last = fraction
 
         if iteration == self.total - 1:
-            sys.stderr.write('\n')
-            sys.stderr.flush()
+            _sys.stderr.write('\n')
+            _sys.stderr.flush()

python/tests/test_Orientation.py

@@ -224,11 +224,11 @@ class TestOrientation:
 
     @pytest.mark.parametrize('family',crystal_families)
     def test_reduced_corner_cases(self,family):
-        # test whether there is always a sym-eq rotation that falls into the FZ
+        # test whether there is always exactly one sym-eq rotation that falls into the FZ
         N = np.random.randint(10,40)
         size = np.ones(3)*np.pi**(2./3.)
         grid = grid_filters.coordinates0_node([N+1,N+1,N+1],size,-size*.5)
-        evenly_distributed = Orientation.from_cubochoric(x=grid[:-2,:-2,:-2],family=family)
+        evenly_distributed = Orientation.from_cubochoric(x=grid,family=family)
         assert evenly_distributed.shape == evenly_distributed.reduced.shape
 
     @pytest.mark.parametrize('family',crystal_families)

python/tests/test_Rotation.py

@@ -301,14 +301,13 @@ def ro2ho(ro):
 #---------- Homochoric vector----------
 def ho2ax(ho):
     """Homochoric vector to axis angle pair."""
-    tfit = np.array([+1.0000000000018852,      -0.5000000002194847,
-                     -0.024999992127593126,    -0.003928701544781374,
-                     -0.0008152701535450438,   -0.0002009500426119712,
-                     -0.00002397986776071756,  -0.00008202868926605841,
-                     +0.00012448715042090092,  -0.0001749114214822577,
-                     +0.0001703481934140054,   -0.00012062065004116828,
-                     +0.000059719705868660826, -0.00001980756723965647,
-                     +0.000003953714684212874, -0.00000036555001439719544])
+    tfit = np.array([+0.9999999999999968,     -0.49999999999986866,     -0.025000000000632055,
+                     -0.003928571496460683,   -0.0008164666077062752,   -0.00019411896443261646,
+                     -0.00004985822229871769, -0.000014164962366386031, -1.9000248160936107e-6,
+                     -5.72184549898506e-6,    +7.772149920658778e-6,    -0.00001053483452909705,
+                     +9.528014229335313e-6,   -5.660288876265125e-6,    +1.2844901692764126e-6,
+                     +1.1255185726258763e-6,  -1.3834391419956455e-6,   +7.513691751164847e-7,
+                     -2.401996891720091e-7,   +4.386887017466388e-8,    -3.5917775353564864e-9])
     # normalize h and store the magnitude
     hmag_squared = np.sum(ho**2.)
     if iszero(hmag_squared):

src/CLI.f90

@@ -18,7 +18,7 @@ module CLI
   use parallelization
   use system_routines
 
-  implicit none
+  implicit none(type,external)
   private
   integer,                       public, protected :: &
     CLI_restartInc = 0                                                                              !< Increment at which calculation starts

src/HDF5_utilities.f90

@@ -18,7 +18,11 @@ module HDF5_utilities
   use prec
   use parallelization
 
+#if (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR>14) && !defined(PETSC_HAVE_MPI_F90MODULE_VISIBILITY)
+  implicit none(type,external)
+#else
   implicit none
+#endif
   private
 
 !--------------------------------------------------------------------------------------------------

src/IO.f90

@@ -12,7 +12,7 @@ module IO
 
   use prec
 
-  implicit none
+  implicit none(type,external)
   private
 
   character(len=*), parameter, public :: &
@@ -24,11 +24,6 @@ module IO
   character, parameter :: &
     CR = achar(13), &
     LF = IO_EOL
-  character(len=*), parameter :: &
-    IO_DIVIDER = '───────────────────'//&
-                 '───────────────────'//&
-                 '───────────────────'//&
-                 '────────────'
 
   public :: &
     IO_init, &
@@ -54,11 +49,11 @@ contains
 !--------------------------------------------------------------------------------------------------
 !> @brief Do self test.
 !--------------------------------------------------------------------------------------------------
-subroutine IO_init
+subroutine IO_init()
 
   print'(/,1x,a)', '<<<+-  IO init  -+>>>'; flush(IO_STDOUT)
 
-  call selfTest
+  call selfTest()
 
 end subroutine IO_init
 
@@ -95,7 +90,7 @@ function IO_readlines(fileName) result(fileContent)
     if (endPos - startPos > pStringLen-1) then
       line = rawData(startPos:startPos+pStringLen-1)
       if (.not. warned) then
-        call IO_warning(207,ext_msg=trim(fileName),el=l)
+        call IO_warning(207,trim(fileName),label1='line',ID1=l)
         warned = .true.
       end if
     else
@@ -129,7 +124,7 @@ function IO_read(fileName) result(fileContent)
   inquire(file = fileName, size=fileLength)
   open(newunit=fileUnit, file=fileName, access='stream',&
        status='old', position='rewind', action='read',iostat=myStat)
-  if (myStat /= 0) call IO_error(100,ext_msg=trim(fileName))
+  if (myStat /= 0) call IO_error(100,trim(fileName))
   allocate(character(len=fileLength)::fileContent)
   if (fileLength==0) then
     close(fileUnit)
@@ -137,7 +132,7 @@ function IO_read(fileName) result(fileContent)
   end if
 
   read(fileUnit,iostat=myStat) fileContent
-  if (myStat /= 0) call IO_error(102,ext_msg=trim(fileName))
+  if (myStat /= 0) call IO_error(102,trim(fileName))
   close(fileUnit)
 
   if (scan(fileContent(:index(fileContent,LF)),CR//LF) /= 0) fileContent = CRLF2LF(fileContent)
@@ -206,7 +201,7 @@ function IO_stringValue(string,chunkPos,myChunk)
 
   validChunk: if (myChunk > chunkPos(1) .or. myChunk < 1) then
     IO_stringValue = ''
-    call IO_error(110,el=myChunk,ext_msg='IO_stringValue: "'//trim(string)//'"')
+    call IO_error(110,'IO_stringValue: "'//trim(string)//'"',label1='chunk',ID1=myChunk)
   else validChunk
     IO_stringValue = string(chunkPos(myChunk*2):chunkPos(myChunk*2+1))
   end if validChunk
@@ -303,10 +298,10 @@ integer function IO_stringAsInt(string)
 
   valid: if (verify(string,VALIDCHARS) == 0) then
     read(string,*,iostat=readStatus) IO_stringAsInt
-    if (readStatus /= 0) call IO_error(111,ext_msg=string)
+    if (readStatus /= 0) call IO_error(111,string)
   else valid
     IO_stringAsInt = 0
-    call IO_error(111,ext_msg=string)
+    call IO_error(111,string)
   end if valid
 
 end function IO_stringAsInt
@@ -325,10 +320,10 @@ real(pReal) function IO_stringAsFloat(string)
 
   valid: if (verify(string,VALIDCHARS) == 0) then
     read(string,*,iostat=readStatus) IO_stringAsFloat
-    if (readStatus /= 0) call IO_error(112,ext_msg=string)
+    if (readStatus /= 0) call IO_error(112,string)
   else valid
     IO_stringAsFloat = 0.0_pReal
-    call IO_error(112,ext_msg=string)
+    call IO_error(112,string)
   end if valid
 
 end function IO_stringAsFloat
@@ -348,33 +343,27 @@ logical function IO_stringAsBool(string)
     IO_stringAsBool = .false.
   else
     IO_stringAsBool = .false.
-    call IO_error(113,ext_msg=string)
+    call IO_error(113,string)
   end if
 
 end function IO_stringAsBool
 
 
+
 !--------------------------------------------------------------------------------------------------
-!> @brief Write error statements to standard out and terminate the run with exit #9xxx
+!> @brief Write error statements and terminate the run with exit #9xxx.
 !--------------------------------------------------------------------------------------------------
-subroutine IO_error(error_ID,el,ip,g,instance,ext_msg)
+subroutine IO_error(error_ID,ext_msg,label1,ID1,label2,ID2)
 
   integer,                    intent(in) :: error_ID
-  integer,          optional, intent(in) :: el,ip,g,instance
-  character(len=*), optional, intent(in) :: ext_msg
+  character(len=*), optional, intent(in) :: ext_msg,label1,label2
+  integer,          optional, intent(in) :: ID1,ID2
 
   external                      :: quit
   character(len=:), allocatable :: msg
-  character(len=pStringLen)     :: formatString
-
 
   select case (error_ID)
 
-!--------------------------------------------------------------------------------------------------
-! internal errors
-    case (0)
-      msg = 'internal check failed:'
-
 !--------------------------------------------------------------------------------------------------
 ! file handling errors
     case (100)
@@ -446,7 +435,7 @@ subroutine IO_error(error_ID,el,ip,g,instance,ext_msg)
     case (190)
       msg = 'unknown element type:'
     case (191)
-      msg = 'mesh consists of more than one element type'
+      msg = 'mesh contains more than one element type'
 
 !--------------------------------------------------------------------------------------------------
 ! plasticity error messages
@@ -483,27 +472,27 @@ subroutine IO_error(error_ID,el,ip,g,instance,ext_msg)
 !------------------------------------------------------------------------------------------------
 ! errors related to YAML data
     case (701)
-      msg = 'Incorrect indent/Null value not allowed'
+      msg = 'incorrect indent/Null value not allowed'
     case (702)
-      msg = 'Invalid use of flow YAML'
+      msg = 'invalid use of flow YAML'
     case (703)
-      msg = 'Invalid YAML'
+      msg = 'invalid YAML'
     case (704)
-      msg = 'Space expected after a colon for <key>: <value> pair'
+      msg = 'space expected after a colon for <key>: <value> pair'
     case (705)
-      msg = 'Unsupported feature'
+      msg = 'unsupported feature'
     case (706)
-      msg = 'Type mismatch in YAML data node'
+      msg = 'type mismatch in YAML data node'
     case (707)
-      msg = 'Abrupt end of file'
+      msg = 'abrupt end of file'
     case (708)
-      msg = '--- expected after YAML file header'
+      msg = '"---" expected after YAML file header'
     case (709)
-      msg = 'Length mismatch'
+      msg = 'length mismatch'
     case (710)
-      msg = 'Closing quotation mark missing in string'
+      msg = 'closing quotation mark missing in string'
     case (711)
-      msg = 'Incorrect type'
+      msg = 'incorrect type'
 
 !-------------------------------------------------------------------------------------------------
 ! errors related to the mesh solver
@@ -540,58 +529,35 @@ subroutine IO_error(error_ID,el,ip,g,instance,ext_msg)
     case (950)
       msg = 'max number of cut back exceeded, terminating'
 
-!-------------------------------------------------------------------------------------------------
-! general error messages
     case default
-      msg = 'unknown error number...'
+      error stop 'invalid error number'
 
   end select
 
-  !$OMP CRITICAL (write2out)
-  write(IO_STDERR,'(/,a)')                ' ┌'//IO_DIVIDER//'┐'
-  write(IO_STDERR,'(a,24x,a,40x,a)')      ' │','error',                                             '│'
-  write(IO_STDERR,'(a,24x,i3,42x,a)')     ' │',error_ID,                                            '│'
-  write(IO_STDERR,'(a)')                  ' ├'//IO_DIVIDER//'┤'
-  write(formatString,'(a,i6.6,a,i6.6,a)') '(1x,a4,a',max(1,len_trim(msg)),',',&
-                                                     max(1,72-len_trim(msg)-4),'x,a)'
-  write(IO_STDERR,formatString)            '│ ',trim(msg),                                          '│'
-  if (present(ext_msg)) then
-    write(formatString,'(a,i6.6,a,i6.6,a)') '(1x,a4,a',max(1,len_trim(ext_msg)),',',&
-                                                       max(1,72-len_trim(ext_msg)-4),'x,a)'
-    write(IO_STDERR,formatString)          '│ ',trim(ext_msg),                                      '│'
-  endif
-  if (present(el)) &
-    write(IO_STDERR,'(a19,1x,i9,44x,a3)') ' │ at element    ',el,                                   '│'
-  if (present(ip)) &
-    write(IO_STDERR,'(a19,1x,i9,44x,a3)') ' │ at IP         ',ip,                                   '│'
-  if (present(g)) &
-    write(IO_STDERR,'(a19,1x,i9,44x,a3)') ' │ at constituent',g,                                    '│'
-  if (present(instance)) &
-    write(IO_STDERR,'(a19,1x,i9,44x,a3)') ' │ at instance   ',instance,                             '│'
-  write(IO_STDERR,'(a,69x,a)')            ' │',                                                     '│'
-  write(IO_STDERR,'(a)')                  ' └'//IO_DIVIDER//'┘'
-  flush(IO_STDERR)
+  call panel('error',error_ID,msg, &
+                 ext_msg=ext_msg, &
+                 label1=label1,ID1=ID1, &
+                 label2=label2,ID2=ID2)
   call quit(9000+error_ID)
-  !$OMP END CRITICAL (write2out)
 
 end subroutine IO_error
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief Write warning statement to standard out.
+!> @brief Write warning statements.
 !--------------------------------------------------------------------------------------------------
-subroutine IO_warning(warning_ID,el,ip,g,ext_msg)
+subroutine IO_warning(warning_ID,ext_msg,label1,ID1,label2,ID2)
 
   integer,                    intent(in) :: warning_ID
-  integer,          optional, intent(in) :: el,ip,g
-  character(len=*), optional, intent(in) :: ext_msg
+  character(len=*), optional, intent(in) :: ext_msg,label1,label2
+  integer,          optional, intent(in) :: ID1,ID2
 
   character(len=:), allocatable :: msg
-  character(len=pStringLen)     :: formatString
+
 
   select case (warning_ID)
     case (47)
-      msg = 'no valid parameter for FFTW, using FFTW_PATIENT'
+      msg = 'invalid parameter for FFTW'
     case (207)
       msg = 'line truncated'
     case (600)
@@ -600,33 +566,15 @@ subroutine IO_warning(warning_ID,el,ip,g,ext_msg)
       msg = 'stiffness close to zero'
     case (709)
       msg = 'read only the first document'
+
     case default
-      msg = 'unknown warning number'
+      error stop 'invalid warning number'
   end select
 
-  !$OMP CRITICAL (write2out)
-  write(IO_STDERR,'(/,a)')                ' ┌'//IO_DIVIDER//'┐'
-  write(IO_STDERR,'(a,24x,a,38x,a)')      ' │','warning',                                           '│'
-  write(IO_STDERR,'(a,24x,i3,42x,a)')     ' │',warning_ID,                                          '│'
-  write(IO_STDERR,'(a)')                  ' ├'//IO_DIVIDER//'┤'
-  write(formatString,'(a,i6.6,a,i6.6,a)') '(1x,a4,a',max(1,len_trim(msg)),',',&
-                                                     max(1,72-len_trim(msg)-4),'x,a)'
-  write(IO_STDERR,formatString)            '│ ',trim(msg),                                          '│'
-  if (present(ext_msg)) then
-    write(formatString,'(a,i6.6,a,i6.6,a)') '(1x,a4,a',max(1,len_trim(ext_msg)),',',&
-                                                       max(1,72-len_trim(ext_msg)-4),'x,a)'
-    write(IO_STDERR,formatString)          '│ ',trim(ext_msg),                                      '│'
-  endif
-  if (present(el)) &
-    write(IO_STDERR,'(a19,1x,i9,44x,a3)') ' │ at element    ',el,                                   '│'
-  if (present(ip)) &
-    write(IO_STDERR,'(a19,1x,i9,44x,a3)') ' │ at IP         ',ip,                                   '│'
-  if (present(g)) &
-    write(IO_STDERR,'(a19,1x,i9,44x,a3)') ' │ at constituent',g,                                    '│'
-  write(IO_STDERR,'(a,69x,a)')            ' │',                                                     '│'
-  write(IO_STDERR,'(a)')                  ' └'//IO_DIVIDER//'┘'
-  flush(IO_STDERR)
-  !$OMP END CRITICAL (write2out)
+  call panel('warning',warning_ID,msg, &
+                 ext_msg=ext_msg, &
+                 label1=label1,ID1=ID1, &
+                 label2=label2,ID2=ID2)
 
 end subroutine IO_warning
 
@@ -654,7 +602,61 @@ pure function CRLF2LF(string)
 
   CRLF2LF = CRLF2LF(:c-n)
 
-end function
+end function CRLF2LF
+
+
+!--------------------------------------------------------------------------------------------------
+!> @brief Write statements to standard error.
+!--------------------------------------------------------------------------------------------------
+subroutine panel(paneltype,ID,msg,ext_msg,label1,ID1,label2,ID2)
+
+  character(len=*),           intent(in) :: paneltype,msg
+  character(len=*), optional, intent(in) :: ext_msg,label1,label2
+  integer,                    intent(in) :: ID
+  integer,          optional, intent(in) :: ID1,ID2
+
+  character(len=pStringLen)              :: formatString
+  integer, parameter                     :: panelwidth = 69
+  character(len=*), parameter            :: DIVIDER = repeat('─',panelwidth)
+
+
+  if (.not. present(label1) .and.       present(ID1)) error stop 'missing label for value 1'
+  if (.not. present(label2) .and.       present(ID2)) error stop 'missing label for value 2'
+  if (      present(label1) .and. .not. present(ID1)) error stop 'missing value for label 1'
+  if (      present(label2) .and. .not. present(ID2)) error stop 'missing value for label 2'
+
+  !$OMP CRITICAL (write2out)
+  write(IO_STDERR,'(/,a)')                ' ┌'//DIVIDER//'┐'
+  write(formatString,'(a,i2,a)') '(a,24x,a,',max(1,panelwidth-24-len_trim(paneltype)),'x,a)'
+  write(IO_STDERR,formatString)          ' │',trim(paneltype),                                      '│'
+  write(formatString,'(a,i2,a)') '(a,24x,i3,',max(1,panelwidth-24-3),'x,a)'
+  write(IO_STDERR,formatString)          ' │',ID,                                                   '│'
+  write(IO_STDERR,'(a)')                  ' ├'//DIVIDER//'┤'
+  write(formatString,'(a,i3.3,a,i3.3,a)') '(1x,a4,a',max(1,len_trim(msg)),',',&
+                                                     max(1,panelwidth+3-len_trim(msg)-4),'x,a)'
+  write(IO_STDERR,formatString)            '│ ',trim(msg),                                          '│'
+  if (present(ext_msg)) then
+    write(formatString,'(a,i3.3,a,i3.3,a)') '(1x,a4,a',max(1,len_trim(ext_msg)),',',&
+                                                       max(1,panelwidth+3-len_trim(ext_msg)-4),'x,a)'
+    write(IO_STDERR,formatString)          '│ ',trim(ext_msg),                                      '│'
+  end if
+  if (present(label1)) then
+    write(formatString,'(a,i3.3,a,i3.3,a)') '(1x,a7,a',max(1,len_trim(label1)),',i9,',&
+                                                       max(1,panelwidth+3-len_trim(label1)-9-7),'x,a)'
+    write(IO_STDERR,formatString)          '│ at ',trim(label1),ID1,                                '│'
+  end if
+  if (present(label2)) then
+    write(formatString,'(a,i3.3,a,i3.3,a)') '(1x,a7,a',max(1,len_trim(label2)),',i9,',&
+                                                       max(1,panelwidth+3-len_trim(label2)-9-7),'x,a)'
+    write(IO_STDERR,formatString)          '│ at ',trim(label2),ID2,                                '│'
+  end if
+  write(formatString,'(a,i2.2,a)') '(a,',max(1,panelwidth),'x,a)'
+  write(IO_STDERR,formatString)          ' │',                                                     '│'
+  write(IO_STDERR,'(a)')                  ' └'//DIVIDER//'┘'
+  flush(IO_STDERR)
+  !$OMP END CRITICAL (write2out)
+
+end subroutine panel
 
 
 !--------------------------------------------------------------------------------------------------
@@ -665,6 +667,7 @@ subroutine selfTest()
   integer, dimension(:), allocatable :: chunkPos
   character(len=:),      allocatable :: str
 
+
   if(dNeq(1.0_pReal, IO_stringAsFloat('1.0')))      error stop 'IO_stringAsFloat'
   if(dNeq(1.0_pReal, IO_stringAsFloat('1e0')))      error stop 'IO_stringAsFloat'
   if(dNeq(0.1_pReal, IO_stringAsFloat('1e-1')))     error stop 'IO_stringAsFloat'

src/LAPACK_interface.f90

@@ -8,6 +8,8 @@ module LAPACK_interface
 
     pure subroutine dgeev(jobvl,jobvr,n,a,lda,wr,wi,vl,ldvl,vr,ldvr,work,lwork,info)
       use prec
+      implicit none(type,external)
+
       character,   intent(in)                             :: jobvl,jobvr
       integer,     intent(in)                             :: n,lda,ldvl,ldvr,lwork
       real(pReal), intent(inout), dimension(lda,n)        :: a
@@ -20,6 +22,8 @@ module LAPACK_interface
 
     pure subroutine dgesv(n,nrhs,a,lda,ipiv,b,ldb,info)
       use prec
+      implicit none(type,external)
+
       integer,     intent(in)                             :: n,nrhs,lda,ldb
       real(pReal), intent(inout), dimension(lda,n)        :: a
       integer,     intent(out),   dimension(n)            :: ipiv
@@ -29,6 +33,8 @@ module LAPACK_interface
 
     pure subroutine dgetrf(m,n,a,lda,ipiv,info)
       use prec
+      implicit none(type,external)
+
       integer,     intent(in)                             :: m,n,lda
       real(pReal), intent(inout), dimension(lda,n)        :: a
       integer,     intent(out),   dimension(min(m,n))     :: ipiv
@@ -37,6 +43,8 @@ module LAPACK_interface
 
     pure subroutine dgetri(n,a,lda,ipiv,work,lwork,info)
       use prec
+      implicit none(type,external)
+
       integer,     intent(in)                             :: n,lda,lwork
       real(pReal), intent(inout), dimension(lda,n)        :: a
       integer,     intent(in),    dimension(n)            :: ipiv
@@ -46,6 +54,8 @@ module LAPACK_interface
 
     pure subroutine dsyev(jobz,uplo,n,a,lda,w,work,lwork,info)
       use prec
+      implicit none(type,external)
+
       character,   intent(in)                             :: jobz,uplo
       integer,     intent(in)                             :: n,lda,lwork
       real(pReal), intent(inout), dimension(lda,n)        :: a

src/Marc/DAMASK_Marc.f90

@@ -25,7 +25,7 @@ module DAMASK_interface
   use ifport, only: &
     CHDIR
 
-  implicit none
+  implicit none(type,external)
   private
 
   logical,          protected, public :: symmetricSolver
@@ -210,8 +210,8 @@ subroutine hypela2(d,g,e,de,s,t,dt,ngens,m,nn,kcus,matus,ndi,nshear,disp, &
   use materialpoint_Marc
   use OMP_LIB
 
-  implicit none
-  integer,                               intent(in) :: &                                            ! according to MSC.Marc 2012 Manual D
+  implicit none(type,external)
+  integer(pI64),                         intent(in) :: &                                            ! according to MSC.Marc 2012 Manual D
     ngens, &                                                                                        !< size of stress-strain law
     nn, &                                                                                           !< integration point number
     ndi, &                                                                                          !< number of direct components
@@ -224,7 +224,7 @@ subroutine hypela2(d,g,e,de,s,t,dt,ngens,m,nn,kcus,matus,ndi,nshear,disp, &
     jtype, &                                                                                        !< element type
     ifr, &                                                                                          !< set to 1 if R has been calculated
     ifu                                                                                             !< set to 1 if stretch has been calculated
-  integer, dimension(2),                 intent(in) :: &                                            ! according to MSC.Marc 2012 Manual D
+  integer(pI64), dimension(2),           intent(in) :: &                                            ! according to MSC.Marc 2012 Manual D
     m, &                                                                                            !< (1) user element number, (2) internal element number
     matus, &                                                                                        !< (1) user material identification number, (2) internal material identification number
     kcus, &                                                                                         !< (1) layer number, (2) internal layer number
@@ -362,7 +362,7 @@ subroutine hypela2(d,g,e,de,s,t,dt,ngens,m,nn,kcus,matus,ndi,nshear,disp, &
   endif
   lastLovl = lovl
 
-  call materialpoint_general(computationMode,ffn,ffn1,t(1),timinc,m(1),nn,stress,ddsdde)
+  call materialpoint_general(computationMode,ffn,ffn1,t(1),timinc,int(m(1)),int(nn),stress,ddsdde)
 
   d = ddsdde(1:ngens,1:ngens)
   s = stress(1:ndi+nshear)
@@ -382,17 +382,18 @@ subroutine flux(f,ts,n,time)
   use homogenization
   use discretization_Marc
 
-  implicit none
+  implicit none(type,external)
   real(pReal),   dimension(6),  intent(in) :: &
     ts
-  integer,     dimension(10),          intent(in) :: &
+  integer(pI64), dimension(10), intent(in) :: &
     n
   real(pReal),                  intent(in) :: &
     time
   real(pReal),   dimension(2),  intent(out) :: &
     f
 
-  f(1) = homogenization_f_T(discretization_Marc_FEM2DAMASK_cell(n(3),n(1)))
+
+  f(1) = homogenization_f_T(discretization_Marc_FEM2DAMASK_cell(int(n(3)),int(n(1))))
   f(2) = 0.0_pReal
 
  end subroutine flux
@@ -409,8 +410,9 @@ subroutine uedinc(inc,incsub)
   use materialpoint_Marc
   use discretization_Marc
 
-  implicit none
-  integer, intent(in) :: inc, incsub
+  implicit none(type,external)
+  integer(pI64), intent(in) :: inc, incsub
+
   integer :: n, nqncomp, nqdatatype
   integer, save :: inc_written
   real(pReal), allocatable, dimension(:,:) :: d_n
@@ -427,9 +429,9 @@ subroutine uedinc(inc,incsub)
     enddo
 
     call discretization_Marc_UpdateNodeAndIpCoords(d_n)
-    call materialpoint_results(inc,cptim)
+    call materialpoint_results(int(inc),cptim)
 
-    inc_written = inc
+    inc_written = int(inc)
   endif
 
 end subroutine uedinc

src/Marc/discretization_Marc.f90

@@ -17,7 +17,7 @@ module discretization_Marc
   use geometry_plastic_nonlocal
   use results
 
-  implicit none
+  implicit none(type,external)
   private
 
   real(pReal),                         public, protected :: &
@@ -80,13 +80,13 @@ subroutine discretization_Marc_init
 
   num_commercialFEM => config_numerics%get('commercialFEM',defaultVal = emptyDict)
   mesh_unitlength = num_commercialFEM%get_asFloat('unitlength',defaultVal=1.0_pReal)                ! set physical extent of a length unit in mesh
-  if (mesh_unitlength <= 0.0_pReal) call IO_error(301,ext_msg='unitlength')
+  if (mesh_unitlength <= 0.0_pReal) call IO_error(301,'unitlength')
 
   call inputRead(elem,node0_elem,connectivity_elem,materialAt)
   nElems = size(connectivity_elem,2)
 
-  if (debug_e < 1 .or. debug_e > nElems)    call IO_error(602,ext_msg='element')
-  if (debug_i < 1 .or. debug_i > elem%nIPs) call IO_error(602,ext_msg='IP')
+  if (debug_e < 1 .or. debug_e > nElems)    call IO_error(602,'element')
+  if (debug_i < 1 .or. debug_i > elem%nIPs) call IO_error(602,'IP')
 
   allocate(cellNodeDefinition(elem%nNodes-1))
   allocate(connectivity_cell(elem%NcellNodesPerCell,elem%nIPs,nElems))
@@ -579,7 +579,7 @@ subroutine inputRead_elemType(elem, &
   character(len=*), dimension(:), intent(in)  :: fileContent                                        !< file content, separated per lines
 
   integer, allocatable, dimension(:) :: chunkPos
-  integer :: i,j,t,l,remainingChunks
+  integer :: i,j,t,t_,l,remainingChunks
 
 
   t = -1
@@ -594,7 +594,8 @@ subroutine inputRead_elemType(elem, &
           t = mapElemtype(IO_stringValue(fileContent(l+1+i+j),chunkPos,2))
           call elem%init(t)
         else
-          if (t /= mapElemtype(IO_stringValue(fileContent(l+1+i+j),chunkPos,2))) call IO_error(191,el=t,ip=i)
+          t_ = mapElemtype(IO_stringValue(fileContent(l+1+i+j),chunkPos,2))
+          if (t /= t_) call IO_error(191,IO_stringValue(fileContent(l+1+i+j),chunkPos,2),label1='type',ID1=t)
         endif
         remainingChunks = elem%nNodes - (chunkPos(1) - 2)
         do while(remainingChunks > 0)
@@ -616,7 +617,8 @@ subroutine inputRead_elemType(elem, &
 
    character(len=*), intent(in) :: what
 
-   select case (IO_lc(what))
+
+   select case (what)
       case (   '6')
         mapElemtype = 1            ! Two-dimensional Plane Strain Triangle
       case ( '125')                ! 155, 128 (need test)
@@ -644,7 +646,7 @@ subroutine inputRead_elemType(elem, &
       case ( '21')
         mapElemtype = 13           ! Three-dimensional Arbitrarily Distorted quadratic hexahedral
       case default
-        call IO_error(error_ID=190,ext_msg=IO_lc(what))
+        call IO_error(190,what)
    end select
 
   end function mapElemtype

src/Marc/element.f90

@@ -5,7 +5,7 @@
 module element
   use IO
 
-  implicit none
+  implicit none(type,external)
   private
 
 !---------------------------------------------------------------------------------------------------
@@ -714,7 +714,7 @@ subroutine tElement_init(self,elemType)
     case(13)
       self%cellNodeParentNodeWeights = CELLNODEPARENTNODEWEIGHTS13
     case default
-      call IO_error(0,ext_msg='invalid element type')
+      error stop 'invalid element type'
   end select
 
 

src/Marc/materialpoint_Marc.f90

@@ -23,7 +23,7 @@ module materialpoint_Marc
   use discretization
   use discretization_Marc
 
-  implicit none
+  implicit none(type,external)
   private
 
   real(pReal), dimension (:,:,:),   allocatable, private :: &
@@ -240,7 +240,8 @@ subroutine materialpoint_general(mode, ffn, ffn1, temperature_inp, dt, elFE, ip,
 
   endif
 
-  if (all(abs(materialpoint_dcsdE(1:6,1:6,ip,elCP)) < 1e-10_pReal)) call IO_warning(601,elCP,ip)
+  if (all(abs(materialpoint_dcsdE(1:6,1:6,ip,elCP)) < 1e-10_pReal)) &
+    call IO_warning(601,label1='element (CP)',ID1=elCP,label2='IP',ID2=ip)
 
   cauchyStress = materialpoint_cs   (1:6,    ip,elCP)
   jacobian     = materialpoint_dcsdE(1:6,1:6,ip,elCP)

src/YAML_parse.f90

@@ -12,7 +12,7 @@ module YAML_parse
   use system_routines
 #endif
 
-  implicit none
+  implicit none(type,external)
   private
 
   public :: &
@@ -24,6 +24,7 @@ module YAML_parse
 
     subroutine to_flow_C(flow,length_flow,mixed) bind(C)
       use, intrinsic :: ISO_C_Binding, only: C_INT, C_CHAR, C_PTR
+      implicit none(type,external)
 
       type(C_PTR), intent(out) :: flow
       integer(C_INT), intent(out) :: length_flow

src/YAML_types.f90

@@ -11,7 +11,7 @@ module YAML_types
   use IO
   use prec
 
-  implicit none
+  implicit none(type,external)
   private
 
   type, abstract, public :: tNode

src/config.f90

@@ -9,7 +9,7 @@ module config
   use results
   use parallelization
 
-  implicit none
+  implicit none(type,external)
   private
 
   class(tNode), pointer, public :: &

src/constants.f90

@@ -5,7 +5,7 @@
 module constants
   use prec
 
-  implicit none
+  implicit none(type,external)
   public
 
   real(pReal), parameter :: &

src/discretization.f90

@@ -7,7 +7,7 @@ module discretization
   use prec
   use results
 
-  implicit none
+  implicit none(type,external)
   private
 
   integer,     public, protected :: &

src/geometry_plastic_nonlocal.f90

@@ -9,7 +9,7 @@ module geometry_plastic_nonlocal
   use prec
   use results
 
-  implicit none
+  implicit none(type,external)
   public
 
   integer, protected :: &

src/grid/DAMASK_grid.f90

@@ -30,7 +30,11 @@ program DAMASK_grid
   use grid_thermal_spectral
   use results
 
+#if (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR>14) && !defined(PETSC_HAVE_MPI_F90MODULE_VISIBILITY)
+  implicit none(type,external)
+#else
   implicit none
+#endif
 
   type :: tLoadCase
     type(tRotation)          :: rot                                                                !< rotation of BC
@@ -272,7 +276,7 @@ program DAMASK_grid
         write(IO_STDOUT,'(2x,a,/,3(3(3x,f12.7,1x)/))',advance='no') 'R:',&
                  transpose(loadCases(l)%rot%asMatrix())
 
-      if (loadCases(l)%r <= 0.0)       errorID = 833
+      if (loadCases(l)%r <= 0.0_pReal) errorID = 833
       if (loadCases(l)%t < 0.0_pReal)  errorID = 834
       if (loadCases(l)%N < 1)          errorID = 835
       if (loadCases(l)%f_out < 1)      errorID = 836
@@ -290,7 +294,7 @@ program DAMASK_grid
       if (loadCases(l)%f_restart < huge(0)) &
         print'(2x,a,1x,i0)',   'f_restart:', loadCases(l)%f_restart
 
-      if (errorID > 0) call IO_error(error_ID = errorID, el = l)
+      if (errorID > 0) call IO_error(errorID,label1='line',ID1=l)
 
     endif reportAndCheck
   enddo
@@ -501,7 +505,7 @@ subroutine getMaskedTensor(values,mask,tensor)
   integer :: i,j
 
 
-  values = 0.0
+  values = 0.0_pReal
   do i = 1,3
     row => tensor%get(i)
     do j = 1,3

src/grid/FFTW.f90

@@ -0,0 +1,13 @@
+!--------------------------------------------------------------------------------------------------
+!> @author Martin Diehl, KU Leuven
+!> @brief Wrap FFTW3 into a module.
+!--------------------------------------------------------------------------------------------------
+module FFTW3
+  use, intrinsic :: ISO_C_binding
+
+  implicit none(type,external)
+  public
+
+  include 'fftw3-mpi.f03'
+
+end module FFTW3

src/grid/VTI.f90

@@ -8,7 +8,7 @@ module VTI
   use base64
   use IO
 
-  implicit none
+  implicit none(type,external)
   private
 
   public :: &
@@ -151,7 +151,7 @@ subroutine VTI_readCellsSizeOrigin(cells,geomSize,origin, &
   character(len=*),          intent(in) :: &
     fileContent
 
-  character(len=:), allocatable :: dataType, headerType
+  character(len=:), allocatable :: headerType
   logical :: inFile, inImage, compressed
   integer(pI64) :: &
     startPos, endPos

src/grid/base64.f90

@@ -7,7 +7,7 @@ module base64
   use prec
   use IO
 
-  implicit none
+  implicit none(type,external)
   private
 
   character(len=*), parameter :: &

src/grid/discretization_grid.f90

@@ -10,6 +10,7 @@ module discretization_grid
 #if (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR>14) && !defined(PETSC_HAVE_MPI_F90MODULE_VISIBILITY)
   use MPI_f08
 #endif
+  use FFTW3
 
   use prec
   use parallelization
@@ -22,7 +23,11 @@ module discretization_grid
   use discretization
   use geometry_plastic_nonlocal
 
+#if (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR>14) && !defined(PETSC_HAVE_MPI_F90MODULE_VISIBILITY)
+  implicit none(type,external)
+#else
   implicit none
+#endif
   private
 
   integer,     dimension(3), public, protected :: &
@@ -50,7 +55,6 @@ subroutine discretization_grid_init(restart)
 
   logical, intent(in) :: restart
 
-  include 'fftw3-mpi.f03'
   real(pReal), dimension(3) :: &
     mySize, &                                                                                       !< domain size of this process
     origin                                                                                          !< (global) distance to origin
@@ -107,10 +111,8 @@ subroutine discretization_grid_init(restart)
 
   if (worldsize>cells(3)) call IO_error(894, ext_msg='number of processes exceeds cells(3)')
 
-  call fftw_mpi_init
-  devNull = fftw_mpi_local_size_3d(int(cells(3),C_INTPTR_T), &
-                                   int(cells(2),C_INTPTR_T), &
-                                   int(cells(1),C_INTPTR_T)/2+1, &
+  call fftw_mpi_init()
+  devNull = fftw_mpi_local_size_3d(int(cells(3),C_INTPTR_T),int(cells(2),C_INTPTR_T),int(cells(1)/2+1,C_INTPTR_T), &
                                    PETSC_COMM_WORLD, &
                                    z, &                                                             ! domain cells size along z
                                    z_offset)                                                        ! domain cells offset along z
@@ -231,9 +233,9 @@ pure function cellSurfaceArea(geomSize,cells)
   real(pReal), dimension(6,1,product(cells)) :: cellSurfaceArea
 
 
-  cellSurfaceArea(1:2,1,:) = geomSize(2)/real(cells(2)) * geomSize(3)/real(cells(3))
-  cellSurfaceArea(3:4,1,:) = geomSize(3)/real(cells(3)) * geomSize(1)/real(cells(1))
-  cellSurfaceArea(5:6,1,:) = geomSize(1)/real(cells(1)) * geomSize(2)/real(cells(2))
+  cellSurfaceArea(1:2,1,:) = geomSize(2)/real(cells(2),pReal) * geomSize(3)/real(cells(3),pReal)
+  cellSurfaceArea(3:4,1,:) = geomSize(3)/real(cells(3),pReal) * geomSize(1)/real(cells(1),pReal)
+  cellSurfaceArea(5:6,1,:) = geomSize(1)/real(cells(1),pReal) * geomSize(2)/real(cells(2),pReal)
 
 end function cellSurfaceArea
 

src/grid/grid_damage_spectral.f90

@@ -22,7 +22,11 @@ module grid_damage_spectral
   use YAML_types
   use config
 
+#if (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR>14) && !defined(PETSC_HAVE_MPI_F90MODULE_VISIBILITY)
+  implicit none(type,external)
+#else
   implicit none
+#endif
   private
 
   type :: tNumerics

src/grid/grid_mech_FEM.f90

@@ -27,7 +27,12 @@ module grid_mechanical_FEM
   use discretization
   use discretization_grid
 
+
+#if (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR>14) && !defined(PETSC_HAVE_MPI_F90MODULE_VISIBILITY)
+  implicit none(type,external)
+#else
   implicit none
+#endif
   private
 
   type(tSolutionParams) :: params
@@ -115,6 +120,8 @@ subroutine grid_mechanical_FEM_init
   class(tNode), pointer :: &
     num_grid, &
     debug_grid
+  character(len=pStringLen) :: &
+    extmsg = ''
 
   print'(/,1x,a)', '<<<+-  grid_mechanical_FEM init  -+>>>'; flush(IO_STDOUT)
 
@@ -134,12 +141,14 @@ subroutine grid_mechanical_FEM_init
   num%itmin           = num_grid%get_asInt  ('itmin',defaultVal=1)
   num%itmax           = num_grid%get_asInt  ('itmax',defaultVal=250)
 
-  if (num%eps_div_atol <= 0.0_pReal)             call IO_error(301,ext_msg='eps_div_atol')
-  if (num%eps_div_rtol < 0.0_pReal)              call IO_error(301,ext_msg='eps_div_rtol')
-  if (num%eps_stress_atol <= 0.0_pReal)          call IO_error(301,ext_msg='eps_stress_atol')
-  if (num%eps_stress_rtol < 0.0_pReal)           call IO_error(301,ext_msg='eps_stress_rtol')
-  if (num%itmax <= 1)                            call IO_error(301,ext_msg='itmax')
-  if (num%itmin > num%itmax .or. num%itmin < 1)  call IO_error(301,ext_msg='itmin')
+  if (num%eps_div_atol <= 0.0_pReal)             extmsg = trim(extmsg)//' eps_div_atol'
+  if (num%eps_div_rtol < 0.0_pReal)              extmsg = trim(extmsg)//' eps_div_rtol'
+  if (num%eps_stress_atol <= 0.0_pReal)          extmsg = trim(extmsg)//' eps_stress_atol'
+  if (num%eps_stress_rtol < 0.0_pReal)           extmsg = trim(extmsg)//' eps_stress_rtol'
+  if (num%itmax <= 1)                            extmsg = trim(extmsg)//' itmax'
+  if (num%itmin > num%itmax .or. num%itmin < 1)  extmsg = trim(extmsg)//' itmin'
+
+  if (extmsg /= '') call IO_error(301,ext_msg=trim(extmsg))
 
 !--------------------------------------------------------------------------------------------------
 ! set default and user defined options for PETSc
@@ -217,14 +226,14 @@ subroutine grid_mechanical_FEM_init
   delta = geomSize/real(cells,pReal)                                                                ! grid spacing
   detJ = product(delta)                                                                             ! cell volume
 
-  BMat = reshape(real([-1.0_pReal/delta(1),-1.0_pReal/delta(2),-1.0_pReal/delta(3), &
-                        1.0_pReal/delta(1),-1.0_pReal/delta(2),-1.0_pReal/delta(3), &
-                       -1.0_pReal/delta(1), 1.0_pReal/delta(2),-1.0_pReal/delta(3), &
-                        1.0_pReal/delta(1), 1.0_pReal/delta(2),-1.0_pReal/delta(3), &
-                       -1.0_pReal/delta(1),-1.0_pReal/delta(2), 1.0_pReal/delta(3), &
-                        1.0_pReal/delta(1),-1.0_pReal/delta(2), 1.0_pReal/delta(3), &
-                       -1.0_pReal/delta(1), 1.0_pReal/delta(2), 1.0_pReal/delta(3), &
-                        1.0_pReal/delta(1), 1.0_pReal/delta(2), 1.0_pReal/delta(3)],pReal), [3,8])/4.0_pReal ! shape function derivative matrix
+  BMat = reshape(real([-delta(1)**(-1),-delta(2)**(-1),-delta(3)**(-1), &
+                        delta(1)**(-1),-delta(2)**(-1),-delta(3)**(-1), &
+                       -delta(1)**(-1), delta(2)**(-1),-delta(3)**(-1), &
+                        delta(1)**(-1), delta(2)**(-1),-delta(3)**(-1), &
+                       -delta(1)**(-1),-delta(2)**(-1), delta(3)**(-1), &
+                        delta(1)**(-1),-delta(2)**(-1), delta(3)**(-1), &
+                       -delta(1)**(-1), delta(2)**(-1), delta(3)**(-1), &
+                        delta(1)**(-1), delta(2)**(-1), delta(3)**(-1)],pReal), [3,8])/4.0_pReal    ! shape function derivative matrix
 
   HGMat = matmul(transpose(HGcomp),HGcomp) &
         * HGCoeff*(delta(1)*delta(2) + delta(2)*delta(3) + delta(3)*delta(1))/16.0_pReal            ! hourglass stabilization matrix
@@ -652,7 +661,7 @@ subroutine formJacobian(da_local,x_local,Jac_pre,Jac,dummy,err_PETSc)
   MatNullSpace                         :: matnull
   PetscErrorCode                       :: err_PETSc
 
-  BMatFull = 0.0
+  BMatFull = 0.0_pReal
   BMatFull(1:3,1 :8 ) = BMat
   BMatFull(4:6,9 :16) = BMat
   BMatFull(7:9,17:24) = BMat
@@ -682,7 +691,7 @@ subroutine formJacobian(da_local,x_local,Jac_pre,Jac,dummy,err_PETSc)
     enddo; enddo; enddo
     row = col
     ce = ce + 1
-    K_ele = 0.0
+    K_ele = 0.0_pReal
     K_ele(1 :8 ,1 :8 ) = HGMat*(homogenization_dPdF(1,1,1,1,ce) + &
                                 homogenization_dPdF(2,2,2,2,ce) + &
                                 homogenization_dPdF(3,3,3,3,ce))/3.0_pReal

src/grid/grid_mech_spectral_basic.f90

@@ -26,7 +26,11 @@ module grid_mechanical_spectral_basic
   use homogenization
   use discretization_grid
 
+#if (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR>14) && !defined(PETSC_HAVE_MPI_F90MODULE_VISIBILITY)
+  implicit none(type,external)
+#else
   implicit none
+#endif
   private
 
   type(tSolutionParams) :: params
@@ -117,6 +121,8 @@ subroutine grid_mechanical_spectral_basic_init
   class (tNode), pointer :: &
     num_grid, &
     debug_grid
+  character(len=pStringLen) :: &
+    extmsg = ''
 
   print'(/,1x,a)', '<<<+-  grid_mechanical_spectral_basic init  -+>>>'; flush(IO_STDOUT)
 
@@ -143,12 +149,14 @@ subroutine grid_mechanical_spectral_basic_init
   num%itmin           = num_grid%get_asInt  ('itmin',defaultVal=1)
   num%itmax           = num_grid%get_asInt  ('itmax',defaultVal=250)
 
-  if (num%eps_div_atol <= 0.0_pReal)             call IO_error(301,ext_msg='eps_div_atol')
-  if (num%eps_div_rtol < 0.0_pReal)              call IO_error(301,ext_msg='eps_div_rtol')
-  if (num%eps_stress_atol <= 0.0_pReal)          call IO_error(301,ext_msg='eps_stress_atol')
-  if (num%eps_stress_rtol < 0.0_pReal)           call IO_error(301,ext_msg='eps_stress_rtol')
-  if (num%itmax <= 1)                            call IO_error(301,ext_msg='itmax')
-  if (num%itmin > num%itmax .or. num%itmin < 1)  call IO_error(301,ext_msg='itmin')
+  if (num%eps_div_atol <= 0.0_pReal)             extmsg = trim(extmsg)//' eps_div_atol'
+  if (num%eps_div_rtol < 0.0_pReal)              extmsg = trim(extmsg)//' eps_div_rtol'
+  if (num%eps_stress_atol <= 0.0_pReal)          extmsg = trim(extmsg)//' eps_stress_atol'
+  if (num%eps_stress_rtol < 0.0_pReal)           extmsg = trim(extmsg)//' eps_stress_rtol'
+  if (num%itmax <= 1)                            extmsg = trim(extmsg)//' itmax'
+  if (num%itmin > num%itmax .or. num%itmin < 1)  extmsg = trim(extmsg)//' itmin'
+
+  if (extmsg /= '') call IO_error(301,ext_msg=trim(extmsg))
 
 !--------------------------------------------------------------------------------------------------
 ! set default and user defined options for PETSc

src/grid/grid_mech_spectral_polarisation.f90

@@ -26,7 +26,11 @@ module grid_mechanical_spectral_polarisation
   use homogenization
   use discretization_grid
 
+#if (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR>14) && !defined(PETSC_HAVE_MPI_F90MODULE_VISIBILITY)
+  implicit none(type,external)
+#else
   implicit none
+#endif
   private
 
   type(tSolutionParams) :: params
@@ -130,6 +134,8 @@ subroutine grid_mechanical_spectral_polarisation_init
   class (tNode), pointer :: &
     num_grid, &
     debug_grid
+  character(len=pStringLen) :: &
+    extmsg = ''
 
   print'(/,1x,a)', '<<<+-  grid_mechanical_spectral_polarization init  -+>>>'; flush(IO_STDOUT)
 
@@ -157,16 +163,18 @@ subroutine grid_mechanical_spectral_polarisation_init
   num%alpha           = num_grid%get_asFloat('alpha',          defaultVal=1.0_pReal)
   num%beta            = num_grid%get_asFloat('beta',           defaultVal=1.0_pReal)
 
-  if (num%eps_div_atol <= 0.0_pReal)                      call IO_error(301,ext_msg='eps_div_atol')
-  if (num%eps_div_rtol < 0.0_pReal)                       call IO_error(301,ext_msg='eps_div_rtol')
-  if (num%eps_curl_atol <= 0.0_pReal)                     call IO_error(301,ext_msg='eps_curl_atol')
-  if (num%eps_curl_rtol < 0.0_pReal)                      call IO_error(301,ext_msg='eps_curl_rtol')
-  if (num%eps_stress_atol <= 0.0_pReal)                   call IO_error(301,ext_msg='eps_stress_atol')
-  if (num%eps_stress_rtol < 0.0_pReal)                    call IO_error(301,ext_msg='eps_stress_rtol')
-  if (num%itmax <= 1)                                     call IO_error(301,ext_msg='itmax')
-  if (num%itmin > num%itmax .or. num%itmin < 1)           call IO_error(301,ext_msg='itmin')
-  if (num%alpha <= 0.0_pReal .or. num%alpha >  2.0_pReal) call IO_error(301,ext_msg='alpha')
-  if (num%beta < 0.0_pReal .or. num%beta > 2.0_pReal)     call IO_error(301,ext_msg='beta')
+  if (num%eps_div_atol <= 0.0_pReal)                      extmsg = trim(extmsg)//' eps_div_atol'
+  if (num%eps_div_rtol < 0.0_pReal)                       extmsg = trim(extmsg)//' eps_div_rtol'
+  if (num%eps_curl_atol <= 0.0_pReal)                     extmsg = trim(extmsg)//' eps_curl_atol'
+  if (num%eps_curl_rtol < 0.0_pReal)                      extmsg = trim(extmsg)//' eps_curl_rtol'
+  if (num%eps_stress_atol <= 0.0_pReal)                   extmsg = trim(extmsg)//' eps_stress_atol'
+  if (num%eps_stress_rtol < 0.0_pReal)                    extmsg = trim(extmsg)//' eps_stress_rtol'
+  if (num%itmax <= 1)                                     extmsg = trim(extmsg)//' itmax'
+  if (num%itmin > num%itmax .or. num%itmin < 1)           extmsg = trim(extmsg)//' itmin'
+  if (num%alpha <= 0.0_pReal .or. num%alpha >  2.0_pReal) extmsg = trim(extmsg)//' alpha'
+  if (num%beta < 0.0_pReal .or. num%beta > 2.0_pReal)     extmsg = trim(extmsg)//' beta'
+
+  if (extmsg /= '') call IO_error(301,ext_msg=trim(extmsg))
 
 !--------------------------------------------------------------------------------------------------
 ! set default and user defined options for PETSc

src/grid/grid_thermal_spectral.f90

@@ -25,7 +25,11 @@ module grid_thermal_spectral
   use YAML_types
   use config
 
+#if (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR>14) && !defined(PETSC_HAVE_MPI_F90MODULE_VISIBILITY)
+  implicit none(type,external)
+#else
   implicit none
+#endif
   private
 
   type :: tNumerics

src/grid/spectral_utilities.f90

@@ -4,13 +4,12 @@
 !> @brief Utilities used by the different spectral solver variants
 !--------------------------------------------------------------------------------------------------
 module spectral_utilities
-  use, intrinsic :: iso_c_binding
-
 #include <petsc/finclude/petscsys.h>
   use PETScSys
 #if (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR>14) && !defined(PETSC_HAVE_MPI_F90MODULE_VISIBILITY)
   use MPI_f08
 #endif
+  use FFTW3
 
   use prec
   use CLI
@@ -23,26 +22,32 @@ module spectral_utilities
   use discretization
   use homogenization
 
-  implicit none
-  private
 
-  include 'fftw3-mpi.f03'
+#if (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR>14) && !defined(PETSC_HAVE_MPI_F90MODULE_VISIBILITY)
+  implicit none(type,external)
+#else
+  implicit none
+#endif
+  private
 
 !--------------------------------------------------------------------------------------------------
 ! grid related information
   real(pReal), protected, public               :: wgt                                               !< weighting factor 1/Nelems
-  integer,     protected,  public                :: cells1Red                                       !< cells(1)/2
   real(pReal), protected, public, dimension(3) :: scaledGeomSize                                    !< scaled geometry size for calculation of divergence
+  integer :: &
+    cells1Red, &                                                                                    !< cells(1)/2+1
+    cells2, &                                                                                       !< (local) cells in 2nd direction
+    cells2Offset                                                                                    !< (local) cells offset in 2nd direction
 
 !--------------------------------------------------------------------------------------------------
 ! variables storing information for spectral method and FFTW
 
-  real   (C_DOUBLE),        public,  dimension(:,:,:,:,:),     pointer     :: tensorField_real      !< real representation (some stress or deformation) of field_fourier
-  complex(C_DOUBLE_COMPLEX),public,  dimension(:,:,:,:,:),     pointer     :: tensorField_fourier   !< field on which the Fourier transform operates
-  real(C_DOUBLE),           public,  dimension(:,:,:,:),       pointer     :: vectorField_real      !< vector field real representation for fftw
-  complex(C_DOUBLE_COMPLEX),public,  dimension(:,:,:,:),       pointer     :: vectorField_fourier   !< vector field fourier representation for fftw
-  real(C_DOUBLE),           public,  dimension(:,:,:),         pointer     :: scalarField_real      !< scalar field real representation for fftw
-  complex(C_DOUBLE_COMPLEX),public,  dimension(:,:,:),         pointer     :: scalarField_fourier   !< scalar field fourier representation for fftw
+  real(C_DOUBLE),         public, dimension(:,:,:,:,:),     pointer     :: tensorField_real         !< tensor field in real space
+  real(C_DOUBLE),         public, dimension(:,:,:,:),       pointer     :: vectorField_real         !< vector field in real space
+  real(C_DOUBLE),         public, dimension(:,:,:),         pointer     :: scalarField_real         !< scalar field in real space
+  complex(C_DOUBLE_COMPLEX),      dimension(:,:,:,:,:),     pointer     :: tensorField_fourier      !< tensor field in Fourier space
+  complex(C_DOUBLE_COMPLEX),      dimension(:,:,:,:),       pointer     :: vectorField_fourier      !< vector field in Fourier space
+  complex(C_DOUBLE_COMPLEX),      dimension(:,:,:),         pointer     :: scalarField_fourier      !< scalar field in Fourier space
   complex(pReal),                 dimension(:,:,:,:,:,:,:), allocatable :: gamma_hat                !< gamma operator (field) for spectral method
   complex(pReal),                 dimension(:,:,:,:),       allocatable :: xi1st                    !< wave vector field for first derivatives
   complex(pReal),                 dimension(:,:,:,:),       allocatable :: xi2nd                    !< wave vector field for second derivatives
@@ -149,13 +154,16 @@ subroutine spectral_utilities_init()
     FFTW_planner_flag
   integer, dimension(3) :: k_s
   type(C_PTR) :: &
-    tensorField, &                                                                                  !< field containing data for FFTW in real and fourier space (in place)
-    vectorField, &                                                                                  !< field containing data for FFTW in real space when debugging FFTW (no in place)
-    scalarField                                                                                     !< field containing data for FFTW in real space when debugging FFTW (no in place)
-  integer(C_INTPTR_T), dimension(3) :: gridFFTW
-  integer(C_INTPTR_T) :: alloc_local, local_K, local_K_offset
+    tensorField, &                                                                                  !< tensor data for FFTW in real and Fourier space (in-place)
+    vectorField, &                                                                                  !< vector data for FFTW in real and Fourier space (in-place)
+    scalarField                                                                                     !< scalar data for FFTW in real and Fourier space (in-place)
+  integer(C_INTPTR_T), dimension(3) :: cellsFFTW
+  integer(C_INTPTR_T) :: N, &
+    cells3FFTW, &                                                                                   !< # of cells in 3. dim on current process in real space
+    cells3_offset, &                                                                                !< offset for cells in 3. dim on current process in real space
+    cells2FFTW, &                                                                                   !< # of cells in 2. dim on current process in Fourier space
+    cells2_offset                                                                                   !< offset for cells in 2. dim on curren process in Fourier space
   integer(C_INTPTR_T), parameter :: &
-    scalarSize = 1_C_INTPTR_T, &
     vectorSize = 3_C_INTPTR_T, &
     tensorSize = 9_C_INTPTR_T
   character(len=*), parameter :: &
@@ -202,7 +210,7 @@ subroutine spectral_utilities_init()
   CHKERRQ(err_PETSc)
 
   cells1Red = cells(1)/2 + 1
-  wgt = 1.0/real(product(cells),pReal)
+  wgt = real(product(cells),pReal)**(-1)
 
   num%memory_efficient      = num_grid%get_asInt('memory_efficient',      defaultVal=1) > 0         ! ToDo: should be logical in YAML file
   num%divergence_correction = num_grid%get_asInt('divergence_correction', defaultVal=2)
@@ -249,7 +257,7 @@ subroutine spectral_utilities_init()
     case('fftw_exhaustive')
       FFTW_planner_flag = FFTW_EXHAUSTIVE
     case default
-      call IO_warning(warning_ID=47,ext_msg=trim(IO_lc(num_grid%get_asString('fftw_plan_mode'))))
+      call IO_warning(47,'using default FFTW_MEASURE instead of "'//trim(num_grid%get_asString('fftw_plan_mode'))//'"')
       FFTW_planner_flag = FFTW_MEASURE
   end select
 
@@ -260,94 +268,107 @@ subroutine spectral_utilities_init()
 
   print'(/,1x,a)', 'FFTW initialized'; flush(IO_STDOUT)
 
-!--------------------------------------------------------------------------------------------------
-! MPI allocation
-  gridFFTW = int(cells,C_INTPTR_T)
-  alloc_local = fftw_mpi_local_size_3d(gridFFTW(3), gridFFTW(2), gridFFTW(1)/2 +1, &
-                                       PETSC_COMM_WORLD, local_K, local_K_offset)
-  allocate (xi1st (3,cells1Red,cells(2),cells3),source = cmplx(0.0_pReal,0.0_pReal,pReal))          ! frequencies for first derivatives, only half the size for first dimension
-  allocate (xi2nd (3,cells1Red,cells(2),cells3),source = cmplx(0.0_pReal,0.0_pReal,pReal))          ! frequencies for second derivatives, only half the size for first dimension
+  cellsFFTW = int(cells,C_INTPTR_T)
 
-  tensorField = fftw_alloc_complex(tensorSize*alloc_local)
-  call c_f_pointer(tensorField, tensorField_real,    [3_C_INTPTR_T,3_C_INTPTR_T, &
-                   2_C_INTPTR_T*(gridFFTW(1)/2_C_INTPTR_T + 1_C_INTPTR_T),gridFFTW(2),local_K])     ! place a pointer for a real tensor representation
-  call c_f_pointer(tensorField, tensorField_fourier, [3_C_INTPTR_T,3_C_INTPTR_T, &
-                   gridFFTW(1)/2_C_INTPTR_T + 1_C_INTPTR_T ,              gridFFTW(2),local_K])     ! place a pointer for a fourier tensor representation
+  N = fftw_mpi_local_size_many_transposed(3,[cellsFFTW(3),cellsFFTW(2),int(cells1Red,C_INTPTR_T)], &
+                                          tensorSize,FFTW_MPI_DEFAULT_BLOCK,FFTW_MPI_DEFAULT_BLOCK,PETSC_COMM_WORLD, &
+                                          cells3FFTW,cells3_offset,cells2FFTW,cells2_offset)
+  cells2 = int(cells2FFTW)
+  cells2Offset = int(cells2_offset)
+  if (int(cells3FFTW) /= cells3) error stop 'domain decomposition mismatch (tensor, real space)'
+  tensorField = fftw_alloc_complex(N)
+  call c_f_pointer(tensorField,tensorField_real, &
+                   [3_C_INTPTR_T,3_C_INTPTR_T,int(cells1Red*2,C_INTPTR_T),cellsFFTW(2),cells3FFTW])
+  call c_f_pointer(tensorField,tensorField_fourier, &
+                   [3_C_INTPTR_T,3_C_INTPTR_T,int(cells1Red,  C_INTPTR_T),cellsFFTW(3),cells2FFTW])
 
-  vectorField = fftw_alloc_complex(vectorSize*alloc_local)
-  call c_f_pointer(vectorField, vectorField_real,   [3_C_INTPTR_T,&
-                   2_C_INTPTR_T*(gridFFTW(1)/2_C_INTPTR_T + 1_C_INTPTR_T),gridFFTW(2),local_K])     ! place a pointer for a real vector representation
-  call c_f_pointer(vectorField, vectorField_fourier,[3_C_INTPTR_T,&
-                   gridFFTW(1)/2_C_INTPTR_T + 1_C_INTPTR_T,               gridFFTW(2),local_K])     ! place a pointer for a fourier vector representation
 
-  scalarField = fftw_alloc_complex(scalarSize*alloc_local)                                          ! allocate data for real representation (no in place transform)
+  N = fftw_mpi_local_size_many_transposed(3,[cellsFFTW(3),cellsFFTW(2),int(cells1Red,C_INTPTR_T)], &
+                                          vectorSize,FFTW_MPI_DEFAULT_BLOCK,FFTW_MPI_DEFAULT_BLOCK,PETSC_COMM_WORLD, &
+                                          cells3FFTW,cells3_offset,cells2FFTW,cells2_offset)
+  if (int(cells3FFTW) /= cells3) error stop 'domain decomposition mismatch (vector, real space)'
+  if (int(cells2FFTW) /= cells2) error stop 'domain decomposition mismatch (vector, Fourier space)'
+  vectorField = fftw_alloc_complex(N)
+  call c_f_pointer(vectorField,vectorField_real, &
+                   [3_C_INTPTR_T,int(cells1Red*2,C_INTPTR_T),cellsFFTW(2),cells3FFTW])
+  call c_f_pointer(vectorField,vectorField_fourier, &
+                   [3_C_INTPTR_T,int(cells1Red,  C_INTPTR_T),cellsFFTW(3),cells2FFTW])
+
+  N = fftw_mpi_local_size_3d_transposed(cellsFFTW(3),cellsFFTW(2),int(cells1Red,C_INTPTR_T), &
+                                        PETSC_COMM_WORLD,cells3FFTW,cells3_offset,cells2FFTW,cells2_offset)
+  if (int(cells3FFTW) /= cells3) error stop 'domain decomposition mismatch (scalar, real space)'
+  if (int(cells2FFTW) /= cells2) error stop 'domain decomposition mismatch (scalar, Fourier space)'
+  scalarField = fftw_alloc_complex(N)
   call c_f_pointer(scalarField,scalarField_real, &
-                   [2_C_INTPTR_T*(gridFFTW(1)/2_C_INTPTR_T + 1),gridFFTW(2),local_K])               ! place a pointer for a real scalar representation
+                   [int(cells1Red*2,C_INTPTR_T),cellsFFTW(2),cells3FFTW])
   call c_f_pointer(scalarField,scalarField_fourier, &
-                    [             gridFFTW(1)/2_C_INTPTR_T + 1 ,gridFFTW(2),local_K])               ! place a pointer for a fourier scarlar representation
+                   [int(cells1Red,  C_INTPTR_T),cellsFFTW(3),cells2FFTW])
+
+!--------------------------------------------------------------------------------------------------
+! allocation
+  allocate (xi1st (3,cells1Red,cells(3),cells2),source = cmplx(0.0_pReal,0.0_pReal,pReal))            ! frequencies for first derivatives, only half the size for first dimension
+  allocate (xi2nd (3,cells1Red,cells(3),cells2),source = cmplx(0.0_pReal,0.0_pReal,pReal))            ! frequencies for second derivatives, only half the size for first dimension
 
 !--------------------------------------------------------------------------------------------------
 ! tensor MPI fftw plans
-  planTensorForth = fftw_mpi_plan_many_dft_r2c(3,gridFFTW(3:1:-1),tensorSize, &
+  planTensorForth = fftw_mpi_plan_many_dft_r2c(3,cellsFFTW(3:1:-1),tensorSize, &
                                                FFTW_MPI_DEFAULT_BLOCK,FFTW_MPI_DEFAULT_BLOCK, &
                                                tensorField_real,tensorField_fourier, &
-                                               PETSC_COMM_WORLD,FFTW_planner_flag)
-  if (.not. c_associated(planTensorForth)) error stop 'FFTW error'
-  planTensorBack  = fftw_mpi_plan_many_dft_c2r(3,gridFFTW(3:1:-1),tensorSize, &
+                                               PETSC_COMM_WORLD,FFTW_planner_flag+FFTW_MPI_TRANSPOSED_OUT)
+  if (.not. c_associated(planTensorForth)) error stop 'FFTW error r2c tensor'
+  planTensorBack  = fftw_mpi_plan_many_dft_c2r(3,cellsFFTW(3:1:-1),tensorSize, &
                                                FFTW_MPI_DEFAULT_BLOCK,FFTW_MPI_DEFAULT_BLOCK, &
                                                tensorField_fourier,tensorField_real, &
-                                               PETSC_COMM_WORLD, FFTW_planner_flag)
-  if (.not. c_associated(planTensorBack)) error stop 'FFTW error'
+                                               PETSC_COMM_WORLD,FFTW_planner_flag+FFTW_MPI_TRANSPOSED_IN)
+  if (.not. c_associated(planTensorBack))  error stop 'FFTW error c2r tensor'
 
 !--------------------------------------------------------------------------------------------------
 ! vector MPI fftw plans
-  planVectorForth = fftw_mpi_plan_many_dft_r2c(3,gridFFTW(3:1:-1),vectorSize, &
+  planVectorForth = fftw_mpi_plan_many_dft_r2c(3,cellsFFTW(3:1:-1),vectorSize, &
                                                FFTW_MPI_DEFAULT_BLOCK,FFTW_MPI_DEFAULT_BLOCK, &
                                                vectorField_real,vectorField_fourier, &
-                                               PETSC_COMM_WORLD,FFTW_planner_flag)
-  if (.not. c_associated(planVectorForth)) error stop 'FFTW error'
-  planVectorBack  = fftw_mpi_plan_many_dft_c2r(3,gridFFTW(3:1:-1),vectorSize, &
+                                               PETSC_COMM_WORLD,FFTW_planner_flag+FFTW_MPI_TRANSPOSED_OUT)
+  if (.not. c_associated(planVectorForth)) error stop 'FFTW error r2c vector'
+  planVectorBack  = fftw_mpi_plan_many_dft_c2r(3,cellsFFTW(3:1:-1),vectorSize, &
                                                FFTW_MPI_DEFAULT_BLOCK,FFTW_MPI_DEFAULT_BLOCK, &
                                                vectorField_fourier,vectorField_real, &
-                                               PETSC_COMM_WORLD, FFTW_planner_flag)
-  if (.not. c_associated(planVectorBack)) error stop 'FFTW error'
+                                               PETSC_COMM_WORLD,FFTW_planner_flag+FFTW_MPI_TRANSPOSED_IN)
+  if (.not. c_associated(planVectorBack))  error stop 'FFTW error c2r vector'
 
 !--------------------------------------------------------------------------------------------------
 ! scalar MPI fftw plans
-  planScalarForth = fftw_mpi_plan_many_dft_r2c(3,gridFFTW(3:1:-1),scalarSize, &
-                                               FFTW_MPI_DEFAULT_BLOCK,FFTW_MPI_DEFAULT_BLOCK, &
+  planScalarForth = fftw_mpi_plan_dft_r2c_3d(cellsFFTW(3),cellsFFTW(2),cellsFFTW(1), &
                                              scalarField_real,scalarField_fourier, &
-                                               PETSC_COMM_WORLD,FFTW_planner_flag)
-  if (.not. c_associated(planScalarForth)) error stop 'FFTW error'
-  planScalarBack  = fftw_mpi_plan_many_dft_c2r(3,gridFFTW(3:1:-1),scalarSize, &
-                                               FFTW_MPI_DEFAULT_BLOCK, FFTW_MPI_DEFAULT_BLOCK, &
+                                             PETSC_COMM_WORLD,FFTW_planner_flag+FFTW_MPI_TRANSPOSED_OUT)
+  if (.not. c_associated(planScalarForth)) error stop 'FFTW error r2c scalar'
+  planScalarBack  = fftw_mpi_plan_dft_c2r_3d(cellsFFTW(3),cellsFFTW(2),cellsFFTW(1), &
                                              scalarField_fourier,scalarField_real, &
-                                               PETSC_COMM_WORLD, FFTW_planner_flag)
-  if (.not. c_associated(planScalarBack)) error stop 'FFTW error'
+                                             PETSC_COMM_WORLD,FFTW_planner_flag+FFTW_MPI_TRANSPOSED_IN)
+  if (.not. c_associated(planScalarBack))  error stop 'FFTW error c2r scalar'
 
 !--------------------------------------------------------------------------------------------------
 ! calculation of discrete angular frequencies, ordered as in FFTW (wrap around)
-  do k = cells3Offset+1, cells3Offset+cells3
-    k_s(3) = k - 1
-    if (k > cells(3)/2 + 1) k_s(3) = k_s(3) - cells(3)                                              ! running from 0,1,...,N/2,N/2+1,-N/2,-N/2+1,...,-1
-      do j = 1, cells(2)
+  do j = cells2Offset+1, cells2Offset+cells2
     k_s(2) = j - 1
     if (j > cells(2)/2 + 1) k_s(2) = k_s(2) - cells(2)                                              ! running from 0,1,...,N/2,N/2+1,-N/2,-N/2+1,...,-1
+      do k = 1, cells(3)
+        k_s(3) = k - 1
+        if (k > cells(3)/2 + 1) k_s(3) = k_s(3) - cells(3)                                          ! running from 0,1,...,N/2,N/2+1,-N/2,-N/2+1,...,-1
           do i = 1, cells1Red
             k_s(1) = i - 1                                                                          ! symmetry, junst running from 0,1,...,N/2,N/2+1
-            xi2nd(1:3,i,j,k-cells3Offset) = utilities_getFreqDerivative(k_s)
+            xi2nd(1:3,i,k,j-cells2Offset) = utilities_getFreqDerivative(k_s)
             where(mod(cells,2)==0 .and. [i,j,k] == cells/2+1 .and. &
                   spectral_derivative_ID == DERIVATIVE_CONTINUOUS_ID)                               ! for even grids, set the Nyquist Freq component to 0.0
-              xi1st(1:3,i,j,k-cells3Offset) = cmplx(0.0_pReal,0.0_pReal,pReal)
+              xi1st(1:3,i,k,j-cells2Offset) = cmplx(0.0_pReal,0.0_pReal,pReal)
             elsewhere
-              xi1st(1:3,i,j,k-cells3Offset) = xi2nd(1:3,i,j,k-cells3Offset)
+              xi1st(1:3,i,k,j-cells2Offset) = xi2nd(1:3,i,k,j-cells2Offset)
             endwhere
   end do; end do; end do
 
   if (num%memory_efficient) then                                                                    ! allocate just single fourth order tensor
     allocate (gamma_hat(3,3,3,3,1,1,1), source = cmplx(0.0_pReal,0.0_pReal,pReal))
   else                                                                                              ! precalculation of gamma_hat field
-    allocate (gamma_hat(3,3,3,3,cells1Red,cells(2),cells3), source = cmplx(0.0_pReal,0.0_pReal,pReal))
+    allocate (gamma_hat(3,3,3,3,cells1Red,cells(3),cells2), source = cmplx(0.0_pReal,0.0_pReal,pReal))
   end if
 
   call selfTest()
@@ -376,33 +397,33 @@ subroutine utilities_updateGamma(C)
   if (.not. num%memory_efficient) then
     gamma_hat = cmplx(0.0_pReal,0.0_pReal,pReal)                                                    ! for the singular point and any non invertible A
     !$OMP PARALLEL DO PRIVATE(l,m,n,o,temp33_cmplx,xiDyad_cmplx,A,A_inv,err)
-    do k = cells3Offset+1, cells3Offset+cells3; do j = 1, cells(2); do i = 1, cells1Red
+    do j = cells2Offset+1, cells2Offset+cells2; do k = 1, cells(3); do i = 1, cells1Red
       if (any([i,j,k] /= 1)) then                                                                   ! singular point at xi=(0.0,0.0,0.0) i.e. i=j=k=1
 #ifndef __INTEL_COMPILER
         do concurrent(l = 1:3, m = 1:3)
-          xiDyad_cmplx(l,m) = conjg(-xi1st(l,i,j,k-cells3Offset))*xi1st(m,i,j,k-cells3Offset)
+          xiDyad_cmplx(l,m) = conjg(-xi1st(l,i,k,j-cells2Offset))*xi1st(m,i,k,j-cells2Offset)
         end do
         do concurrent(l = 1:3, m = 1:3)
-          temp33_cmplx(l,m) = sum(cmplx(C_ref(l,1:3,m,1:3),0.0_pReal)*xiDyad_cmplx)
+          temp33_cmplx(l,m) = sum(cmplx(C_ref(l,1:3,m,1:3),0.0_pReal,pReal)*xiDyad_cmplx)
         end do
 #else
         forall(l = 1:3, m = 1:3) &
-          xiDyad_cmplx(l,m) = conjg(-xi1st(l,i,j,k-cells3Offset))*xi1st(m,i,j,k-cells3Offset)
+          xiDyad_cmplx(l,m) = conjg(-xi1st(l,i,k,j-cells2Offset))*xi1st(m,i,k,j-cells2Offset)
         forall(l = 1:3, m = 1:3) &
-          temp33_cmplx(l,m) = sum(cmplx(C_ref(l,1:3,m,1:3),0.0_pReal)*xiDyad_cmplx)
+          temp33_cmplx(l,m) = sum(cmplx(C_ref(l,1:3,m,1:3),0.0_pReal,pReal)*xiDyad_cmplx)
 #endif
         A(1:3,1:3) = temp33_cmplx%re; A(4:6,4:6) =  temp33_cmplx%re
         A(1:3,4:6) = temp33_cmplx%im; A(4:6,1:3) = -temp33_cmplx%im
-        if (abs(math_det33(A(1:3,1:3))) > 1e-16) then
+        if (abs(math_det33(A(1:3,1:3))) > 1.e-16_pReal) then
           call math_invert(A_inv, err, A)
           temp33_cmplx = cmplx(A_inv(1:3,1:3),A_inv(1:3,4:6),pReal)
 #ifndef __INTEL_COMPILER
           do concurrent(l=1:3, m=1:3, n=1:3, o=1:3)
-            gamma_hat(l,m,n,o,i,j,k-cells3Offset) = temp33_cmplx(l,n) * xiDyad_cmplx(o,m)
+            gamma_hat(l,m,n,o,i,k,j-cells2Offset) = temp33_cmplx(l,n) * xiDyad_cmplx(o,m)
           end do
 #else
           forall(l=1:3, m=1:3, n=1:3, o=1:3) &
-            gamma_hat(l,m,n,o,i,j,k-cells3Offset) = temp33_cmplx(l,n) * xiDyad_cmplx(o,m)
+            gamma_hat(l,m,n,o,i,k,j-cells2Offset) = temp33_cmplx(l,n) * xiDyad_cmplx(o,m)
 #endif
         end if
       end if
@@ -509,24 +530,24 @@ subroutine utilities_fourierGammaConvolution(fieldAim)
 ! do the actual spectral method calculation (mechanical equilibrium)
   memoryEfficient: if (num%memory_efficient) then
     !$OMP PARALLEL DO PRIVATE(l,m,n,o,temp33_cmplx,xiDyad_cmplx,A,A_inv,err,gamma_hat)
-    do k = 1, cells3; do j = 1, cells(2); do i = 1, cells1Red
-      if (any([i,j,k+cells3Offset] /= 1)) then                                                      ! singular point at xi=(0.0,0.0,0.0) i.e. i=j=k=1
+    do j = 1, cells2; do k = 1, cells(3); do i = 1, cells1Red
+      if (any([i,j+cells2Offset,k] /= 1)) then                                                      ! singular point at xi=(0.0,0.0,0.0) i.e. i=j=k=1
 #ifndef __INTEL_COMPILER
         do concurrent(l = 1:3, m = 1:3)
-          xiDyad_cmplx(l,m) = conjg(-xi1st(l,i,j,k))*xi1st(m,i,j,k)
+          xiDyad_cmplx(l,m) = conjg(-xi1st(l,i,k,j))*xi1st(m,i,k,j)
         end do
         do concurrent(l = 1:3, m = 1:3)
-          temp33_cmplx(l,m) = sum(cmplx(C_ref(l,1:3,m,1:3),0.0_pReal)*xiDyad_cmplx)
+          temp33_cmplx(l,m) = sum(cmplx(C_ref(l,1:3,m,1:3),0.0_pReal,pReal)*xiDyad_cmplx)
         end do
 #else
         forall(l = 1:3, m = 1:3) &
-          xiDyad_cmplx(l,m) = conjg(-xi1st(l,i,j,k))*xi1st(m,i,j,k)
+          xiDyad_cmplx(l,m) = conjg(-xi1st(l,i,k,j))*xi1st(m,i,k,j)
         forall(l = 1:3, m = 1:3) &
-          temp33_cmplx(l,m) = sum(cmplx(C_ref(l,1:3,m,1:3),0.0_pReal)*xiDyad_cmplx)
+          temp33_cmplx(l,m) = sum(cmplx(C_ref(l,1:3,m,1:3),0.0_pReal,pReal)*xiDyad_cmplx)
 #endif
         A(1:3,1:3) = temp33_cmplx%re; A(4:6,4:6) =  temp33_cmplx%re
         A(1:3,4:6) = temp33_cmplx%im; A(4:6,1:3) = -temp33_cmplx%im
-        if (abs(math_det33(A(1:3,1:3))) > 1e-16) then
+        if (abs(math_det33(A(1:3,1:3))) > 1.e-16_pReal) then
           call math_invert(A_inv, err, A)
           temp33_cmplx = cmplx(A_inv(1:3,1:3),A_inv(1:3,4:6),pReal)
 #ifndef __INTEL_COMPILER
@@ -534,33 +555,33 @@ subroutine utilities_fourierGammaConvolution(fieldAim)
             gamma_hat(l,m,n,o,1,1,1) = temp33_cmplx(l,n)*xiDyad_cmplx(o,m)
           end do
           do concurrent(l = 1:3, m = 1:3)
-            temp33_cmplx(l,m) = sum(gamma_hat(l,m,1:3,1:3,1,1,1)*tensorField_fourier(1:3,1:3,i,j,k))
+            temp33_cmplx(l,m) = sum(gamma_hat(l,m,1:3,1:3,1,1,1)*tensorField_fourier(1:3,1:3,i,k,j))
           end do
 #else
           forall(l=1:3, m=1:3, n=1:3, o=1:3) &
             gamma_hat(l,m,n,o,1,1,1) = temp33_cmplx(l,n)*xiDyad_cmplx(o,m)
           forall(l = 1:3, m = 1:3) &
-            temp33_cmplx(l,m) = sum(gamma_hat(l,m,1:3,1:3,1,1,1)*tensorField_fourier(1:3,1:3,i,j,k))
+            temp33_cmplx(l,m) = sum(gamma_hat(l,m,1:3,1:3,1,1,1)*tensorField_fourier(1:3,1:3,i,k,j))
 #endif
-          tensorField_fourier(1:3,1:3,i,j,k) = temp33_cmplx
+          tensorField_fourier(1:3,1:3,i,k,j) = temp33_cmplx
         else
-          tensorField_fourier(1:3,1:3,i,j,k) = cmplx(0.0_pReal,0.0_pReal,pReal)
+          tensorField_fourier(1:3,1:3,i,k,j) = cmplx(0.0_pReal,0.0_pReal,pReal)
         end if
       end if
     end do; end do; end do
     !$OMP END PARALLEL DO
   else memoryEfficient
     !$OMP PARALLEL DO PRIVATE(l,m,temp33_cmplx)
-    do k = 1, cells3;  do j = 1, cells(2);  do i = 1,cells1Red
+    do j = 1, cells2;  do k = 1, cells(3);  do i = 1,cells1Red
 #ifndef __INTEL_COMPILER
       do concurrent(l = 1:3, m = 1:3)
-        temp33_cmplx(l,m) = sum(gamma_hat(l,m,1:3,1:3,i,j,k)*tensorField_fourier(1:3,1:3,i,j,k))
+        temp33_cmplx(l,m) = sum(gamma_hat(l,m,1:3,1:3,i,k,j)*tensorField_fourier(1:3,1:3,i,k,j))
       end do
 #else
       forall(l = 1:3, m = 1:3) &
-        temp33_cmplx(l,m) = sum(gamma_hat(l,m,1:3,1:3,i,j,k)*tensorField_fourier(1:3,1:3,i,j,k))
+        temp33_cmplx(l,m) = sum(gamma_hat(l,m,1:3,1:3,i,k,j)*tensorField_fourier(1:3,1:3,i,k,j))
 #endif
-      tensorField_fourier(1:3,1:3,i,j,k) = temp33_cmplx
+      tensorField_fourier(1:3,1:3,i,k,j) = temp33_cmplx
     end do; end do; end do
     !$OMP END PARALLEL DO
   end if memoryEfficient
@@ -583,11 +604,11 @@ subroutine utilities_fourierGreenConvolution(D_ref, mu_ref, Delta_t)
 !--------------------------------------------------------------------------------------------------
 ! do the actual spectral method calculation
   !$OMP PARALLEL DO PRIVATE(GreenOp_hat)
-  do k = 1, cells3; do j = 1, cells(2) ;do i = 1, cells1Red
+  do j = 1, cells2; do k = 1, cells(3); do i = 1, cells1Red
     GreenOp_hat = cmplx(1.0_pReal,0.0_pReal,pReal) &
-                / (cmplx(mu_ref,0.0_pReal,pReal) + cmplx(Delta_t,0.0_pReal) &
-                   * sum(conjg(xi1st(1:3,i,j,k))* matmul(cmplx(D_ref,0.0_pReal),xi1st(1:3,i,j,k))))
-    scalarField_fourier(i,j,k) = scalarField_fourier(i,j,k)*GreenOp_hat
+                / (cmplx(mu_ref,0.0_pReal,pReal) + cmplx(Delta_t,0.0_pReal,pReal) &
+                   * sum(conjg(xi1st(1:3,i,k,j))* matmul(cmplx(D_ref,0.0_pReal,pReal),xi1st(1:3,i,k,j))))
+    scalarField_fourier(i,k,j) = scalarField_fourier(i,k,j)*GreenOp_hat
   end do; end do; end do
   !$OMP END PARALLEL DO
 
@@ -603,31 +624,32 @@ real(pReal) function utilities_divergenceRMS()
   integer(MPI_INTEGER_KIND) :: err_MPI
   complex(pReal), dimension(3) :: rescaledGeom
 
+
   print'(/,1x,a)', '... calculating divergence ................................................'
   flush(IO_STDOUT)
 
-  rescaledGeom = cmplx(geomSize/scaledGeomSize,0.0_pReal)
+  rescaledGeom = cmplx(geomSize/scaledGeomSize,0.0_pReal,pReal)
 
 !--------------------------------------------------------------------------------------------------
 ! calculating RMS divergence criterion in Fourier space
   utilities_divergenceRMS = 0.0_pReal
-  do k = 1, cells3; do j = 1, cells(2)
+  do j = 1, cells2; do k = 1, cells(3)
     do i = 2, cells1Red -1                                                                          ! Has somewhere a conj. complex counterpart. Therefore count it twice.
       utilities_divergenceRMS = utilities_divergenceRMS &
-            + 2.0_pReal*(sum (real(matmul(tensorField_fourier(1:3,1:3,i,j,k), &                     ! (sqrt(real(a)**2 + aimag(a)**2))**2 = real(a)**2 + aimag(a)**2, i.e. do not take square root and square again
-                                          conjg(-xi1st(1:3,i,j,k))*rescaledGeom))**2) &             ! --> sum squared L_2 norm of vector
-                        +sum(aimag(matmul(tensorField_fourier(1:3,1:3,i,j,k),&
-                                          conjg(-xi1st(1:3,i,j,k))*rescaledGeom))**2))
+            + 2.0_pReal*(sum (real(matmul(tensorField_fourier(1:3,1:3,i,k,j), &                     ! (sqrt(real(a)**2 + aimag(a)**2))**2 = real(a)**2 + aimag(a)**2, i.e. do not take square root and square again
+                                          conjg(-xi1st(1:3,i,k,j))*rescaledGeom))**2) &             ! --> sum squared L_2 norm of vector
+                        +sum(aimag(matmul(tensorField_fourier(1:3,1:3,i,k,j),&
+                                          conjg(-xi1st(1:3,i,k,j))*rescaledGeom))**2))
     end do
     utilities_divergenceRMS = utilities_divergenceRMS &                                             ! these two layers (DC and Nyquist) do not have a conjugate complex counterpart (if cells(1) /= 1)
-               + sum( real(matmul(tensorField_fourier(1:3,1:3,1       ,j,k), &
-                                  conjg(-xi1st(1:3,1,j,k))*rescaledGeom))**2) &
-               + sum(aimag(matmul(tensorField_fourier(1:3,1:3,1       ,j,k), &
-                                  conjg(-xi1st(1:3,1,j,k))*rescaledGeom))**2) &
-               + sum( real(matmul(tensorField_fourier(1:3,1:3,cells1Red,j,k), &
-                                  conjg(-xi1st(1:3,cells1Red,j,k))*rescaledGeom))**2) &
-               + sum(aimag(matmul(tensorField_fourier(1:3,1:3,cells1Red,j,k), &
-                                  conjg(-xi1st(1:3,cells1Red,j,k))*rescaledGeom))**2)
+               + sum( real(matmul(tensorField_fourier(1:3,1:3,1       ,k,j), &
+                                  conjg(-xi1st(1:3,1,k,j))*rescaledGeom))**2) &
+               + sum(aimag(matmul(tensorField_fourier(1:3,1:3,1       ,k,j), &
+                                  conjg(-xi1st(1:3,1,k,j))*rescaledGeom))**2) &
+               + sum( real(matmul(tensorField_fourier(1:3,1:3,cells1Red,k,j), &
+                                  conjg(-xi1st(1:3,cells1Red,k,j))*rescaledGeom))**2) &
+               + sum(aimag(matmul(tensorField_fourier(1:3,1:3,cells1Red,k,j), &
+                                  conjg(-xi1st(1:3,cells1Red,k,j))*rescaledGeom))**2)
   end do; end do
   if (cells(1) == 1) utilities_divergenceRMS = utilities_divergenceRMS * 0.5_pReal                   ! counted twice in case of cells(1) == 1
   call MPI_Allreduce(MPI_IN_PLACE,utilities_divergenceRMS,1_MPI_INTEGER_KIND,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD,err_MPI)
@@ -650,42 +672,42 @@ real(pReal) function utilities_curlRMS()
   print'(/,1x,a)', '... calculating curl ......................................................'
   flush(IO_STDOUT)
 
-  rescaledGeom = cmplx(geomSize/scaledGeomSize,0.0_pReal)
+  rescaledGeom = cmplx(geomSize/scaledGeomSize,0.0_pReal,pReal)
 
 !--------------------------------------------------------------------------------------------------
 ! calculating max curl criterion in Fourier space
   utilities_curlRMS = 0.0_pReal
 
-  do k = 1, cells3; do j = 1, cells(2);
+  do j = 1, cells2; do k = 1, cells(3);
     do i = 2, cells1Red - 1
       do l = 1, 3
-        curl_fourier(l,1) = (+tensorField_fourier(l,3,i,j,k)*xi1st(2,i,j,k)*rescaledGeom(2) &
-                             -tensorField_fourier(l,2,i,j,k)*xi1st(3,i,j,k)*rescaledGeom(3))
-        curl_fourier(l,2) = (+tensorField_fourier(l,1,i,j,k)*xi1st(3,i,j,k)*rescaledGeom(3) &
-                             -tensorField_fourier(l,3,i,j,k)*xi1st(1,i,j,k)*rescaledGeom(1))
-        curl_fourier(l,3) = (+tensorField_fourier(l,2,i,j,k)*xi1st(1,i,j,k)*rescaledGeom(1) &
-                             -tensorField_fourier(l,1,i,j,k)*xi1st(2,i,j,k)*rescaledGeom(2))
+        curl_fourier(l,1) = (+tensorField_fourier(l,3,i,k,j)*xi1st(2,i,k,j)*rescaledGeom(2) &
+                             -tensorField_fourier(l,2,i,k,j)*xi1st(3,i,k,j)*rescaledGeom(3))
+        curl_fourier(l,2) = (+tensorField_fourier(l,1,i,k,j)*xi1st(3,i,k,j)*rescaledGeom(3) &
+                             -tensorField_fourier(l,3,i,k,j)*xi1st(1,i,k,j)*rescaledGeom(1))
+        curl_fourier(l,3) = (+tensorField_fourier(l,2,i,k,j)*xi1st(1,i,k,j)*rescaledGeom(1) &
+                             -tensorField_fourier(l,1,i,k,j)*xi1st(2,i,k,j)*rescaledGeom(2))
       end do
       utilities_curlRMS = utilities_curlRMS &
                         +2.0_pReal*sum(curl_fourier%re**2+curl_fourier%im**2)                       ! Has somewhere a conj. complex counterpart. Therefore count it twice.
     end do
     do l = 1, 3
-       curl_fourier = (+tensorField_fourier(l,3,1,j,k)*xi1st(2,1,j,k)*rescaledGeom(2) &
-                       -tensorField_fourier(l,2,1,j,k)*xi1st(3,1,j,k)*rescaledGeom(3))
-       curl_fourier = (+tensorField_fourier(l,1,1,j,k)*xi1st(3,1,j,k)*rescaledGeom(3) &
-                       -tensorField_fourier(l,3,1,j,k)*xi1st(1,1,j,k)*rescaledGeom(1))
-       curl_fourier = (+tensorField_fourier(l,2,1,j,k)*xi1st(1,1,j,k)*rescaledGeom(1) &
-                       -tensorField_fourier(l,1,1,j,k)*xi1st(2,1,j,k)*rescaledGeom(2))
+       curl_fourier = (+tensorField_fourier(l,3,1,k,j)*xi1st(2,1,k,j)*rescaledGeom(2) &
+                       -tensorField_fourier(l,2,1,k,j)*xi1st(3,1,k,j)*rescaledGeom(3))
+       curl_fourier = (+tensorField_fourier(l,1,1,k,j)*xi1st(3,1,k,j)*rescaledGeom(3) &
+                       -tensorField_fourier(l,3,1,k,j)*xi1st(1,1,k,j)*rescaledGeom(1))
+       curl_fourier = (+tensorField_fourier(l,2,1,k,j)*xi1st(1,1,k,j)*rescaledGeom(1) &
+                       -tensorField_fourier(l,1,1,k,j)*xi1st(2,1,k,j)*rescaledGeom(2))
     end do
     utilities_curlRMS = utilities_curlRMS &
                       + sum(curl_fourier%re**2 + curl_fourier%im**2)                                ! this layer (DC) does not have a conjugate complex counterpart (if cells(1) /= 1)
     do l = 1, 3
-      curl_fourier = (+tensorField_fourier(l,3,cells1Red,j,k)*xi1st(2,cells1Red,j,k)*rescaledGeom(2) &
-                      -tensorField_fourier(l,2,cells1Red,j,k)*xi1st(3,cells1Red,j,k)*rescaledGeom(3))
-      curl_fourier = (+tensorField_fourier(l,1,cells1Red,j,k)*xi1st(3,cells1Red,j,k)*rescaledGeom(3) &
-                      -tensorField_fourier(l,3,cells1Red,j,k)*xi1st(1,cells1Red,j,k)*rescaledGeom(1))
-      curl_fourier = (+tensorField_fourier(l,2,cells1Red,j,k)*xi1st(1,cells1Red,j,k)*rescaledGeom(1) &
-                      -tensorField_fourier(l,1,cells1Red,j,k)*xi1st(2,cells1Red,j,k)*rescaledGeom(2))
+      curl_fourier = (+tensorField_fourier(l,3,cells1Red,k,j)*xi1st(2,cells1Red,k,j)*rescaledGeom(2) &
+                      -tensorField_fourier(l,2,cells1Red,k,j)*xi1st(3,cells1Red,k,j)*rescaledGeom(3))
+      curl_fourier = (+tensorField_fourier(l,1,cells1Red,k,j)*xi1st(3,cells1Red,k,j)*rescaledGeom(3) &
+                      -tensorField_fourier(l,3,cells1Red,k,j)*xi1st(1,cells1Red,k,j)*rescaledGeom(1))
+      curl_fourier = (+tensorField_fourier(l,2,cells1Red,k,j)*xi1st(1,cells1Red,k,j)*rescaledGeom(1) &
+                      -tensorField_fourier(l,1,cells1Red,k,j)*xi1st(2,cells1Red,k,j)*rescaledGeom(2))
     end do
     utilities_curlRMS = utilities_curlRMS &
                       + sum(curl_fourier%re**2 + curl_fourier%im**2)                                ! this layer (Nyquist) does not have a conjugate complex counterpart (if cells(1) /= 1)
@@ -777,8 +799,8 @@ subroutine utilities_fourierScalarGradient()
   integer :: i, j, k
 
 
-  do k = 1, cells3;  do j = 1, cells(2);  do i = 1,cells1Red
-    vectorField_fourier(1:3,i,j,k) = scalarField_fourier(i,j,k)*xi1st(1:3,i,j,k)                    ! ToDo: no -conjg?
+  do j = 1, cells2;  do k = 1, cells(3);  do i = 1,cells1Red
+    vectorField_fourier(1:3,i,k,j) = scalarField_fourier(i,k,j)*xi1st(1:3,i,k,j)                    ! ToDo: no -conjg?
   end do; end do; end do
 
 end subroutine utilities_fourierScalarGradient
@@ -789,8 +811,7 @@ end subroutine utilities_fourierScalarGradient
 !--------------------------------------------------------------------------------------------------
 subroutine utilities_fourierVectorDivergence()
 
-
-  scalarField_fourier(1:cells1Red,1:cells(2),1:cells3) = sum(vectorField_fourier(1:3,1:cells1Red,1:cells(2),1:cells3) &
+  scalarField_fourier(1:cells1Red,1:cells(3),1:cells2) = sum(vectorField_fourier(1:3,1:cells1Red,1:cells(3),1:cells2) &
                                                              *conjg(-xi1st),1)
 
 end subroutine utilities_fourierVectorDivergence
@@ -803,10 +824,9 @@ subroutine utilities_fourierVectorGradient()
 
   integer :: i, j, k, m, n
 
-
-  do k = 1, cells3;  do j = 1, cells(2);  do i = 1,cells1Red
+  do j = 1, cells2;  do k = 1, cells(3);  do i = 1,cells1Red
     do m = 1, 3; do n = 1, 3
-      tensorField_fourier(m,n,i,j,k) = vectorField_fourier(m,i,j,k)*xi1st(n,i,j,k)
+      tensorField_fourier(m,n,i,k,j) = vectorField_fourier(m,i,k,j)*xi1st(n,i,k,j)
     end do; end do
   end do; end do; end do
 
@@ -820,9 +840,8 @@ subroutine utilities_fourierTensorDivergence()
 
   integer :: i, j, k
 
-
-  do k = 1, cells3;  do j = 1, cells(2);  do i = 1,cells1Red
-    vectorField_fourier(:,i,j,k) = matmul(tensorField_fourier(:,:,i,j,k),conjg(-xi1st(:,i,j,k)))
+  do j = 1, cells2;  do k = 1, cells(3);  do i = 1,cells1Red
+    vectorField_fourier(:,i,k,j) = matmul(tensorField_fourier(:,:,i,k,j),conjg(-xi1st(:,i,k,j)))
   end do; end do; end do
 
 end subroutine utilities_fourierTensorDivergence
@@ -946,19 +965,21 @@ function utilities_forwardField(Delta_t,field_lastInc,rate,aim)
     rate                                                                                            !< rate by which to forward
   real(pReal), intent(in), optional, dimension(3,3) :: &
     aim                                                                                             !< average field value aim
+
   real(pReal),                       dimension(3,3,cells(1),cells(2),cells3) :: &
     utilities_forwardField
   real(pReal),                       dimension(3,3) :: fieldDiff                                    !< <a + adot*t> - aim
   integer(MPI_INTEGER_KIND) :: err_MPI
 
+
   utilities_forwardField = field_lastInc + rate*Delta_t
   if (present(aim)) then                                                                            !< correct to match average
     fieldDiff = sum(sum(sum(utilities_forwardField,dim=5),dim=4),dim=3)*wgt
     call MPI_Allreduce(MPI_IN_PLACE,fieldDiff,9_MPI_INTEGER_KIND,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD,err_MPI)
     if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error'
     fieldDiff = fieldDiff - aim
-    utilities_forwardField = utilities_forwardField - &
-                     spread(spread(spread(fieldDiff,3,cells(1)),4,cells(2)),5,cells3)
+    utilities_forwardField = utilities_forwardField &
+                           - spread(spread(spread(fieldDiff,3,cells(1)),4,cells(2)),5,cells3)
   end if
 
 end function utilities_forwardField
@@ -972,8 +993,10 @@ end function utilities_forwardField
 pure function utilities_getFreqDerivative(k_s)
 
   integer, intent(in),  dimension(3) :: k_s                                                         !< indices of frequency
+
   complex(pReal),       dimension(3) :: utilities_getFreqDerivative
 
+
   select case (spectral_derivative_ID)
     case (DERIVATIVE_CONTINUOUS_ID)
       utilities_getFreqDerivative = cmplx(0.0_pReal, TAU*real(k_s,pReal)/geomSize,pReal)
@@ -1059,12 +1082,12 @@ subroutine utilities_updateCoords(F)
   call utilities_FFTtensorForward()
 
   !$OMP PARALLEL DO
-  do k = 1, cells3; do j = 1, cells(2); do i = 1, cells1Red
-    if (any([i,j,k+cells3Offset] /= 1)) then
-      vectorField_fourier(1:3,i,j,k) = matmul(tensorField_fourier(1:3,1:3,i,j,k),xi2nd(1:3,i,j,k)) &
-                                     / sum(conjg(-xi2nd(1:3,i,j,k))*xi2nd(1:3,i,j,k)) * cmplx(wgt,0.0,pReal)
+  do j = 1, cells2; do k = 1, cells(3); do i = 1, cells1Red
+    if (any([i,j+cells2Offset,k] /= 1)) then
+      vectorField_fourier(1:3,i,k,j) = matmul(tensorField_fourier(1:3,1:3,i,k,j),xi2nd(1:3,i,k,j)) &
+                                     / sum(conjg(-xi2nd(1:3,i,k,j))*xi2nd(1:3,i,k,j)) * cmplx(wgt,0.0,pReal)
     else
-      vectorField_fourier(1:3,i,j,k) = cmplx(0.0,0.0,pReal)
+      vectorField_fourier(1:3,i,k,j) = cmplx(0.0,0.0,pReal)
     end if
   end do; end do; end do
   !$OMP END PARALLEL DO
@@ -1073,7 +1096,7 @@ subroutine utilities_updateCoords(F)
 
  !--------------------------------------------------------------------------------------------------
  ! average F
-  if (cells3Offset == 0) Favg = real(tensorField_fourier(1:3,1:3,1,1,1),pReal)*wgt
+  if (cells3Offset == 0) Favg = tensorField_fourier(1:3,1:3,1,1,1)%re*wgt
   call MPI_Bcast(Favg,9_MPI_INTEGER_KIND,MPI_DOUBLE,0_MPI_INTEGER_KIND,MPI_COMM_WORLD,err_MPI)
   if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error'
 
@@ -1107,7 +1130,7 @@ subroutine utilities_updateCoords(F)
  !--------------------------------------------------------------------------------------------------
  ! calculate nodal displacements
   nodeCoords = 0.0_pReal
-  do k = 0,cells3; do j = 0,cells(2); do i = 0,cells(1)
+  do j = 0,cells(2); do k = 0,cells3; do i = 0,cells(1)
     nodeCoords(1:3,i+1,j+1,k+1) = matmul(Favg,step*(real([i,j,k+cells3Offset],pReal)))
     averageFluct: do n = 1,8
       me = [i+neighbor(1,n),j+neighbor(2,n),k+neighbor(3,n)]
@@ -1137,6 +1160,7 @@ subroutine utilities_saveReferenceStiffness
   integer :: &
     fileUnit,ierr
 
+
   if (worldrank == 0) then
     print'(/,1x,a)', '... writing reference stiffness data required for restart to file .........'; flush(IO_STDOUT)
     open(newunit=fileUnit, file=getSolverJobName()//'.C_ref',&
@@ -1163,6 +1187,10 @@ subroutine selfTest()
   tensorField_real(1:3,1:3,cells(1)+1:cells1Red*2,:,:) = 0.0_pReal
   tensorField_real_ = tensorField_real
   call utilities_FFTtensorForward()
+  if (worldsize==1) then
+    if (any(dNeq(sum(sum(sum(tensorField_real_,dim=5),dim=4),dim=3)/tensorField_fourier(:,:,1,1,1)%re,1.0_pReal,1.0e-12_pReal))) &
+      error stop 'tensorField avg'
+  endif
   call utilities_FFTtensorBackward()
   tensorField_real(1:3,1:3,cells(1)+1:cells1Red*2,:,:) = 0.0_pReal
   if (maxval(abs(tensorField_real_ - tensorField_real))>5.0e-15_pReal) error stop 'tensorField'
@@ -1171,6 +1199,10 @@ subroutine selfTest()
   vectorField_real(1:3,cells(1)+1:cells1Red*2,:,:) = 0.0_pReal
   vectorField_real_ = vectorField_real
   call utilities_FFTvectorForward()
+  if (worldsize==1) then
+    if (any(dNeq(sum(sum(sum(vectorField_real_,dim=4),dim=3),dim=2)/vectorField_fourier(:,1,1,1)%re,1.0_pReal,1.0e-12_pReal))) &
+      error stop 'vector avg'
+  endif
   call utilities_FFTvectorBackward()
   vectorField_real(1:3,cells(1)+1:cells1Red*2,:,:) = 0.0_pReal
   if (maxval(abs(vectorField_real_ - vectorField_real))>5.0e-15_pReal) error stop 'vectorField'
@@ -1179,6 +1211,10 @@ subroutine selfTest()
   scalarField_real(cells(1)+1:cells1Red*2,:,:) = 0.0_pReal
   scalarField_real_ = scalarField_real
   call utilities_FFTscalarForward()
+  if (worldsize==1) then
+    if (dNeq(sum(sum(sum(scalarField_real_,dim=3),dim=2),dim=1)/scalarField_fourier(1,1,1)%re,1.0_pReal,1.0e-12_pReal)) &
+      error stop 'scalar avg'
+  endif
   call utilities_FFTscalarBackward()
   scalarField_real(cells(1)+1:cells1Red*2,:,:) = 0.0_pReal
   if (maxval(abs(scalarField_real_ - scalarField_real))>5.0e-15_pReal) error stop 'scalarField'

src/grid/zlib.f90

@@ -5,7 +5,7 @@
 module zlib
   use prec
 
-  implicit none
+  implicit none(type,external)
   private
 
   public :: &
@@ -14,13 +14,12 @@ module zlib
   interface
 
     subroutine inflate_C(s_deflated,s_inflated,deflated,inflated) bind(C)
-    use, intrinsic :: ISO_C_Binding, only: &
-      C_SIGNED_CHAR, C_INT64_T
+      use, intrinsic :: ISO_C_Binding, only: C_SIGNED_CHAR, C_INT64_T
+      implicit none(type,external)
 
       integer(C_INT64_T),                            intent(in)  :: s_deflated,s_inflated
       integer(C_SIGNED_CHAR), dimension(s_deflated), intent(in)  :: deflated
       integer(C_SIGNED_CHAR), dimension(s_inflated), intent(out) :: inflated
-
     end subroutine inflate_C
 
   end interface
@@ -37,6 +36,7 @@ function zlib_inflate(deflated,size_inflated)
 
   integer(C_SIGNED_CHAR), dimension(size_inflated) :: zlib_inflate
 
+
   call inflate_C(size(deflated,kind=C_INT64_T),int(size_inflated,C_INT64_T),deflated,zlib_inflate)
 
 end function zlib_inflate

src/homogenization.f90

@@ -18,7 +18,7 @@ module homogenization
   use results
   use lattice
 
-  implicit none
+  implicit none(type,external)
   private
 
   type :: tState

src/homogenization_damage.f90

@@ -3,8 +3,6 @@
 !--------------------------------------------------------------------------------------------------
 submodule(homogenization) damage
 
-  use lattice
-
   interface
 
     module subroutine pass_init
@@ -79,9 +77,10 @@ end subroutine damage_init
 !--------------------------------------------------------------------------------------------------
 module subroutine damage_partition(ce)
 
-  real(pReal) :: phi
   integer, intent(in) :: ce
 
+  real(pReal) :: phi
+
 
   if(damageState_h(material_homogenizationID(ce))%sizeState < 1) return
   phi = damagestate_h(material_homogenizationID(ce))%state(1,material_homogenizationEntry(ce))
@@ -91,7 +90,7 @@ end subroutine damage_partition
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief Homogenized damage viscosity.
+!> @brief Homogenize damage viscosity.
 !--------------------------------------------------------------------------------------------------
 module function homogenization_mu_phi(ce) result(mu)
 
@@ -105,7 +104,7 @@ end function homogenization_mu_phi
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief Homogenized damage conductivity/diffusivity in reference configuration.
+!> @brief Homogenize damage conductivity.
 !--------------------------------------------------------------------------------------------------
 module function homogenization_K_phi(ce) result(K)
 
@@ -119,13 +118,12 @@ end function homogenization_K_phi
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief Homogenized damage driving force.
+!> @brief Homogenize damage driving force.
 !--------------------------------------------------------------------------------------------------
 module function homogenization_f_phi(phi,ce) result(f)
 
   integer, intent(in) :: ce
-  real(pReal), intent(in) :: &
-    phi
+  real(pReal), intent(in) :: phi
   real(pReal) :: f
 
 
@@ -140,8 +138,7 @@ end function homogenization_f_phi
 module subroutine homogenization_set_phi(phi,ce)
 
   integer, intent(in) :: ce
-  real(pReal),   intent(in) :: &
-    phi
+  real(pReal), intent(in) :: phi
 
   integer :: &
     ho, &
@@ -166,6 +163,7 @@ module subroutine damage_results(ho,group)
 
   integer :: o
 
+
   associate(prm => param(ho))
       outputsLoop: do o = 1,size(prm%output)
         select case(prm%output(o))

src/homogenization_mechanical_RGC.f90

@@ -601,9 +601,9 @@ module function RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHappy)
     ! calculate the stress and penalty due to volume discrepancy
     vPen      = 0.0_pReal
     do i = 1,nGrain
-      vPen(:,:,i) = -1.0_pReal/real(nGrain,pReal)*num%volDiscrMod*num%volDiscrPow/num%maxVolDiscr* &
-                         sign((abs(vDiscrep)/num%maxVolDiscr)**(num%volDiscrPow - 1.0),vDiscrep)* &
-                         gVol(i)*transpose(math_inv33(fDef(:,:,i)))
+      vPen(:,:,i) = -real(nGrain,pReal)**(-1)*num%volDiscrMod*num%volDiscrPow/num%maxVolDiscr &
+                  * sign((abs(vDiscrep)/num%maxVolDiscr)**(num%volDiscrPow - 1.0_pReal),vDiscrep) &
+                  * gVol(i)*transpose(math_inv33(fDef(:,:,i)))
     end do
 
   end subroutine volumePenalty

src/homogenization_thermal.f90

@@ -3,8 +3,6 @@
 !--------------------------------------------------------------------------------------------------
 submodule(homogenization) thermal
 
-  use lattice
-
   interface
 
     module subroutine pass_init
@@ -105,7 +103,7 @@ end subroutine thermal_partition
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief Homogenized thermal viscosity.
+!> @brief Homogenize thermal viscosity.
 !--------------------------------------------------------------------------------------------------
 module function homogenization_mu_T(ce) result(mu)
 
@@ -124,7 +122,7 @@ end function homogenization_mu_T
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief Homogenized thermal conductivity in reference configuration.
+!> @brief Homogenize thermal conductivity.
 !--------------------------------------------------------------------------------------------------
 module function homogenization_K_T(ce) result(K)
 
@@ -143,7 +141,7 @@ end function homogenization_K_T
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief Homogenized heat generation rate.
+!> @brief Homogenize heat generation rate.
 !--------------------------------------------------------------------------------------------------
 module function homogenization_f_T(ce) result(f)
 
@@ -187,6 +185,7 @@ module subroutine thermal_results(ho,group)
 
   integer :: o
 
+
   associate(prm => param(ho))
     outputsLoop: do o = 1,size(prm%output)
       select case(trim(prm%output(o)))

src/lattice.f90

@@ -13,7 +13,7 @@ module lattice
   use math
   use rotations
 
-  implicit none
+  implicit none(type,external)
   private
 
 !--------------------------------------------------------------------------------------------------
@@ -576,7 +576,7 @@ function lattice_C66_trans(Ntrans,C_parent66,lattice_target, &
 
   do i = 1,6
     if (abs(C_target_unrotated66(i,i))<tol_math_check) &
-    call IO_error(135,el=i,ext_msg='matrix diagonal "el"ement in transformation')
+    call IO_error(135,'matrix diagonal in transformation',label1='entry',ID1=i)
   end do
 
   call buildTransformationSystem(Q,S,Ntrans,cOverA_trans,a_cF,a_cI)
@@ -1431,7 +1431,7 @@ function lattice_SchmidMatrix_slip(Nslip,lattice,cOverA) result(SchmidMatrix)
   do i = 1, sum(Nslip)
     SchmidMatrix(1:3,1:3,i) = math_outer(coordinateSystem(1:3,1,i),coordinateSystem(1:3,2,i))
     if (abs(math_trace33(SchmidMatrix(1:3,1:3,i))) > tol_math_check) &
-      call IO_error(0,i,ext_msg = 'dilatational Schmid matrix for slip')
+      error stop 'dilatational Schmid matrix for slip'
   end do
 
 end function lattice_SchmidMatrix_slip
@@ -1478,7 +1478,7 @@ function lattice_SchmidMatrix_twin(Ntwin,lattice,cOverA) result(SchmidMatrix)
   do i = 1, sum(Ntwin)
     SchmidMatrix(1:3,1:3,i) = math_outer(coordinateSystem(1:3,1,i),coordinateSystem(1:3,2,i))
     if (abs(math_trace33(SchmidMatrix(1:3,1:3,i))) > tol_math_check) &
-      call IO_error(0,i,ext_msg = 'dilatational Schmid matrix for twin')
+      error stop 'dilatational Schmid matrix for twin'
   end do
 
 end function lattice_SchmidMatrix_twin
@@ -2008,7 +2008,7 @@ subroutine buildTransformationSystem(Q,S,Ntrans,cOverA,a_cF,a_cI)
       ],pReal),shape(CFTOHP_SYSTEMTRANS))
 
   real(pReal), dimension(4,cF_Ntrans), parameter :: &
-    CFTOCI_SYSTEMTRANS = reshape([&
+    CFTOCI_SYSTEMTRANS = real(reshape([&
       0.0, 1.0, 0.0,     10.26, &                                                                   ! Pitsch OR (Ma & Hartmaier 2014, Table 3)
       0.0,-1.0, 0.0,     10.26, &
       0.0, 0.0, 1.0,     10.26, &
@@ -2021,7 +2021,7 @@ subroutine buildTransformationSystem(Q,S,Ntrans,cOverA,a_cF,a_cI)
      -1.0, 0.0, 0.0,     10.26, &
       0.0, 1.0, 0.0,     10.26, &
       0.0,-1.0, 0.0,     10.26  &
-      ],shape(CFTOCI_SYSTEMTRANS))
+      ],shape(CFTOCI_SYSTEMTRANS)),pReal)
 
   integer, dimension(9,cF_Ntrans), parameter :: &
     CFTOCI_BAINVARIANT = reshape( [&
@@ -2040,7 +2040,7 @@ subroutine buildTransformationSystem(Q,S,Ntrans,cOverA,a_cF,a_cI)
       ],shape(CFTOCI_BAINVARIANT))
 
   real(pReal), dimension(4,cF_Ntrans), parameter :: &
-    CFTOCI_BAINROT = reshape([&
+    CFTOCI_BAINROT = real(reshape([&
       1.0, 0.0, 0.0,     45.0, &                                                                    ! Rotate cF austensite to bain variant
       1.0, 0.0, 0.0,     45.0, &
       1.0, 0.0, 0.0,     45.0, &
@@ -2053,7 +2053,7 @@ subroutine buildTransformationSystem(Q,S,Ntrans,cOverA,a_cF,a_cI)
       0.0, 0.0, 1.0,     45.0, &
       0.0, 0.0, 1.0,     45.0, &
       0.0, 0.0, 1.0,     45.0  &
-      ],shape(CFTOCI_BAINROT))
+      ],shape(CFTOCI_BAINROT)),pReal)
 
   if (present(a_cI) .and. present(a_cF)) then
     do i = 1,sum(Ntrans)

src/material.f90

@@ -14,7 +14,7 @@ module material
   use discretization
   use YAML_types
 
-  implicit none
+  implicit none(type,external)
   private
 
   type, public :: tRotationContainer

src/materialpoint.f90

@@ -31,7 +31,7 @@ module materialpoint
   use discretization_grid
 #endif
 
-  implicit none
+  implicit none(type,external)
   public
 
 contains

src/math.f90

@@ -12,7 +12,7 @@ module math
   use YAML_types
   use LAPACK_interface
 
-  implicit none
+  implicit none(type,external)
   public
 #if __INTEL_COMPILER >= 1900
   ! do not make use of associated entities available to other modules
@@ -1110,7 +1110,7 @@ pure function math_rotationalPart(F) result(R)
 
   C = matmul(transpose(F),F)
   I_C = math_invariantsSym33(C)
-  I_F = [math_trace33(F), 0.5*(math_trace33(F)**2 - math_trace33(matmul(F,F)))]
+  I_F = [math_trace33(F), 0.5_pReal*(math_trace33(F)**2 - math_trace33(matmul(F,F)))]
 
   x = math_clip(I_C(1)**2 -3.0_pReal*I_C(2),0.0_pReal)**(3.0_pReal/2.0_pReal)
   if (dNeq0(x)) then
@@ -1129,7 +1129,7 @@ pure function math_rotationalPart(F) result(R)
     - I_U(1)*I_F(1) * transpose(F) &
     + I_U(1) * transpose(matmul(F,F)) &
     - matmul(F,C)
-  R = R /(I_U(1)*I_U(2)-I_U(3))
+  R = R*math_det33(R)**(-1.0_pReal/3.0_pReal)
 
 end function math_rotationalPart
 
@@ -1386,7 +1386,7 @@ subroutine selfTest()
   call random_number(v3_3)
   call random_number(v3_4)
 
-  if (dNeq(abs(dot_product(math_cross(v3_1-v3_4,v3_2-v3_4),v3_3-v3_4))/6.0, &
+  if (dNeq(abs(dot_product(math_cross(v3_1-v3_4,v3_2-v3_4),v3_3-v3_4))/6.0_pReal, &
           math_volTetrahedron(v3_1,v3_2,v3_3,v3_4),tol=1.0e-12_pReal)) &
     error stop 'math_volTetrahedron'
 
@@ -1422,7 +1422,9 @@ subroutine selfTest()
   t33_2 = math_rotationalPart(transpose(t33))
   t33   = math_rotationalPart(t33)
   if (any(dNeq0(matmul(t33_2,t33) - math_I3,tol=1.0e-10_pReal))) &
-    error stop 'math_rotationalPart'
+    error stop 'math_rotationalPart (forward-backward)'
+  if (dNeq(1.0_pReal,math_det33(math_rotationalPart(t33)),tol=1.0e-10_pReal)) &
+    error stop 'math_rotationalPart (determinant)'
 
   call random_number(r)
   d = int(r*5.0_pReal) + 1

src/mesh/DAMASK_mesh.f90

@@ -20,7 +20,7 @@ program DAMASK_mesh
   use FEM_Utilities
   use mesh_mechanical_FEM
 
-  implicit none
+  implicit none(type,external)
 
   type :: tLoadCase
     real(pReal)  :: time                   = 0.0_pReal                                              !< length of increment

src/mesh/FEM_quadrature.f90

@@ -5,7 +5,7 @@
 module FEM_quadrature
   use prec
 
-  implicit none
+  implicit none(type,external)
   private
 
   integer, parameter :: &

src/mesh/FEM_utilities.f90

@@ -21,7 +21,11 @@ module FEM_utilities
   use homogenization
   use FEM_quadrature
 
+#if (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR>14) && !defined(PETSC_HAVE_MPI_F90MODULE_VISIBILITY)
+  implicit none(type,external)
+#else
   implicit none
+#endif
   private
 
   logical,     public             :: cutBack = .false.                                              !< cut back of BVP solver in case convergence is not achieved or a material point is terminally ill
@@ -65,6 +69,11 @@ module FEM_utilities
     type(tComponentBC), allocatable, dimension(:) :: componentBC
   end type tFieldBC
 
+  external :: &                                                                                     ! ToDo: write interfaces
+    PetscSectionGetFieldComponents, &
+    PetscSectionGetFieldDof, &
+    PetscSectionGetFieldOffset
+
   public :: &
     FEM_utilities_init, &
     utilities_constitutiveResponse, &
@@ -131,7 +140,7 @@ subroutine FEM_utilities_init
   call PetscOptionsInsertString(PETSC_NULL_OPTIONS,trim(petsc_optionsOrder),err_PETSc)
   CHKERRQ(err_PETSc)
 
-  wgt = 1.0/real(mesh_maxNips*mesh_NcpElemsGlobal,pReal)
+  wgt = real(mesh_maxNips*mesh_NcpElemsGlobal,pReal)**(-1)
 
 
 end subroutine FEM_utilities_init

src/mesh/discretization_mesh.f90

@@ -25,7 +25,11 @@ module discretization_mesh
   use YAML_types
   use prec
 
+#if (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR>14) && !defined(PETSC_HAVE_MPI_F90MODULE_VISIBILITY)
+  implicit none(type,external)
+#else
   implicit none
+#endif
   private
 
   integer, public, protected :: &
@@ -52,6 +56,11 @@ module discretization_mesh
   real(pReal), dimension(:,:,:), allocatable :: &
     mesh_ipCoordinates                                                                              !< IP x,y,z coordinates (after deformation!)
 
+  external :: &
+#ifdef PETSC_USE_64BIT_INDICES
+    DMDestroy, &
+#endif
+    DMView                                                                                          ! ToDo: write interface
   public :: &
     discretization_mesh_init, &
     mesh_FEM_build_ipVolumes, &
@@ -242,12 +251,12 @@ subroutine mesh_FEM_build_ipCoordinates(dimPlex,qPoints)
     call DMPlexComputeCellGeometryAffineFEM(geomMesh,cell,pV0,pCellJ,pInvcellJ,detJ,err_PETSc)
     CHKERRQ(err_PETSc)
     qOffset = 0
-    do qPt = 1, mesh_maxNips
-      do dirI = 1, dimPlex
+    do qPt = 1_pPETSCINT, mesh_maxNips
+      do dirI = 1_pPETSCINT, dimPlex
         mesh_ipCoordinates(dirI,qPt,cell+1) = pV0(dirI)
-        do dirJ = 1, dimPlex
+        do dirJ = 1_pPETSCINT, dimPlex
           mesh_ipCoordinates(dirI,qPt,cell+1) = mesh_ipCoordinates(dirI,qPt,cell+1) + &
-                                                pCellJ((dirI-1)*dimPlex+dirJ)*(qPoints(qOffset+dirJ) + 1.0)
+                                                pCellJ((dirI-1)*dimPlex+dirJ)*(qPoints(qOffset+dirJ) + 1.0_pReal)
         enddo
       enddo
       qOffset = qOffset + dimPlex

src/mesh/mesh_mech_FEM.f90

@@ -26,7 +26,11 @@ module mesh_mechanical_FEM
   use homogenization
   use math
 
+#if (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR>14) && !defined(PETSC_HAVE_MPI_F90MODULE_VISIBILITY)
+  implicit none(type,external)
+#else
   implicit none
+#endif
   private
 
 !--------------------------------------------------------------------------------------------------
@@ -67,6 +71,18 @@ module mesh_mechanical_FEM
   logical :: ForwardData
   real(pReal), parameter :: eps = 1.0e-18_pReal
 
+  external :: &                                                                                     ! ToDo: write interfaces
+#ifdef PETSC_USE_64BIT_INDICES
+    ISDestroy, &
+#endif
+    PetscSectionGetNumFields, &
+    PetscFESetQuadrature, &
+    PetscFEGetDimension, &
+    PetscFEDestroy, &
+    PetscSectionGetDof, &
+    PetscFEGetDualSpace, &
+    PetscDualSpaceGetFunctional
+
   public :: &
     FEM_mechanical_init, &
     FEM_mechanical_solution, &
@@ -230,14 +246,14 @@ subroutine FEM_mechanical_init(fieldBC)
   CHKERRQ(err_PETSc)
   call SNESSetConvergenceTest(mechanical_snes,FEM_mechanical_converged,PETSC_NULL_VEC,PETSC_NULL_FUNCTION,err_PETSc)
   CHKERRQ(err_PETSc)
-  call SNESSetTolerances(mechanical_snes,1.0,0.0,0.0,num%itmax,num%itmax,err_PETSc)
+  call SNESSetTolerances(mechanical_snes,1.0_pReal,0.0_pReal,0.0_pReal,num%itmax,num%itmax,err_PETSc)
   CHKERRQ(err_PETSc)
   call SNESSetFromOptions(mechanical_snes,err_PETSc); CHKERRQ(err_PETSc)
 
 !--------------------------------------------------------------------------------------------------
 ! init fields
-  call VecSet(solution        ,0.0,err_PETSc); CHKERRQ(err_PETSc)
-  call VecSet(solution_rate   ,0.0,err_PETSc); CHKERRQ(err_PETSc)
+  call VecSet(solution     ,0.0_pReal,err_PETSc); CHKERRQ(err_PETSc)
+  call VecSet(solution_rate,0.0_pReal,err_PETSc); CHKERRQ(err_PETSc)
   allocate(x_scal(cellDof))
   allocate(nodalWeightsP(1))
   allocate(nodalPointsP(dimPlex))
@@ -338,10 +354,9 @@ subroutine FEM_mechanical_formResidual(dm_local,xx_local,f_local,dummy,err_PETSc
   PetscInt                           :: cellStart, cellEnd, cell, field, face, &
                                         qPt, basis, comp, cidx, &
                                         numFields, &
-                                        bcSize,m
+                                        bcSize,m,i
   PetscReal                          :: detFAvg, detJ
   PetscReal, dimension(dimPlex*dimPlex,cellDof) :: BMat
-
   IS :: bcPoints
 
 
@@ -358,9 +373,9 @@ subroutine FEM_mechanical_formResidual(dm_local,xx_local,f_local,dummy,err_PETSc
   CHKERRQ(err_PETSc)
   call DMGetLocalVector(dm_local,x_local,err_PETSc)
   CHKERRQ(err_PETSc)
-  call VecWAXPY(x_local,1.0,xx_local,solution_local,err_PETSc)
+  call VecWAXPY(x_local,1.0_pReal,xx_local,solution_local,err_PETSc)
   CHKERRQ(err_PETSc)
-  do field = 1, dimPlex; do face = 1, mesh_Nboundaries
+  do field = 1_pPETSCINT, dimPlex; do face = 1_pPETSCINT, mesh_Nboundaries
     if (params%fieldBC%componentBC(field)%Mask(face)) then
       call DMGetStratumSize(dm_local,'Face Sets',mesh_boundaries(face),bcSize,err_PETSc)
       if (bcSize > 0) then
@@ -375,7 +390,7 @@ subroutine FEM_mechanical_formResidual(dm_local,xx_local,f_local,dummy,err_PETSc
 
 !--------------------------------------------------------------------------------------------------
 ! evaluate field derivatives
-  do cell = cellStart, cellEnd-1                                                                    !< loop over all elements
+  do cell = cellStart, cellEnd-1_pPETSCINT                                                          !< loop over all elements
 
     call PetscSectionGetNumFields(section,numFields,err_PETSc)
     CHKERRQ(err_PETSc)
@@ -384,25 +399,25 @@ subroutine FEM_mechanical_formResidual(dm_local,xx_local,f_local,dummy,err_PETSc
     call  DMPlexComputeCellGeometryAffineFEM(dm_local,cell,pV0,pCellJ,pInvcellJ,detJ,err_PETSc)
     CHKERRQ(err_PETSc)
     IcellJMat = reshape(pInvcellJ,shape=[dimPlex,dimPlex])
-    do qPt = 0, nQuadrature-1
-      m = cell*nQuadrature + qPt+1
-      BMat = 0.0
-      do basis = 0, nBasis-1
-        do comp = 0, dimPlex-1
+    do qPt = 0_pPETSCINT, nQuadrature-1_pPETSCINT
+      m = cell*nQuadrature + qPt+1_pPETSCINT
+      BMat = 0.0_pReal
+      do basis = 0_pPETSCINT, nBasis-1_pPETSCINT
+        do comp = 0_pPETSCINT, dimPlex-1_pPETSCINT
           cidx = basis*dimPlex+comp
-          BMat(comp*dimPlex+1:(comp+1)*dimPlex,basis*dimPlex+comp+1) = &
-            matmul(IcellJMat,basisFieldDer((((qPt*nBasis + basis)*dimPlex + comp)*dimPlex+comp  )*dimPlex+1: &
-                                           (((qPt*nBasis + basis)*dimPlex + comp)*dimPlex+comp+1)*dimPlex))
+          i = ((qPt*nBasis + basis)*dimPlex + comp)*dimPlex+comp
+          BMat(comp*dimPlex+1_pPETSCINT:(comp+1_pPETSCINT)*dimPlex,basis*dimPlex+comp+1_pPETSCINT) = &
+            matmul(IcellJMat,basisFieldDer(i*dimPlex+1_pPETSCINT:(i+1_pPETSCINT)*dimPlex))
         enddo
       enddo
       homogenization_F(1:dimPlex,1:dimPlex,m) = reshape(matmul(BMat,x_scal),shape=[dimPlex,dimPlex], order=[2,1])
     enddo
     if (num%BBarStabilisation) then
-      detFAvg = math_det33(sum(homogenization_F(1:3,1:3,cell*nQuadrature+1:(cell+1)*nQuadrature),dim=3)/real(nQuadrature))
+      detFAvg = math_det33(sum(homogenization_F(1:3,1:3,cell*nQuadrature+1:(cell+1)*nQuadrature),dim=3)/real(nQuadrature,pReal))
       do qPt = 0, nQuadrature-1
         m = cell*nQuadrature + qPt+1
         homogenization_F(1:dimPlex,1:dimPlex,m) = homogenization_F(1:dimPlex,1:dimPlex,m) &
-                                                * (detFAvg/math_det33(homogenization_F(1:3,1:3,m)))**(1.0/real(dimPlex))
+                                                * (detFAvg/math_det33(homogenization_F(1:3,1:3,m)))**(1.0_pReal/real(dimPlex,pReal))
 
       enddo
     endif
@@ -425,22 +440,22 @@ subroutine FEM_mechanical_formResidual(dm_local,xx_local,f_local,dummy,err_PETSc
     call  DMPlexComputeCellGeometryAffineFEM(dm_local,cell,pV0,pCellJ,pInvcellJ,detJ,err_PETSc)
     CHKERRQ(err_PETSc)
     IcellJMat = reshape(pInvcellJ,shape=[dimPlex,dimPlex])
-    f_scal = 0.0
-    do qPt = 0, nQuadrature-1
-      m = cell*nQuadrature + qPt+1
-      BMat = 0.0
-      do basis = 0, nBasis-1
-        do comp = 0, dimPlex-1
+    f_scal = 0.0_pReal
+    do qPt = 0_pPETSCINT, nQuadrature-1_pPETSCINT
+      m = cell*nQuadrature + qPt+1_pPETSCINT
+      BMat = 0.0_pReal
+      do basis = 0_pPETSCINT, nBasis-1_pPETSCINT
+        do comp = 0_pPETSCINT, dimPlex-1_pPETSCINT
           cidx = basis*dimPlex+comp
-          BMat(comp*dimPlex+1:(comp+1)*dimPlex,basis*dimPlex+comp+1) = &
-            matmul(IcellJMat,basisFieldDer((((qPt*nBasis + basis)*dimPlex + comp)*dimPlex+comp  )*dimPlex+1: &
-                                           (((qPt*nBasis + basis)*dimPlex + comp)*dimPlex+comp+1)*dimPlex))
+          i = ((qPt*nBasis + basis)*dimPlex + comp)*dimPlex+comp
+          BMat(comp*dimPlex+1_pPETSCINT:(comp+1_pPETSCINT)*dimPlex,basis*dimPlex+comp+1_pPETSCINT) = &
+            matmul(IcellJMat,basisFieldDer(i*dimPlex+1_pPETSCINT:(i+1_pPETSCINT)*dimPlex))
         enddo
       enddo
-      f_scal = f_scal + &
-               matmul(transpose(BMat), &
+      f_scal = f_scal &
+             +  matmul(transpose(BMat), &
                        reshape(transpose(homogenization_P(1:dimPlex,1:dimPlex,m)), &
-                              shape=[dimPlex*dimPlex]))*qWeights(qPt+1)
+                               shape=[dimPlex*dimPlex]))*qWeights(qPt+1_pPETSCINT)
     enddo
     f_scal = f_scal*abs(detJ)
     pf_scal => f_scal
@@ -467,7 +482,6 @@ subroutine FEM_mechanical_formJacobian(dm_local,xx_local,Jac_pre,Jac,dummy,err_P
 
   PetscDS      :: prob
   Vec          :: x_local, xx_local
-
   PetscSection :: section, gSection
 
   PetscReal, dimension(1,         cellDof)  :: MatB
@@ -480,12 +494,10 @@ subroutine FEM_mechanical_formJacobian(dm_local,xx_local,Jac_pre,Jac,dummy,err_P
   PetscScalar, dimension(:),   pointer :: pK_e, x_scal
 
   PetscScalar,dimension(cellDOF,cellDOF),  target :: K_e
-  PetscScalar,dimension(cellDOF,cellDOF)  :: K_eA  , &
-                                          K_eB
+  PetscScalar,dimension(cellDOF,cellDOF) :: K_eA, K_eB
 
   PetscInt :: cellStart, cellEnd, cell, field, face, &
-                                          qPt, basis, comp, cidx,bcSize, m
-
+              qPt, basis, comp, cidx,bcSize, m, i
   IS :: bcPoints
 
 
@@ -530,30 +542,29 @@ subroutine FEM_mechanical_formJacobian(dm_local,xx_local,Jac_pre,Jac,dummy,err_P
     CHKERRQ(err_PETSc)
     call  DMPlexComputeCellGeometryAffineFEM(dm_local,cell,pV0,pCellJ,pInvcellJ,detJ,err_PETSc)
     CHKERRQ(err_PETSc)
-    K_eA = 0.0
-    K_eB = 0.0
-    MatB = 0.0
-    FAvg = 0.0
-    BMatAvg = 0.0
-    do qPt = 0, nQuadrature-1
-      m = cell*nQuadrature + qPt + 1
-      BMat = 0.0
-      do basis = 0, nBasis-1
-        do comp = 0, dimPlex-1
+    K_eA = 0.0_pReal
+    K_eB = 0.0_pReal
+    MatB = 0.0_pReal
+    FAvg = 0.0_pReal
+    BMatAvg = 0.0_pReal
+    do qPt = 0_pPETSCINT, nQuadrature-1_pPETSCINT
+      m = cell*nQuadrature + qPt + 1_pPETSCINT
+      BMat = 0.0_pReal
+      do basis = 0_pPETSCINT, nBasis-1_pPETSCINT
+        do comp = 0_pPETSCINT, dimPlex-1_pPETSCINT
           cidx = basis*dimPlex+comp
-          BMat(comp*dimPlex+1:(comp+1)*dimPlex,basis*dimPlex+comp+1) = &
-            matmul( reshape(pInvcellJ, shape = [dimPlex,dimPlex]),&
-                             basisFieldDer((((qPt*nBasis + basis)*dimPlex + comp)*dimPlex+comp  )*dimPlex+1: &
-                                           (((qPt*nBasis + basis)*dimPlex + comp)*dimPlex+comp+1)*dimPlex))
+          i = ((qPt*nBasis + basis)*dimPlex + comp)*dimPlex+comp
+          BMat(comp*dimPlex+1_pPETSCINT:(comp+1_pPETSCINT)*dimPlex,basis*dimPlex+comp+1_pPETSCINT) = &
+            matmul(reshape(pInvcellJ,[dimPlex,dimPlex]),basisFieldDer(i*dimPlex+1_pPETSCINT:(i+1_pPETSCINT)*dimPlex))
         enddo
       enddo
       MatA = matmul(reshape(reshape(homogenization_dPdF(1:dimPlex,1:dimPlex,1:dimPlex,1:dimPlex,m), &
                                     shape=[dimPlex,dimPlex,dimPlex,dimPlex], order=[2,1,4,3]), &
-                            shape=[dimPlex*dimPlex,dimPlex*dimPlex]),BMat)*qWeights(qPt+1)
+                            shape=[dimPlex*dimPlex,dimPlex*dimPlex]),BMat)*qWeights(qPt+1_pPETSCINT)
       if (num%BBarStabilisation) then
         F(1:dimPlex,1:dimPlex) = reshape(matmul(BMat,x_scal),shape=[dimPlex,dimPlex])
         FInv = math_inv33(F)
-        K_eA = K_eA + matmul(transpose(BMat),MatA)*math_det33(FInv)**(1.0/real(dimPlex))
+        K_eA = K_eA + matmul(transpose(BMat),MatA)*math_det33(FInv)**(1.0_pReal/real(dimPlex,pReal))
         K_eB = K_eB - &
                matmul(transpose(matmul(reshape(homogenization_F(1:dimPlex,1:dimPlex,m),shape=[dimPlex**2,1_pPETSCINT]), &
                                        matmul(reshape(FInv(1:dimPlex,1:dimPlex), &
@@ -568,10 +579,10 @@ subroutine FEM_mechanical_formJacobian(dm_local,xx_local,Jac_pre,Jac,dummy,err_P
     enddo
     if (num%BBarStabilisation) then
       FInv = math_inv33(FAvg)
-      K_e = K_eA*math_det33(FAvg/real(nQuadrature))**(1.0/real(dimPlex)) + &
+      K_e = K_eA*math_det33(FAvg/real(nQuadrature,pReal))**(1.0_pReal/real(dimPlex,pReal)) + &
             (matmul(matmul(transpose(BMatAvg), &
                            reshape(FInv(1:dimPlex,1:dimPlex),shape=[dimPlex**2,1_pPETSCINT],order=[2,1])),MatB) + &
-             K_eB)/real(dimPlex)
+             K_eB)/real(dimPlex,pReal)
     else
       K_e = K_eA
     endif
@@ -641,7 +652,7 @@ subroutine FEM_mechanical_forward(guess,timeinc,timeinc_old,fieldBC)
     CHKERRQ(err_PETSc)
     call DMGlobalToLocalEnd(dm_local,solution,INSERT_VALUES,x_local,err_PETSc)
     CHKERRQ(err_PETSc)
-    call VecAXPY(solution_local,1.0,x_local,err_PETSc); CHKERRQ(err_PETSc)
+    call VecAXPY(solution_local,1.0_pReal,x_local,err_PETSc); CHKERRQ(err_PETSc)
     do field = 1, dimPlex; do face = 1, mesh_Nboundaries
       if (fieldBC%componentBC(field)%Mask(face)) then
         call DMGetStratumSize(dm_local,'Face Sets',mesh_boundaries(face),bcSize,err_PETSc)
@@ -659,7 +670,7 @@ subroutine FEM_mechanical_forward(guess,timeinc,timeinc_old,fieldBC)
 !--------------------------------------------------------------------------------------------------
 ! update rate and forward last inc
     call VecCopy(solution,solution_rate,err_PETSc); CHKERRQ(err_PETSc)
-    call VecScale(solution_rate,1.0/timeinc_old,err_PETSc); CHKERRQ(err_PETSc)
+    call VecScale(solution_rate,timeinc_old**(-1),err_PETSc); CHKERRQ(err_PETSc)
   endif
   call VecCopy(solution_rate,solution,err_PETSc); CHKERRQ(err_PETSc)
   call VecScale(solution,timeinc,err_PETSc); CHKERRQ(err_PETSc)
@@ -685,9 +696,8 @@ subroutine FEM_mechanical_converged(snes_local,PETScIter,xnorm,snorm,fnorm,reaso
   call SNESConvergedDefault(snes_local,PETScIter,xnorm,snorm,fnorm/divTol,reason,dummy,err_PETSc)
   CHKERRQ(err_PETSc)
   if (terminallyIll) reason = SNES_DIVERGED_FUNCTION_DOMAIN
-  print'(/,1x,a,a,i0,a,i0,f0.3)', trim(incInfo), &
-                  ' @ Iteration ',PETScIter,' mechanical residual norm = ', &
-                                                  int(fnorm/divTol),fnorm/divTol-int(fnorm/divTol)  ! ToDo: int casting?
+  print'(/,1x,a,a,i0,a,f0.3)', trim(incInfo), &
+                  ' @ Iteration ',PETScIter,' mechanical residual norm = ',fnorm/divTol
   print'(/,1x,a,/,2(3(2x,f12.4,1x)/),3(2x,f12.4,1x))', &
     'Piola--Kirchhoff stress / MPa =',transpose(P_av)*1.e-6_pReal
   flush(IO_STDOUT)
@@ -747,7 +757,7 @@ subroutine FEM_mechanical_updateCoords()
   call PetscDSGetTabulation(mechQuad,0_pPETSCINT,basisField,basisFieldDer,err_PETSc)
   CHKERRQ(err_PETSc)
   allocate(ipCoords(3,nQuadrature,mesh_NcpElems),source=0.0_pReal)
-  do c=cellStart,cellEnd-1
+  do c=cellStart,cellEnd-1_pPETSCINT
     qOffset=0
     call DMPlexVecGetClosure(dm_local,section,x_local,c,x_scal,err_PETSc)                           !< get nodal coordinates of each element
     CHKERRQ(err_PETSc)

src/parallelization.f90

@@ -18,7 +18,11 @@ module parallelization
 
   use prec
 
+#if (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR>14) && !defined(PETSC_HAVE_MPI_F90MODULE_VISIBILITY)
+  implicit none(type,external)
+#else
   implicit none
+#endif
   private
 
 #ifndef PETSC
@@ -90,7 +94,7 @@ subroutine parallelization_init
 #ifdef LOGFILE
   write(rank_str,'(i4.4)') worldrank
   open(OUTPUT_UNIT,file='out.'//rank_str,status='replace',encoding='UTF-8')
-  open(ERROR_UNIT,file='error.'//rank_str,status='replace',encoding='UTF-8')
+  open(ERROR_UNIT,file='err.'//rank_str,status='replace',encoding='UTF-8')
 #else
   if (worldrank /= 0) then
     close(OUTPUT_UNIT)                                                                              ! disable output

src/phase.f90

@@ -19,7 +19,7 @@ module phase
   use HDF5
   use HDF5_utilities
 
-  implicit none
+  implicit none(type,external)
   private
 
   type :: tState
@@ -539,7 +539,8 @@ subroutine crystallite_init()
   class(tNode), pointer :: &
     num_crystallite, &
     phases
-
+  character(len=pStringLen) :: &
+    extmsg = ''
 
   num_crystallite => config_numerics%get('crystallite',defaultVal=emptyDict)
 
@@ -555,22 +556,19 @@ subroutine crystallite_init()
   num%nState                 = num_crystallite%get_asInt   ('nState',           defaultVal=20)
   num%nStress                = num_crystallite%get_asInt   ('nStress',          defaultVal=40)
 
-  if (num%subStepMinCryst   <= 0.0_pReal)      call IO_error(301,ext_msg='subStepMinCryst')
-  if (num%subStepSizeCryst  <= 0.0_pReal)      call IO_error(301,ext_msg='subStepSizeCryst')
-  if (num%stepIncreaseCryst <= 0.0_pReal)      call IO_error(301,ext_msg='stepIncreaseCryst')
-
-  if (num%subStepSizeLp <= 0.0_pReal)          call IO_error(301,ext_msg='subStepSizeLp')
-  if (num%subStepSizeLi <= 0.0_pReal)          call IO_error(301,ext_msg='subStepSizeLi')
-
-  if (num%rtol_crystalliteState  <= 0.0_pReal) call IO_error(301,ext_msg='rtol_crystalliteState')
-  if (num%rtol_crystalliteStress <= 0.0_pReal) call IO_error(301,ext_msg='rtol_crystalliteStress')
-  if (num%atol_crystalliteStress <= 0.0_pReal) call IO_error(301,ext_msg='atol_crystalliteStress')
-
-  if (num%iJacoLpresiduum < 1)                 call IO_error(301,ext_msg='iJacoLpresiduum')
-
-  if (num%nState < 1)                          call IO_error(301,ext_msg='nState')
-  if (num%nStress< 1)                          call IO_error(301,ext_msg='nStress')
+  if (num%subStepMinCryst   <= 0.0_pReal)      extmsg = trim(extmsg)//' subStepMinCryst'
+  if (num%subStepSizeCryst  <= 0.0_pReal)      extmsg = trim(extmsg)//' subStepSizeCryst'
+  if (num%stepIncreaseCryst <= 0.0_pReal)      extmsg = trim(extmsg)//' stepIncreaseCryst'
+  if (num%subStepSizeLp <= 0.0_pReal)          extmsg = trim(extmsg)//' subStepSizeLp'
+  if (num%subStepSizeLi <= 0.0_pReal)          extmsg = trim(extmsg)//' subStepSizeLi'
+  if (num%rtol_crystalliteState  <= 0.0_pReal) extmsg = trim(extmsg)//' rtol_crystalliteState'
+  if (num%rtol_crystalliteStress <= 0.0_pReal) extmsg = trim(extmsg)//' rtol_crystalliteStress'
+  if (num%atol_crystalliteStress <= 0.0_pReal) extmsg = trim(extmsg)//' atol_crystalliteStress'
+  if (num%iJacoLpresiduum < 1)                 extmsg = trim(extmsg)//' iJacoLpresiduum'
+  if (num%nState  < 1)                         extmsg = trim(extmsg)//' nState'
+  if (num%nStress < 1)                         extmsg = trim(extmsg)//' nStress'
 
+  if (extmsg /= '') call IO_error(301,ext_msg=trim(extmsg))
 
   phases => config_material%get('phase')
 
@@ -631,6 +629,7 @@ function crystallite_push33ToRef(co,ce, tensor33)
   real(pReal), dimension(3,3) :: T
   integer :: ph, en
 
+
   ph = material_phaseID(co,ce)
   en = material_phaseEntry(co,ce)
   T = matmul(phase_O_0(ph)%data(en)%asMatrix(),transpose(math_inv33(phase_F(co,ce))))               ! ToDo: initial orientation correct?

src/phase_damage.f90

@@ -4,8 +4,9 @@
 submodule(phase) damage
 
   type :: tDamageParameters
-    real(pReal) ::                 mu = 0.0_pReal                                                   !< viscosity
-    real(pReal), dimension(3,3) :: D  = 0.0_pReal                                                   !< conductivity/diffusivity
+    real(pReal) :: &
+      mu = 0.0_pReal, &                                                                             !< viscosity
+      l_c = 0.0_pReal                                                                               !< characteristic length
   end type tDamageParameters
 
   enum, bind(c); enumerator :: &
@@ -105,9 +106,7 @@ module subroutine damage_init
       damage_active = .true.
       source => sources%get(1)
       param(ph)%mu = source%get_asFloat('mu')
-      param(ph)%D(1,1) = source%get_asFloat('D_11')
-      if (any(phase_lattice(ph) == ['hP','tI'])) param(ph)%D(3,3) = source%get_asFloat('D_33')
-      param(ph)%D = lattice_symmetrize_33(param(ph)%D,phase_lattice(ph))
+      param(ph)%l_c = source%get_asFloat('l_c')
     end if
 
   end do
@@ -385,9 +384,9 @@ module function phase_K_phi(co,ce) result(K)
 
   integer, intent(in) :: co, ce
   real(pReal), dimension(3,3) :: K
-  real(pReal), parameter :: l = 1.0_pReal
 
-  K = crystallite_push33ToRef(co,ce,param(material_phaseID(co,ce))%D) * l**2
+
+  K = crystallite_push33ToRef(co,ce,param(material_phaseID(co,ce))%l_c**2*math_I3)
 
 end function phase_K_phi
 
@@ -403,6 +402,7 @@ function phase_damage_deltaState(Fe, ph, en) result(broken)
     en
   real(pReal),   intent(in), dimension(3,3) :: &
     Fe                                                                                              !< elastic deformation gradient
+
   integer :: &
     myOffset, &
     mySize

src/phase_damage_anisobrittle.f90

@@ -157,7 +157,7 @@ module subroutine anisobrittle_results(phase,group)
     outputsLoop: do o = 1,size(prm%output)
       select case(trim(prm%output(o)))
         case ('f_phi')
-          call results_writeDataset(stt,group,trim(prm%output(o)),'driving force','J/m³')
+          call results_writeDataset(stt,group,trim(prm%output(o)),'driving force','-')
       end select
     end do outputsLoop
   end associate

src/phase_damage_isobrittle.f90

@@ -63,7 +63,7 @@ module function isobrittle_init() result(mySources)
       associate(prm => param(ph), dlt => deltaState(ph), stt => state(ph))
         src => sources%get(1)
 
-        prm%W_crit = src%get_asFloat('W_crit')
+        prm%W_crit = src%get_asFloat('G_crit')/src%get_asFloat('l_c')
 
 #if defined (__GFORTRAN__)
         prm%output = output_as1dString(src)
@@ -75,7 +75,7 @@ module function isobrittle_init() result(mySources)
         if (prm%W_crit <= 0.0_pReal) extmsg = trim(extmsg)//' W_crit'
 
         Nmembers = count(material_phaseID==ph)
-        call phase_allocateState(damageState(ph),Nmembers,1,1,1)
+        call phase_allocateState(damageState(ph),Nmembers,1,0,1)
         damageState(ph)%atol = src%get_asFloat('atol_phi',defaultVal=1.0e-9_pReal)
         if (any(damageState(ph)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' atol_phi'
 
@@ -139,7 +139,7 @@ module subroutine isobrittle_results(phase,group)
     outputsLoop: do o = 1,size(prm%output)
       select case(trim(prm%output(o)))
         case ('f_phi')
-          call results_writeDataset(stt,group,trim(prm%output(o)),'driving force','J/m³') ! Wrong, this is dimensionless
+          call results_writeDataset(stt,group,trim(prm%output(o)),'driving force','-')
       end select
     end do outputsLoop
 

src/phase_mechanical.f90

@@ -593,7 +593,7 @@ function integrateStateFPI(F_0,F,subFp0,subFi0,subState0,Delta_t,ph,en) result(b
 
   iteration: do NiterationState = 1, num%nState
 
-    dotState_last(1:sizeDotState,2) = merge(dotState_last(1:sizeDotState,1),0.0, nIterationState > 1)
+    dotState_last(1:sizeDotState,2) = merge(dotState_last(1:sizeDotState,1),0.0_pReal, nIterationState > 1)
     dotState_last(1:sizeDotState,1) = dotState
 
     broken = integrateStress(F,subFp0,subFi0,Delta_t,ph,en)
@@ -756,7 +756,7 @@ function integrateStateRK4(F_0,F,subFp0,subFi0,subState0,Delta_t,ph,en) result(b
   real(pReal), dimension(3), parameter :: &
     C = [0.5_pReal, 0.5_pReal, 1.0_pReal]
   real(pReal), dimension(4), parameter :: &
-    B = [1.0_pReal/6.0_pReal, 1.0_pReal/3.0_pReal, 1.0_pReal/3.0_pReal, 1.0_pReal/6.0_pReal]
+    B = [6.0_pReal, 3.0_pReal, 3.0_pReal, 6.0_pReal]**(-1)
 
 
   broken = integrateStateRK(F_0,F,subFp0,subFi0,subState0,Delta_t,ph,en,A,B,C)
@@ -1167,8 +1167,8 @@ module function phase_mechanical_dPdF(Delta_t,co,ce) result(dPdF)
     end do; end do
     call math_invert(temp_99,error,math_3333to99(lhs_3333))
     if (error) then
-      call IO_warning(warning_ID=600, &
-                      ext_msg='inversion error in analytic tangent calculation')
+      call IO_warning(600,'inversion error in analytic tangent calculation', &
+                      label1='phase',ID1=ph,label2='entry',ID2=en)
       dFidS = 0.0_pReal
     else
       dFidS = math_mul3333xx3333(math_99to3333(temp_99),rhs_3333)
@@ -1201,8 +1201,8 @@ module function phase_mechanical_dPdF(Delta_t,co,ce) result(dPdF)
 
   call math_invert(temp_99,error,math_eye(9)+math_3333to99(lhs_3333))
   if (error) then
-    call IO_warning(warning_ID=600, &
-                    ext_msg='inversion error in analytic tangent calculation')
+    call IO_warning(600,'inversion error in analytic tangent calculation', &
+                    label1='phase',ID1=ph,label2='entry',ID2=en)
     dSdF = rhs_3333
   else
     dSdF = math_mul3333xx3333(math_99to3333(temp_99),rhs_3333)
@@ -1263,8 +1263,6 @@ module subroutine mechanical_restartRead(groupHandle,ph)
   integer(HID_T), intent(in) :: groupHandle
   integer, intent(in) :: ph
 
-  integer :: en
-
 
   call HDF5_read(plasticState(ph)%state0,groupHandle,'omega_plastic')
   call HDF5_read(phase_mechanical_S0(ph)%data,groupHandle,'S')

src/phase_mechanical_plastic_dislotungsten.f90

@@ -252,7 +252,7 @@ module function plastic_dislotungsten_init() result(myPlasticity)
 
 !--------------------------------------------------------------------------------------------------
 !  exit if any parameter is out of range
-    if (extmsg /= '') call IO_error(211,ext_msg=trim(extmsg)//'(dislotungsten)')
+    if (extmsg /= '') call IO_error(211,ext_msg=trim(extmsg))
 
   end do
 

src/phase_mechanical_plastic_dislotwin.f90

@@ -9,6 +9,7 @@
 !--------------------------------------------------------------------------------------------------
 submodule(phase:plastic) dislotwin
 
+  real(pReal), parameter :: gamma_char_tr = sqrt(0.125_pReal)                                       !< Characteristic shear for transformation
   type :: tParameters
     real(pReal) :: &
       Q_cl       = 1.0_pReal, &                                                                     !< activation energy for dislocation climb
@@ -22,11 +23,10 @@ submodule(phase:plastic) dislotwin
       L_tr       = 1.0_pReal, &                                                                     !< length of trans nuclei
       x_c        = 1.0_pReal, &                                                                     !< critical distance for formation of twin/trans nucleus
       V_cs       = 1.0_pReal, &                                                                     !< cross slip volume
-      xi_sb               = 1.0_pReal, &                                                            !< value for shearband resistance
-      v_sb                = 1.0_pReal, &                                                            !< value for shearband velocity_0
+      tau_sb     = 1.0_pReal, &                                                                     !< value for shearband resistance
+      gamma_0_sb = 1.0_pReal, &                                                                     !< value for shearband velocity_0
       E_sb       = 1.0_pReal, &                                                                     !< activation energy for shear bands
       h          = 1.0_pReal, &                                                                     !< stack height of hex nucleus
-      gamma_char_tr       = sqrt(0.125_pReal), &                                                    !< Characteristic shear for transformation
       a_cF       = 1.0_pReal, &
       cOverA_hP  = 1.0_pReal, &
       V_mol      = 1.0_pReal, &
@@ -331,15 +331,15 @@ module function plastic_dislotwin_init() result(myPlasticity)
 
 !--------------------------------------------------------------------------------------------------
 ! shearband related parameters
-    prm%v_sb = pl%get_asFloat('v_sb',defaultVal=0.0_pReal)
-    if (prm%v_sb > 0.0_pReal) then
-      prm%xi_sb        = pl%get_asFloat('xi_sb')
+    prm%gamma_0_sb = pl%get_asFloat('gamma_0_sb',defaultVal=0.0_pReal)
+    if (prm%gamma_0_sb > 0.0_pReal) then
+      prm%tau_sb       = pl%get_asFloat('tau_sb')
       prm%E_sb         = pl%get_asFloat('Q_sb')
       prm%p_sb         = pl%get_asFloat('p_sb')
       prm%q_sb         = pl%get_asFloat('q_sb')
 
       ! sanity checks
-      if (prm%xi_sb         <  0.0_pReal) extmsg = trim(extmsg)//' xi_sb'
+      if (prm%tau_sb <  0.0_pReal) extmsg = trim(extmsg)//' tau_sb'
       if (prm%E_sb   <  0.0_pReal) extmsg = trim(extmsg)//' Q_sb'
       if (prm%p_sb   <= 0.0_pReal) extmsg = trim(extmsg)//' p_sb'
       if (prm%q_sb   <= 0.0_pReal) extmsg = trim(extmsg)//' q_sb'
@@ -430,7 +430,7 @@ module function plastic_dislotwin_init() result(myPlasticity)
 
 !--------------------------------------------------------------------------------------------------
 !  exit if any parameter is out of range
-    if (extmsg /= '') call IO_error(211,ext_msg=trim(extmsg)//'(dislotwin)')
+    if (extmsg /= '') call IO_error(211,ext_msg=trim(extmsg))
 
   end do
 
@@ -569,7 +569,7 @@ module subroutine dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,en)
     Lp      = Lp      * f_matrix
     dLp_dMp = dLp_dMp * f_matrix
 
-    shearBandingContribution: if (dNeq0(prm%v_sb)) then
+    shearBandingContribution: if (dNeq0(prm%gamma_0_sb)) then
 
       E_kB_T = prm%E_sb/(K_B*T)
       call math_eigh33(eigValues,eigVectors,Mp)                                                     ! is Mp symmetric by design?
@@ -580,10 +580,10 @@ module subroutine dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,en)
         tau = math_tensordot(Mp,P_sb)
 
         significantShearBandStress: if (abs(tau) > tol_math_check) then
-          StressRatio_p = (abs(tau)/prm%xi_sb)**prm%p_sb
-          dot_gamma_sb = sign(prm%v_sb*exp(-E_kB_T*(1-StressRatio_p)**prm%q_sb), tau)
-          ddot_gamma_dtau = abs(dot_gamma_sb)*E_kB_T*prm%p_sb*prm%q_sb/prm%xi_sb &
-                          * (abs(tau)/prm%xi_sb)**(prm%p_sb-1.0_pReal) &
+          StressRatio_p = (abs(tau)/prm%tau_sb)**prm%p_sb
+          dot_gamma_sb = sign(prm%gamma_0_sb*exp(-E_kB_T*(1-StressRatio_p)**prm%q_sb), tau)
+          ddot_gamma_dtau = abs(dot_gamma_sb)*E_kB_T*prm%p_sb*prm%q_sb/prm%tau_sb &
+                          * (abs(tau)/prm%tau_sb)**(prm%p_sb-1.0_pReal) &
                           * (1.0_pReal-StressRatio_p)**(prm%q_sb-1.0_pReal)
 
           Lp = Lp + dot_gamma_sb * P_sb
@@ -697,7 +697,7 @@ module function dislotwin_dotState(Mp,ph,en) result(dotState)
     dot_f_tw = f_matrix*dot_gamma_tw/prm%gamma_char_tw
 
     if (prm%sum_N_tr > 0) call kinetics_tr(Mp,T,abs_dot_gamma_sl,ph,en,dot_gamma_tr)
-    dot_f_tr = f_matrix*dot_gamma_tr/prm%gamma_char_tr
+    dot_f_tr = f_matrix*dot_gamma_tr/gamma_char_tr
 
   end associate
 
@@ -912,7 +912,7 @@ pure subroutine kinetics_tw(Mp,T,abs_dot_gamma_sl,ph,en,&
   real(pReal), dimension(param(ph)%sum_N_tw), optional, intent(out) :: &
     ddot_gamma_dtau_tw
 
-  real :: &
+  real(pReal) :: &
     tau, tau_r, tau_hat, &
     dot_N_0, &
     x0, V, &
@@ -988,7 +988,7 @@ pure subroutine kinetics_tr(Mp,T,abs_dot_gamma_sl,ph,en,&
   real(pReal), dimension(param(ph)%sum_N_tr), optional, intent(out) :: &
     ddot_gamma_dtau_tr
 
-  real :: &
+  real(pReal) :: &
     tau, tau_r, tau_hat, &
     dot_N_0, &
     x0, V, &
@@ -1026,9 +1026,9 @@ pure subroutine kinetics_tr(Mp,T,abs_dot_gamma_sl,ph,en,&
         dP_ncs_dtau = prm%V_cs / (K_B * T) * (P_ncs - 1.0_pReal)
 
         V = PI/4.0_pReal*dst%Lambda_tr(i,en)**2*prm%t_tr(i)
-        dot_gamma_tr(i) = V*dot_N_0*P_ncs*P*prm%gamma_char_tr
+        dot_gamma_tr(i) = V*dot_N_0*P_ncs*P*gamma_char_tr
         if (present(ddot_gamma_dtau_tr)) &
-          ddot_gamma_dtau_tr(i) = V*dot_N_0*(P*dP_ncs_dtau + P_ncs*dP_dtau)*prm%gamma_char_tr
+          ddot_gamma_dtau_tr(i) = V*dot_N_0*(P*dP_ncs_dtau + P_ncs*dP_dtau)*gamma_char_tr
       else
         dot_gamma_tr(i) = 0.0_pReal
         if (present(ddot_gamma_dtau_tr)) ddot_gamma_dtau_tr(i) = 0.0_pReal

src/phase_mechanical_plastic_isotropic.f90

@@ -69,8 +69,8 @@ module function plastic_isotropic_init() result(myPlasticity)
   print'(/,1x,a)', '<<<+-  phase:mechanical:plastic:isotropic init  -+>>>'
   print'(/,a,i0)', ' # phases: ',count(myPlasticity); flush(IO_STDOUT)
 
-  print'(/,a)', 'T. Maiti and P. Eisenlohr, Scripta Materialia 145:37–40, 2018'
-  print'(/,a)', 'https://doi.org/10.1016/j.scriptamat.2017.09.047'
+  print'(/,1x,a)', 'T. Maiti and P. Eisenlohr, Scripta Materialia 145:37–40, 2018'
+  print'(  1x,a)', 'https://doi.org/10.1016/j.scriptamat.2017.09.047'
 
   phases => config_material%get('phase')
   allocate(param(phases%length))
@@ -135,7 +135,7 @@ module function plastic_isotropic_init() result(myPlasticity)
 
 !--------------------------------------------------------------------------------------------------
 !  exit if any parameter is out of range
-    if (extmsg /= '') call IO_error(211,ext_msg=trim(extmsg)//'(isotropic)')
+    if (extmsg /= '') call IO_error(211,ext_msg=trim(extmsg))
 
   end do
 

src/phase_mechanical_plastic_kinehardening.f90

@@ -224,7 +224,7 @@ module function plastic_kinehardening_init() result(myPlasticity)
 
 !--------------------------------------------------------------------------------------------------
 !  exit if any parameter is out of range
-    if (extmsg /= '') call IO_error(211,ext_msg=trim(extmsg)//'(kinehardening)')
+    if (extmsg /= '') call IO_error(211,ext_msg=trim(extmsg))
 
   end do
 

src/phase_mechanical_plastic_nonlocal.f90

@@ -504,7 +504,7 @@ module function plastic_nonlocal_init() result(myPlasticity)
 
 !--------------------------------------------------------------------------------------------------
 !  exit if any parameter is out of range
-    if (extmsg /= '') call IO_error(211,ext_msg=trim(extmsg)//'(nonlocal)')
+    if (extmsg /= '') call IO_error(211,ext_msg=trim(extmsg))
 
   end do
 

src/phase_mechanical_plastic_phenopowerlaw.f90

@@ -269,7 +269,7 @@ module function plastic_phenopowerlaw_init() result(myPlasticity)
 
 !--------------------------------------------------------------------------------------------------
 !  exit if any parameter is out of range
-    if (extmsg /= '') call IO_error(211,ext_msg=trim(extmsg)//'(phenopowerlaw)')
+    if (extmsg /= '') call IO_error(211,ext_msg=trim(extmsg))
 
   end do
 

src/phase_thermal.f90

@@ -72,7 +72,7 @@ submodule(phase) thermal
 contains
 
 !----------------------------------------------------------------------------------------------
-!< @brief initializes thermal sources and kinematics mechanism
+!< @brief Initializes thermal sources and kinematics mechanism.
 !----------------------------------------------------------------------------------------------
 module subroutine thermal_init(phases)
 
@@ -146,7 +146,7 @@ end subroutine thermal_init
 
 
 !----------------------------------------------------------------------------------------------
-!< @brief calculates thermal dissipation rate
+!< @brief Calculate thermal source.
 !----------------------------------------------------------------------------------------------
 module function phase_f_T(ph,en) result(f)
 
@@ -176,7 +176,7 @@ end function phase_f_T
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief contains the constitutive equation for calculating the rate of change of microstructure
+!> @brief tbd.
 !--------------------------------------------------------------------------------------------------
 function phase_thermal_collectDotState(ph,en) result(broken)
 
@@ -216,7 +216,7 @@ end function phase_mu_T
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief Thermal conductivity/diffusivity in reference configuration.
+!> @brief Thermal conductivity in reference configuration.
 !--------------------------------------------------------------------------------------------------
 module function phase_K_T(co,ce) result(K)
 
@@ -255,6 +255,7 @@ function integrateThermalState(Delta_t, ph,en) result(broken)
     so, &
     sizeDotState
 
+
   broken = phase_thermal_collectDotState(ph,en)
   if (broken) return
 

src/polynomials.f90

@@ -8,7 +8,7 @@ module polynomials
   use YAML_parse
   use YAML_types
 
-  implicit none
+  implicit none(type,external)
   private
 
   type, public :: tPolynomial

src/prec.f90

@@ -15,7 +15,7 @@ module prec
   use PETScSys
 #endif
 
-  implicit none
+  implicit none(type,external)
   public
 
   ! https://stevelionel.com/drfortran/2017/03/27/doctor-fortran-in-it-takes-all-kinds

src/quit.f90

@@ -12,8 +12,14 @@ subroutine quit(stop_id)
 #endif
   use HDF5
 
+#if (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR>14) && !defined(PETSC_HAVE_MPI_F90MODULE_VISIBILITY)
+  implicit none(type,external)
+#else
   implicit none
+#endif
+
   integer, intent(in) :: stop_id
+
   integer, dimension(8) :: dateAndTime
   integer :: err_HDF5
   integer(MPI_INTEGER_KIND) :: err_MPI

src/results.f90

@@ -21,7 +21,11 @@ module results
   use DAMASK_interface
 #endif
 
+#if (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR>14) && !defined(PETSC_HAVE_MPI_F90MODULE_VISIBILITY)
+  implicit none(type,external)
+#else
   implicit none
+#endif
   private
 
   integer(HID_T) :: resultsFile
@@ -95,11 +99,11 @@ subroutine results_init(restart)
     call results_openJobFile
     call get_command(commandLine)
     call results_addAttribute('call (restart at '//date//')',trim(commandLine))
-    call h5gmove_f(resultsFile,'setup','tmp',hdferr)
+    call H5Gmove_f(resultsFile,'setup','tmp',hdferr)
     call results_addAttribute('description','input data used to run the simulation up to restart at '//date,'tmp')
     call results_closeGroup(results_addGroup('setup'))
     call results_addAttribute('description','input data used to run the simulation','setup')
-    call h5gmove_f(resultsFile,'tmp','setup/previous',hdferr)
+    call H5Gmove_f(resultsFile,'tmp','setup/previous',hdferr)
   end if
 
   call results_closeJobFile
@@ -333,8 +337,8 @@ subroutine results_removeLink(link)
   integer                      :: hdferr
 
 
-  call h5ldelete_f(resultsFile,link, hdferr)
-  if (hdferr < 0) call IO_error(1,ext_msg = 'results_removeLink: h5ldelete_soft_f ('//trim(link)//')')
+  call H5Ldelete_f(resultsFile,link, hdferr)
+  if (hdferr < 0) call IO_error(1,ext_msg = 'results_removeLink: H5Ldelete_soft_f ('//trim(link)//')')
 
 end subroutine results_removeLink
 
@@ -522,7 +526,7 @@ subroutine results_mapping_phase(ID,entry,label)
   writeSize = 0
   writeSize(worldrank) = size(entry(1,:))                                                           ! total number of entries of this process
 
-  call h5pcreate_f(H5P_DATASET_XFER_F, plist_id, hdferr)
+  call H5Pcreate_f(H5P_DATASET_XFER_F, plist_id, hdferr)
   if(hdferr < 0) error stop 'HDF5 error'
 
 #ifndef PETSC
@@ -530,7 +534,7 @@ subroutine results_mapping_phase(ID,entry,label)
 #else
 !--------------------------------------------------------------------------------------------------
 ! MPI settings and communication
-  call h5pset_dxpl_mpio_f(plist_id, H5FD_MPIO_COLLECTIVE_F, hdferr)
+  call H5Pset_dxpl_mpio_f(plist_id, H5FD_MPIO_COLLECTIVE_F, hdferr)
   if(hdferr < 0) error stop 'HDF5 error'
 
   call MPI_Allreduce(MPI_IN_PLACE,writeSize,worldsize,MPI_INTEGER,MPI_SUM,MPI_COMM_WORLD,err_MPI)   ! get output at each process
@@ -558,82 +562,82 @@ subroutine results_mapping_phase(ID,entry,label)
 
 !---------------------------------------------------------------------------------------------------
 ! compound type: label(ID) + entry
-  call h5tcopy_f(H5T_NATIVE_CHARACTER, dt_id, hdferr)
+  call H5Tcopy_f(H5T_NATIVE_CHARACTER, dt_id, hdferr)
   if(hdferr < 0) error stop 'HDF5 error'
-  call h5tset_size_f(dt_id, int(len(label(1)),SIZE_T), hdferr)
+  call H5Tset_size_f(dt_id, int(len(label(1)),SIZE_T), hdferr)
   if(hdferr < 0) error stop 'HDF5 error'
-  call h5tget_size_f(dt_id, type_size_string, hdferr)
+  call H5Tget_size_f(dt_id, type_size_string, hdferr)
   if(hdferr < 0) error stop 'HDF5 error'
 
   pI64_t = h5kind_to_type(kind(entryGlobal),H5_INTEGER_KIND)
-  call h5tget_size_f(pI64_t, type_size_int, hdferr)
+  call H5Tget_size_f(pI64_t, type_size_int, hdferr)
   if(hdferr < 0) error stop 'HDF5 error'
 
-  call h5tcreate_f(H5T_COMPOUND_F, type_size_string + type_size_int, dtype_id, hdferr)
+  call H5Tcreate_f(H5T_COMPOUND_F, type_size_string + type_size_int, dtype_id, hdferr)
   if(hdferr < 0) error stop 'HDF5 error'
-  call h5tinsert_f(dtype_id, 'label', 0_SIZE_T, dt_id,hdferr)
+  call H5Tinsert_f(dtype_id, 'label', 0_SIZE_T, dt_id,hdferr)
   if(hdferr < 0) error stop 'HDF5 error'
-  call h5tinsert_f(dtype_id, 'entry', type_size_string, pI64_t, hdferr)
+  call H5Tinsert_f(dtype_id, 'entry', type_size_string, pI64_t, hdferr)
   if(hdferr < 0) error stop 'HDF5 error'
 
 !--------------------------------------------------------------------------------------------------
 ! create memory types for each component of the compound type
-  call h5tcreate_f(H5T_COMPOUND_F, type_size_string, label_id, hdferr)
+  call H5Tcreate_f(H5T_COMPOUND_F, type_size_string, label_id, hdferr)
   if(hdferr < 0) error stop 'HDF5 error'
-  call h5tinsert_f(label_id, 'label', 0_SIZE_T, dt_id, hdferr)
+  call H5Tinsert_f(label_id, 'label', 0_SIZE_T, dt_id, hdferr)
   if(hdferr < 0) error stop 'HDF5 error'
 
-  call h5tcreate_f(H5T_COMPOUND_F, type_size_int, entry_id, hdferr)
+  call H5Tcreate_f(H5T_COMPOUND_F, type_size_int, entry_id, hdferr)
   if(hdferr < 0) error stop 'HDF5 error'
-  call h5tinsert_f(entry_id, 'entry', 0_SIZE_T, pI64_t, hdferr)
+  call H5Tinsert_f(entry_id, 'entry', 0_SIZE_T, pI64_t, hdferr)
   if(hdferr < 0) error stop 'HDF5 error'
 
-  call h5tclose_f(dt_id, hdferr)
+  call H5Tclose_f(dt_id, hdferr)
   if(hdferr < 0) error stop 'HDF5 error'
 
 !--------------------------------------------------------------------------------------------------
 ! create dataspace in memory (local shape = hyperslab) and in file (global shape)
-  call h5screate_simple_f(2,myShape,memspace_id,hdferr,myShape)
+  call H5Screate_simple_f(2,myShape,memspace_id,hdferr,myShape)
   if(hdferr < 0) error stop 'HDF5 error'
 
-  call h5screate_simple_f(2,totalShape,filespace_id,hdferr,totalShape)
+  call H5Screate_simple_f(2,totalShape,filespace_id,hdferr,totalShape)
   if(hdferr < 0) error stop 'HDF5 error'
 
-  call h5sselect_hyperslab_f(filespace_id, H5S_SELECT_SET_F, myOffset, myShape, hdferr)
+  call H5Sselect_hyperslab_f(filespace_id, H5S_SELECT_SET_F, myOffset, myShape, hdferr)
   if(hdferr < 0) error stop 'HDF5 error'
 
 !--------------------------------------------------------------------------------------------------
 ! write the components of the compound type individually
-  call h5pset_preserve_f(plist_id, .true., hdferr)
+  call H5Pset_preserve_f(plist_id, .true., hdferr)
   if(hdferr < 0) error stop 'HDF5 error'
 
   loc_id = results_openGroup('/cell_to')
-  call h5dcreate_f(loc_id, 'phase', dtype_id, filespace_id, dset_id, hdferr)
+  call H5Dcreate_f(loc_id, 'phase', dtype_id, filespace_id, dset_id, hdferr)
   if(hdferr < 0) error stop 'HDF5 error'
 
-  call h5dwrite_f(dset_id, label_id, reshape(label(pack(ID,.true.)),myShape), &
+  call H5Dwrite_f(dset_id, label_id, reshape(label(pack(ID,.true.)),myShape), &
                   myShape, hdferr, file_space_id = filespace_id, mem_space_id = memspace_id, xfer_prp = plist_id)
   if(hdferr < 0) error stop 'HDF5 error'
-  call h5dwrite_f(dset_id, entry_id, reshape(pack(entryGlobal,.true.),myShape), &
+  call H5Dwrite_f(dset_id, entry_id, reshape(pack(entryGlobal,.true.),myShape), &
                   myShape, hdferr, file_space_id = filespace_id, mem_space_id = memspace_id, xfer_prp = plist_id)
   if(hdferr < 0) error stop 'HDF5 error'
 
 !--------------------------------------------------------------------------------------------------
 ! close all
   call HDF5_closeGroup(loc_id)
-  call h5pclose_f(plist_id, hdferr)
+  call H5Pclose_f(plist_id, hdferr)
   if(hdferr < 0) error stop 'HDF5 error'
-  call h5sclose_f(filespace_id, hdferr)
+  call H5Sclose_f(filespace_id, hdferr)
   if(hdferr < 0) error stop 'HDF5 error'
-  call h5sclose_f(memspace_id, hdferr)
+  call H5Sclose_f(memspace_id, hdferr)
   if(hdferr < 0) error stop 'HDF5 error'
-  call h5dclose_f(dset_id, hdferr)
+  call H5Dclose_f(dset_id, hdferr)
   if(hdferr < 0) error stop 'HDF5 error'
-  call h5tclose_f(dtype_id, hdferr)
+  call H5Tclose_f(dtype_id, hdferr)
   if(hdferr < 0) error stop 'HDF5 error'
-  call h5tclose_f(label_id, hdferr)
+  call H5Tclose_f(label_id, hdferr)
   if(hdferr < 0) error stop 'HDF5 error'
-  call h5tclose_f(entry_id, hdferr)
+  call H5Tclose_f(entry_id, hdferr)
 
   call executionStamp('cell_to/phase','cell ID and constituent ID to phase results')
 
@@ -678,7 +682,7 @@ subroutine results_mapping_homogenization(ID,entry,label)
   writeSize = 0
   writeSize(worldrank) = size(entry)                                                                ! total number of entries of this process
 
-  call h5pcreate_f(H5P_DATASET_XFER_F, plist_id, hdferr)
+  call H5Pcreate_f(H5P_DATASET_XFER_F, plist_id, hdferr)
   if(hdferr < 0) error stop 'HDF5 error'
 
 #ifndef PETSC
@@ -686,7 +690,7 @@ subroutine results_mapping_homogenization(ID,entry,label)
 #else
 !--------------------------------------------------------------------------------------------------
 ! MPI settings and communication
-  call h5pset_dxpl_mpio_f(plist_id, H5FD_MPIO_COLLECTIVE_F, hdferr)
+  call H5Pset_dxpl_mpio_f(plist_id, H5FD_MPIO_COLLECTIVE_F, hdferr)
   if(hdferr < 0) error stop 'HDF5 error'
 
   call MPI_Allreduce(MPI_IN_PLACE,writeSize,worldsize,MPI_INTEGER,MPI_SUM,MPI_COMM_WORLD,err_MPI)   ! get output at each process
@@ -710,82 +714,82 @@ subroutine results_mapping_homogenization(ID,entry,label)
 
 !---------------------------------------------------------------------------------------------------
 ! compound type: label(ID) + entry
-  call h5tcopy_f(H5T_NATIVE_CHARACTER, dt_id, hdferr)
+  call H5Tcopy_f(H5T_NATIVE_CHARACTER, dt_id, hdferr)
   if(hdferr < 0) error stop 'HDF5 error'
-  call h5tset_size_f(dt_id, int(len(label(1)),SIZE_T), hdferr)
+  call H5Tset_size_f(dt_id, int(len(label(1)),SIZE_T), hdferr)
   if(hdferr < 0) error stop 'HDF5 error'
-  call h5tget_size_f(dt_id, type_size_string, hdferr)
+  call H5Tget_size_f(dt_id, type_size_string, hdferr)
   if(hdferr < 0) error stop 'HDF5 error'
 
   pI64_t = h5kind_to_type(kind(entryGlobal),H5_INTEGER_KIND)
-  call h5tget_size_f(pI64_t, type_size_int, hdferr)
+  call H5Tget_size_f(pI64_t, type_size_int, hdferr)
   if(hdferr < 0) error stop 'HDF5 error'
 
-  call h5tcreate_f(H5T_COMPOUND_F, type_size_string + type_size_int, dtype_id, hdferr)
+  call H5Tcreate_f(H5T_COMPOUND_F, type_size_string + type_size_int, dtype_id, hdferr)
   if(hdferr < 0) error stop 'HDF5 error'
-  call h5tinsert_f(dtype_id, 'label', 0_SIZE_T, dt_id,hdferr)
+  call H5Tinsert_f(dtype_id, 'label', 0_SIZE_T, dt_id,hdferr)
   if(hdferr < 0) error stop 'HDF5 error'
-  call h5tinsert_f(dtype_id, 'entry', type_size_string, pI64_t, hdferr)
+  call H5Tinsert_f(dtype_id, 'entry', type_size_string, pI64_t, hdferr)
   if(hdferr < 0) error stop 'HDF5 error'
 
 !--------------------------------------------------------------------------------------------------
 ! create memory types for each component of the compound type
-  call h5tcreate_f(H5T_COMPOUND_F, type_size_string, label_id, hdferr)
+  call H5Tcreate_f(H5T_COMPOUND_F, type_size_string, label_id, hdferr)
   if(hdferr < 0) error stop 'HDF5 error'
-  call h5tinsert_f(label_id, 'label', 0_SIZE_T, dt_id, hdferr)
+  call H5Tinsert_f(label_id, 'label', 0_SIZE_T, dt_id, hdferr)
   if(hdferr < 0) error stop 'HDF5 error'
 
-  call h5tcreate_f(H5T_COMPOUND_F, type_size_int, entry_id, hdferr)
+  call H5Tcreate_f(H5T_COMPOUND_F, type_size_int, entry_id, hdferr)
   if(hdferr < 0) error stop 'HDF5 error'
-  call h5tinsert_f(entry_id, 'entry', 0_SIZE_T, pI64_t, hdferr)
+  call H5Tinsert_f(entry_id, 'entry', 0_SIZE_T, pI64_t, hdferr)
   if(hdferr < 0) error stop 'HDF5 error'
 
-  call h5tclose_f(dt_id, hdferr)
+  call H5Tclose_f(dt_id, hdferr)
   if(hdferr < 0) error stop 'HDF5 error'
 
 !--------------------------------------------------------------------------------------------------
 ! create dataspace in memory (local shape = hyperslab) and in file (global shape)
-  call h5screate_simple_f(1,myShape,memspace_id,hdferr,myShape)
+  call H5Screate_simple_f(1,myShape,memspace_id,hdferr,myShape)
   if(hdferr < 0) error stop 'HDF5 error'
 
-  call h5screate_simple_f(1,totalShape,filespace_id,hdferr,totalShape)
+  call H5Screate_simple_f(1,totalShape,filespace_id,hdferr,totalShape)
   if(hdferr < 0) error stop 'HDF5 error'
 
-  call h5sselect_hyperslab_f(filespace_id, H5S_SELECT_SET_F, myOffset, myShape, hdferr)
+  call H5Sselect_hyperslab_f(filespace_id, H5S_SELECT_SET_F, myOffset, myShape, hdferr)
   if(hdferr < 0) error stop 'HDF5 error'
 
 !--------------------------------------------------------------------------------------------------
 ! write the components of the compound type individually
-  call h5pset_preserve_f(plist_id, .true., hdferr)
+  call H5Pset_preserve_f(plist_id, .true., hdferr)
   if(hdferr < 0) error stop 'HDF5 error'
 
   loc_id = results_openGroup('/cell_to')
-  call h5dcreate_f(loc_id, 'homogenization', dtype_id, filespace_id, dset_id, hdferr)
+  call H5Dcreate_f(loc_id, 'homogenization', dtype_id, filespace_id, dset_id, hdferr)
   if(hdferr < 0) error stop 'HDF5 error'
 
-  call h5dwrite_f(dset_id, label_id, reshape(label(pack(ID,.true.)),myShape), &
+  call H5Dwrite_f(dset_id, label_id, reshape(label(pack(ID,.true.)),myShape), &
                   myShape, hdferr, file_space_id = filespace_id, mem_space_id = memspace_id, xfer_prp = plist_id)
   if(hdferr < 0) error stop 'HDF5 error'
-  call h5dwrite_f(dset_id, entry_id, reshape(pack(entryGlobal,.true.),myShape), &
+  call H5Dwrite_f(dset_id, entry_id, reshape(pack(entryGlobal,.true.),myShape), &
                   myShape, hdferr, file_space_id = filespace_id, mem_space_id = memspace_id, xfer_prp = plist_id)
   if(hdferr < 0) error stop 'HDF5 error'
 
 !--------------------------------------------------------------------------------------------------
 ! close all
   call HDF5_closeGroup(loc_id)
-  call h5pclose_f(plist_id, hdferr)
+  call H5Pclose_f(plist_id, hdferr)
   if(hdferr < 0) error stop 'HDF5 error'
-  call h5sclose_f(filespace_id, hdferr)
+  call H5Sclose_f(filespace_id, hdferr)
   if(hdferr < 0) error stop 'HDF5 error'
-  call h5sclose_f(memspace_id, hdferr)
+  call H5Sclose_f(memspace_id, hdferr)
   if(hdferr < 0) error stop 'HDF5 error'
-  call h5dclose_f(dset_id, hdferr)
+  call H5Dclose_f(dset_id, hdferr)
   if(hdferr < 0) error stop 'HDF5 error'
-  call h5tclose_f(dtype_id, hdferr)
+  call H5Tclose_f(dtype_id, hdferr)
   if(hdferr < 0) error stop 'HDF5 error'
-  call h5tclose_f(label_id, hdferr)
+  call H5Tclose_f(label_id, hdferr)
   if(hdferr < 0) error stop 'HDF5 error'
-  call h5tclose_f(entry_id, hdferr)
+  call H5Tclose_f(entry_id, hdferr)
   if(hdferr < 0) error stop 'HDF5 error'
 
   call executionStamp('cell_to/homogenization','cell ID to homogenization results')

src/signals.f90

@@ -6,7 +6,7 @@ module signals
   use prec
   use system_routines
 
-  implicit none
+  implicit none(type,external)
   private
 
   logical, volatile, public, protected :: &

src/system_routines.f90

@@ -7,7 +7,7 @@ module system_routines
 
   use prec
 
-  implicit none
+  implicit none(type,external)
   private
 
   public :: &
@@ -24,8 +24,10 @@ module system_routines
 
   interface
 
+
     function setCWD_C(cwd) bind(C)
       use, intrinsic :: ISO_C_Binding, only: C_INT, C_CHAR
+      implicit none(type,external)
 
       integer(C_INT) :: setCWD_C
       character(kind=C_CHAR), dimension(*), intent(in) :: cwd
@@ -34,6 +36,7 @@ module system_routines
     subroutine getCWD_C(cwd, stat) bind(C)
       use, intrinsic :: ISO_C_Binding, only: C_INT, C_CHAR
       use prec
+      implicit none(type,external)
 
       character(kind=C_CHAR), dimension(pPathLen+1), intent(out) :: cwd                             ! NULL-terminated array
       integer(C_INT),                                intent(out) :: stat
@@ -42,6 +45,7 @@ module system_routines
     subroutine getHostName_C(hostname, stat) bind(C)
       use, intrinsic :: ISO_C_Binding, only: C_INT, C_CHAR
       use prec
+      implicit none(type,external)
 
       character(kind=C_CHAR), dimension(pStringLen+1), intent(out) :: hostname                        ! NULL-terminated array
       integer(C_INT),                                  intent(out) :: stat
@@ -50,6 +54,7 @@ module system_routines
     subroutine getUserName_C(username, stat) bind(C)
       use, intrinsic :: ISO_C_Binding, only: C_INT, C_CHAR
       use prec
+      implicit none(type,external)
 
       character(kind=C_CHAR), dimension(pStringLen+1), intent(out) :: username                        ! NULL-terminated array
       integer(C_INT),                                  intent(out) :: stat
@@ -57,27 +62,31 @@ module system_routines
 
     subroutine signalint_C(handler) bind(C)
       use, intrinsic :: ISO_C_Binding, only: C_FUNPTR
+      implicit none(type,external)
 
       type(C_FUNPTR), intent(in), value :: handler
     end subroutine signalint_C
 
     subroutine signalusr1_C(handler) bind(C)
       use, intrinsic :: ISO_C_Binding, only: C_FUNPTR
+      implicit none(type,external)
 
       type(C_FUNPTR), intent(in), value :: handler
     end subroutine signalusr1_C
 
     subroutine signalusr2_C(handler) bind(C)
       use, intrinsic :: ISO_C_Binding, only: C_FUNPTR
+      implicit none(type,external)
 
       type(C_FUNPTR), intent(in), value :: handler
     end subroutine signalusr2_C
 
     subroutine free_C(ptr) bind(C,name='free')
-      import c_ptr
-      type(c_ptr), value :: ptr
-    end subroutine free_C
+      use, intrinsic :: ISO_C_Binding, only: C_PTR
+      implicit none(type,external)
 
+      type(C_PTR), value :: ptr
+    end subroutine free_C
 
   end interface