Merge branch 'development' into 44-column-major-access-to-interaction-matrices

.gitignore

@@ -1,6 +1,7 @@
 *.pyc
 *.mod
 *.o
+*.hdf5
 *.exe
 *.bak
 *~

.gitlab-ci.yml

@@ -343,6 +343,15 @@ Spectral_MPI:
     - master
     - release
 
+SpectralAll_restartMPI:
+  stage: spectral
+  script:
+    - module load $IntelCompiler $MPICH_Intel $PETSc_MPICH_Intel
+    - SpectralAll_restartMPI/test.py
+  except:
+    - master
+    - release
+
 Plasticity_DetectChanges:
   stage: spectral
   script: Plasticity_DetectChanges/test.py

VERSION

@@ -1 +1 @@
-v2.0.2-1277-g53bc24cc
+v2.0.2-1389-g070952db

installation/mods_Abaqus/abaqus_v6.env

@@ -12,7 +12,15 @@
 #
 import os, re, glob, driverUtils
 from damask import version as DAMASKVERSION
+from damask import Environment
+myEnv = Environment()
 
+if myEnv.options['DAMASK_HDF5'] == 'ON':
+  # use hdf5 compiler wrapper in $PATH
+  fortCmd   = os.popen('h5fc -shlib -show').read().replace('\n','')                                 # complicated way needed to pass in DAMASKVERSION string
+  link_sl  += fortCmd.split()[1:]
+  fortCmd  +=" -DDAMASKHDF5"
+else:
   # Use the version in $PATH
   fortCmd = "ifort"
 
@@ -50,4 +58,6 @@ ask_delete=OFF
 
 # Remove the temporary names from the namespace
 del fortCmd
+del Environment
+del myEnv
 del DAMASKVERSION

installation/mods_Abaqus/abaqus_v6_debug.env

@@ -12,7 +12,15 @@
 #
 import os, re, glob, driverUtils
 from damask import version as DAMASKVERSION
+from damask import Environment
+myEnv = Environment()
 
+if myEnv.options['DAMASK_HDF5'] == 'ON':
+  # use hdf5 compiler wrapper in $PATH
+  fortCmd   = os.popen('h5fc -shlib -show').read().replace('\n','')                                 # complicated way needed to pass in DAMASKVERSION string
+  link_sl  += fortCmd.split()[1:]
+  fortCmd  +=" -DDAMASKHDF5"
+else:
   # Use the version in $PATH
   fortCmd = "ifort"
 
@@ -55,4 +63,6 @@ ask_delete=OFF
 
 # Remove the temporary names from the namespace
 del fortCmd
+del Environment
+del myEnv
 del DAMASKVERSION

installation/mods_MarcMentat/2018.1/Marc_tools/include_linux64

@@ -63,7 +63,6 @@ else
  INTEGER_PATH=/$MARC_INTEGER_SIZE
 fi
 
-FCOMP=ifort
 INTELPATH="/opt/intel/compilers_and_libraries_2017/linux"
 
 # find the root directory of the compiler installation:
@@ -99,6 +98,16 @@ else
  FCOMPROOT=
 fi
 
+# DAMASK uses the HDF5 compiler wrapper around the Intel compiler
+if test "$DAMASK_HDF5" = "ON";then
+  H5FC="$(h5fc -shlib -show)"
+  HDF5_LIB=${H5FC//ifort/}
+  FCOMP="$H5FC -DDAMASKHDF5"
+  echo $FCOMP
+else
+  FCOMP=ifort
+fi
+
 # AEM
 if test "$MARCDLLOUTDIR" = ""; then
  DLLOUTDIR="$MARC_LIB"
@@ -535,23 +544,17 @@ else
   DAMASKVERSION="'N/A'"
 fi
 
-if test "$DAMASK_HDF5" = "ON";then
-  DFCOMP="$(h5fc -show) -DDAMASKHDF5"
-else
-  DFCOMP=$FCOMP
-fi
 
-#
 # DAMASK compiler calls: additional flags are in line 2 OpenMP flags in line 3
-DFORTLOWMP="$DFCOMP -c -implicitnone -stand f08 -standard-semantics -assume nostd_mod_proc_name -safe_cray_ptr $PROFILE -zero -mp1 -WB  -O0 $I8FFLAGS -I$MARC_SOURCE/common \
+DFORTLOWMP="$FCOMP -c -implicitnone -stand f08 -standard-semantics -assume nostd_mod_proc_name -safe_cray_ptr $PROFILE -zero -mp1 -WB  -O0 $I8FFLAGS -I$MARC_SOURCE/common \
  -fpp -ftz -diag-disable 5268 -warn declarations -warn general -warn usage -warn interfaces -warn ignore_loc -warn alignments -DMarc4DAMASK=2018.1 -DDAMASKVERSION=$DAMASKVERSION \
  -qopenmp -qopenmp-threadprivate=compat\
   $MUMPS_INCLUDE $I8DEFINES -DLinux -DLINUX -DLinux_intel $FDEFINES $DDM $SOLVERFLAGS -I$KDTREE2_MOD"
-DFORTRANMP="$DFCOMP -c -implicitnone -stand f08 -standard-semantics -assume nostd_mod_proc_name -safe_cray_ptr $PROFILE -zero -mp1 -WB  -O1 $I8FFLAGS -I$MARC_SOURCE/common \
+DFORTRANMP="$FCOMP -c -implicitnone -stand f08 -standard-semantics -assume nostd_mod_proc_name -safe_cray_ptr $PROFILE -zero -mp1 -WB  -O1 $I8FFLAGS -I$MARC_SOURCE/common \
  -fpp -ftz -diag-disable 5268 -warn declarations -warn general -warn usage -warn interfaces -warn ignore_loc -warn alignments -DMarc4DAMASK=2018.1 -DDAMASKVERSION=$DAMASKVERSION \
  -qopenmp -qopenmp-threadprivate=compat\
   $MUMPS_INCLUDE $I8DEFINES -DLinux -DLINUX -DLinux_intel $FDEFINES $DDM $SOLVERFLAGS -I$KDTREE2_MOD"
-DFORTHIGHMP="$DFCOMP -c -implicitnone -stand f08 -standard-semantics -assume nostd_mod_proc_name -safe_cray_ptr $PROFILE -zero -mp1 -WB  -fno-alias -O2 $I8FFLAGS -I$MARC_SOURCE/common \
+DFORTHIGHMP="$FCOMP -c -implicitnone -stand f08 -standard-semantics -assume nostd_mod_proc_name -safe_cray_ptr $PROFILE -zero -mp1 -WB  -fno-alias -O2 $I8FFLAGS -I$MARC_SOURCE/common \
  -fpp -ftz -diag-disable 5268 -warn declarations -warn general -warn usage -warn interfaces -warn ignore_loc -warn alignments -DMarc4DAMASK=2018.1 -DDAMASKVERSION=$DAMASKVERSION \
  -qopenmp -qopenmp-threadprivate=compat\
   $MUMPS_INCLUDE $I8DEFINES -DLinux -DLINUX -DLinux_intel $FDEFINES $DDM $SOLVERFLAGS -I$KDTREE2_MOD"
@@ -570,15 +573,15 @@ then
 fi
 
       # DAMASK compiler calls: additional flags are in line 2 OpenMP flags in line 3
-  DFORTLOWMP="$DFCOMP -c -implicitnone -stand f08 -standard-semantics -assume nostd_mod_proc_name -safe_cray_ptr $PROFILE -zero -mp1 -WB  $I8FFLAGS -I$MARC_SOURCE/common \
+  DFORTLOWMP="$FCOMP -c -implicitnone -stand f08 -standard-semantics -assume nostd_mod_proc_name -safe_cray_ptr $PROFILE -zero -mp1 -WB  $I8FFLAGS -I$MARC_SOURCE/common \
    -fpp -ftz -diag-disable 5268 -warn declarations -warn general -warn usage -warn interfaces -warn ignore_loc -warn alignments -DMarc4DAMASK=2018.1 -DDAMASKVERSION=$DAMASKVERSION \
    -qopenmp -qopenmp-threadprivate=compat\
     $MUMPS_INCLUDE $I8DEFINES -DLinux -DLINUX -DLinux_intel $FDEFINES $DDM $SOLVERFLAGS -I$KDTREE2_MOD"
-  DFORTRANMP="$DFCOMP -c -implicitnone -stand f08 -standard-semantics -assume nostd_mod_proc_name -safe_cray_ptr $PROFILE -zero -mp1 -WB  $I8FFLAGS -I$MARC_SOURCE/common \
+  DFORTRANMP="$FCOMP -c -implicitnone -stand f08 -standard-semantics -assume nostd_mod_proc_name -safe_cray_ptr $PROFILE -zero -mp1 -WB  $I8FFLAGS -I$MARC_SOURCE/common \
    -fpp -ftz -diag-disable 5268 -warn declarations -warn general -warn usage -warn interfaces -warn ignore_loc -warn alignments -DMarc4DAMASK=2018.1 -DDAMASKVERSION=$DAMASKVERSION \
    -qopenmp -qopenmp-threadprivate=compat\
     $MUMPS_INCLUDE $I8DEFINES  -DLinux -DLINUX -DLinux_intel $FDEFINES $DDM $SOLVERFLAGS -I$KDTREE2_MOD"
-  DFORTHIGHMP="$DFCOMP -c -implicitnone -stand f08 -standard-semantics -assume nostd_mod_proc_name -safe_cray_ptr $PROFILE -zero -mp1 -WB  -fno-alias $I8FFLAGS -I$MARC_SOURCE/common \
+  DFORTHIGHMP="$FCOMP -c -implicitnone -stand f08 -standard-semantics -assume nostd_mod_proc_name -safe_cray_ptr $PROFILE -zero -mp1 -WB  -fno-alias $I8FFLAGS -I$MARC_SOURCE/common \
    -fpp -ftz -diag-disable 5268 -warn declarations -warn general -warn usage -warn interfaces -warn ignore_loc -warn alignments -DMarc4DAMASK=2018.1 -DDAMASKVERSION=$DAMASKVERSION \
    -qopenmp -qopenmp-threadprivate=compat\
     $MUMPS_INCLUDE $I8DEFINES -DLinux -DLINUX -DLinux_intel $FDEFINES $DDM $SOLVERFLAGS -I$KDTREE2_MOD"
@@ -744,7 +747,7 @@ SECLIBS="-L$MARC_LIB -llapi"
 
 SOLVERLIBS="${BCSSOLVERLIBS} ${VKISOLVERLIBS} ${CASISOLVERLIBS} ${MF2SOLVERLIBS} \
   $MKLLIB -L$MARC_MKL -liomp5 \
-  $MARC_LIB/blas_src.a ${ACSI_LIB}/ACSI_MarcLib.a $KDTREE2_LIB/kdtree2.a "
+  $MARC_LIB/blas_src.a ${ACSI_LIB}/ACSI_MarcLib.a $KDTREE2_LIB/kdtree2.a $HDF5_LIB "
 
 SOLVERLIBS_DLL=${SOLVERLIBS}
 if test "$AEM_DLL" -eq 1

processing/post/ascii2hdf5.py

@@ -1,142 +0,0 @@
-#!/usr/bin/env python2.7
-# -*- coding: UTF-8 no BOM -*-
-
-# ------------------------------------------------------------------- #
-# NOTE:                                                               #
-#   1. Not all output is defined in the DS_HDF5.xml, please add new   #
-#      new one to the system wide definition file                     #
-#             <DAMASK_ROOT>/lib/damask/DS_HDF5.xml                    #
-#      or specify your own when initializing HDF5 class               #
-#   2. Somehow the point cloud structure cannot be properly handled   #
-#      by Xdmf, which is a descriptive wrapper for visualizing HDF5   #
-#      using Paraview. The current solution is using cell structured  #
-#      HDF5 so that Xdmf can describe the data shape as a rectangular #
-#      mesh rather than polyvertex.                                   #
-# TODO:                                                               #
-#   1. remove the <ASCII_TABLE>._tmp file, basically need a way to    #
-#      just load data from ASCII table.                               #
-#   2. a progress monitor when transferring data from ASCII table     #
-#      to HDF5.                                                       #
-#   3. a more flexible way handle the data structure rather than a    #
-#      xml file.                                                      #
-# ------------------------------------------------------------------- #
-
-import os
-import damask
-import numpy as np
-from optparse import OptionParser
-
-
-scriptName = os.path.splitext(os.path.basename(__file__))[0]
-scriptID = ' '.join([scriptName, damask.version])
-
-
-# ----- helper function ----- #
-def get_rectMshVectors(xyz_array, posNum):
-        """Get Vx, Vy, Vz for rectLinear grid"""
-        # need some improvement, and only works for rectangular grid
-        v = sorted(list(set(xyz_array[:, posNum])))
-        v_interval = (v[2]+v[1])/2.0 - (v[1]+v[0])/2.0
-        v_start = (v[1]+v[0])/2.0 - v_interval
-        v_end = (v[-1]+v[-2])/2.0 + v_interval
-        V = np.linspace(v_start, v_end, len(v)+1)
-        return V
-
-
-# ----- MAIN ---- #
-desp_msg = "Convert DAMASK ascii table to HDF5 file"
-parser = OptionParser(option_class=damask.extendableOption,
-                      usage='%prog options [file[s]]',
-                      description=desp_msg,
-                      version=scriptID)
-parser.add_option('-D', '--DefinitionFile',
-                  dest='storage definition file',
-                  type='string',
-                  metavar='string',
-                  help='definition file for H5 data storage')
-parser.add_option('-p', '--pos', '--position',
-                  dest='pos',
-                  type='string', metavar='string',
-                  help='label of coordinates [%default]')
-
-parser.set_defaults(DefinitionFile='default',
-                    pos='pos')
-
-(options, filenames) = parser.parse_args()
-
-filename = filenames[0]
-
-if options.DefinitionFile == 'default':
-        defFile = None
-else:
-        defFile = options.DefinitionFile
-
-# ----- read in data using DAMASK ASCII table class ----- #
-asciiTable = damask.ASCIItable(name=filename, buffered=False)
-asciiTable.head_read()
-asciiTable.data_readArray()
-incNum = int(asciiTable.data[asciiTable.label_index('inc'), 0])
-fullTable = np.copy(asciiTable.data)  # deep copy all data, just to be safe
-labels = asciiTable.labels()
-labels_idx = [asciiTable.label_index(label) for label in labels]
-featuresDim = [labels_idx[i+1] - labels_idx[i] for i in range(len(labels)-1)]
-featuresDim.append(fullTable.shape[1] - labels_idx[-1])
-
-# ----- figure out size and grid ----- #
-pos_idx = asciiTable.label_index('pos')
-xyz_array = asciiTable.data[:, pos_idx:pos_idx+3]
-Vx = get_rectMshVectors(xyz_array, 0)
-Vy = get_rectMshVectors(xyz_array, 1)
-Vz = get_rectMshVectors(xyz_array, 2)
-# use the dimension of the rectangular grid to reshape all other data
-mshGridDim = [len(Vx)-1, len(Vy)-1, len(Vz)-1]
-
-# ----- compose cmd log ----- #
-cmd_log = " ".join([scriptID, filename])
-
-# ----- create a new HDF5 file and save the data -----#
-# force remove existing HDF5 file
-h5fName = filename.replace(".txt", ".h5")
-try:
-        os.remove(h5fName)
-except OSError:
-        pass
-h5f = damask.H5Table(h5fName,
-                     new_file=True,
-                     dsXMLFile=defFile)
-# adding increment number as root level attributes
-h5f.add_attr('inc', incNum)
-# add the mesh grid data now
-h5f.add_data("Vx", Vx, cmd_log=cmd_log)
-h5f.add_data("Vy", Vy, cmd_log=cmd_log)
-h5f.add_data("Vz", Vz, cmd_log=cmd_log)
-
-# add the rest of data from table
-labelsProcessed = ['inc']
-for fi in range(len(labels)):
-        featureName = labels[fi]
-        # remove trouble maker "("" and ")" from label/feature name
-        if "(" in featureName:
-            featureName = featureName.replace("(", "")
-        if ")" in featureName:
-            featureName = featureName.replace(")", "")
-        # skip increment and duplicated columns in the ASCII table
-        if featureName in labelsProcessed:
-            continue
-
-        featureIdx = labels_idx[fi]
-        featureDim = featuresDim[fi]
-        # grab the data hook
-        dataset = fullTable[:, featureIdx:featureIdx+featureDim]
-        # mapping 2D data onto a 3D rectangular mesh to get 4D data
-        # WARNING: In paraview, the data for a recmesh is mapped as:
-        # -->  len(z), len(y), len(x), size(data)
-        # dataset = dataset.reshape((mshGridDim[0],
-        #                            mshGridDim[1],
-        #                            mshGridDim[2],
-        #                            dataset.shape[1]))
-        # write out data
-        print("adding {}...".format(featureName))
-        h5f.add_data(featureName, dataset, cmd_log=cmd_log)
-        # write down the processed label
-        labelsProcessed.append(featureName)

src/CMakeLists.txt

@@ -29,24 +29,28 @@ add_library(IO OBJECT "IO.f90")
 add_dependencies(IO DAMASK_INTERFACE)
 list(APPEND OBJECTFILES $<TARGET_OBJECTS:IO>)
 
-add_library(HDF5_UTILITIES OBJECT "HDF5_utilities.f90")
-add_dependencies(HDF5_UTILITIES IO)
-list(APPEND OBJECTFILES $<TARGET_OBJECTS:HDF5_UTILITIES>)
-
 add_library(NUMERICS OBJECT "numerics.f90")
-add_dependencies(NUMERICS HDF5_UTILITIES)
+add_dependencies(NUMERICS IO)
 list(APPEND OBJECTFILES $<TARGET_OBJECTS:NUMERICS>)
 
 add_library(DEBUG OBJECT "debug.f90")
 add_dependencies(DEBUG NUMERICS)
 list(APPEND OBJECTFILES $<TARGET_OBJECTS:DEBUG>)
 
-add_library(CONFIG OBJECT "config.f90")
+add_library(DAMASK_CONFIG OBJECT "config.f90")
-add_dependencies(CONFIG DEBUG)
+add_dependencies(DAMASK_CONFIG DEBUG)
-list(APPEND OBJECTFILES $<TARGET_OBJECTS:CONFIG>)
+list(APPEND OBJECTFILES $<TARGET_OBJECTS:DAMASK_CONFIG>)
+
+add_library(HDF5_UTILITIES OBJECT "HDF5_utilities.f90")
+add_dependencies(HDF5_UTILITIES DAMASK_CONFIG)
+list(APPEND OBJECTFILES $<TARGET_OBJECTS:HDF5_UTILITIES>)
+
+add_library(RESULTS OBJECT "results.f90")
+add_dependencies(RESULTS HDF5_UTILITIES)
+list(APPEND OBJECTFILES $<TARGET_OBJECTS:RESULTS>)
 
 add_library(FEsolving OBJECT "FEsolving.f90")
-add_dependencies(FEsolving DEBUG)
+add_dependencies(FEsolving RESULTS)
 list(APPEND OBJECTFILES $<TARGET_OBJECTS:FEsolving>)
 
 add_library(DAMASK_MATH OBJECT "math.f90")
@@ -68,7 +72,7 @@ elseif (PROJECT_NAME STREQUAL "DAMASK_FEM")
 endif()
 
 add_library(MATERIAL OBJECT "material.f90")
-add_dependencies(MATERIAL MESH CONFIG)
+add_dependencies(MATERIAL MESH DAMASK_CONFIG)
 list(APPEND OBJECTFILES $<TARGET_OBJECTS:MATERIAL>)
 
 add_library(DAMASK_HELPERS OBJECT "lattice.f90")

src/CPFEM2.f90

@@ -10,8 +10,8 @@ module CPFEM2
 
  public :: &
    CPFEM_age, &
-   CPFEM_initAll
+   CPFEM_initAll, &
-
+   CPFEM_results
 contains
 
 
@@ -20,8 +20,7 @@ contains
 !--------------------------------------------------------------------------------------------------
 subroutine CPFEM_initAll()
  use prec, only: &
-   pInt
+   pInt, &
- use prec, only: &
    prec_init
  use numerics, only: &
    numerics_init
@@ -39,6 +38,8 @@ subroutine CPFEM_initAll()
    material_init
  use HDF5_utilities, only: &
    HDF5_utilities_init
+ use results, only: &
+   results_init
  use lattice, only: &
    lattice_init
  use constitutive, only: &
@@ -73,6 +74,7 @@ subroutine CPFEM_initAll()
  call lattice_init
  call material_init
  call HDF5_utilities_init
+ call results_init
  call constitutive_init
  call crystallite_init
  call homogenization_init
@@ -105,8 +107,7 @@ subroutine CPFEM_init
    debug_levelBasic, &
    debug_levelExtensive
  use FEsolving, only: &
-   restartRead, &
+   restartRead
-   modelName
  use material, only: &
    material_phase, &
    homogState, &
@@ -125,23 +126,22 @@ subroutine CPFEM_init
  use hdf5
  use HDF5_utilities, only: &
    HDF5_openFile, &
-   HDF5_openGroup2, &
+   HDF5_closeFile, &
+   HDF5_openGroup, &
+   HDF5_closeGroup, &
    HDF5_read
  use DAMASK_interface, only: &
    getSolverJobName
 
  implicit none
- integer(pInt) :: k,l,m,ph,homog
+ integer(pInt) :: ph,homog
  character(len=1024) :: rankStr, PlasticItem, HomogItem
- integer(HID_T) :: fileReadID, groupPlasticID, groupHomogID
+ integer(HID_T) :: fileHandle, groupPlasticID, groupHomogID
- integer        :: hdferr
 
- mainProcess: if (worldrank == 0) then 
  write(6,'(/,a)')   ' <<<+-  CPFEM init  -+>>>'
  write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
 #include "compilation_info.f90"
  flush(6)
- endif mainProcess
 
  ! *** restore the last converged values of each essential variable from the binary file
  if (restartRead) then
@@ -152,34 +152,39 @@ subroutine CPFEM_init
    
    write(rankStr,'(a1,i0)')'_',worldrank
 
-   fileReadID = HDF5_openFile(trim(getSolverJobName())//trim(rankStr)//'.hdf5')
+   fileHandle = HDF5_openFile(trim(getSolverJobName())//trim(rankStr)//'.hdf5')
   
-   call HDF5_read(material_phase,      fileReadID,'recordedPhase')
+   call HDF5_read(fileHandle,material_phase,      'recordedPhase')
-   call HDF5_read(crystallite_F0,      fileReadID,'convergedF')
+   call HDF5_read(fileHandle,crystallite_F0,      'convergedF')
-   call HDF5_read(crystallite_Fp0,     fileReadID,'convergedFp')
+   call HDF5_read(fileHandle,crystallite_Fp0,     'convergedFp')
-   call HDF5_read(crystallite_Fi0,     fileReadID,'convergedFi')
+   call HDF5_read(fileHandle,crystallite_Fi0,     'convergedFi')
-   call HDF5_read(crystallite_Lp0,     fileReadID,'convergedLp')
+   call HDF5_read(fileHandle,crystallite_Lp0,     'convergedLp')
-   call HDF5_read(crystallite_Li0,     fileReadID,'convergedLi')
+   call HDF5_read(fileHandle,crystallite_Li0,     'convergedLi')
-   call HDF5_read(crystallite_dPdF0,   fileReadID,'convergeddPdF')
+   call HDF5_read(fileHandle,crystallite_dPdF0,   'convergeddPdF')
-   call HDF5_read(crystallite_Tstar0_v,fileReadID,'convergedTstar')
+   call HDF5_read(fileHandle,crystallite_Tstar0_v,'convergedTstar')
    
-   groupPlasticID = HDF5_openGroup2(fileReadID,'PlasticPhases')
+   groupPlasticID = HDF5_openGroup(fileHandle,'PlasticPhases')
    do ph = 1_pInt,size(phase_plasticity)
     write(PlasticItem,*) ph,'_'
-    call HDF5_read(plasticState(ph)%state0,groupPlasticID,trim(PlasticItem)//'convergedStateConst')
+    call HDF5_read(groupPlasticID,plasticState(ph)%state0,trim(PlasticItem)//'convergedStateConst')
    enddo
+   call HDF5_closeGroup(groupPlasticID)
    
-   groupHomogID = HDF5_openGroup2(fileReadID,'HomogStates')
+   groupHomogID = HDF5_openGroup(fileHandle,'HomogStates')
    do homog = 1_pInt, material_Nhomogenization
     write(HomogItem,*) homog,'_'
-    call HDF5_read(homogState(homog)%state0, groupHomogID,trim(HomogItem)//'convergedStateHomog')
+    call HDF5_read(groupHomogID,homogState(homog)%state0, trim(HomogItem)//'convergedStateHomog')
    enddo
+   call HDF5_closeGroup(groupHomogID)
+
+   call HDF5_closeFile(fileHandle)
    
    restartRead = .false.
  endif
 
 end subroutine CPFEM_init
 
+
 !--------------------------------------------------------------------------------------------------
 !> @brief forwards data after successful increment
 !--------------------------------------------------------------------------------------------------
@@ -229,31 +234,28 @@ subroutine CPFEM_age()
  use HDF5_utilities, only: &
    HDF5_openFile, &
    HDF5_closeFile, &
+   HDF5_addGroup, &
    HDF5_closeGroup, &
-   HDF5_addGroup2, &
    HDF5_write
  use hdf5
  use DAMASK_interface, only: &
    getSolverJobName
 
  implicit none
-
+ integer(pInt)    ::  i, ph, homog, mySource
- integer(pInt) ::                                    i, k, l, m, ph, homog, mySource
  character(len=32) :: rankStr, PlasticItem, HomogItem
  integer(HID_T) :: fileHandle, groupPlastic, groupHomog
- integer        :: hdferr
- integer(HSIZE_T)  :: hdfsize
 
  if (iand(debug_level(debug_CPFEM), debug_levelBasic) /= 0_pInt) &
    write(6,'(a)') '<< CPFEM >> aging states'
 
-  crystallite_F0  = crystallite_partionedF                                                    ! crystallite deformation (_subF is perturbed...)
+ crystallite_F0  = crystallite_partionedF
-  crystallite_Fp0 = crystallite_Fp                                                            ! crystallite plastic deformation
+ crystallite_Fp0 = crystallite_Fp
-  crystallite_Lp0 = crystallite_Lp                                                            ! crystallite plastic velocity
+ crystallite_Lp0 = crystallite_Lp
-  crystallite_Fi0 = crystallite_Fi                                                            ! crystallite intermediate deformation
+ crystallite_Fi0 = crystallite_Fi
-  crystallite_Li0 = crystallite_Li                                                            ! crystallite intermediate velocity
+ crystallite_Li0 = crystallite_Li
-  crystallite_dPdF0 = crystallite_dPdF                                                        ! crystallite stiffness
+ crystallite_dPdF0 = crystallite_dPdF
-  crystallite_Tstar0_v = crystallite_Tstar_v                                                  ! crystallite 2nd Piola Kirchhoff stress
+ crystallite_Tstar0_v = crystallite_Tstar_v
  
  forall (i = 1:size(plasticState)) plasticState(i)%state0 = plasticState(i)%state            ! copy state in this lengthy way because: A component cannot be an array if the encompassing structure is an array
  
@@ -271,30 +273,30 @@ if (iand(debug_level(debug_CPFEM), debug_levelBasic) /= 0_pInt) &
  if (restartWrite) then
    if (iand(debug_level(debug_CPFEM), debug_levelBasic) /= 0_pInt) &
      write(6,'(a)') '<< CPFEM >> writing restart variables of last converged step to hdf5 file'
-  write(rankStr,'(a1,i0)')'_',worldrank
    
+   write(rankStr,'(a1,i0)')'_',worldrank
    fileHandle = HDF5_openFile(trim(getSolverJobName())//trim(rankStr)//'.hdf5','w')
    
-  call HDF5_write(material_phase,      fileHandle,'recordedPhase')
+   call HDF5_write(fileHandle,material_phase,      'recordedPhase')
-  call HDF5_write(crystallite_F0,      fileHandle,'convergedF')
+   call HDF5_write(fileHandle,crystallite_F0,      'convergedF')
-  call HDF5_write(crystallite_Fp0,     fileHandle,'convergedFp')
+   call HDF5_write(fileHandle,crystallite_Fp0,     'convergedFp')
-  call HDF5_write(crystallite_Fi0,     fileHandle,'convergedFi')
+   call HDF5_write(fileHandle,crystallite_Fi0,     'convergedFi')
-  call HDF5_write(crystallite_Lp0,     fileHandle,'convergedLp')
+   call HDF5_write(fileHandle,crystallite_Lp0,     'convergedLp')
-  call HDF5_write(crystallite_Li0,     fileHandle,'convergedLi')
+   call HDF5_write(fileHandle,crystallite_Li0,     'convergedLi')
-  call HDF5_write(crystallite_dPdF0,   fileHandle,'convergeddPdF')
+   call HDF5_write(fileHandle,crystallite_dPdF0,   'convergeddPdF')
-  call HDF5_write(crystallite_Tstar0_v,fileHandle,'convergedTstar')
+   call HDF5_write(fileHandle,crystallite_Tstar0_v,'convergedTstar')
    
-  groupPlastic = HDF5_addGroup2(fileHandle,'PlasticPhases')
+   groupPlastic = HDF5_addGroup(fileHandle,'PlasticPhases')
    do ph = 1_pInt,size(phase_plasticity)
      write(PlasticItem,*) ph,'_'
-    call HDF5_write(plasticState(ph)%state0,groupPlastic,trim(PlasticItem)//'convergedStateConst')
+     call HDF5_write(groupPlastic,plasticState(ph)%state0,trim(PlasticItem)//'convergedStateConst')
    enddo
    call HDF5_closeGroup(groupPlastic)
 
-  groupHomog = HDF5_addGroup2(fileHandle,'HomogStates')
+   groupHomog = HDF5_addGroup(fileHandle,'HomogStates')
    do homog = 1_pInt, material_Nhomogenization
      write(HomogItem,*) homog,'_'
-    call HDF5_write(homogState(homog)%state0,groupHomog,trim(HomogItem)//'convergedStateHomog')
+     call HDF5_write(groupHomog,homogState(homog)%state0,trim(HomogItem)//'convergedStateHomog')
    enddo
    call HDF5_closeGroup(groupHomog)
    
@@ -307,4 +309,27 @@ if (iand(debug_level(debug_CPFEM), debug_levelBasic) /= 0_pInt) &
 
 end subroutine CPFEM_age
 
+!--------------------------------------------------------------------------------------------------
+!> @brief triggers writing of the results
+!--------------------------------------------------------------------------------------------------
+subroutine CPFEM_results(inc,time)
+ use prec, only: &
+   pInt
+ use results
+ use HDF5_utilities
+ use constitutive, only: &
+   constitutive_results
+
+ implicit none
+ integer(pInt), intent(in) :: inc
+ real(pReal),   intent(in) :: time
+
+ call results_openJobFile
+ call results_addIncrement(inc,time)
+ call constitutive_results()
+ call results_removeLink('current') ! ToDo: put this into closeJobFile
+ call results_closeJobFile
+
+end subroutine CPFEM_results
+
 end module CPFEM2

src/DAMASK_spectral.f90

@@ -46,7 +46,8 @@ program DAMASK_spectral
    grid, &
    geomSize
  use CPFEM2, only: &
-   CPFEM_initAll
+   CPFEM_initAll, &
+   CPFEM_results
  use FEsolving, only: &
    restartWrite, &
    restartInc
@@ -80,6 +81,7 @@ program DAMASK_spectral
  use spectral_mech_Polarisation
  use spectral_damage
  use spectral_thermal
+ use results
 
  implicit none
 
@@ -157,6 +159,9 @@ program DAMASK_spectral
  write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
 #include "compilation_info.f90"
 
+
+ call results_openJobFile()
+ call results_closeJobFile()
 !--------------------------------------------------------------------------------------------------
 ! initialize field solver information
  nActiveFields = 1
@@ -420,6 +425,7 @@ program DAMASK_spectral
 
  writeUndeformed: if (interface_restartInc < 1_pInt) then
    write(6,'(1/,a)') ' ... writing initial configuration to file ........................'
+   call CPFEM_results(0_pInt,0.0_pReal)
    do i = 1, size(materialpoint_results,3)/(maxByteOut/(materialpoint_sizeResults*pReal))+1         ! slice the output of my process in chunks not exceeding the limit for one output
      outputIndex = int([(i-1_pInt)*((maxRealOut)/materialpoint_sizeResults)+1_pInt, &               ! QUESTION: why not starting i at 0 instead of murky 1?
                              min(i*((maxRealOut)/materialpoint_sizeResults),size(materialpoint_results,3))],pLongInt)
@@ -596,6 +602,7 @@ program DAMASK_spectral
            if(ierr /=0_pInt) call IO_error(894_pInt, ext_msg='MPI_file_write')
          enddo
          fileOffset = fileOffset + sum(outputSize)                                                  ! forward to current file position
+         call CPFEM_results(totalIncsCounter,time)
        endif
        if (              loadCases(currentLoadCase)%restartFrequency > 0_pInt &                     ! writing of restart info requested ...
            .and. mod(inc,loadCases(currentLoadCase)%restartFrequency) == 0_pInt) then               ! ... and at frequency of writing restart information

src/commercialFEM_fileList.f90

@@ -4,12 +4,12 @@
 !> @details List of files needed by MSC.Marc, Abaqus/Explicit, and Abaqus/Standard
 !--------------------------------------------------------------------------------------------------
 #include "IO.f90"
-#ifdef DAMASKHDF5
-#include "HDF5_utilities.f90"
-#endif
 #include "numerics.f90"
 #include "debug.f90"
 #include "config.f90"
+#ifdef DAMASKHDF5
+#include "HDF5_utilities.f90"
+#endif
 #include "math.f90"
 #include "FEsolving.f90"
 #include "mesh.f90"

src/config.f90

@@ -1,4 +1,4 @@
-!--------------------------------------------------------------------------------------------------
+!-------------------------------------------------------------------------------------------------
 !> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
 !> @brief Reads in the material configuration from file
 !> @details Reads the material configuration file, where solverJobName.materialConfig takes
@@ -220,7 +220,7 @@ subroutine parseFile(sectionNames,part,line, &
  partPosition = [partPosition, i]                                                                   ! needed when actually storing content
 
  do i = 1_pInt, size(partPosition) -1_pInt
-   sectionNames(i) = trim(adjustl(fileContent(partPosition(i))))
+   sectionNames(i) = trim(adjustl(IO_getTag(fileContent(partPosition(i)),'[',']')))
    do j = partPosition(i) + 1_pInt,  partPosition(i+1) -1_pInt
      call part(i)%add(trim(adjustl(fileContent(j))))
    enddo

src/constitutive.f90

@@ -25,7 +25,8 @@ module constitutive
    constitutive_SandItsTangents, &
    constitutive_collectDotState, &
    constitutive_collectDeltaState, &
-   constitutive_postResults
+   constitutive_postResults, &
+   constitutive_results
 
  private :: &
    constitutive_hooke_SandItsTangents
@@ -149,7 +150,7 @@ subroutine constitutive_init()
  if (any(phase_plasticity == PLASTICITY_NONE_ID))          call plastic_none_init
  if (any(phase_plasticity == PLASTICITY_ISOTROPIC_ID))     call plastic_isotropic_init
  if (any(phase_plasticity == PLASTICITY_PHENOPOWERLAW_ID)) call plastic_phenopowerlaw_init
- if (any(phase_plasticity == PLASTICITY_KINEHARDENING_ID)) call plastic_kinehardening_init(FILEUNIT)
+ if (any(phase_plasticity == PLASTICITY_KINEHARDENING_ID)) call plastic_kinehardening_init
  if (any(phase_plasticity == PLASTICITY_DISLOTWIN_ID))     call plastic_dislotwin_init
  if (any(phase_plasticity == PLASTICITY_DISLOUCLA_ID))     call plastic_disloucla_init
  if (any(phase_plasticity == PLASTICITY_NONLOCAL_ID)) then
@@ -489,8 +490,9 @@ subroutine constitutive_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, S6, Fi, ipc, ip, e
      call plastic_phenopowerlaw_LpAndItsTangent   (Lp,dLp_dMp,Mp,instance,of)
 
    case (PLASTICITY_KINEHARDENING_ID) plasticityType
-     call plastic_kinehardening_LpAndItsTangent   (Lp,dLp_dMp99, math_Mandel33to6(Mp),ipc,ip,el)
+     of = phasememberAt(ipc,ip,el)
-     dLp_dMp = math_Plain99to3333(dLp_dMp99)                                                        ! ToDo: We revert here the last statement in plastic_xx_LpAndItsTanget
+     instance = phase_plasticityInstance(material_phase(ipc,ip,el))
+     call plastic_kinehardening_LpAndItsTangent   (Lp,dLp_dMp, Mp,instance,of)
 
    case (PLASTICITY_NONLOCAL_ID) plasticityType
      call plastic_nonlocal_LpAndItsTangent        (Lp,dLp_dMp99, math_Mandel33to6(Mp), &
@@ -873,7 +875,9 @@ subroutine constitutive_collectDotState(S6, FeArray, Fi, FpArray, subdt, subfrac
      call plastic_phenopowerlaw_dotState(Mp,instance,of)
 
    case (PLASTICITY_KINEHARDENING_ID) plasticityType
-     call plastic_kinehardening_dotState(math_Mandel33to6(Mp),ipc,ip,el)
+     of = phasememberAt(ipc,ip,el)
+     instance = phase_plasticityInstance(material_phase(ipc,ip,el))
+     call plastic_kinehardening_dotState(Mp,instance,of)
 
    case (PLASTICITY_DISLOTWIN_ID) plasticityType
      of = phasememberAt(ipc,ip,el)
@@ -929,6 +933,8 @@ subroutine constitutive_collectDeltaState(S6, Fe, Fi, ipc, ip, el)
    math_Mandel33to6, &
    math_mul33x33
  use material, only: &
+   phasememberAt, &
+   phase_plasticityInstance, &
    phase_plasticity, &
    phase_source, &
    phase_Nsources, &
@@ -954,19 +960,22 @@ subroutine constitutive_collectDeltaState(S6, Fe, Fi, ipc, ip, el)
    Fe, &                                                                                            !< elastic deformation gradient
    Fi                                                                                               !< intermediate deformation gradient
  real(pReal),               dimension(3,3) :: &
-   Mstar
+   Mp
  integer(pInt) :: &
-   s                                                                                                !< counter in source loop
+   s, &                                                                                                !< counter in source loop
+   instance, of
 
- Mstar  = math_mul33x33(math_mul33x33(transpose(Fi),Fi),math_Mandel6to33(S6))
+ Mp  = math_mul33x33(math_mul33x33(transpose(Fi),Fi),math_Mandel6to33(S6))
 
  plasticityType: select case (phase_plasticity(material_phase(ipc,ip,el)))
 
    case (PLASTICITY_KINEHARDENING_ID) plasticityType
-     call plastic_kinehardening_deltaState(math_Mandel33to6(Mstar),ipc,ip,el)
+     of = phasememberAt(ipc,ip,el)
+     instance = phase_plasticityInstance(material_phase(ipc,ip,el))
+     call plastic_kinehardening_deltaState(Mp,instance,of)
 
    case (PLASTICITY_NONLOCAL_ID) plasticityType
-     call plastic_nonlocal_deltaState(math_Mandel33to6(Mstar),ip,el)
+     call plastic_nonlocal_deltaState(math_Mandel33to6(Mp),ip,el)
 
  end select plasticityType
 
@@ -1089,8 +1098,10 @@ function constitutive_postResults(S6, Fi, FeArray, ipc, ip, el)
        plastic_phenopowerlaw_postResults(Mp,instance,of)
 
    case (PLASTICITY_KINEHARDENING_ID) plasticityType
+     of = phasememberAt(ipc,ip,el)
+     instance = phase_plasticityInstance(material_phase(ipc,ip,el))
      constitutive_postResults(startPos:endPos) = &
-       plastic_kinehardening_postResults(S6,ipc,ip,el)
+       plastic_kinehardening_postResults(Mp,instance,of)
 
    case (PLASTICITY_DISLOTWIN_ID) plasticityType
      of = phasememberAt(ipc,ip,el)
@@ -1126,4 +1137,43 @@ function constitutive_postResults(S6, Fi, FeArray, ipc, ip, el)
 
 end function constitutive_postResults
 
+
+!--------------------------------------------------------------------------------------------------
+!> @brief writes constitutive results to HDF5 output file
+!--------------------------------------------------------------------------------------------------
+subroutine constitutive_results()
+ use material, only: &
+   PLASTICITY_ISOTROPIC_ID, &
+   PLASTICITY_PHENOPOWERLAW_ID, &
+   PLASTICITY_KINEHARDENING_ID, &
+   PLASTICITY_DISLOTWIN_ID, &
+   PLASTICITY_DISLOUCLA_ID, &
+   PLASTICITY_NONLOCAL_ID
+#if defined(PETSc) || defined(DAMASKHDF5)
+ use results
+ use HDF5_utilities
+ use config, only: &
+   config_name_phase => phase_name                                                                  ! anticipate logical name
+   
+ use material, only: &
+   phase_plasticityInstance, &
+   material_phase_plasticity_type => phase_plasticity
+ use plastic_phenopowerlaw, only: &
+   plastic_phenopowerlaw_results
+ 
+ implicit none
+ integer(pInt) :: p  
+ call HDF5_closeGroup(results_addGroup('current/phase'))                                              
+ do p=1,size(config_name_phase)                                                                           
+   call HDF5_closeGroup(results_addGroup('current/phase/'//trim(config_name_phase(p))))
+   if (material_phase_plasticity_type(p) == PLASTICITY_PHENOPOWERLAW_ID) then
+     call plastic_phenopowerlaw_results(phase_plasticityInstance(p),'current/phase/'//trim(config_name_phase(p)))
+   endif
+ enddo      
+
+#endif
+
+
+end subroutine constitutive_results
+
 end module constitutive

src/plastic_disloUCLA.f90

@@ -63,8 +63,7 @@ module plastic_disloUCLA
      interaction_SlipSlip, &                                                                        !< slip resistance from slip activity
      forestProjectionEdge
    real(pReal),                 allocatable, dimension(:,:,:) :: &
-     Schmid_slip, &
+     Schmid, &
-     Schmid_twin, &
      nonSchmid_pos, &
      nonSchmid_neg
    integer(pInt) :: &
@@ -77,13 +76,11 @@ module plastic_disloUCLA
      dipoleformation
  end type                                                                                           !< container type for internal constitutive parameters
 
- type(tParameters), dimension(:), allocatable, private :: param                                     !< containers of constitutive parameters (len Ninstance)
-
  type, private :: tDisloUCLAState
    real(pReal), pointer, dimension(:,:) :: &
      rhoEdge, &
      rhoEdgeDip, &
-     accshear_slip
+     accshear
  end type
 
  type, private :: tDisloUCLAdependentState
@@ -93,6 +90,8 @@ module plastic_disloUCLA
      threshold_stress
  end type tDisloUCLAdependentState
 
+
+ type(tParameters), dimension(:), allocatable, private :: param                                     !< containers of constitutive parameters (len Ninstance)
  type(tDisloUCLAState ), allocatable, dimension(:), private :: &
    dotState, &
    state
@@ -110,6 +109,7 @@ module plastic_disloUCLA
 
 contains
 
+
 !--------------------------------------------------------------------------------------------------
 !> @brief module initialization
 !> @details reads in material parameters, allocates arrays, and does sanity checks
@@ -152,7 +152,7 @@ subroutine plastic_disloUCLA_init()
    f,j,k,o, &
    Ninstance, &
    p, i, &
-   NipcMyPhase, outputSize, &
+   NipcMyPhase, &
    sizeState, sizeDotState, &
    startIndex, endIndex
 
@@ -213,7 +213,7 @@ subroutine plastic_disloUCLA_init()
    prm%Nslip      = config%getInts('nslip',defaultVal=emptyIntArray)
    prm%totalNslip = sum(prm%Nslip)
    slipActive: if (prm%totalNslip > 0_pInt) then
-     prm%Schmid_slip          = lattice_SchmidMatrix_slip(prm%Nslip,structure(1:3),&
+     prm%Schmid          = lattice_SchmidMatrix_slip(prm%Nslip,structure(1:3),&
                                                      config%getFloat('c/a',defaultVal=0.0_pReal))
      if(structure=='bcc') then
        prm%nonSchmidCoeff     = config%getFloats('nonschmid_coefficients',&
@@ -221,8 +221,8 @@ subroutine plastic_disloUCLA_init()
        prm%nonSchmid_pos      = lattice_nonSchmidMatrix(prm%Nslip,prm%nonSchmidCoeff,+1_pInt)
        prm%nonSchmid_neg      = lattice_nonSchmidMatrix(prm%Nslip,prm%nonSchmidCoeff,-1_pInt)
      else
-       prm%nonSchmid_pos      = prm%Schmid_slip
+       prm%nonSchmid_pos      = prm%Schmid
-       prm%nonSchmid_neg      = prm%Schmid_slip
+       prm%nonSchmid_neg      = prm%Schmid
      endif
      prm%interaction_SlipSlip = lattice_interaction_SlipSlip(prm%Nslip, &
                                                              config%getFloats('interaction_slipslip'), &
@@ -272,17 +272,17 @@ subroutine plastic_disloUCLA_init()
      prm%tau0 = prm%tau_peierls + prm%SolidSolutionStrength
 
      ! sanity checks
+     if (    prm%D0             <= 0.0_pReal)  extmsg = trim(extmsg)//' d0'
+     if (    prm%Qsd            <= 0.0_pReal)  extmsg = trim(extmsg)//' qsd'
      if (any(prm%rho0           <  0.0_pReal)) extmsg = trim(extmsg)//' rhoedge0'
      if (any(prm%rhoDip0        <  0.0_pReal)) extmsg = trim(extmsg)//' rhoedgedip0'
      if (any(prm%v0             <  0.0_pReal)) extmsg = trim(extmsg)//' v0'
      if (any(prm%burgers        <= 0.0_pReal)) extmsg = trim(extmsg)//' slipburgers'
      if (any(prm%H0kp           <= 0.0_pReal)) extmsg = trim(extmsg)//' qedge'
      if (any(prm%tau_peierls    <  0.0_pReal)) extmsg = trim(extmsg)//' tau_peierls'
-     if (    prm%D0           <= 0.0_pReal)  extmsg = trim(extmsg)//' d0'
+     if (any(prm%minDipDistance <= 0.0_pReal)) extmsg = trim(extmsg)//' cedgedipmindistance or slipburgers'
-     if (    prm%Qsd          <= 0.0_pReal)  extmsg = trim(extmsg)//' qsd'
+     if (any(prm%atomicVolume   <= 0.0_pReal)) extmsg = trim(extmsg)//' catomicvolume or slipburgers'
 
-     !if (plastic_disloUCLA_CAtomicVolume(instance) <= 0.0_pReal) &
-     !  call IO_error(211_pInt,el=instance,ext_msg='cAtomicVolume ('//PLASTICITY_DISLOUCLA_label//')')
    else slipActive
      allocate(prm%rho0(0))
      allocate(prm%rhoDip0(0))
@@ -299,7 +299,6 @@ subroutine plastic_disloUCLA_init()
    allocate(prm%outputID(0))
    do i=1_pInt, size(outputs)
      outputID = undefined_ID
-     outputSize = prm%totalNslip
      select case(trim(outputs(i)))
 
        case ('edge_density')
@@ -321,7 +320,7 @@ subroutine plastic_disloUCLA_init()
 
      if (outputID /= undefined_ID) then
        plastic_disloUCLA_output(i,phase_plasticityInstance(p)) = outputs(i)
-       plastic_disloUCLA_sizePostResult(i,phase_plasticityInstance(p)) = outputSize
+       plastic_disloUCLA_sizePostResult(i,phase_plasticityInstance(p)) = prm%totalNslip
        prm%outputID = [prm%outputID, outputID]
      endif
 
@@ -373,9 +372,9 @@ subroutine plastic_disloUCLA_init()
 
    startIndex = endIndex + 1_pInt
    endIndex   = endIndex + prm%totalNslip
-   stt%accshear_slip=>plasticState(p)%state(startIndex:endIndex,:)
+   stt%accshear=>plasticState(p)%state(startIndex:endIndex,:)
-   dot%accshear_slip=>plasticState(p)%dotState(startIndex:endIndex,:)
+   dot%accshear=>plasticState(p)%dotState(startIndex:endIndex,:)
-   plasticState(p)%aTolState(startIndex:endIndex) = 1e6_pReal
+   plasticState(p)%aTolState(startIndex:endIndex) = 1e6_pReal                  !ToDo: better make optional parameter
    ! global alias
    plasticState(p)%slipRate        => plasticState(p)%dotState(startIndex:endIndex,:)
    plasticState(p)%accumulatedSlip => plasticState(p)%state(startIndex:endIndex,:)
@@ -393,36 +392,6 @@ subroutine plastic_disloUCLA_init()
 end subroutine plastic_disloUCLA_init
 
 
-!--------------------------------------------------------------------------------------------------
-!> @brief calculates derived quantities from state
-!--------------------------------------------------------------------------------------------------
-subroutine plastic_disloUCLA_dependentState(instance,of)
-
- implicit none
- integer(pInt), intent(in) :: instance, of
-
- integer(pInt) :: &
-   i
-
- associate(prm => param(instance), stt => state(instance),dst => dependentState(instance))
-
- forall (i = 1_pInt:prm%totalNslip)
-   dst%dislocationSpacing(i,of) = sqrt(dot_product(stt%rhoEdge(:,of)+stt%rhoEdgeDip(:,of), &
-                                       prm%forestProjectionEdge(:,i)))
-   dst%threshold_stress(i,of) = prm%mu*prm%burgers(i) &
-                              * sqrt(dot_product(stt%rhoEdge(:,of)+stt%rhoEdgeDip(:,of), &
-                                                 prm%interaction_SlipSlip(i,:)))
- end forall
-
- dst%mfp(:,of) = prm%grainSize/(1.0_pReal+prm%grainSize*dst%dislocationSpacing(:,of)/prm%Clambda)
- dst%dislocationSpacing(:,of) = dst%mfp(:,of)                                                       ! ToDo: Hack to recover wrong behavior for the moment
- 
- end associate
-
-
-end subroutine plastic_disloUCLA_dependentState
-
-
 !--------------------------------------------------------------------------------------------------
 !> @brief calculates plastic velocity gradient and its tangent
 !--------------------------------------------------------------------------------------------------
@@ -445,22 +414,22 @@ pure subroutine plastic_disloUCLA_LpAndItsTangent(Lp,dLp_dMp,Mp,Temperature,inst
  integer(pInt) :: &
    i,k,l,m,n
  real(pReal), dimension(param(instance)%totalNslip) :: &
-   dgdot_dtauslip_pos,dgdot_dtauslip_neg, &
+   gdot_pos,gdot_neg, &
-   gdot_slip_pos,gdot_slip_neg
+   dgdot_dtau_pos,dgdot_dtau_neg
 
  Lp = 0.0_pReal
  dLp_dMp = 0.0_pReal
 
  associate(prm => param(instance))
 
- call kinetics(Mp,Temperature,instance,of,gdot_slip_pos,gdot_slip_neg,dgdot_dtauslip_pos,dgdot_dtauslip_neg)
+ call kinetics(Mp,Temperature,instance,of,gdot_pos,gdot_neg,dgdot_dtau_pos,dgdot_dtau_neg)
- slipSystems: do i = 1_pInt, prm%totalNslip
+ do i = 1_pInt, prm%totalNslip
-   Lp = Lp + (gdot_slip_pos(i)+gdot_slip_neg(i))*prm%Schmid_slip(1:3,1:3,i)
+   Lp = Lp + (gdot_pos(i)+gdot_neg(i))*prm%Schmid(1:3,1:3,i)
    forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) &
      dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) &
-                      + dgdot_dtauslip_pos(i) * prm%Schmid_slip(k,l,i) * prm%nonSchmid_pos(m,n,i) &
+                      + dgdot_dtau_pos(i) * prm%Schmid(k,l,i) * prm%nonSchmid_pos(m,n,i) &
-                      + dgdot_dtauslip_neg(i) * prm%Schmid_slip(k,l,i) * prm%nonSchmid_neg(m,n,i)
+                      + dgdot_dtau_neg(i) * prm%Schmid(k,l,i) * prm%nonSchmid_neg(m,n,i)
- enddo slipSystems
+ enddo
 
  end associate
 
@@ -484,34 +453,35 @@ subroutine plastic_disloUCLA_dotState(Mp,Temperature,instance,of)
  real(pReal),                  intent(in) :: &
    temperature                                                                                      !< temperature
  integer(pInt),                intent(in) :: &
-   instance, of
+   instance, &
+   of
 
  real(pReal) :: &
    VacancyDiffusion
  real(pReal), dimension(param(instance)%totalNslip) :: &
-   gdot_slip_pos, gdot_slip_neg,&
+   gdot_pos, gdot_neg,&
-   tau_slip_pos,&
+   tau_pos,&
-   tau_slip_neg, &
+   tau_neg, &
    DotRhoDipFormation, ClimbVelocity, EdgeDipDistance, &
    DotRhoEdgeDipClimb
 
  associate(prm => param(instance), stt => state(instance),dot => dotState(instance), dst => dependentState(instance))
 
  call kinetics(Mp,Temperature,instance,of,&
-               gdot_slip_pos,gdot_slip_neg, &
+               gdot_pos,gdot_neg, &
-               tau_slip_pos1 = tau_slip_pos,tau_slip_neg1 = tau_slip_neg)
+               tau_pos1 = tau_pos,tau_neg1 = tau_neg)
 
- dot%accshear_slip(:,of) = (gdot_slip_pos+gdot_slip_neg)                                            ! ToDo: needs to be abs
+ dot%accshear(:,of) = (gdot_pos+gdot_neg)                                                      ! ToDo: needs to be abs
  VacancyDiffusion = prm%D0*exp(-prm%Qsd/(kB*Temperature))
 
- where(dEq0(tau_slip_pos))                                                                          ! ToDo: use avg of pos and neg
+ where(dEq0(tau_pos))                                                                          ! ToDo: use avg of pos and neg
    DotRhoDipFormation = 0.0_pReal
    DotRhoEdgeDipClimb = 0.0_pReal
  else where
-   EdgeDipDistance = math_clip((3.0_pReal*prm%mu*prm%burgers)/(16.0_pReal*PI*abs(tau_slip_pos)), &
+   EdgeDipDistance = math_clip((3.0_pReal*prm%mu*prm%burgers)/(16.0_pReal*PI*abs(tau_pos)), &
                                prm%minDipDistance, &                                                ! lower limit
                                dst%mfp(:,of))                                                       ! upper limit
-   DotRhoDipFormation = merge(((2.0_pReal*EdgeDipDistance)/prm%burgers)* stt%rhoEdge(:,of)*abs(dot%accshear_slip(:,of)), & ! ToDo: ignore region of spontaneous annihilation
+   DotRhoDipFormation = merge(((2.0_pReal*EdgeDipDistance)/prm%burgers)* stt%rhoEdge(:,of)*abs(dot%accshear(:,of)), & ! ToDo: ignore region of spontaneous annihilation
                               0.0_pReal, &
                               prm%dipoleformation)
    ClimbVelocity = (3.0_pReal*prm%mu*VacancyDiffusion*prm%atomicVolume/(2.0_pReal*pi*kB*Temperature)) &
@@ -519,11 +489,11 @@ subroutine plastic_disloUCLA_dotState(Mp,Temperature,instance,of)
    DotRhoEdgeDipClimb = (4.0_pReal*ClimbVelocity*stt%rhoEdgeDip(:,of))/(EdgeDipDistance-prm%minDipDistance) ! ToDo: Discuss with Franz: Stress dependency?
  end where
 
- dot%rhoEdge(:,of) = abs(dot%accshear_slip(:,of))/(prm%burgers*dst%mfp(:,of)) &                     ! multiplication
+ dot%rhoEdge(:,of) = abs(dot%accshear(:,of))/(prm%burgers*dst%mfp(:,of)) &                     ! multiplication
                    - DotRhoDipFormation &
-                   - (2.0_pReal*prm%minDipDistance)/prm%burgers*stt%rhoEdge(:,of)*abs(dot%accshear_slip(:,of)) !* Spontaneous annihilation of 2 single edge dislocations
+                   - (2.0_pReal*prm%minDipDistance)/prm%burgers*stt%rhoEdge(:,of)*abs(dot%accshear(:,of)) !* Spontaneous annihilation of 2 single edge dislocations
  dot%rhoEdgeDip(:,of) = DotRhoDipFormation &
-                      - (2.0_pReal*prm%minDipDistance)/prm%burgers*stt%rhoEdgeDip(:,of)*abs(dot%accshear_slip(:,of)) & !* Spontaneous annihilation of a single edge dislocation with a dipole constituent
+                      - (2.0_pReal*prm%minDipDistance)/prm%burgers*stt%rhoEdgeDip(:,of)*abs(dot%accshear(:,of)) & !* Spontaneous annihilation of a single edge dislocation with a dipole constituent
                       - DotRhoEdgeDipClimb
 
  end associate
@@ -531,6 +501,37 @@ subroutine plastic_disloUCLA_dotState(Mp,Temperature,instance,of)
 end subroutine plastic_disloUCLA_dotState
 
 
+!--------------------------------------------------------------------------------------------------
+!> @brief calculates derived quantities from state
+!--------------------------------------------------------------------------------------------------
+subroutine plastic_disloUCLA_dependentState(instance,of)
+
+ implicit none
+ integer(pInt),                intent(in) :: &
+   instance, &
+   of
+
+ integer(pInt) :: &
+   i
+
+ associate(prm => param(instance), stt => state(instance),dst => dependentState(instance))
+
+ forall (i = 1_pInt:prm%totalNslip)
+   dst%dislocationSpacing(i,of) = sqrt(dot_product(stt%rhoEdge(:,of)+stt%rhoEdgeDip(:,of), &
+                                       prm%forestProjectionEdge(:,i)))
+   dst%threshold_stress(i,of) = prm%mu*prm%burgers(i) &
+                              * sqrt(dot_product(stt%rhoEdge(:,of)+stt%rhoEdgeDip(:,of), &
+                                                 prm%interaction_SlipSlip(i,:)))
+ end forall
+
+ dst%mfp(:,of) = prm%grainSize/(1.0_pReal+prm%grainSize*dst%dislocationSpacing(:,of)/prm%Clambda)
+ dst%dislocationSpacing(:,of) = dst%mfp(:,of)                                                       ! ToDo: Hack to recover wrong behavior for the moment
+
+ end associate
+
+end subroutine plastic_disloUCLA_dependentState
+
+
 !--------------------------------------------------------------------------------------------------
 !> @brief return array of constitutive results
 !--------------------------------------------------------------------------------------------------
@@ -556,7 +557,7 @@ function plastic_disloUCLA_postResults(Mp,Temperature,instance,of) result(postRe
  integer(pInt) :: &
    o,c,i
  real(pReal), dimension(param(instance)%totalNslip) :: &
-   gdot_slip_pos,gdot_slip_neg
+   gdot_pos,gdot_neg
 
  c = 0_pInt
 
@@ -570,10 +571,10 @@ function plastic_disloUCLA_postResults(Mp,Temperature,instance,of) result(postRe
      case (rhoDip_ID)
        postResults(c+1_pInt:c+prm%totalNslip) = stt%rhoEdgeDip(1_pInt:prm%totalNslip,of)
      case (shearrate_ID)
-       call kinetics(Mp,Temperature,instance,of,gdot_slip_pos,gdot_slip_neg)
+       call kinetics(Mp,Temperature,instance,of,gdot_pos,gdot_neg)
-       postResults(c+1:c+prm%totalNslip) = gdot_slip_pos + gdot_slip_neg
+       postResults(c+1:c+prm%totalNslip) = gdot_pos + gdot_neg
      case (accumulatedshear_ID)
-       postResults(c+1_pInt:c+prm%totalNslip) = stt%accshear_slip(1_pInt:prm%totalNslip, of)
+       postResults(c+1_pInt:c+prm%totalNslip) = stt%accshear(1_pInt:prm%totalNslip, of)
      case (mfp_ID)
        postResults(c+1_pInt:c+prm%totalNslip) = dst%mfp(1_pInt:prm%totalNslip, of)
      case (thresholdstress_ID)
@@ -608,7 +609,7 @@ end function plastic_disloUCLA_postResults
 ! have the optional arguments at the end
 !--------------------------------------------------------------------------------------------------
 pure subroutine kinetics(Mp,Temperature,instance,of, &
-                 gdot_slip_pos,gdot_slip_neg,dgdot_dtauslip_pos,dgdot_dtauslip_neg,tau_slip_pos1,tau_slip_neg1)
+                 gdot_pos,gdot_neg,dgdot_dtau_pos,dgdot_dtau_neg,tau_pos1,tau_neg1)
  use prec, only: &
    tol_math_check, &
    dEq, dNeq0
@@ -626,119 +627,119 @@ pure subroutine kinetics(Mp,Temperature,instance,of, &
    of
 
  real(pReal),                  intent(out), dimension(param(instance)%totalNslip) :: &
-   gdot_slip_pos, &
+   gdot_pos, &
-   gdot_slip_neg
+   gdot_neg
  real(pReal),                  intent(out), optional, dimension(param(instance)%totalNslip) :: &
-   dgdot_dtauslip_pos, &
+   dgdot_dtau_pos, &
-   dgdot_dtauslip_neg, &
+   dgdot_dtau_neg, &
-   tau_slip_pos1, &
+   tau_pos1, &
-   tau_slip_neg1
+   tau_neg1
  real(pReal), dimension(param(instance)%totalNslip) :: &
    StressRatio, &
    StressRatio_p,StressRatio_pminus1, &
-   dvel_slip, vel_slip, &
+   dvel, vel, &
-   tau_slip_pos,tau_slip_neg, &
+   tau_pos,tau_neg, &
    needsGoodName                                                                                    ! ToDo: @Karo: any idea?
  integer(pInt) :: j
 
  associate(prm => param(instance), stt => state(instance), dst => dependentState(instance))
 
  do j = 1_pInt, prm%totalNslip
-   tau_slip_pos(j) = math_mul33xx33(Mp,prm%nonSchmid_pos(1:3,1:3,j))
+   tau_pos(j) = math_mul33xx33(Mp,prm%nonSchmid_pos(1:3,1:3,j))
-   tau_slip_neg(j) = math_mul33xx33(Mp,prm%nonSchmid_neg(1:3,1:3,j))
+   tau_neg(j) = math_mul33xx33(Mp,prm%nonSchmid_neg(1:3,1:3,j))
  enddo
 
 
- if (present(tau_slip_pos1)) tau_slip_pos1 = tau_slip_pos
+ if (present(tau_pos1)) tau_pos1 = tau_pos
- if (present(tau_slip_neg1)) tau_slip_neg1 = tau_slip_neg
+ if (present(tau_neg1)) tau_neg1 = tau_neg
 
  associate(BoltzmannRatio  => prm%H0kp/(kB*Temperature), &
            DotGamma0       => stt%rhoEdge(:,of)*prm%burgers*prm%v0, &
            effectiveLength => dst%mfp(:,of) - prm%w)
 
- significantPositiveTau: where(abs(tau_slip_pos)-dst%threshold_stress(:,of) > tol_math_check)
+ significantPositiveTau: where(abs(tau_pos)-dst%threshold_stress(:,of) > tol_math_check)
-   StressRatio = (abs(tau_slip_pos)-dst%threshold_stress(:,of))/prm%tau0
+   StressRatio = (abs(tau_pos)-dst%threshold_stress(:,of))/prm%tau0
    StressRatio_p       = StressRatio** prm%p
    StressRatio_pminus1 = StressRatio**(prm%p-1.0_pReal)
    needsGoodName       = exp(-BoltzmannRatio*(1-StressRatio_p) ** prm%q)
 
-   vel_slip = 2.0_pReal*prm%burgers * prm%kink_height * prm%omega  &
+   vel = 2.0_pReal*prm%burgers * prm%kink_height * prm%omega  &
-            * effectiveLength * tau_slip_pos * needsGoodName  &
+            * effectiveLength * tau_pos * needsGoodName  &
-          / (  2.0_pReal*(prm%burgers**2.0_pReal)*tau_slip_pos &
+          / (  2.0_pReal*(prm%burgers**2.0_pReal)*tau_pos &
              + prm%omega * prm%B * effectiveLength**2.0_pReal* needsGoodName  &
             )
 
-   gdot_slip_pos = DotGamma0 * sign(vel_slip,tau_slip_pos) * 0.5_pReal
+   gdot_pos = DotGamma0 * sign(vel,tau_pos) * 0.5_pReal
  else where significantPositiveTau
-   gdot_slip_pos = 0.0_pReal
+   gdot_pos = 0.0_pReal
  end where significantPositiveTau
 
- if (present(dgdot_dtauslip_pos)) then
+ if (present(dgdot_dtau_pos)) then
- significantPositiveTau2: where(abs(tau_slip_pos)-dst%threshold_stress(:,of) > tol_math_check)
+ significantPositiveTau2: where(abs(tau_pos)-dst%threshold_stress(:,of) > tol_math_check)
-   dvel_slip = 2.0_pReal*prm%burgers * prm%kink_height * prm%omega* effectiveLength &
+   dvel = 2.0_pReal*prm%burgers * prm%kink_height * prm%omega* effectiveLength &
              * ( &
-                 (needsGoodName + tau_slip_pos * abs(needsGoodName)*BoltzmannRatio*prm%p &
+                 (needsGoodName + tau_pos * abs(needsGoodName)*BoltzmannRatio*prm%p &
                                 * prm%q/prm%tau0 &
                                 * StressRatio_pminus1*(1-StressRatio_p)**(prm%q-1.0_pReal)   &
                  ) &
-               * (  2.0_pReal*(prm%burgers**2.0_pReal)*tau_slip_pos &
+               * (  2.0_pReal*(prm%burgers**2.0_pReal)*tau_pos &
                   +  prm%omega * prm%B* effectiveLength **2.0_pReal* needsGoodName &
                  ) &
-               - tau_slip_pos * needsGoodName * (2.0_pReal*prm%burgers**2.0_pReal &
+               - tau_pos * needsGoodName * (2.0_pReal*prm%burgers**2.0_pReal &
                                                  +  prm%omega * prm%B *effectiveLength **2.0_pReal&
                                                   * (abs(needsGoodName)*BoltzmannRatio*prm%p *prm%q/prm%tau0 &
                                                 *StressRatio_pminus1*(1-StressRatio_p)**(prm%q-1.0_pReal)  )&
                                                ) &
               )  &
-            /(2.0_pReal*prm%burgers**2.0_pReal*tau_slip_pos &
+            /(2.0_pReal*prm%burgers**2.0_pReal*tau_pos &
               + prm%omega * prm%B* effectiveLength**2.0_pReal* needsGoodName )**2.0_pReal
 
-   dgdot_dtauslip_pos = DotGamma0 * dvel_slip* 0.5_pReal
+   dgdot_dtau_pos = DotGamma0 * dvel* 0.5_pReal
  else where significantPositiveTau2
-   dgdot_dtauslip_pos = 0.0_pReal
+   dgdot_dtau_pos = 0.0_pReal
  end where significantPositiveTau2
  endif
 
- significantNegativeTau: where(abs(tau_slip_neg)-dst%threshold_stress(:,of) > tol_math_check)
+ significantNegativeTau: where(abs(tau_neg)-dst%threshold_stress(:,of) > tol_math_check)
-   StressRatio = (abs(tau_slip_neg)-dst%threshold_stress(:,of))/prm%tau0
+   StressRatio = (abs(tau_neg)-dst%threshold_stress(:,of))/prm%tau0
    StressRatio_p       = StressRatio** prm%p
    StressRatio_pminus1 = StressRatio**(prm%p-1.0_pReal)
    needsGoodName       = exp(-BoltzmannRatio*(1-StressRatio_p) ** prm%q)
 
-   vel_slip = 2.0_pReal*prm%burgers * prm%kink_height * prm%omega  &
+   vel = 2.0_pReal*prm%burgers * prm%kink_height * prm%omega  &
-            * effectiveLength * tau_slip_neg * needsGoodName  &
+            * effectiveLength * tau_neg * needsGoodName  &
-          / (  2.0_pReal*(prm%burgers**2.0_pReal)*tau_slip_neg &
+          / (  2.0_pReal*(prm%burgers**2.0_pReal)*tau_neg &
              + prm%omega * prm%B * effectiveLength**2.0_pReal* needsGoodName  &
             )
 
-   gdot_slip_neg = DotGamma0 * sign(vel_slip,tau_slip_neg) * 0.5_pReal
+   gdot_neg = DotGamma0 * sign(vel,tau_neg) * 0.5_pReal
  else where significantNegativeTau
-   gdot_slip_neg = 0.0_pReal
+   gdot_neg = 0.0_pReal
  end where significantNegativeTau
 
- if (present(dgdot_dtauslip_neg)) then
+ if (present(dgdot_dtau_neg)) then
- significantNegativeTau2: where(abs(tau_slip_neg)-dst%threshold_stress(:,of) > tol_math_check)
+ significantNegativeTau2: where(abs(tau_neg)-dst%threshold_stress(:,of) > tol_math_check)
-   dvel_slip = 2.0_pReal*prm%burgers * prm%kink_height * prm%omega* effectiveLength &
+   dvel = 2.0_pReal*prm%burgers * prm%kink_height * prm%omega* effectiveLength &
              * ( &
-                 (needsGoodName + tau_slip_neg * abs(needsGoodName)*BoltzmannRatio*prm%p &
+                 (needsGoodName + tau_neg * abs(needsGoodName)*BoltzmannRatio*prm%p &
                                 * prm%q/prm%tau0 &
                                 * StressRatio_pminus1*(1-StressRatio_p)**(prm%q-1.0_pReal)   &
                  ) &
-               * (  2.0_pReal*(prm%burgers**2.0_pReal)*tau_slip_neg &
+               * (  2.0_pReal*(prm%burgers**2.0_pReal)*tau_neg &
                   +  prm%omega * prm%B* effectiveLength **2.0_pReal* needsGoodName &
                  ) &
-               - tau_slip_neg * needsGoodName * (2.0_pReal*prm%burgers**2.0_pReal &
+               - tau_neg * needsGoodName * (2.0_pReal*prm%burgers**2.0_pReal &
                                                  +  prm%omega * prm%B *effectiveLength **2.0_pReal&
                                                   * (abs(needsGoodName)*BoltzmannRatio*prm%p *prm%q/prm%tau0 &
                                                 *StressRatio_pminus1*(1-StressRatio_p)**(prm%q-1.0_pReal)  )&
                                                ) &
               )  &
-            /(2.0_pReal*prm%burgers**2.0_pReal*tau_slip_neg &
+            /(2.0_pReal*prm%burgers**2.0_pReal*tau_neg &
               + prm%omega * prm%B* effectiveLength**2.0_pReal* needsGoodName )**2.0_pReal
 
-   dgdot_dtauslip_neg = DotGamma0 * dvel_slip * 0.5_pReal
+   dgdot_dtau_neg = DotGamma0 * dvel * 0.5_pReal
  else where significantNegativeTau2
-   dgdot_dtauslip_neg = 0.0_pReal
+   dgdot_dtau_neg = 0.0_pReal
  end where significantNegativeTau2
  end if
  end associate

src/plastic_isotropic.f90

@@ -50,14 +50,14 @@ module plastic_isotropic
      dilatation
  end type
 
- type(tParameters), dimension(:), allocatable, private :: param                                     !< containers of constitutive parameters (len Ninstance)
-
  type, private :: tIsotropicState
    real(pReal), pointer, dimension(:) :: &
      flowstress, &
      accumulatedShear
  end type
 
+
+ type(tParameters), dimension(:), allocatable, private :: param                                     !< containers of constitutive parameters (len Ninstance)
  type(tIsotropicState), allocatable, dimension(:), private :: &
    dotState, &
    state
@@ -231,17 +231,17 @@ subroutine plastic_isotropic_init()
 
 !--------------------------------------------------------------------------------------------------
 ! locally defined state aliases and initialization of state0 and aTolState
-   stt%flowstress  => plasticState(p)%state    (1,1:NipcMyPhase)
+   stt%flowstress  => plasticState(p)%state   (1,:)
    stt%flowstress  = prm%tau0
-   dot%flowstress  => plasticState(p)%dotState (1,1:NipcMyPhase)
+   dot%flowstress  => plasticState(p)%dotState(1,:)
    plasticState(p)%aTolState(1) = prm%aTolFlowstress
 
-   stt%accumulatedShear  => plasticState(p)%state    (2,1:NipcMyPhase)
+   stt%accumulatedShear  => plasticState(p)%state   (2,:)
-   dot%accumulatedShear  => plasticState(p)%dotState (2,1:NipcMyPhase)
+   dot%accumulatedShear  => plasticState(p)%dotState(2,:)
    plasticState(p)%aTolState(2) = prm%aTolShear
    ! global alias
-   plasticState(p)%slipRate           => plasticState(p)%dotState(2:2,1:NipcMyPhase)
+   plasticState(p)%slipRate        => plasticState(p)%dotState(2:2,:)
-   plasticState(p)%accumulatedSlip    => plasticState(p)%state   (2:2,1:NipcMyPhase)
+   plasticState(p)%accumulatedSlip => plasticState(p)%state   (2:2,:)
 
    plasticState(p)%state0 = plasticState(p)%state                                                   ! ToDo: this could be done centrally
 

src/plastic_none.f90

@@ -60,7 +60,6 @@ subroutine plastic_none_init
 !--------------------------------------------------------------------------------------------------
 ! allocate state arrays
    NipcMyPhase = count(material_phase == p)
-
    call material_allocatePlasticState(p,NipcMyPhase,0_pInt,0_pInt,0_pInt, &
                                       0_pInt,0_pInt,0_pInt)
    plasticState(p)%sizePostResults = 0_pInt

src/plastic_phenopowerlaw.f90

@@ -74,8 +74,6 @@ module plastic_phenopowerlaw
      outputID                                                                                       !< ID of each post result output
  end type                                                                                           !< container type for internal constitutive parameters
 
- type(tParameters), dimension(:), allocatable, private :: param                                     !< containers of constitutive parameters (len Ninstance)
-
  type, private :: tPhenopowerlawState
    real(pReal), pointer, dimension(:,:) :: &
      xi_slip, &
@@ -84,6 +82,8 @@ module plastic_phenopowerlaw
      gamma_twin
  end type
 
+
+ type(tParameters), dimension(:), allocatable, private :: param                                     !< containers of constitutive parameters (len Ninstance)
  type(tPhenopowerlawState), allocatable, dimension(:), private :: &
    dotState, &
    state
@@ -92,13 +92,15 @@ module plastic_phenopowerlaw
    plastic_phenopowerlaw_init, &
    plastic_phenopowerlaw_LpAndItsTangent, &
    plastic_phenopowerlaw_dotState, &
-   plastic_phenopowerlaw_postResults
+   plastic_phenopowerlaw_postResults, &
+   plastic_phenopowerlaw_results
  private :: &
    kinetics_slip, &
    kinetics_twin
 
 contains
 
+
 !--------------------------------------------------------------------------------------------------
 !> @brief module initialization
 !> @details reads in material parameters, allocates arrays, and does sanity checks
@@ -341,8 +343,8 @@ subroutine plastic_phenopowerlaw_init
 !--------------------------------------------------------------------------------------------------
 ! allocate state arrays
    NipcMyPhase = count(material_phase == p)
-   sizeDotState = size(['tau_slip  ','gamma_slip']) * prm%TotalNslip &
+   sizeDotState = size(['tau_slip  ','gamma_slip']) * prm%totalNslip &
-                + size(['tau_twin  ','gamma_twin']) * prm%TotalNtwin
+                + size(['tau_twin  ','gamma_twin']) * prm%totalNtwin
    sizeState = sizeDotState
 
    call material_allocatePlasticState(p,NipcMyPhase,sizeState,sizeDotState,0_pInt, &
@@ -391,8 +393,8 @@ end subroutine plastic_phenopowerlaw_init
 
 !--------------------------------------------------------------------------------------------------
 !> @brief calculates plastic velocity gradient and its tangent
-!> @details asumme that deformation by dislocation glide affects twinned and untwinned volume 
+!> @details asummes that deformation by dislocation glide affects twinned and untwinned volume
-!  equally (Taylor assumption). Twinning happens only in untwinned volume (
+!  equally (Taylor assumption). Twinning happens only in untwinned volume
 !--------------------------------------------------------------------------------------------------
 pure subroutine plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,of)
 
@@ -411,8 +413,8 @@ pure subroutine plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,of)
  integer(pInt) :: &
    i,k,l,m,n
  real(pReal), dimension(param(instance)%totalNslip) :: &
-   dgdot_dtauslip_pos,dgdot_dtauslip_neg, &
+   gdot_slip_pos,gdot_slip_neg, &
-   gdot_slip_pos,gdot_slip_neg
+   dgdot_dtauslip_pos,dgdot_dtauslip_neg
  real(pReal), dimension(param(instance)%totalNtwin) :: &
    gdot_twin,dgdot_dtautwin
 
@@ -575,6 +577,35 @@ function plastic_phenopowerlaw_postResults(Mp,instance,of) result(postResults)
 end function plastic_phenopowerlaw_postResults
 
 
+!--------------------------------------------------------------------------------------------------
+!> @brief writes results to HDF5 output file
+!--------------------------------------------------------------------------------------------------
+subroutine plastic_phenopowerlaw_results(instance,group)
+#if defined(PETSc) || defined(DAMASKHDF5)
+ use results
+
+ implicit none
+ integer(pInt), intent(in) :: instance
+ character(len=*) :: group
+ integer(pInt) :: o
+
+ associate(prm => param(instance), stt => state(instance))
+ outputsLoop: do o = 1_pInt,size(prm%outputID)
+   select case(prm%outputID(o))
+     case (resistance_slip_ID)
+       call results_writeVectorDataset(group,stt%xi_slip,'xi_slip','Pa')
+     case (accumulatedshear_slip_ID)
+       call results_writeVectorDataset(group,stt%gamma_slip,'gamma_slip','-')
+   end select
+ enddo outputsLoop
+ end associate
+#else
+ integer(pInt), intent(in) :: instance
+ character(len=*) :: group
+#endif
+end subroutine plastic_phenopowerlaw_results
+
+
 !--------------------------------------------------------------------------------------------------
 !> @brief Shear rates on slip systems and their derivatives with respect to resolved stress
 !> @details Derivatives are calculated only optionally.
@@ -601,6 +632,7 @@ pure subroutine kinetics_slip(Mp,instance,of, &
  real(pReal),                  intent(out), optional, dimension(param(instance)%totalNslip) :: &
    dgdot_dtau_slip_pos, &
    dgdot_dtau_slip_neg
+
  real(pReal), dimension(param(instance)%totalNslip) :: &
    tau_slip_pos, &
    tau_slip_neg
@@ -625,7 +657,7 @@ pure subroutine kinetics_slip(Mp,instance,of, &
  end where
 
  where(dNeq0(tau_slip_neg))
-   gdot_slip_neg = 0.5_pReal*prm%gdot0_slip &
+   gdot_slip_neg = prm%gdot0_slip * 0.5_pReal &                                                     ! only used if non-Schmid active, always 1/2
                  * sign(abs(tau_slip_neg/stt%xi_slip(:,of))**prm%n_slip,  tau_slip_neg)
  else where
    gdot_slip_neg = 0.0_pReal

src/quit.f90

@@ -23,6 +23,8 @@ subroutine quit(stop_id)
  integer(pInt) :: error = 0_pInt
  PetscErrorCode :: ierr = 0
 
+ call h5open_f(hdferr)
+ if (hdferr /= 0) write(6,'(a,i5)') ' Error in h5open_f',hdferr                                     ! prevents error if not opened yet
  call h5close_f(hdferr)
  if (hdferr /= 0) write(6,'(a,i5)') ' Error in h5close_f',hdferr
 

src/results.f90

@@ -0,0 +1,974 @@
+!--------------------------------------------------------------------------------------------------
+!> @author Vitesh Shah, Max-Planck-Institut für Eisenforschung GmbH
+!> @author Yi-Chin Yang, Max-Planck-Institut für Eisenforschung GmbH
+!> @author Jennifer Nastola, Max-Planck-Institut für Eisenforschung GmbH
+!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
+!--------------------------------------------------------------------------------------------------
+module results
+  use prec
+  use IO
+  use HDF5
+  use HDF5_utilities
+#ifdef PETSc
+  use PETSC
+#endif
+
+ implicit none
+ private
+ integer(HID_T), public, protected :: tempCoordinates, tempResults
+ integer(HID_T), private :: resultsFile, currentIncID, plist_id
+
+
+ public :: &
+   results_init, &
+   results_openJobFile, &
+   results_closeJobFile, &
+   results_addIncrement, &
+   results_addGroup, &
+   results_openGroup, &
+   results_writeVectorDataset, &
+   results_setLink, &
+   results_removeLink
+contains
+
+subroutine results_init
+ use, intrinsic :: &
+   iso_fortran_env  ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment)
+ use DAMASK_interface, only: &
+   getSolverJobName
+ implicit none
+
+ write(6,'(/,a)') ' <<<+-  results init  -+>>>'
+#include "compilation_info.f90"
+
+ call HDF5_closeFile(HDF5_openFile(trim(getSolverJobName())//'.hdf5','w',.true.))
+
+end subroutine results_init
+
+
+!--------------------------------------------------------------------------------------------------
+!> @brief opens the results file to append data
+!--------------------------------------------------------------------------------------------------
+subroutine results_openJobFile()
+ use DAMASK_interface, only: &
+   getSolverJobName
+
+ implicit none
+ character(len=pStringLen) :: commandLine
+
+ resultsFile = HDF5_openFile(trim(getSolverJobName())//'.hdf5','a',.true.)
+ call HDF5_addAttribute(resultsFile,'DADF5',0.1_pReal)
+ call HDF5_addAttribute(resultsFile,'DAMASK',DAMASKVERSION)
+ call get_command(commandLine)
+ call HDF5_addAttribute(resultsFile,'call',trim(commandLine))
+ 
+end subroutine results_openJobFile
+
+
+!--------------------------------------------------------------------------------------------------
+!> @brief closes the results file
+!--------------------------------------------------------------------------------------------------
+subroutine results_closeJobFile()
+ implicit none
+
+ call HDF5_closeFile(resultsFile)
+
+end subroutine results_closeJobFile
+
+
+!--------------------------------------------------------------------------------------------------
+!> @brief closes the results file
+!--------------------------------------------------------------------------------------------------
+subroutine results_addIncrement(inc,time)
+ 
+ implicit none
+ integer(pInt), intent(in) :: inc
+ real(pReal),   intent(in) :: time
+ character(len=pStringLen) :: incChar
+
+ write(incChar,*) inc
+ call HDF5_closeGroup(results_addGroup(trim('inc'//trim(adjustl(incChar)))))
+ call results_setLink(trim('inc'//trim(adjustl(incChar))),'current')
+ call HDF5_addAttribute(resultsFile,'time/s',time,trim('inc'//trim(adjustl(incChar))))
+
+end subroutine results_addIncrement
+
+!--------------------------------------------------------------------------------------------------
+!> @brief open a group from the results file
+!--------------------------------------------------------------------------------------------------
+integer(HID_T) function results_openGroup(groupName)
+
+ implicit none
+ character(len=*), intent(in) :: groupName
+ 
+ results_openGroup = HDF5_openGroup(resultsFile,groupName)
+
+end function results_openGroup
+
+
+!--------------------------------------------------------------------------------------------------
+!> @brief adds a new group to the results file
+!--------------------------------------------------------------------------------------------------
+integer(HID_T) function results_addGroup(groupName)
+
+ implicit none
+ character(len=*), intent(in) :: groupName
+ 
+ results_addGroup = HDF5_addGroup(resultsFile,groupName)
+
+end function results_addGroup
+
+
+!--------------------------------------------------------------------------------------------------
+!> @brief set link to object in results file
+!--------------------------------------------------------------------------------------------------
+subroutine results_setLink(path,link)
+ use hdf5_utilities, only: &
+  HDF5_setLink
+
+ implicit none
+ character(len=*), intent(in) :: path, link
+
+ call HDF5_setLink(resultsFile,path,link)
+
+end subroutine results_setLink
+
+
+!--------------------------------------------------------------------------------------------------
+!> @brief remove link to an object
+!--------------------------------------------------------------------------------------------------
+subroutine results_removeLink(link)
+ use hdf5
+
+ implicit none
+ character(len=*), intent(in) :: link
+ integer                      :: hdferr
+
+ call h5ldelete_f(resultsFile,link, hdferr)
+ if (hdferr < 0) call IO_error(1_pInt,ext_msg = 'results_removeLink: h5ldelete_soft_f ('//trim(link)//')')
+
+end subroutine results_removeLink
+
+
+!--------------------------------------------------------------------------------------------------
+!> @brief stores a vector dataset in a group
+!--------------------------------------------------------------------------------------------------
+subroutine results_writeVectorDataset(group,dataset,label,SIunit)
+
+ implicit none
+ character(len=*), intent(in)                 :: SIunit,label,group
+ real(pReal),      intent(inout), dimension(:,:) :: dataset
+ integer(HID_T)        :: groupHandle
+ 
+ groupHandle = results_openGroup(group)
+ call HDF5_write(groupHandle,dataset,label)
+ if (HDF5_objectExists(groupHandle,label)) &
+   call HDF5_addAttribute(groupHandle,'Unit',SIunit,label)
+ call HDF5_closeGroup(groupHandle)
+
+end subroutine results_writeVectorDataset
+
+
+!--------------------------------------------------------------------------------------------------
+!> @brief adds the unique mapping from spatial position and constituent ID to results
+!--------------------------------------------------------------------------------------------------
+subroutine HDF5_mappingPhase(mapping,mapping2,Nconstituents,material_phase,phase_name,dataspace_size,mpiOffset,mpiOffset_phase)
+ use hdf5
+
+ implicit none
+ integer(pInt),    intent(in)                   :: Nconstituents, dataspace_size, mpiOffset
+ integer(pInt),    intent(in), dimension(:)     :: mapping, mapping2
+ character(len=*), intent(in), dimension(:)     :: phase_name
+ integer(pInt),    intent(in), dimension(:)     :: mpiOffset_phase
+ integer(pInt),    intent(in), dimension(:,:,:) :: material_phase
+
+ character(len=len(phase_name(1))), dimension(:), allocatable :: namesNA
+ character(len=len(phase_name(1)))                            :: a
+ character(len=*),                  parameter                 :: n = "NULL"
+
+ integer(pInt)   :: hdferr, NmatPoints, i, j, k
+ integer(HID_T)  :: mapping_id, dtype_id, dset_id, space_id, name_id, position_id, plist_id, memspace
+ integer(HID_T)  :: dt5_id      ! Memory datatype identifier
+ integer(SIZE_T) :: typesize, type_sizec, type_sizei, type_size
+
+ integer(HSIZE_T),  dimension(2)                :: counter
+ integer(HSSIZE_T), dimension(2)                :: fileOffset
+ integer(pInt),     dimension(:,:), allocatable :: arrOffset
+
+ a = n
+ allocate(namesNA(0:size(phase_name)),source=[a,phase_name])
+ NmatPoints = size(mapping,1)/Nconstituents
+ mapping_ID = results_openGroup("current/mapGeometry")
+
+ allocate(arrOffset(Nconstituents,NmatPoints))
+ do i=1_pInt, NmatPoints
+   do k=1_pInt, Nconstituents
+     do j=1_pInt, size(phase_name)
+       if(material_phase(k,1,i) == j) &
+         arrOffset(k,i) = mpiOffset_phase(j)
+     enddo
+   enddo
+ enddo
+
+!--------------------------------------------------------------------------------------------------
+! create dataspace
+ call h5screate_simple_f(2, int([Nconstituents,dataspace_size],HSIZE_T), space_id, hdferr, &
+                            int([Nconstituents,dataspace_size],HSIZE_T))
+ if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_writeMapping')
+
+!--------------------------------------------------------------------------------------------------
+! compound type
+     ! First calculate total size by calculating sizes of each member
+     !
+     CALL h5tcopy_f(H5T_NATIVE_CHARACTER, dt5_id, hdferr)
+     typesize = len(phase_name(1))
+     CALL h5tset_size_f(dt5_id, typesize, hdferr)
+     CALL h5tget_size_f(dt5_id, type_sizec, hdferr)
+     CALL h5tget_size_f(H5T_STD_I32LE,type_sizei, hdferr)
+     type_size = type_sizec + type_sizei
+ call h5tcreate_f(H5T_COMPOUND_F, type_size, dtype_id, hdferr)
+ if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_writeMapping: h5tcreate_f dtype_id')
+
+ call h5tinsert_f(dtype_id, "Name",          0_SIZE_T, dt5_id, hdferr)
+ if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingPhase: h5tinsert_f 0')
+ call h5tinsert_f(dtype_id, "Position",  type_sizec, H5T_STD_I32LE, hdferr)
+ if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingPhase: h5tinsert_f 2')
+
+!--------------------------------------------------------------------------------------------------
+! create Dataset
+ call h5dcreate_f(mapping_id, 'constitutive', dtype_id, space_id, dset_id, hdferr)
+ if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingPhase')
+
+!--------------------------------------------------------------------------------------------------
+! Create memory types (one compound datatype for each member)
+ call h5tcreate_f(H5T_COMPOUND_F, int(type_sizec,SIZE_T), name_id, hdferr)
+ if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingPhase: h5tcreate_f instance_id')
+ call h5tinsert_f(name_id, "Name",        0_SIZE_T, dt5_id, hdferr)
+ if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingPhase: h5tinsert_f instance_id')
+
+ call h5tcreate_f(H5T_COMPOUND_F, int(pInt,SIZE_T), position_id, hdferr)
+ if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingPhase: h5tcreate_f position_id')
+ call h5tinsert_f(position_id, "Position", 0_SIZE_T, H5T_STD_I32LE, hdferr)
+ if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingPhase: h5tinsert_f position_id')
+
+!--------------------------------------------------------------------------------------------------
+! Define and select hyperslabs
+ counter(1)    = Nconstituents                  ! how big i am
+ counter(2)    = NmatPoints
+ fileOffset(1) = 0                              ! where i start to write my data
+ fileOffset(2) = mpiOffset
+
+ call h5screate_simple_f(2, counter, memspace, hdferr)
+ if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingPhase: h5screate_simple_f')
+ call h5dget_space_f(dset_id, space_id, hdferr)
+ if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingPhase: h5dget_space_f')
+ call h5sselect_hyperslab_f(space_id, H5S_SELECT_SET_F, fileOffset, counter, hdferr)
+ if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingPhase: h5sselect_hyperslab_f')
+
+!--------------------------------------------------------------------------------------------------
+ ! Create property list for collective dataset write
+#ifdef PETSc
+ call h5pcreate_f(H5P_DATASET_XFER_F, plist_id, hdferr)
+ if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingPhase: h5pcreate_f')
+ call h5pset_dxpl_mpio_f(plist_id, H5FD_MPIO_COLLECTIVE_F, hdferr)
+ if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingPhase: h5pset_dxpl_mpio_f')
+#endif
+
+!--------------------------------------------------------------------------------------------------
+! write data by fields in the datatype. Fields order is not important.
+ call h5dwrite_f(dset_id, name_id, reshape(namesNA(mapping),[Nconstituents,NmatPoints]), &
+                          int([Nconstituents,  dataspace_size],HSIZE_T), hdferr, &
+                          file_space_id = space_id, mem_space_id = memspace, xfer_prp = plist_id)
+ if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingPhase: h5dwrite_f position_id')
+
+ call h5dwrite_f(dset_id, position_id, reshape(mapping2-1_pInt,[Nconstituents,NmatPoints])+arrOffset, &
+                          int([Nconstituents,  dataspace_size],HSIZE_T), hdferr, &
+                          file_space_id = space_id, mem_space_id = memspace, xfer_prp = plist_id)
+ if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingPhase: h5dwrite_f instance_id')
+
+!--------------------------------------------------------------------------------------------------
+! close types, dataspaces
+ call h5tclose_f(dtype_id, hdferr)
+ if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingPhase: h5tclose_f dtype_id')
+ call h5tclose_f(position_id, hdferr)
+ if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingPhase: h5tclose_f position_id')
+ call h5tclose_f(name_id, hdferr)
+ if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingPhase: h5tclose_f name_id  ')
+ call h5tclose_f(dt5_id, hdferr)
+ if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingPhase: h5tclose_f dt5_id')
+ call h5dclose_f(dset_id, hdferr)
+ if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingPhase: h5dclose_f')
+ call h5sclose_f(space_id, hdferr)
+ if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingPhase: h5sclose_f space_id')
+ call h5sclose_f(memspace, hdferr)
+ if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingPhase: h5sclose_f memspace')
+ call h5pclose_f(plist_id, hdferr)
+ if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingPhase: h5pclose_f')
+ call HDF5_closeGroup(mapping_ID)
+
+end subroutine HDF5_mappingPhase
+
+!--------------------------------------------------------------------------------------------------
+!> @brief adds the backward mapping from spatial position and constituent ID to results
+!--------------------------------------------------------------------------------------------------
+subroutine HDF5_backwardMappingPhase(material_phase,phasememberat,phase_name,dataspace_size,mpiOffset,mpiOffset_phase)
+ use hdf5
+
+ implicit none
+ integer(pInt),    intent(in), dimension(:,:,:) :: material_phase, phasememberat
+ character(len=*), intent(in), dimension(:)     :: phase_name
+ integer(pInt),    intent(in), dimension(:)     :: dataspace_size, mpiOffset_phase
+ integer(pInt),    intent(in)                   :: mpiOffset
+
+ integer(pInt)   :: hdferr, NmatPoints, Nconstituents, i, j
+ integer(HID_T)  :: mapping_id, dtype_id, dset_id, space_id, position_id, plist_id, memspace
+ integer(SIZE_T) :: type_size
+
+ integer(pInt), dimension(:,:), allocatable :: arr
+
+ integer(HSIZE_T),  dimension(1) :: counter
+ integer(HSSIZE_T), dimension(1) :: fileOffset
+
+ character(len=64) :: phaseID
+
+ Nconstituents = size(phasememberat,1)
+ NmatPoints = count(material_phase /=0_pInt)/Nconstituents
+
+ allocate(arr(2,NmatPoints*Nconstituents))
+
+ do i=1_pInt, NmatPoints
+   do j=Nconstituents-1_pInt, 0_pInt, -1_pInt
+     arr(1,Nconstituents*i-j) = i-1_pInt
+   enddo
+ enddo
+ arr(2,:) = pack(material_phase,material_phase/=0_pInt)
+
+ do i=1_pInt, size(phase_name)
+   write(phaseID, '(i0)') i
+   mapping_ID = results_openGroup('/current/constitutive/'//trim(phaseID)//'_'//phase_name(i))
+   NmatPoints = count(material_phase == i)
+
+!--------------------------------------------------------------------------------------------------
+  ! create dataspace
+   call h5screate_simple_f(1, int([dataspace_size(i)],HSIZE_T), space_id, hdferr, &
+                              int([dataspace_size(i)],HSIZE_T))
+   if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_writeBackwardMapping')
+
+!--------------------------------------------------------------------------------------------------
+  ! compound type
+   call h5tget_size_f(H5T_STD_I32LE, type_size, hdferr)
+   call h5tcreate_f(H5T_COMPOUND_F, type_size, dtype_id, hdferr)
+   if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_writeBackwardMapping: h5tcreate_f dtype_id')
+
+   call h5tinsert_f(dtype_id, "Position",          0_SIZE_T, H5T_STD_I32LE, hdferr)
+   if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_backwardMappingPhase: h5tinsert_f 0')
+
+!--------------------------------------------------------------------------------------------------
+  ! create Dataset
+   call h5dcreate_f(mapping_id, 'mapGeometry', dtype_id, space_id, dset_id, hdferr)
+   if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_backwardMappingPhase')
+
+!--------------------------------------------------------------------------------------------------
+  ! Create memory types (one compound datatype for each member)
+   call h5tcreate_f(H5T_COMPOUND_F, int(pInt,SIZE_T), position_id, hdferr)
+   if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_backwardMappingPhase: h5tcreate_f position_id')
+   call h5tinsert_f(position_id, "Position", 0_SIZE_T, H5T_STD_I32LE, hdferr)
+   if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_backwardMappingPhase: h5tinsert_f position_id')
+
+!--------------------------------------------------------------------------------------------------
+  ! Define and select hyperslabs
+   counter = NmatPoints                       ! how big i am
+   fileOffset = mpiOffset_phase(i)            ! where i start to write my data
+
+   call h5screate_simple_f(1, counter, memspace, hdferr)
+   if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_backwardMappingPhase: h5screate_simple_f')
+   call h5dget_space_f(dset_id, space_id, hdferr)
+   if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_backwardMappingPhase: h5dget_space_f')
+   call h5sselect_hyperslab_f(space_id, H5S_SELECT_SET_F, fileOffset, counter, hdferr)
+   if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_backwardMappingPhase: h5sselect_hyperslab_f')
+
+!--------------------------------------------------------------------------------------------------
+ ! Create property list for collective dataset write
+#ifdef PETSc
+   call h5pcreate_f(H5P_DATASET_XFER_F, plist_id, hdferr)
+   if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_backwardMappingPhase: h5pcreate_f')
+   call h5pset_dxpl_mpio_f(plist_id, H5FD_MPIO_COLLECTIVE_F, hdferr)
+   if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_backwardMappingPhase: h5pset_dxpl_mpio_f')
+#endif
+
+!--------------------------------------------------------------------------------------------------
+  ! write data by fields in the datatype. Fields order is not important.
+   call h5dwrite_f(dset_id, position_id, pack(arr(1,:),arr(2,:)==i)+mpiOffset, int([dataspace_size(i)],HSIZE_T),&
+                              hdferr, file_space_id = space_id, mem_space_id = memspace, xfer_prp = plist_id)
+   if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_backwardMappingPhase: h5dwrite_f instance_id')
+
+!--------------------------------------------------------------------------------------------------
+  !close types, dataspaces
+   call h5tclose_f(dtype_id, hdferr)
+   if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_backwardMappingPhase: h5tclose_f dtype_id')
+   call h5tclose_f(position_id, hdferr)
+   if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_backwardMappingPhase: h5tclose_f position_id')
+   call h5dclose_f(dset_id, hdferr)
+   if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_backwardMappingPhase: h5dclose_f')
+   call h5sclose_f(space_id, hdferr)
+   if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_backwardMappingPhase: h5sclose_f space_id')
+   call h5sclose_f(memspace, hdferr)
+   if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_backwardMappingPhase: h5sclose_f memspace')
+   call h5pclose_f(plist_id, hdferr)
+   if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_backwardMappingPhase: h5pclose_f')
+   call HDF5_closeGroup(mapping_ID)
+
+ enddo
+
+end subroutine HDF5_backwardMappingPhase
+
+!--------------------------------------------------------------------------------------------------
+!> @brief adds the unique mapping from spatial position and constituent ID to results
+!--------------------------------------------------------------------------------------------------
+subroutine HDF5_mappingHomog(material_homog,homogmemberat,homogenization_name,dataspace_size,mpiOffset,mpiOffset_homog)
+ use hdf5
+
+ implicit none
+ integer(pInt),    intent(in), dimension(:,:) :: material_homog, homogmemberat
+ character(len=*), intent(in), dimension(:)   :: homogenization_name
+ integer(pInt),    intent(in), dimension(:)   :: mpiOffset_homog
+ integer(pInt),    intent(in)                 :: dataspace_size, mpiOffset
+
+ integer(pInt)   :: hdferr, NmatPoints, i, j
+ integer(HID_T)  :: mapping_id, dtype_id, dset_id, space_id, name_id, position_id, plist_id, memspace
+
+ integer(HID_T)  :: dt5_id      ! Memory datatype identifier
+ integer(SIZE_T) :: typesize, type_sizec, type_sizei, type_size
+
+ integer(HSIZE_T),  dimension(1)              :: counter
+ integer(HSSIZE_T), dimension(1)              :: fileOffset
+ integer(pInt),     dimension(:), allocatable :: arrOffset
+
+ NmatPoints = count(material_homog /=0_pInt)
+ mapping_ID = results_openGroup("current/mapGeometry")
+
+ allocate(arrOffset(NmatPoints))
+ do i=1_pInt, NmatPoints
+   do j=1_pInt, size(homogenization_name)
+     if(material_homog(1,i) == j) &
+       arrOffset(i) = mpiOffset_homog(j)
+   enddo
+ enddo
+
+!--------------------------------------------------------------------------------------------------
+! create dataspace
+ call h5screate_simple_f(1, int([dataspace_size],HSIZE_T), space_id, hdferr, &
+                            int([dataspace_size],HSIZE_T))
+ if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_writeMapping')
+
+!--------------------------------------------------------------------------------------------------
+! compound type
+     ! First calculate total size by calculating sizes of each member
+     !
+     CALL h5tcopy_f(H5T_NATIVE_CHARACTER, dt5_id, hdferr)
+     typesize = len(homogenization_name(1))
+     CALL h5tset_size_f(dt5_id, typesize, hdferr)
+     CALL h5tget_size_f(dt5_id, type_sizec, hdferr)
+     CALL h5tget_size_f(H5T_STD_I32LE,type_sizei, hdferr)
+     type_size = type_sizec + type_sizei
+ call h5tcreate_f(H5T_COMPOUND_F, type_size, dtype_id, hdferr)
+ if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_writeMapping: h5tcreate_f dtype_id')
+
+ call h5tinsert_f(dtype_id, "Name",          0_SIZE_T, dt5_id, hdferr)
+ if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingHomog: h5tinsert_f 0')
+ call h5tinsert_f(dtype_id, "Position",  type_sizec, H5T_STD_I32LE, hdferr)
+ if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingHomog: h5tinsert_f 2')
+
+!--------------------------------------------------------------------------------------------------
+! create Dataset
+ call h5dcreate_f(mapping_id, 'homogenization', dtype_id, space_id, dset_id, hdferr)
+ if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingHomog')
+
+!--------------------------------------------------------------------------------------------------
+! Create memory types (one compound datatype for each member)
+ call h5tcreate_f(H5T_COMPOUND_F, int(type_sizec,SIZE_T), name_id, hdferr)
+ if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingHomog: h5tcreate_f instance_id')
+ call h5tinsert_f(name_id, "Name",        0_SIZE_T, dt5_id, hdferr)
+ if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingHomog: h5tinsert_f instance_id')
+
+ call h5tcreate_f(H5T_COMPOUND_F, int(pInt,SIZE_T), position_id, hdferr)
+ if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingHomog: h5tcreate_f position_id')
+ call h5tinsert_f(position_id, "Position", 0_SIZE_T, H5T_STD_I32LE, hdferr)
+ if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingHomog: h5tinsert_f position_id')
+
+!--------------------------------------------------------------------------------------------------
+! Define and select hyperslabs
+ counter = NmatPoints                    ! how big i am
+ fileOffset = mpiOffset                  ! where i start to write my data
+
+ call h5screate_simple_f(1, counter, memspace, hdferr)
+ if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingHomog: h5screate_simple_f')
+ call h5dget_space_f(dset_id, space_id, hdferr)
+ if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingHomog: h5dget_space_f')
+ call h5sselect_hyperslab_f(space_id, H5S_SELECT_SET_F, fileOffset, counter, hdferr)
+ if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingHomog: h5sselect_hyperslab_f')
+
+!--------------------------------------------------------------------------------------------------
+! Create property list for collective dataset write
+#ifdef PETSc
+ call h5pcreate_f(H5P_DATASET_XFER_F, plist_id, hdferr)
+ if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingHomog: h5pcreate_f')
+ call h5pset_dxpl_mpio_f(plist_id, H5FD_MPIO_COLLECTIVE_F, hdferr)
+ if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingHomog: h5pset_dxpl_mpio_f')
+#endif
+
+!--------------------------------------------------------------------------------------------------
+! write data by fields in the datatype. Fields order is not important.
+ call h5dwrite_f(dset_id, name_id, homogenization_name(pack(material_homog,material_homog/=0_pInt)), &
+                          int([dataspace_size],HSIZE_T), hdferr, file_space_id = space_id, &
+                          mem_space_id = memspace, xfer_prp = plist_id)
+ if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingHomog: h5dwrite_f position_id')
+
+ call h5dwrite_f(dset_id, position_id, pack(homogmemberat-1_pInt,homogmemberat/=0_pInt) + arrOffset, &
+                          int([dataspace_size],HSIZE_T), hdferr, file_space_id = space_id, &
+                          mem_space_id = memspace, xfer_prp = plist_id)
+ if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingHomog: h5dwrite_f instance_id')
+
+!--------------------------------------------------------------------------------------------------
+!close types, dataspaces
+call h5tclose_f(dtype_id, hdferr)
+if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingHomog: h5tclose_f dtype_id')
+call h5tclose_f(position_id, hdferr)
+if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingHomog: h5tclose_f position_id')
+call h5tclose_f(name_id, hdferr)
+if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingHomog: h5tclose_f name_id  ')
+call h5tclose_f(dt5_id, hdferr)
+if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingHomog: h5tclose_f dt5_id')
+call h5dclose_f(dset_id, hdferr)
+if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingHomog: h5dclose_f')
+call h5sclose_f(space_id, hdferr)
+if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingHomog: h5sclose_f space_id')
+call h5sclose_f(memspace, hdferr)
+if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingHomog: h5sclose_f memspace')
+call h5pclose_f(plist_id, hdferr)
+if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingHomog: h5pclose_f')
+call HDF5_closeGroup(mapping_ID)
+
+end subroutine HDF5_mappingHomog
+
+!--------------------------------------------------------------------------------------------------
+!> @brief adds the backward mapping from spatial position and constituent ID to results
+!--------------------------------------------------------------------------------------------------
+subroutine HDF5_backwardMappingHomog(material_homog,homogmemberat,homogenization_name,dataspace_size,mpiOffset,mpiOffset_homog)
+ use hdf5
+
+ implicit none
+ integer(pInt),    intent(in), dimension(:,:) :: material_homog, homogmemberat
+ character(len=*), intent(in), dimension(:)   :: homogenization_name
+ integer(pInt),    intent(in), dimension(:)   :: dataspace_size, mpiOffset_homog
+ integer(pInt),    intent(in)                 :: mpiOffset
+
+ integer(pInt)   :: hdferr, NmatPoints, i
+ integer(HID_T)  :: mapping_id, dtype_id, dset_id, space_id, position_id, plist_id, memspace
+ integer(SIZE_T) :: type_size
+
+ integer(pInt),     dimension(:,:), allocatable :: arr
+
+ integer(HSIZE_T),  dimension(1) :: counter
+ integer(HSSIZE_T), dimension(1) :: fileOffset
+
+ character(len=64) :: homogID
+
+ NmatPoints = count(material_homog /=0_pInt)
+ allocate(arr(2,NmatPoints))
+
+ arr(1,:) = (/(i, i=0_pint,NmatPoints-1_pInt)/)
+ arr(2,:) = pack(material_homog,material_homog/=0_pInt)
+
+ do i=1_pInt, size(homogenization_name)
+   write(homogID, '(i0)') i
+   mapping_ID = results_openGroup('/current/homogenization/'//trim(homogID)//'_'//homogenization_name(i))
+
+!--------------------------------------------------------------------------------------------------
+  ! create dataspace
+   call h5screate_simple_f(1, int([dataspace_size(i)],HSIZE_T), space_id, hdferr, &
+                              int([dataspace_size(i)],HSIZE_T))
+   if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_writeBackwardMapping')
+
+!--------------------------------------------------------------------------------------------------
+  ! compound type
+   call h5tget_size_f(H5T_STD_I32LE, type_size, hdferr)
+   call h5tcreate_f(H5T_COMPOUND_F, type_size, dtype_id, hdferr)
+   if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_writeBackwardMapping: h5tcreate_f dtype_id')
+
+   call h5tinsert_f(dtype_id, "Position",          0_SIZE_T, H5T_STD_I32LE, hdferr)
+   if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_backwardMappingHomog: h5tinsert_f 0')
+
+!--------------------------------------------------------------------------------------------------
+  ! create Dataset
+   call h5dcreate_f(mapping_id, 'mapGeometry', dtype_id, space_id, dset_id, hdferr)
+   if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_backwardMappingHomog')
+
+!--------------------------------------------------------------------------------------------------
+  ! Create memory types (one compound datatype for each member)
+   call h5tcreate_f(H5T_COMPOUND_F, int(pInt,SIZE_T), position_id, hdferr)
+   if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_backwardMappingHomog: h5tcreate_f position_id')
+   call h5tinsert_f(position_id, "Position", 0_SIZE_T, H5T_STD_I32LE, hdferr)
+   if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_backwardMappingHomog: h5tinsert_f position_id')
+
+!--------------------------------------------------------------------------------------------------
+  ! Define and select hyperslabs
+   counter = NmatPoints                           ! how big i am
+   fileOffset = mpiOffset_homog(i)                ! where i start to write my data
+
+   call h5screate_simple_f(1, counter, memspace, hdferr)
+   if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_backwardMappingHomog: h5screate_simple_f')
+   call h5dget_space_f(dset_id, space_id, hdferr)
+   if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_backwardMappingHomog: h5dget_space_f')
+   call h5sselect_hyperslab_f(space_id, H5S_SELECT_SET_F, fileOffset, counter, hdferr)
+   if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_backwardMappingHomog: h5sselect_hyperslab_f')
+
+!--------------------------------------------------------------------------------------------------
+ ! Create property list for collective dataset write
+#ifdef PETSc
+   call h5pcreate_f(H5P_DATASET_XFER_F, plist_id, hdferr)
+   if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_backwardMappingHomog: h5pcreate_f')
+   call h5pset_dxpl_mpio_f(plist_id, H5FD_MPIO_COLLECTIVE_F, hdferr)
+   if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_backwardMappingHomog: h5pset_dxpl_mpio_f')
+#endif
+
+!--------------------------------------------------------------------------------------------------
+  ! write data by fields in the datatype. Fields order is not important.
+   call h5dwrite_f(dset_id, position_id, pack(arr(1,:),arr(2,:)==i)+mpiOffset,int([dataspace_size(i)],HSIZE_T),&
+                   hdferr, file_space_id = space_id, mem_space_id = memspace, xfer_prp = plist_id)
+   if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_backwardMappingHomog: h5dwrite_f instance_id')
+
+!--------------------------------------------------------------------------------------------------
+  !close types, dataspaces
+   call h5tclose_f(dtype_id, hdferr)
+   if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_backwardMappingHomog: h5tclose_f dtype_id')
+   call h5tclose_f(position_id, hdferr)
+   if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_backwardMappingHomog: h5tclose_f position_id')
+   call h5dclose_f(dset_id, hdferr)
+   if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_backwardMappingHomog: h5dclose_f')
+   call h5sclose_f(space_id, hdferr)
+   if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_backwardMappingHomog: h5sclose_f space_id')
+   call h5sclose_f(memspace, hdferr)
+   if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_backwardMappingHomog: h5sclose_f memspace')
+   call h5pclose_f(plist_id, hdferr)
+   if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_backwardMappingHomog: h5pclose_f')
+   call HDF5_closeGroup(mapping_ID)
+
+ enddo
+
+end subroutine HDF5_backwardMappingHomog
+
+!--------------------------------------------------------------------------------------------------
+!> @brief adds the unique mapping from spatial position and constituent ID to results
+!--------------------------------------------------------------------------------------------------
+subroutine HDF5_mappingCrystallite(crystalliteAt,crystmemberAt,crystallite_name,dataspace_size,mpiOffset,mpiOffset_cryst)
+ use hdf5
+
+ implicit none
+ integer(pInt),    intent(in), dimension(:,:)   :: crystalliteAt
+ integer(pInt),    intent(in), dimension(:,:,:) :: crystmemberAt
+ character(len=*), intent(in), dimension(:)     :: crystallite_name
+ integer(pInt),    intent(in), dimension(:)     :: mpiOffset_cryst
+ integer(pInt),    intent(in)                   :: dataspace_size, mpiOffset
+
+ integer         :: hdferr
+ integer(pInt)   :: NmatPoints, Nconstituents, i, j
+ integer(HID_T)  :: mapping_id, dtype_id, dset_id, space_id, name_id, plist_id, memspace
+
+ integer(HID_T), dimension(:), allocatable :: position_id
+
+ integer(HID_T)  :: dt5_id      ! Memory datatype identifier
+ integer(SIZE_T) :: typesize, type_sizec, type_sizei, type_size
+
+ integer(HSIZE_T),  dimension(1)              :: counter
+ integer(HSSIZE_T), dimension(1)              :: fileOffset
+ integer(pInt),     dimension(:), allocatable :: arrOffset
+
+ character(len=64) :: m
+
+ Nconstituents = size(crystmemberAt,1)
+ NmatPoints = count(crystalliteAt /=0_pInt)
+ mapping_ID = results_openGroup("current/mapGeometry")
+
+ allocate(position_id(Nconstituents))
+
+ allocate(arrOffset(NmatPoints))
+ do i=1_pInt, NmatPoints
+   do j=1_pInt, size(crystallite_name)
+     if(crystalliteAt(1,i) == j) &
+       arrOffset(i) = Nconstituents*mpiOffset_cryst(j)
+   enddo
+ enddo
+
+!--------------------------------------------------------------------------------------------------
+! create dataspace
+ call h5screate_simple_f(1, int([dataspace_size],HSIZE_T), space_id, hdferr, &
+                            int([dataspace_size],HSIZE_T))
+ if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_writeMapping')
+
+!--------------------------------------------------------------------------------------------------
+! compound type
+     ! First calculate total size by calculating sizes of each member
+     !
+     CALL h5tcopy_f(H5T_NATIVE_CHARACTER, dt5_id, hdferr)
+     typesize = len(crystallite_name(1))
+     CALL h5tset_size_f(dt5_id, typesize, hdferr)
+     CALL h5tget_size_f(dt5_id, type_sizec, hdferr)
+     CALL h5tget_size_f(H5T_STD_I32LE, type_sizei, hdferr)
+     type_size = type_sizec + type_sizei*Nconstituents
+ call h5tcreate_f(H5T_COMPOUND_F, type_size, dtype_id, hdferr)
+ if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_writeMapping: h5tcreate_f dtype_id')
+
+ call h5tinsert_f(dtype_id, "Name",          0_SIZE_T, dt5_id, hdferr)
+ if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingCrystallite: h5tinsert_f 0')
+ do i=1_pInt, Nconstituents
+   write(m, '(i0)') i
+   call h5tinsert_f(dtype_id, "Position "//trim(m),  type_sizec+(i-1)*type_sizei, H5T_STD_I32LE, hdferr)
+   if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingCrystallite: h5tinsert_f 2 '//trim(m))
+ enddo
+
+!--------------------------------------------------------------------------------------------------
+! create Dataset
+ call h5dcreate_f(mapping_id, 'crystallite', dtype_id, space_id, dset_id, hdferr)
+ if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingCrystallite')
+
+!--------------------------------------------------------------------------------------------------
+! Create memory types (one compound datatype for each member)
+ call h5tcreate_f(H5T_COMPOUND_F, int(type_sizec,SIZE_T), name_id, hdferr)
+ if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingCrystallite: h5tcreate_f instance_id')
+ call h5tinsert_f(name_id, "Name",        0_SIZE_T, dt5_id, hdferr)
+ if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingCrystallite: h5tinsert_f instance_id')
+
+ do i=1_pInt, Nconstituents
+   write(m, '(i0)') i
+   call h5tcreate_f(H5T_COMPOUND_F, int(pInt,SIZE_T), position_id(i), hdferr)
+   if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingCrystallite: h5tcreate_f position_id')
+   call h5tinsert_f(position_id(i), "Position "//trim(m), 0_SIZE_T, H5T_STD_I32LE, hdferr)
+   if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingCrystallite: h5tinsert_f position_id')
+ enddo
+
+!--------------------------------------------------------------------------------------------------
+! Define and select hyperslabs
+ counter = NmatPoints                    ! how big i am
+ fileOffset = mpiOffset                  ! where i start to write my data
+
+ call h5screate_simple_f(1, counter, memspace, hdferr)
+ if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingCrystallite: h5screate_simple_f')
+ call h5dget_space_f(dset_id, space_id, hdferr)
+ if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingCrystallite: h5dget_space_f')
+ call h5sselect_hyperslab_f(space_id, H5S_SELECT_SET_F, fileOffset, counter, hdferr)
+ if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingCrystallite: h5sselect_hyperslab_f')
+
+!--------------------------------------------------------------------------------------------------
+ ! Create property list for collective dataset write
+#ifdef PETSc
+ call h5pcreate_f(H5P_DATASET_XFER_F, plist_id, hdferr)
+ if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingCrystallite: h5pcreate_f')
+ call h5pset_dxpl_mpio_f(plist_id, H5FD_MPIO_COLLECTIVE_F, hdferr)
+ if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingCrystallite: h5pset_dxpl_mpio_f')
+#endif
+
+!--------------------------------------------------------------------------------------------------
+! write data by fields in the datatype. Fields order is not important.
+ call h5dwrite_f(dset_id, name_id, crystallite_name(pack(crystalliteAt,crystalliteAt/=0_pInt)), &
+                          int([dataspace_size],HSIZE_T), hdferr, file_space_id = space_id, &
+                          mem_space_id = memspace, xfer_prp = plist_id)
+ if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingCrystallite: h5dwrite_f position_id')
+
+ do i=1_pInt, Nconstituents
+   call h5dwrite_f(dset_id, position_id(i), pack(crystmemberAt(i,:,:)-1_pInt,crystmemberAt(i,:,:)/=0_pInt)+arrOffset,&
+                            int([dataspace_size],HSIZE_T), hdferr, file_space_id = space_id, &
+                            mem_space_id = memspace, xfer_prp = plist_id)
+   if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingCrystallite: h5dwrite_f instance_id')
+ enddo
+
+!--------------------------------------------------------------------------------------------------
+!close types, dataspaces
+ call h5tclose_f(dtype_id, hdferr)
+ if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingCrystallite: h5tclose_f dtype_id')
+ do i=1_pInt, Nconstituents
+   call h5tclose_f(position_id(i), hdferr)
+ enddo
+ if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingCrystallite: h5tclose_f position_id')
+ call h5tclose_f(name_id, hdferr)
+ if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingCrystallite: h5tclose_f name_id')
+ call h5tclose_f(dt5_id, hdferr)
+ if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingCrystallite: h5tclose_f dt5_id')
+ call h5dclose_f(dset_id, hdferr)
+ if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingCrystallite: h5dclose_f')
+ call h5sclose_f(space_id, hdferr)
+ if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingCrystallite: h5sclose_f space_id')
+ call h5sclose_f(memspace, hdferr)
+ if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingCrystallite: h5sclose_f memspace')
+ call h5pclose_f(plist_id, hdferr)
+ if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingCrystallite: h5pclose_f')
+ call HDF5_closeGroup(mapping_ID)
+
+end subroutine HDF5_mappingCrystallite
+
+
+!--------------------------------------------------------------------------------------------------
+!> @brief adds the backward mapping from spatial position and constituent ID to results
+!--------------------------------------------------------------------------------------------------
+subroutine HDF5_backwardMappingCrystallite(crystalliteAt,crystmemberAt,crystallite_name,dataspace_size,mpiOffset,mpiOffset_cryst)
+ use hdf5
+
+ implicit none
+ integer(pInt),    intent(in), dimension(:,:)   :: crystalliteAt
+ integer(pInt),    intent(in), dimension(:,:,:) :: crystmemberAt
+ character(len=*), intent(in), dimension(:)     :: crystallite_name
+ integer(pInt),    intent(in), dimension(:)     :: dataspace_size, mpiOffset_cryst
+ integer(pInt),    intent(in)                   :: mpiOffset
+
+ integer         :: hdferr
+ integer(pInt)   :: NmatPoints, Nconstituents, i, j
+ integer(HID_T)  :: mapping_id, dtype_id, dset_id, space_id, position_id, plist_id, memspace
+ integer(SIZE_T) :: type_size
+
+ integer(pInt), dimension(:,:), allocatable     :: h_arr, arr
+
+ integer(HSIZE_T),  dimension(1) :: counter
+ integer(HSSIZE_T), dimension(1) :: fileOffset
+
+ character(len=64) :: crystallID
+
+ Nconstituents = size(crystmemberAt,1)
+ NmatPoints = count(crystalliteAt /=0_pInt)
+
+ allocate(h_arr(2,NmatPoints))
+ allocate(arr(2,Nconstituents*NmatPoints))
+
+ h_arr(1,:) = (/(i, i=0_pInt,NmatPoints-1_pInt)/)
+ h_arr(2,:) = pack(crystalliteAt,crystalliteAt/=0_pInt)
+
+ do i=1_pInt, NmatPoints
+   do j=Nconstituents-1_pInt, 0_pInt, -1_pInt
+     arr(1,Nconstituents*i-j) = h_arr(1,i)
+     arr(2,Nconstituents*i-j) = h_arr(2,i)
+   enddo
+ enddo
+
+ do i=1_pInt, size(crystallite_name)
+   if (crystallite_name(i) == 'none') cycle
+   write(crystallID, '(i0)') i
+   mapping_ID = results_openGroup('/current/crystallite/'//trim(crystallID)//'_'//crystallite_name(i))
+   NmatPoints = count(crystalliteAt == i)
+
+!--------------------------------------------------------------------------------------------------
+  ! create dataspace
+   call h5screate_simple_f(1, int([Nconstituents*dataspace_size(i)],HSIZE_T), space_id, hdferr, &
+                              int([Nconstituents*dataspace_size(i)],HSIZE_T))
+   if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_writeBackwardMapping')
+
+!--------------------------------------------------------------------------------------------------
+  ! compound type
+   call h5tget_size_f(H5T_STD_I32LE, type_size, hdferr)
+   call h5tcreate_f(H5T_COMPOUND_F, type_size, dtype_id, hdferr)
+   if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_writeBackwardMapping: h5tcreate_f dtype_id')
+
+   call h5tinsert_f(dtype_id, "Position",          0_SIZE_T, H5T_STD_I32LE, hdferr)
+   if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_backwardMappingCrystallite: h5tinsert_f 0')
+
+!--------------------------------------------------------------------------------------------------
+  ! create Dataset
+   call h5dcreate_f(mapping_id, 'mapGeometry', dtype_id, space_id, dset_id, hdferr)
+   if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_backwardMappingCrystallite')
+
+!--------------------------------------------------------------------------------------------------
+  ! Create memory types
+   call h5tcreate_f(H5T_COMPOUND_F, int(pInt,SIZE_T), position_id, hdferr)
+   if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_backwardMappingCrystallite: h5tcreate_f position_id')
+   call h5tinsert_f(position_id, "Position", 0_SIZE_T, H5T_STD_I32LE, hdferr)
+   if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_backwardMappingCrystallite: h5tinsert_f position_id')
+
+!--------------------------------------------------------------------------------------------------
+  ! Define and select hyperslabs
+   counter = Nconstituents*NmatPoints                       ! how big i am
+   fileOffset = Nconstituents*mpiOffset_cryst(i)            ! where i start to write my data
+
+   call h5screate_simple_f(1, counter, memspace, hdferr)
+   if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_backwardMappingCrystallite: h5screate_simple_f')
+   call h5dget_space_f(dset_id, space_id, hdferr)
+   if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_backwardMappingCrystallite: h5dget_space_f')
+   call h5sselect_hyperslab_f(space_id, H5S_SELECT_SET_F, fileOffset, counter, hdferr)
+   if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_backwardMappingCrystallite: h5sselect_hyperslab_f')
+
+!--------------------------------------------------------------------------------------------------
+ ! Create property list for collective dataset write
+#ifdef PETSc
+   call h5pcreate_f(H5P_DATASET_XFER_F, plist_id, hdferr)
+   if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_backwardMappingCrystallite: h5pcreate_f')
+   call h5pset_dxpl_mpio_f(plist_id, H5FD_MPIO_COLLECTIVE_F, hdferr)
+   if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_backwardMappingCrystallite: h5pset_dxpl_mpio_f')
+#endif
+
+!--------------------------------------------------------------------------------------------------
+  ! write data by fields in the datatype. Fields order is not important.
+   call h5dwrite_f(dset_id, position_id, pack(arr(1,:),arr(2,:)==i) + mpiOffset,&
+                              int([Nconstituents*dataspace_size(i)],HSIZE_T), hdferr, file_space_id = space_id, &
+                              mem_space_id = memspace, xfer_prp = plist_id)
+   if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_backwardMappingCrystallite: h5dwrite_f instance_id')
+
+!--------------------------------------------------------------------------------------------------
+  !close types, dataspaces
+   call h5tclose_f(dtype_id, hdferr)
+   if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_backwardMappingCrystallite: h5tclose_f dtype_id')
+   call h5tclose_f(position_id, hdferr)
+   if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_backwardMappingCrystallite: h5tclose_f position_id')
+   call h5dclose_f(dset_id, hdferr)
+   if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_backwardMappingCrystallite: h5dclose_f')
+   call h5sclose_f(space_id, hdferr)
+   if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_backwardMappingCrystallite: h5sclose_f space_id')
+   call h5sclose_f(memspace, hdferr)
+   if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_backwardMappingCrystallite: h5sclose_f memspace')
+   call h5pclose_f(plist_id, hdferr)
+   if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_backwardMappingCrystallite: h5pclose_f')
+   call HDF5_closeGroup(mapping_ID)
+
+ enddo
+
+end subroutine HDF5_backwardMappingCrystallite
+
+!--------------------------------------------------------------------------------------------------
+!> @brief adds the unique cell to node mapping
+!--------------------------------------------------------------------------------------------------
+subroutine HDF5_mappingCells(mapping)
+ use hdf5
+
+ implicit none
+ integer(pInt), intent(in), dimension(:) :: mapping
+
+ integer        :: hdferr, Nnodes
+ integer(HID_T) :: mapping_id, dset_id, space_id
+
+ Nnodes=size(mapping)
+ mapping_ID = results_openGroup("mapping")
+
+!--------------------------------------------------------------------------------------------------
+! create dataspace
+ call h5screate_simple_f(1, int([Nnodes],HSIZE_T), space_id, hdferr, &
+                            int([Nnodes],HSIZE_T))
+ if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingCells: h5screate_simple_f')
+
+!--------------------------------------------------------------------------------------------------
+! create Dataset
+ call h5dcreate_f(mapping_id, "Cell",H5T_NATIVE_INTEGER, space_id, dset_id, hdferr)
+ if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingCells')
+
+!--------------------------------------------------------------------------------------------------
+! write data by fields in the datatype. Fields order is not important.
+ call h5dwrite_f(dset_id, H5T_NATIVE_INTEGER, mapping, int([Nnodes],HSIZE_T), hdferr)
+ if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingCells: h5dwrite_f instance_id')
+
+!--------------------------------------------------------------------------------------------------
+!close types, dataspaces
+ call h5dclose_f(dset_id, hdferr)
+ if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingConstitutive: h5dclose_f')
+ call h5sclose_f(space_id, hdferr)
+ if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingConstitutive: h5sclose_f')
+ call HDF5_closeGroup(mapping_ID)
+
+end subroutine HDF5_mappingCells
+
+
+end module results

src/spectral_mech_Basic.f90

@@ -80,6 +80,7 @@ subroutine basic_init
 #endif
  use IO, only: &
    IO_intOut, &
+   IO_error, &
    IO_read_realFile, &
    IO_timeStamp
  use debug, only: &
@@ -173,7 +174,11 @@ subroutine basic_init
    call IO_read_realFile(777,'F_aimDot',trim(getSolverJobName()),size(F_aimDot))
    read (777,rec=1) F_aimDot; close (777)
    F_aim         = reshape(sum(sum(sum(F,dim=4),dim=3),dim=2) * wgt, [3,3])                         ! average of F
+   call MPI_Allreduce(MPI_IN_PLACE,F_aim,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr)
+   if(ierr /=0_pInt) call IO_error(894_pInt, ext_msg='F_aim')
    F_aim_lastInc = sum(sum(sum(F_lastInc,dim=5),dim=4),dim=3) * wgt                                 ! average of F_lastInc 
+   call MPI_Allreduce(MPI_IN_PLACE,F_aim_lastInc,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr)
+   if(ierr /=0_pInt) call IO_error(894_pInt, ext_msg='F_aim_lastInc')
  elseif (restartInc == 0_pInt) then restart
    F_lastInc = spread(spread(spread(math_I3,3,grid(1)),4,grid(2)),5,grid3)                          ! initialize to identity
    F = reshape(F_lastInc,[9,grid(1),grid(2),grid3])

src/spectral_mech_Polarisation.f90

@@ -78,7 +78,6 @@ contains
 
 !--------------------------------------------------------------------------------------------------
 !> @brief allocates all necessary fields and fills them with data, potentially from restart info
-!> @todo use sourced allocation, e.g. allocate(Fdot,source = F_lastInc)
 !--------------------------------------------------------------------------------------------------
 subroutine Polarisation_init
 #if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800
@@ -88,6 +87,7 @@ subroutine Polarisation_init
 #endif
  use IO, only: &
    IO_intOut, &
+   IO_error, &
    IO_read_realFile, &
    IO_timeStamp
  use debug, only: &
@@ -191,7 +191,11 @@ subroutine Polarisation_init
    call IO_read_realFile(777,'F_aimDot',trim(getSolverJobName()),size(F_aimDot))
    read (777,rec=1) F_aimDot; close (777)
    F_aim         = reshape(sum(sum(sum(F,dim=4),dim=3),dim=2) * wgt, [3,3])                         ! average of F
+   call MPI_Allreduce(MPI_IN_PLACE,F_aim,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr)
+   if(ierr /=0_pInt) call IO_error(894_pInt, ext_msg='F_aim')
    F_aim_lastInc = sum(sum(sum(F_lastInc,dim=5),dim=4),dim=3) * wgt                                 ! average of F_lastInc 
+   call MPI_Allreduce(MPI_IN_PLACE,F_aim_lastInc,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr)
+   if(ierr /=0_pInt) call IO_error(894_pInt, ext_msg='F_aim_lastInc')
  elseif (restartInc == 0_pInt) then restart
    F_lastInc = spread(spread(spread(math_I3,3,grid(1)),4,grid(2)),5,grid3)                          ! initialize to identity
    F = reshape(F_lastInc,[9,grid(1),grid(2),grid3])