Merge branch 'development' into improve-Lp-guessing

.gitlab-ci.yml

@@ -213,6 +213,14 @@ Post_OrientationConversion:
     - master
     - release
 
+Post_OrientationAverageMisorientation:
+  stage: postprocessing
+  script:
+    - OrientationAverageMisorientation/test.py
+  except:
+    - master
+    - release
+
 ###################################################################################################
 grid_mech_compile_Intel:
   stage: compilePETSc

PRIVATE

@@ -1 +1 @@
-Subproject commit 76d367e297c054de78f5568074470b047427395b
+Subproject commit 183d7a3a3bafa0a308c3eac858ca03c08fc03d50

VERSION

@@ -1 +1 @@
-v2.0.3-152-gd74599d3
+v2.0.3-252-g0e77be25

processing/post/DADF5_postResults.py

@@ -0,0 +1,83 @@
+#!/usr/bin/env python3
+# -*- coding: UTF-8 no BOM -*-
+
+import os
+import numpy as np
+import argparse
+import damask
+
+scriptName = os.path.splitext(os.path.basename(__file__))[0]
+scriptID   = ' '.join([scriptName,damask.version])
+
+# --------------------------------------------------------------------
+#                                MAIN
+# --------------------------------------------------------------------
+parser = argparse.ArgumentParser()
+
+#ToDo:  We need to decide on a way of handling arguments of variable lentght
+#https://stackoverflow.com/questions/15459997/passing-integer-lists-to-python
+
+#parser.add_argument('--version', action='version', version='%(prog)s {}'.format(scriptID))
+parser.add_argument('filenames', nargs='+',
+                    help='DADF5 files')
+parser.add_argument('-d','--dir', dest='dir',default='postProc',metavar='string',
+                    help='name of subdirectory to hold output')
+
+options = parser.parse_args()
+
+options.labels = ['Fe','Fp','xi_sl']
+
+# --- loop over input files ------------------------------------------------------------------------
+
+for filename in options.filenames:
+  results = damask.DADF5(filename)
+  
+  if not results.structured: continue
+  delta = results.size/results.grid*0.5
+  x, y, z = np.meshgrid(np.linspace(delta[2],results.size[2]-delta[2],results.grid[2]),
+                        np.linspace(delta[1],results.size[1]-delta[1],results.grid[1]),
+                        np.linspace(delta[0],results.size[0]-delta[0],results.grid[0]),
+                        indexing = 'ij')
+
+  coords = np.concatenate((z[:,:,:,None],y[:,:,:,None],x[:,:,:,None]),axis = 3) 
+  
+  for i,inc in enumerate(results.increments):
+    print('Output step {}/{}'.format(i+1,len(results.increments)))
+
+    header = '1 header\n'
+    
+    data = np.array([inc['inc'] for j in range(np.product(results.grid))]).reshape([np.product(results.grid),1])
+    header+= 'inc'
+    
+    data = np.concatenate((data,np.array([j+1 for j in range(np.product(results.grid))]).reshape([np.product(results.grid),1])),1)
+    header+=' node'
+
+    coords = coords.reshape([np.product(results.grid),3])
+    data = np.concatenate((data,coords),1)
+    header+=' 1_pos 2_pos 3_pos'
+        
+    results.active['increments'] = [inc]
+    for label in options.labels:
+      for o in results.c_output_types:
+        results.active['c_output_types'] = [o]
+        for c in results.constituents:
+          results.active['constituents'] = [c]
+          x = results.get_dataset_location(label)
+          if len(x) == 0:
+            continue
+          label = x[0].split('/')[-1]
+          array = results.read_dataset(x,0)
+          d = np.product(np.shape(array)[1:])
+          array = np.reshape(array,[np.product(results.grid),d])
+          data = np.concatenate((data,array),1)
+          
+          header+= ''.join([' {}_{}'.format(j+1,label) for j in range(d)])
+
+
+    dirname  = os.path.abspath(os.path.join(os.path.dirname(filename),options.dir))
+    try:
+      os.mkdir(dirname)
+    except FileExistsError:
+      pass
+    file_out = '{}_inc{:04d}.txt'.format(filename.split('.')[0],inc['inc'])
+    np.savetxt(os.path.join(dirname,file_out),data,header=header,comments='')

processing/post/DADF5_vtk_cells.py

@@ -21,6 +21,8 @@ parser = argparse.ArgumentParser()
 #parser.add_argument('--version', action='version', version='%(prog)s {}'.format(scriptID))
 parser.add_argument('filenames', nargs='+',
                     help='DADF5 files')
+parser.add_argument('-d','--dir', dest='dir',default='postProc',metavar='string',
+                    help='name of subdirectory to hold output')
 
 options = parser.parse_args()
 
@@ -80,12 +82,17 @@ for filename in options.filenames:
     if results.structured:
       writer = vtk.vtkXMLRectilinearGridWriter()
 
+
+    dirname  = os.path.abspath(os.path.join(os.path.dirname(filename),options.dir))
+    try:
+      os.mkdir(dirname)
+    except FileExistsError:
+      pass
+    file_out = '{}_inc{:04d}.{}'.format(filename.split('.')[0],inc['inc'],writer.GetDefaultFileExtension())
+    
     writer.SetCompressorTypeToZLib()
     writer.SetDataModeToBinary()
-    writer.SetFileName(os.path.join(os.path.split(filename)[0],
+    writer.SetFileName(os.path.join(dirname,file_out))
-                                    os.path.splitext(os.path.split(filename)[1])[0] +
-                                    '_inc{:04d}'.format(i) +                                        # ToDo: adjust to length of increments
-                                    '.' + writer.GetDefaultFileExtension()))
     if results.structured:
       writer.SetInputData(rGrid)
     

python/damask/dadf5.py

@@ -20,14 +20,14 @@ class DADF5():
       
     with h5py.File(filename,'r') as f:
       
-      if f.attrs['DADF5-major'] != 0 or f.attrs['DADF5-minor'] != 1:
+      if f.attrs['DADF5-major'] != 0 or f.attrs['DADF5-minor'] != 2:
         raise TypeError('Unsupported DADF5 version {} '.format(f.attrs['DADF5-version']))
     
-      self.structured = 'grid' in f['mapping'].attrs.keys()
+      self.structured = 'grid' in f['geometry'].attrs.keys()
     
       if self.structured:
-        self.grid = f['mapping'].attrs['grid']
+        self.grid = f['geometry'].attrs['grid']
-        self.size = f['mapping'].attrs['size']
+        self.size = f['geometry'].attrs['size']
         
       r=re.compile('inc[0-9]+')
       self.increments = [{'inc':  int(u[3:]),

python/damask/quaternion.py

@@ -0,0 +1,214 @@
+# -*- coding: UTF-8 no BOM -*-
+
+import numpy as np
+
+P = -1                                                                                              # convention (see DOI:10.1088/0965-0393/23/8/083501)
+
+####################################################################################################
+class Quaternion:
+    u"""
+    Quaternion with basic operations
+    
+    q is the real part, p = (x, y, z) are the imaginary parts.
+    Definition of multiplication depends on variable P, P ∈ {-1,1}.
+    """                                                                              
+
+    def __init__(self,
+                 q    = 0.0,
+                 p    = np.zeros(3,dtype=float)):
+      """Initializes to identity unless specified"""
+      self.q = float(q)
+      self.p = np.array(p,dtype=float)
+
+
+    def __copy__(self):
+      """Copy"""
+      return self.__class__(q=self.q,
+                            p=self.p.copy())
+
+    copy = __copy__
+
+
+    def __iter__(self):
+      """Components"""
+      return iter(self.asList())
+
+    def __repr__(self):
+      """Readable string"""
+      return 'Quaternion: (real={q:+.6f}, imag=<{p[0]:+.6f}, {p[1]:+.6f}, {p[2]:+.6f}>)'.format(q=self.q,p=self.p)
+      
+
+    def __add__(self, other):
+      """Addition"""
+      if isinstance(other, Quaternion):
+        return self.__class__(q=self.q + other.q,
+                              p=self.p + other.p)
+      else:
+        return NotImplemented
+
+    def __iadd__(self, other):
+      """In-place addition"""
+      if isinstance(other, Quaternion):
+        self.q += other.q
+        self.p += other.p
+        return self
+      else:
+        return NotImplemented
+      
+    def __pos__(self):                                                                      
+      """Unary positive operator"""                                                                       
+      return self
+      
+      
+    def __sub__(self, other):
+      """Subtraction"""
+      if isinstance(other, Quaternion):
+        return self.__class__(q=self.q - other.q,
+                              p=self.p - other.p)
+      else:
+        return NotImplemented
+
+    def __isub__(self, other):
+      """In-place subtraction"""
+      if isinstance(other, Quaternion):
+        self.q -= other.q
+        self.p -= other.p
+        return self
+      else:
+        return NotImplemented
+        
+    def __neg__(self):                                                                      
+      """Unary negative operator"""  
+      return self * -1.0
+      
+      
+    def __mul__(self, other):
+      """Multiplication with quaternion or scalar"""
+      if isinstance(other, Quaternion):
+        return self.__class__(q=self.q*other.q - np.dot(self.p,other.p),
+                              p=self.q*other.p + other.q*self.p + P * np.cross(self.p,other.p))
+      elif isinstance(other, (int, float)):
+        return self.__class__(q=self.q*other,
+                              p=self.p*other)
+      else:
+        return NotImplemented
+
+    def __imul__(self, other):
+      """In-place multiplication with quaternion or scalar"""
+      if isinstance(other, Quaternion):
+        self.q = self.q*other.q - np.dot(self.p,other.p)
+        self.p = self.q*other.p + other.q*self.p + P * np.cross(self.p,other.p)
+        return self
+      elif isinstance(other, (int, float)):
+        self.q *= other
+        self.p *= other
+        return self
+      else:
+        return NotImplemented
+
+
+    def __truediv__(self, other):
+      """Divsion with quaternion or scalar"""
+      if isinstance(other, Quaternion):
+        s = other.conjugate()/abs(other)**2.
+        return self.__class__(q=self.q * s,
+                              p=self.p * s)
+      elif isinstance(other, (int, float)):
+        return self.__class__(q=self.q / other,
+                              p=self.p / other)
+      else:
+        return NotImplemented
+
+    def __itruediv__(self, other):
+      """In-place divsion with quaternion or scalar"""
+      if isinstance(other, Quaternion):
+        s = other.conjugate()/abs(other)**2.
+        self *= s
+        return self
+      elif isinstance(other, (int, float)):
+        self.q /= other
+        self.p /= other
+        return self
+      else:
+        return NotImplemented
+
+
+    def __pow__(self, exponent):
+      """Power"""
+      if isinstance(exponent, (int, float)):
+        omega = np.acos(self.q)
+        return self.__class__(q=         np.cos(exponent*omega),
+                              p=self.p * np.sin(exponent*omega)/np.sin(omega))
+      else:
+        return NotImplemented
+
+    def __ipow__(self, exponent):
+      """In-place power"""
+      if isinstance(exponent, (int, float)):
+        omega = np.acos(self.q)
+        self.q  = np.cos(exponent*omega)
+        self.p *= np.sin(exponent*omega)/np.sin(omega)
+        return self
+      else:
+        return NotImplemented
+      
+
+    def __abs__(self):
+      """Norm"""
+      return np.sqrt(self.q ** 2 + np.dot(self.p,self.p))
+
+    magnitude = __abs__
+
+
+    def __eq__(self,other):
+      """Equal (sufficiently close) to each other"""
+      return np.isclose(( self-other).magnitude(),0.0) \
+          or np.isclose((-self-other).magnitude(),0.0)
+
+    def __ne__(self,other):
+      """Not equal (sufficiently close) to each other"""
+      return not self.__eq__(other)
+
+
+    def asM(self):
+      """Intermediate representation useful for quaternion averaging (see F. Landis Markley et al.)"""
+      return np.outer(self.asArray(),self.asArray())
+
+    def asArray(self):
+      """As numpy array"""
+      return np.array((self.q,self.p[0],self.p[1],self.p[2]))
+ 
+    def asList(self):
+      """As list"""
+      return [self.q]+list(self.p)
+
+
+    def normalize(self):
+      """Normalizes in-place"""
+      d = self.magnitude()
+      if d > 0.0: self /= d
+      return self
+      
+    def normalized(self):
+      """Returns normalized copy"""
+      return self.copy().normalize()
+
+
+    def conjugate(self):
+      """Conjugates in-place"""
+      self.p = -self.p
+      return self
+      
+    def conjugated(self):
+      """Returns conjugated copy"""
+      return self.copy().conjugate()
+
+
+    def homomorph(self):
+      """Homomorphs in-place"""
+      self *= -1.0
+      return self
+      
+    def homomorphed(self):
+      """Returns homomorphed copy"""
+      return self.copy().homomorph()

src/CPFEM.f90

@@ -47,7 +47,8 @@ module CPFEM
 
  public :: &
    CPFEM_general, &
-   CPFEM_initAll
+   CPFEM_initAll, &
+   CPFEM_results
 
 contains
 
@@ -633,4 +634,35 @@ subroutine CPFEM_general(mode, parallelExecution, ffn, ffn1, temperature_inp, dt
 
 end subroutine CPFEM_general
 
+
+!--------------------------------------------------------------------------------------------------
+!> @brief triggers writing of the results
+!--------------------------------------------------------------------------------------------------
+subroutine CPFEM_results(inc,time)
+ use prec, only: &
+   pInt
+#ifdef DAMASK_HDF5
+ use results
+ use HDF5_utilities
+#endif
+ use constitutive, only: &
+   constitutive_results
+ use crystallite, only: &
+   crystallite_results
+
+ implicit none
+ integer(pInt), intent(in) :: inc
+ real(pReal),   intent(in) :: time
+
+#ifdef DAMASK_HDF5
+ call results_openJobFile
+ call results_addIncrement(inc,time)
+ call constitutive_results
+ call crystallite_results
+ call results_removeLink('current') ! ToDo: put this into closeJobFile
+ call results_closeJobFile
+#endif
+
+end subroutine CPFEM_results
+
 end module CPFEM

src/CPFEM2.f90

@@ -257,7 +257,7 @@ subroutine CPFEM_age()
      write(6,'(a)') '<< CPFEM >> writing restart variables of last converged step to hdf5 file'
    
    write(rankStr,'(a1,i0)')'_',worldrank
-   fileHandle = HDF5_openFile(trim(getSolverJobName())//trim(rankStr)//'.hdf5','w')
+   fileHandle = HDF5_openFile(trim(getSolverJobName())//trim(rankStr)//'.hdf5','a')
    
    call HDF5_write(fileHandle,material_phase,      'recordedPhase')
    call HDF5_write(fileHandle,crystallite_F0,      'convergedF')
@@ -290,6 +290,7 @@ subroutine CPFEM_age()
 
 end subroutine CPFEM_age
 
+
 !--------------------------------------------------------------------------------------------------
 !> @brief triggers writing of the results
 !--------------------------------------------------------------------------------------------------

src/DAMASK_interface.f90

@@ -45,7 +45,7 @@ contains
 !> @brief initializes the solver by interpreting the command line arguments. Also writes
 !! information on computation to screen
 !--------------------------------------------------------------------------------------------------
-subroutine DAMASK_interface_init()
+subroutine DAMASK_interface_init
   use, intrinsic :: &
     iso_fortran_env
   use, intrinsic :: &
@@ -130,11 +130,11 @@ subroutine DAMASK_interface_init()
   call MPI_Init_Thread(MPI_THREAD_FUNNELED,threadLevel,mpi_err)
   if (mpi_err /= 0) call quit(1)
   if (threadLevel<MPI_THREAD_FUNNELED) then
-    write(6,'(a)') ' MPI library does not support OpenMP'
+    write(6,'(/,a)') ' ERROR: MPI library does not support OpenMP'
     call quit(1)
   endif
 #endif
-  call PETScInitialize(PETSC_NULL_CHARACTER,petsc_err)                                              ! according to PETSc manual, that should be the first line in the code
+  call PETScInitializeNoArguments(petsc_err)                                                        ! according to PETSc manual, that should be the first line in the code
   CHKERRQ(petsc_err)                                                                                ! this is a macro definition, it is case sensitive
 
   call MPI_Comm_rank(PETSC_COMM_WORLD,worldrank,mpi_err)
@@ -144,11 +144,11 @@ subroutine DAMASK_interface_init()
 
   mainProcess: if (worldrank == 0) then
     if (output_unit /= 6) then
-      write(output_unit,'(a)') ' STDOUT != 6'
+      write(output_unit,'(/,a)') ' ERROR: STDOUT != 6'
       call quit(1)
     endif
     if (error_unit /= 0) then
-      write(output_unit,'(a)') ' STDERR != 0'
+      write(output_unit,'(/,a)') ' ERROR: STDERR != 0'
       call quit(1)
     endif
   else mainProcess
@@ -254,14 +254,14 @@ subroutine DAMASK_interface_init()
         call get_command_argument(i+1,arg)
         read(arg,*,iostat=stat) interface_restartInc
         if (interface_restartInc < 0 .or. stat /=0) then
-          write(6,'(a)') ' Could not parse restart increment: '//trim(arg)
+          write(6,'(/,a)') ' ERROR: Could not parse restart increment: '//trim(arg)
           call quit(1)
         endif
     end select
   enddo
 
   if (len_trim(loadcaseArg) == 0 .or. len_trim(geometryArg) == 0) then
-    write(6,'(a)') ' Please specify geometry AND load case (-h for help)'
+    write(6,'(/,a)') ' ERROR: Please specify geometry AND load case (-h for help)'
     call quit(1)
   endif
 
@@ -324,7 +324,7 @@ subroutine setWorkingDirectory(workingDirectoryArg)
   workingDirectory = trim(rectifyPath(workingDirectory))
   error = setCWD(trim(workingDirectory))
   if(error) then
-    write(6,'(a20,a,a16)') ' Working directory "',trim(workingDirectory),'" does not exist'
+    write(6,'(/,a)') ' ERROR: Working directory "'//trim(workingDirectory)//'" does not exist'
     call quit(1)
   endif
 
@@ -374,7 +374,7 @@ character(len=1024) function getGeometryFile(geometryParameter)
 
   inquire(file=trim(getGeometryFile), exist=file_exists)
   if (.not. file_exists) then
-    write(6,'(a)') ' Geometry file does not exists ('//trim(getGeometryFile)//')'
+    write(6,'(/,a)') ' ERROR: Geometry file does not exists ('//trim(getGeometryFile)//')'
     call quit(1)
   endif
 
@@ -399,7 +399,7 @@ character(len=1024) function getLoadCaseFile(loadCaseParameter)
 
   inquire(file=trim(getLoadCaseFile), exist=file_exists)
   if (.not. file_exists) then
-    write(6,'(a)') ' Load case file does not exists ('//trim(getLoadCaseFile)//')'
+    write(6,'(/,a)') ' ERROR: Load case file does not exists ('//trim(getLoadCaseFile)//')'
     call quit(1)
   endif
 

src/DAMASK_marc.f90

@@ -16,6 +16,7 @@
 !> @details  Marc subroutines used:
 !> @details   - hypela2
 !> @details   - plotv
+!> @details   - uedinc
 !> @details   - flux
 !> @details   - quit
 !> @details  Marc common blocks included:
@@ -273,26 +274,20 @@ subroutine hypela2(d,g,e,de,s,t,dt,ngens,m,nn,kcus,matus,ndi,nshear,disp, &
        outdatedByNewInc = .false.                                                                   ! no aging of state
        calcMode = .false.                                                                           ! pretend last step was collection
        lastLovl = lovl                                                                              ! pretend that this is NOT the first after a lovl change
-       !$OMP CRITICAL (write2out)
+       write(6,'(a,i6,1x,i2)') '<< HYPELA2 >> start of analysis..! ',m(1),nn
-         write(6,'(a,i6,1x,i2)') '<< HYPELA2 >> start of analysis..! ',m(1),nn
+       flush(6)
-         flush(6)
-       !$OMP END CRITICAL (write2out)
      else if (inc - theInc > 1) then                                                                ! >> restart of broken analysis <<
        lastIncConverged = .false.                                                                   ! no Jacobian backup
        outdatedByNewInc = .false.                                                                   ! no aging of state
        calcMode = .true.                                                                            ! pretend last step was calculation
-       !$OMP CRITICAL (write2out)
+       write(6,'(a,i6,1x,i2)') '<< HYPELA2 >> restart of analysis..! ',m(1),nn
-         write(6,'(a,i6,1x,i2)') '<< HYPELA2 >> restart of analysis..! ',m(1),nn
+       flush(6)
-         flush(6)
-       !$OMP END CRITICAL (write2out)
      else                                                                                           ! >> just the next inc <<
        lastIncConverged = .true.                                                                    ! request Jacobian backup
        outdatedByNewInc = .true.                                                                    ! request aging of state
        calcMode = .true.                                                                            ! assure last step was calculation
-       !$OMP CRITICAL (write2out)
+       write(6,'(a,i6,1x,i2)') '<< HYPELA2 >> new increment..! ',m(1),nn
-         write(6,'(a,i6,1x,i2)') '<< HYPELA2 >> new increment..! ',m(1),nn
+       flush(6)
-         flush(6)
-       !$OMP END CRITICAL (write2out)
      endif
    else if ( timinc < theDelta ) then                                                               ! >> cutBack <<
      lastIncConverged = .false.                                                                     ! no Jacobian backup
@@ -300,10 +295,8 @@ subroutine hypela2(d,g,e,de,s,t,dt,ngens,m,nn,kcus,matus,ndi,nshear,disp, &
      terminallyIll = .false.
      cycleCounter = -1                                                                              ! first calc step increments this to cycle = 0
      calcMode = .true.                                                                              ! pretend last step was calculation
-     !$OMP CRITICAL (write2out)
+     write(6,'(a,i6,1x,i2)') '<< HYPELA2 >> cutback detected..! ',m(1),nn
-       write(6,'(a,i6,1x,i2)') '<< HYPELA2 >> cutback detected..! ',m(1),nn
+     flush(6)
-       flush(6)
-     !$OMP END CRITICAL (write2out)
    endif                                                                                            ! convergence treatment end
 
 
@@ -365,7 +358,6 @@ subroutine hypela2(d,g,e,de,s,t,dt,ngens,m,nn,kcus,matus,ndi,nshear,disp, &
  lastLovl = lovl                                                                                    ! record lovl
 
  call CPFEM_general(computationMode,usePingPong,ffn,ffn1,t(1),timinc,m(1),nn,stress,ddsdde)
-
 !     Mandel: 11, 22, 33, SQRT(2)*12, SQRT(2)*23, SQRT(2)*13
 !     Marc:   11, 22, 33, 12, 23, 13
 !     Marc:   11, 22, 33, 12
@@ -407,6 +399,26 @@ subroutine flux(f,ts,n,time)
  end subroutine flux
 
 
+!--------------------------------------------------------------------------------------------------
+!> @brief sets user defined output variables for Marc
+!> @details select a variable contour plotting (user subroutine).
+!--------------------------------------------------------------------------------------------------
+subroutine uedinc(inc,incsub)
+  use prec, only: &
+    pReal, &
+    pInt
+  use CPFEM, only: &
+    CPFEM_results
+
+  implicit none
+  integer, intent(in) :: inc, incsub
+#include QUOTE(PASTE(./MarcInclude/creeps,Marc4DAMASK))                                             ! creeps is needed for timinc (time increment)
+
+  call CPFEM_results(inc,cptim)
+
+end subroutine uedinc
+
+
 !--------------------------------------------------------------------------------------------------
 !> @brief sets user defined output variables for Marc
 !> @details select a variable contour plotting (user subroutine).

src/grid/DAMASK_grid.f90

@@ -77,6 +77,7 @@ program DAMASK_spectral
  use grid_mech_FEM
  use grid_damage_spectral
  use grid_thermal_spectral
+ use HDF5_utilities
  use results
 
  implicit none
@@ -155,8 +156,6 @@ program DAMASK_spectral
  write(6,'(/,a)') ' Shanthraj et al., Handbook of Mechanics of Materials, 2019'
  write(6,'(a)')   ' https://doi.org/10.1007/978-981-10-6855-3_80'
 
- call results_openJobFile()
- call results_closeJobFile()
 !--------------------------------------------------------------------------------------------------
 ! initialize field solver information
  nActiveFields = 1
@@ -301,7 +300,7 @@ program DAMASK_spectral
 
    reportAndCheck: if (worldrank == 0) then
      write (loadcase_string, '(i6)' ) currentLoadCase
-     write(6,'(1x,a,i6)') 'load case: ', currentLoadCase
+     write(6,'(/,1x,a,i6)') 'load case: ', currentLoadCase
      if (.not. newLoadCase%followFormerTrajectory) write(6,'(2x,a)') 'drop guessing along trajectory'
      if (newLoadCase%deformation%myType == 'l') then
        do j = 1_pInt, 3_pInt
@@ -348,10 +347,10 @@ program DAMASK_spectral
      if (newLoadCase%time < 0.0_pReal) errorID = 834_pInt                                           ! negative time increment
      write(6,'(2x,a,f12.6)') 'time:       ', newLoadCase%time
      if (newLoadCase%incs < 1_pInt)    errorID = 835_pInt                                           ! non-positive incs count
-     write(6,'(2x,a,i5)')    'increments: ', newLoadCase%incs
+     write(6,'(2x,a,i5)')  'increments: ', newLoadCase%incs
      if (newLoadCase%outputfrequency < 1_pInt)  errorID = 836_pInt                                  ! non-positive result frequency
-     write(6,'(2x,a,i5)')    'output  frequency:  ', newLoadCase%outputfrequency
+     write(6,'(2x,a,i5)')  'output  frequency:  ', newLoadCase%outputfrequency
-     write(6,'(2x,a,i5,/)')  'restart frequency:  ', newLoadCase%restartfrequency
+     write(6,'(2x,a,i5)')  'restart frequency:  ', newLoadCase%restartfrequency
      if (errorID > 0_pInt) call IO_error(error_ID = errorID, ext_msg = loadcase_string)             ! exit with error message
    endif reportAndCheck
    loadCases = [loadCases,newLoadCase]                                                              ! load case is ok, append it
@@ -359,8 +358,9 @@ program DAMASK_spectral
  close(fileUnit)
 
  call results_openJobFile
- call results_addAttribute('grid',grid,'mapping')
+ call HDF5_closeGroup(results_addGroup('geometry'))
- call results_addAttribute('size',geomSize,'mapping')
+ call results_addAttribute('grid',grid,'geometry')
+ call results_addAttribute('size',geomSize,'geometry')
  call results_closeJobFile
 
 !--------------------------------------------------------------------------------------------------

src/grid/grid_damage_spectral.f90

@@ -64,7 +64,6 @@ subroutine grid_damage_spectral_init
     worldsize, &
     petsc_options
     
-  implicit none
   PetscInt, dimension(worldsize) :: localK  
   integer :: i, j, k, cell
   DM :: damage_grid
@@ -164,7 +163,6 @@ function grid_damage_spectral_solution(timeinc,timeinc_old,loadCaseTime) result(
   use damage_nonlocal, only: &
     damage_nonlocal_putNonLocalDamage
  
-  implicit none
   real(pReal), intent(in) :: &
     timeinc, &                                                                                      !< increment in time for current solution
     timeinc_old, &                                                                                  !< increment in time of last increment
@@ -236,7 +234,6 @@ subroutine grid_damage_spectral_forward
     damage_nonlocal_getDiffusion33, &
     damage_nonlocal_getMobility
     
-  implicit none
   integer                               :: i, j, k, cell
   DM :: dm_local
   PetscScalar,  dimension(:,:,:), pointer     :: x_scal
@@ -301,7 +298,6 @@ subroutine formResidual(in,x_scal,f_scal,dummy,ierr)
     damage_nonlocal_getDiffusion33, &
     damage_nonlocal_getMobility
  
-  implicit none
   DMDALocalInfo, dimension(DMDA_LOCAL_INFO_SIZE) :: &
     in
   PetscScalar, dimension( &

src/grid/grid_mech_spectral_basic.f90

@@ -7,6 +7,8 @@
 module grid_mech_spectral_basic
 #include <petsc/finclude/petscsnes.h>
 #include <petsc/finclude/petscdmda.h>
+  use DAMASK_interface
+  use HDF5_utilities
   use PETScdmda
   use PETScsnes
   use prec, only: &
@@ -25,8 +27,7 @@ module grid_mech_spectral_basic
   type(tSolutionParams), private :: params
   
   type, private :: tNumerics
-    logical :: &
+    logical :: update_gamma                                                                         !< update gamma operator with current stiffness
-      update_gamma !< update gamma operator with current stiffness
   end type tNumerics
   
   type(tNumerics) :: num                                                                            ! numerics parameters. Better name?
@@ -40,8 +41,8 @@ module grid_mech_spectral_basic
 !--------------------------------------------------------------------------------------------------
 ! common pointwise data
   real(pReal), private, dimension(:,:,:,:,:), allocatable ::  &
-    F_lastInc, &
+    F_lastInc, &                                                                                    !< field of previous compatible deformation gradients
-    Fdot
+    Fdot                                                                                            !< field of assumed rate of compatible deformation gradient
 
 !--------------------------------------------------------------------------------------------------
 ! stress, stiffness and compliance average etc.
@@ -99,23 +100,22 @@ subroutine grid_mech_spectral_basic_init
   use spectral_utilities, only: &
     utilities_constitutiveResponse, &
     utilities_updateGamma, &
-    utilities_updateIPcoords, &
+    utilities_updateIPcoords
-    wgt
   use mesh, only: &
     grid, &
     grid3
   use math, only: &
     math_invSym3333
     
-  implicit none
   real(pReal), dimension(3,3,grid(1),grid(2),grid3) :: P
   real(pReal), dimension(3,3) :: &
     temp33_Real = 0.0_pReal
  
   PetscErrorCode :: ierr
   PetscScalar, pointer, dimension(:,:,:,:) :: &
-   F                                                                                                ! pointer to solution data
+    F                                                                                               ! pointer to solution data
   PetscInt, dimension(worldsize) :: localK  
+  integer(HID_T) :: fileHandle
   integer :: fileUnit
   character(len=1024) :: rankStr
  
@@ -163,7 +163,7 @@ subroutine grid_mech_spectral_basic_init
   call DMcreateGlobalVector(da,solution_vec,ierr); CHKERRQ(ierr)                                    ! global solution vector (grid x 9, i.e. every def grad tensor)
   call DMDASNESsetFunctionLocal(da,INSERT_VALUES,formResidual,PETSC_NULL_SNES,ierr)                 ! residual vector of same shape as solution vector
   CHKERRQ(ierr) 
-  call SNESsetConvergenceTest(snes,converged,PETSC_NULL_SNES,PETSC_NULL_FUNCTION,ierr)! specify custom convergence check function "converged"
+  call SNESsetConvergenceTest(snes,converged,PETSC_NULL_SNES,PETSC_NULL_FUNCTION,ierr)              ! specify custom convergence check function "converged"
   CHKERRQ(ierr)
   call SNESsetFromOptions(snes,ierr); CHKERRQ(ierr)                                                 ! pull it all together with additional CLI arguments
 
@@ -174,19 +174,14 @@ subroutine grid_mech_spectral_basic_init
   restart: if (restartInc > 0) then                                                     
     write(6,'(/,a,'//IO_intOut(restartInc)//',a)') ' reading values of increment ', restartInc, ' from file'
 
-    fileUnit = IO_open_jobFile_binary('F_aim')
-    read(fileUnit) F_aim; close(fileUnit)
-    fileUnit = IO_open_jobFile_binary('F_aim_lastInc')
-    read(fileUnit) F_aim_lastInc; close(fileUnit)
-    fileUnit = IO_open_jobFile_binary('F_aimDot')
-    read(fileUnit) F_aimDot; close(fileUnit)
- 
     write(rankStr,'(a1,i0)')'_',worldrank
+    fileHandle = HDF5_openFile(trim(getSolverJobName())//trim(rankStr)//'.hdf5')
  
-    fileUnit = IO_open_jobFile_binary('F'//trim(rankStr))
+    call HDF5_read(fileHandle,F_aim,        'F_aim')
-    read(fileUnit) F; close (fileUnit)
+    call HDF5_read(fileHandle,F_aim_lastInc,'F_aim_lastInc')
-    fileUnit = IO_open_jobFile_binary('F_lastInc'//trim(rankStr))
+    call HDF5_read(fileHandle,F_aimDot,     'F_aimDot')
-    read(fileUnit) F_lastInc; close (fileUnit)
+    call HDF5_read(fileHandle,F,            'F')
+    call HDF5_read(fileHandle,F_lastInc,    'F_lastInc')
  
   elseif (restartInc == 0) then restart
     F_lastInc = spread(spread(spread(math_I3,3,grid(1)),4,grid(2)),5,grid3)                         ! initialize to identity
@@ -203,15 +198,15 @@ subroutine grid_mech_spectral_basic_init
  
   restartRead: if (restartInc > 0) then
     write(6,'(/,a,'//IO_intOut(restartInc)//',a)') 'reading more values of increment ', restartInc, ' from file'
-    fileUnit = IO_open_jobFile_binary('C_volAvg')
+    call HDF5_read(fileHandle,C_volAvg,       'C_volAvg')
-    read(fileUnit) C_volAvg; close(fileUnit)
+    call HDF5_read(fileHandle,C_volAvgLastInc,'C_volAvgLastInc')
-    fileUnit = IO_open_jobFile_binary('C_volAvgLastInv')
+    call HDF5_closeFile(fileHandle)
-    read(fileUnit) C_volAvgLastInc; close(fileUnit)
+
     fileUnit = IO_open_jobFile_binary('C_ref')
     read(fileUnit) C_minMaxAvg; close(fileUnit)
   endif restartRead
 
- call Utilities_updateGamma(C_minMaxAvg,.true.)
+  call utilities_updateGamma(C_minMaxAvg,.true.)
 
 end subroutine grid_mech_spectral_basic_init
 
@@ -228,8 +223,6 @@ function grid_mech_spectral_basic_solution(incInfoIn,timeinc,timeinc_old,stress_
     restartWrite, &
     terminallyIll
  
-  implicit none
-
 !--------------------------------------------------------------------------------------------------
 ! input data for solution
   character(len=*),            intent(in) :: &
@@ -251,8 +244,8 @@ function grid_mech_spectral_basic_solution(incInfoIn,timeinc,timeinc_old,stress_
 
 !--------------------------------------------------------------------------------------------------
 ! update stiffness (and gamma operator)
-  S = Utilities_maskedCompliance(rotation_BC,stress_BC%maskLogical,C_volAvg)
+  S = utilities_maskedCompliance(rotation_BC,stress_BC%maskLogical,C_volAvg)
-  if (num%update_gamma) call Utilities_updateGamma(C_minMaxAvg,restartWrite)
+  if (num%update_gamma) call utilities_updateGamma(C_minMaxAvg,restartWrite)
  
 !--------------------------------------------------------------------------------------------------
 ! set module wide available data 
@@ -306,7 +299,6 @@ subroutine grid_mech_spectral_basic_forward(guess,timeinc,timeinc_old,loadCaseTi
   use FEsolving, only: &
     restartWrite
 
-  implicit none
   logical,                     intent(in) :: &
     guess
   real(pReal),                 intent(in) :: &
@@ -321,7 +313,7 @@ subroutine grid_mech_spectral_basic_forward(guess,timeinc,timeinc_old,loadCaseTi
   PetscErrorCode :: ierr
   PetscScalar, dimension(:,:,:,:), pointer :: F
   
-  integer :: fileUnit
+  integer(HID_T) :: fileHandle
   character(len=32) :: rankStr
 
   call DMDAVecGetArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr)
@@ -331,29 +323,24 @@ subroutine grid_mech_spectral_basic_forward(guess,timeinc,timeinc_old,loadCaseTi
     C_minMaxAvg = C_minMaxAvgLastInc                                                               ! QUESTION: where is this required?
   else
   !--------------------------------------------------------------------------------------------------
-    ! restart information for spectral solver
+  ! restart information for spectral solver
-    if (restartWrite) then                                                                           ! QUESTION: where is this logical properly set?
+    if (restartWrite) then
-      write(6,'(/,a)') ' writing converged results for restart'
+      write(6,'(/,a)') ' writing converged results for restart';flush(6)
-      flush(6)
+      
-
-      if (worldrank == 0) then
-        fileUnit = IO_open_jobFile_binary('C_volAvg','w')
-        write(fileUnit) C_volAvg; close(fileUnit)
-        fileUnit = IO_open_jobFile_binary('C_volAvgLastInv','w')
-        write(fileUnit) C_volAvgLastInc; close(fileUnit)
-        fileUnit = IO_open_jobFile_binary('F_aim','w')
-        write(fileUnit) F_aim; close(fileUnit)
-        fileUnit = IO_open_jobFile_binary('F_aim_lastInc','w')
-        write(fileUnit) F_aim_lastInc; close(fileUnit)
-        fileUnit = IO_open_jobFile_binary('F_aimDot','w')
-        write(fileUnit) F_aimDot; close(fileUnit)
-      endif
-
       write(rankStr,'(a1,i0)')'_',worldrank
-      fileUnit = IO_open_jobFile_binary('F'//trim(rankStr),'w')
+      fileHandle = HDF5_openFile(trim(getSolverJobName())//trim(rankStr)//'.hdf5','w')
-      write(fileUnit) F; close (fileUnit)
+      
-      fileUnit = IO_open_jobFile_binary('F_lastInc'//trim(rankStr),'w')
+      call HDF5_write(fileHandle,F_aim,        'F_aim')
-      write(fileUnit) F_lastInc; close (fileUnit)
+      call HDF5_write(fileHandle,F_aim_lastInc,'F_aim_lastInc')
+      call HDF5_write(fileHandle,F_aimDot,     'F_aimDot')
+      call HDF5_write(fileHandle,F,            'F')
+      call HDF5_write(fileHandle,F_lastInc,    'F_lastInc')
+
+      call HDF5_write(fileHandle,C_volAvg,       'C_volAvg')
+      call HDF5_write(fileHandle,C_volAvgLastInc,'C_volAvgLastInc')
+      call HDF5_write(fileHandle,C_minMaxAvg,    'C_minMaxAvg')
+
+      call HDF5_closeFile(fileHandle)
     endif
 
     call CPFEM_age                                                                                  ! age state and kinematics
@@ -399,7 +386,7 @@ end subroutine grid_mech_spectral_basic_forward
 !--------------------------------------------------------------------------------------------------
 !> @brief convergence check
 !--------------------------------------------------------------------------------------------------
-subroutine converged(snes_local,PETScIter,xnorm,snorm,fnorm,reason,dummy,ierr)
+subroutine converged(snes_local,PETScIter,devNull1,devNull2,devNull3,reason,dummy,ierr)
   use numerics, only: &
     itmax, &
     itmin, &
@@ -410,13 +397,12 @@ subroutine converged(snes_local,PETScIter,xnorm,snorm,fnorm,reason,dummy,ierr)
   use FEsolving, only: &
     terminallyIll
 
-  implicit none
   SNES :: snes_local
-  PetscInt :: PETScIter
+  PetscInt,  intent(in) :: PETScIter
-  PetscReal :: &
+  PetscReal, intent(in) :: &
-    xnorm, &                                                                                        ! not used
+    devNull1, &
-    snorm, &                                                                                        ! not used
+    devNull2, &
-    fnorm                                                                                           ! not used
+    devNull3 
   SNESConvergedReason :: reason
   PetscObject :: dummy
   PetscErrorCode :: ierr
@@ -452,7 +438,7 @@ end subroutine converged
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief forms the basic residual vector
+!> @brief forms the residual vector
 !--------------------------------------------------------------------------------------------------
 subroutine formResidual(in, F, &
                         residuum, dummy, ierr)
@@ -481,7 +467,6 @@ subroutine formResidual(in, F, &
   use FEsolving, only: &
     terminallyIll
 
-  implicit none
   DMDALocalInfo, dimension(DMDA_LOCAL_INFO_SIZE) :: in                                              !< DMDA info (needs to be named "in" for macros like XRANGE to work)
   PetscScalar, dimension(3,3,XG_RANGE,YG_RANGE,ZG_RANGE), &
     intent(in) :: F                                                                                 !< deformation gradient field
@@ -515,7 +500,7 @@ subroutine formResidual(in, F, &
 
 !--------------------------------------------------------------------------------------------------
 ! evaluate constitutive response
-  call Utilities_constitutiveResponse(residuum, &                                                   ! "residuum" gets field of first PK stress (to save memory)
+  call utilities_constitutiveResponse(residuum, &                                                   ! "residuum" gets field of first PK stress (to save memory)
                                       P_av,C_volAvg,C_minMaxAvg, &
                                       F,params%timeinc,params%rotation_BC)
   call MPI_Allreduce(MPI_IN_PLACE,terminallyIll,1,MPI_LOGICAL,MPI_LOR,PETSC_COMM_WORLD,ierr)

src/grid/grid_thermal_spectral.f90

@@ -69,7 +69,6 @@ subroutine grid_thermal_spectral_init
     worldsize, &
     petsc_options
 
-  implicit none
   PetscInt, dimension(worldsize) :: localK  
   integer :: i, j, k, cell
   DM :: thermal_grid
@@ -167,7 +166,6 @@ function grid_thermal_spectral_solution(timeinc,timeinc_old,loadCaseTime) result
   use thermal_conduction, only: &
     thermal_conduction_putTemperatureAndItsRate
  
-  implicit none
   real(pReal), intent(in) :: &
     timeinc, &                                                                                      !< increment in time for current solution
     timeinc_old, &                                                                                  !< increment in time of last increment
@@ -242,7 +240,6 @@ subroutine grid_thermal_spectral_forward
     thermal_conduction_getMassDensity, &
     thermal_conduction_getSpecificHeat
     
-  implicit none
   integer :: i, j, k, cell
   DM :: dm_local
   PetscScalar,  dimension(:,:,:), pointer :: x_scal
@@ -311,7 +308,6 @@ subroutine formResidual(in,x_scal,f_scal,dummy,ierr)
     thermal_conduction_getMassDensity, &
     thermal_conduction_getSpecificHeat
  
-  implicit none
   DMDALocalInfo, dimension(DMDA_LOCAL_INFO_SIZE) :: &
     in
   PetscScalar, dimension( &

src/material.f90

@@ -183,14 +183,11 @@ module material
 
  integer(pInt), private :: &
    microstructure_maxNconstituents, &                                                               !< max number of constituents in any phase
-   texture_maxNgauss, &                                                                             !< max number of Gauss components in any texture
+   texture_maxNgauss                                                                                !< max number of Gauss components in any texture
-   texture_maxNfiber                                                                                !< max number of Fiber components in any texture
 
  integer(pInt), dimension(:), allocatable, private :: &
    microstructure_Nconstituents, &                                                                  !< number of constituents in each microstructure
-   texture_symmetry, &                                                                              !< number of symmetric orientations per texture
+   texture_Ngauss                                                                                   !< number of Gauss components per texture
-   texture_Ngauss, &                                                                                !< number of Gauss components per texture
-   texture_Nfiber                                                                                   !< number of Fiber components per texture
 
  integer(pInt), dimension(:,:), allocatable, private :: &
    microstructure_phase, &                                                                          !< phase IDs of each microstructure
@@ -200,9 +197,7 @@ module material
    microstructure_fraction                                                                          !< vol fraction of each constituent in microstructure
 
  real(pReal), dimension(:,:,:), allocatable, private :: &
-   material_volume, &                                                                               !< volume of each grain,IP,element
    texture_Gauss, &                                                                                 !< data of each Gauss component
-   texture_Fiber, &                                                                                 !< data of each Fiber component
    texture_transformation                                                                           !< transformation for each texture
 
  logical, dimension(:), allocatable, private :: &
@@ -807,31 +802,27 @@ subroutine material_parseTexture
    math_det33
 
  implicit none
- integer(pInt) :: section, gauss, fiber, j, t, i
+ integer(pInt) :: section, gauss, j, t, i
  character(len=65536), dimension(:), allocatable ::  strings                                     ! Values for given key in material config 
  integer(pInt), dimension(:), allocatable :: chunkPos
 
- allocate(texture_symmetry(size(config_texture)), source=1_pInt)
  allocate(texture_Ngauss(size(config_texture)),   source=0_pInt)
- allocate(texture_Nfiber(size(config_texture)),   source=0_pInt)
 
  do t=1_pInt, size(config_texture)
-   texture_Ngauss(t) =  config_texture(t)%countKeys('(gauss)') &
+   texture_Ngauss(t) =  config_texture(t)%countKeys('(gauss)')
-                     +  config_texture(t)%countKeys('(random)')
+   if (config_texture(t)%keyExists('symmetry')) call IO_error(147,ext_msg='symmetry')
-   texture_Nfiber(t) =  config_texture(t)%countKeys('(fiber)')
+   if (config_texture(t)%keyExists('(random)')) call IO_error(147,ext_msg='(random)')
+   if (config_texture(t)%keyExists('(fiber)'))  call IO_error(147,ext_msg='(fiber)')
  enddo
 
  texture_maxNgauss = maxval(texture_Ngauss)
- texture_maxNfiber = maxval(texture_Nfiber)
  allocate(texture_Gauss (5,texture_maxNgauss,size(config_texture)), source=0.0_pReal)
- allocate(texture_Fiber (6,texture_maxNfiber,size(config_texture)), source=0.0_pReal)
  allocate(texture_transformation(3,3,size(config_texture)),         source=0.0_pReal)
           texture_transformation = spread(math_I3,3,size(config_texture))
 
  do t=1_pInt, size(config_texture)
    section = t
    gauss = 0_pInt
-   fiber = 0_pInt
    
    if (config_texture(t)%keyExists('axes')) then
      strings = config_texture(t)%getStrings('axes')
@@ -856,10 +847,6 @@ subroutine material_parseTexture
      if(dNeq(math_det33(texture_transformation(1:3,1:3,t)),1.0_pReal)) call IO_error(157_pInt,t)
    endif
 
-   if (config_texture(t)%keyExists('symmetry')) call IO_error(147,ext_msg='symmetry')
-   if (config_texture(t)%keyExists('(random)')) call IO_error(147,ext_msg='(random)')
-   if (config_texture(t)%keyExists('(fiber)'))  call IO_error(147,ext_msg='(fiber)')
-   
    if (config_texture(t)%keyExists('(gauss)')) then
      gauss = gauss + 1_pInt
      strings = config_texture(t)%getStrings('(gauss)',raw= .true.)
@@ -873,10 +860,6 @@ subroutine material_parseTexture
                  texture_Gauss(2,gauss,t) = IO_floatValue(strings(i),chunkPos,j+1_pInt)*inRad
              case('phi2')
                  texture_Gauss(3,gauss,t) = IO_floatValue(strings(i),chunkPos,j+1_pInt)*inRad
-             case('scatter')
-                 texture_Gauss(4,gauss,t) = IO_floatValue(strings(i),chunkPos,j+1_pInt)*inRad
-             case('fraction')
-                 texture_Gauss(5,gauss,t) = IO_floatValue(strings(i),chunkPos,j+1_pInt)
           end select
       enddo
      enddo
@@ -984,28 +967,19 @@ end subroutine material_allocateSourceState
 !! calculates the volume of the grains and deals with texture components
 !--------------------------------------------------------------------------------------------------
 subroutine material_populateGrains
- use prec, only: &
-   dEq
  use math, only: &
+   math_EulertoR, &
    math_RtoEuler, &
-   math_EulerToR, &
    math_mul33x33, &
    math_range
  use mesh, only: &
-   theMesh, &
+   theMesh
-   mesh_ipVolume
  use config, only: &
    config_homogenization, &
    config_microstructure, &
-   config_deallocate, &
+   config_deallocate
-   homogenization_name, &
-   microstructure_name 
  use IO, only: &
    IO_error
- use debug, only: &
-   debug_level, &
-   debug_material, &
-   debug_levelBasic
 
  implicit none
  integer(pInt), dimension (:,:), allocatable :: Ngrains
@@ -1015,21 +989,17 @@ subroutine material_populateGrains
    randomOrder
  real(pReal), dimension (microstructure_maxNconstituents)  :: &
    rndArray
- real(pReal), dimension (:),     allocatable :: volumeOfGrain
  real(pReal), dimension (:,:),   allocatable :: orientationOfGrain
  real(pReal), dimension (3)                  :: orientation
- real(pReal), dimension (3,3)                :: symOrientation
  integer(pInt), dimension (:),   allocatable :: phaseOfGrain, textureOfGrain
- integer(pInt) :: t,e,i,g,j,m,c,r,homog,micro,sgn,hme, myDebug, &
+ integer(pInt) :: t,e,i,g,j,m,c,r,homog,micro,sgn,hme, &
                   phaseID,textureID,dGrains,myNgrains,myNorientations,myNconstituents, &
-                  grain,constituentGrain,ipGrain,symExtension, ip
+                  grain,constituentGrain,ipGrain,ip
  real(pReal) :: deviation,extreme,rnd
  integer(pInt),         dimension (:,:), allocatable :: Nelems                                      ! counts number of elements in homog, micro array
  type(group_int), dimension (:,:), allocatable :: elemsOfHomogMicro                                 ! lists element number in homog, micro array
 
- myDebug = debug_level(debug_material)
 
- allocate(material_volume(homogenization_maxNgrains,theMesh%elem%nIPs,theMesh%Nelems),       source=0.0_pReal)
  allocate(material_phase(homogenization_maxNgrains,theMesh%elem%nIPs,theMesh%Nelems),        source=0_pInt)
  allocate(material_texture(homogenization_maxNgrains,theMesh%elem%nIPs,theMesh%Nelems),      source=0_pInt)
  allocate(material_EulerAngles(3,homogenization_maxNgrains,theMesh%elem%nIPs,theMesh%Nelems),source=0.0_pReal)
@@ -1038,6 +1008,25 @@ subroutine material_populateGrains
  allocate(Nelems (size(config_homogenization),size(config_microstructure)),            source=0_pInt)
 
 
+  do e = 1, theMesh%Nelems
+    do i = 1, theMesh%elem%nIPs
+      homog = theMesh%homogenizationAt(e)
+      micro = theMesh%microstructureAt(e)
+      do c = 1, homogenization_Ngrains(homog)
+        material_phase(c,i,e)   = microstructure_phase(c,micro)
+        material_texture(c,i,e) = microstructure_texture(c,micro)
+        material_EulerAngles(1:3,c,i,e) = texture_Gauss(1:3,1,material_texture(c,i,e))
+        material_EulerAngles(1:3,c,i,e) = math_RtoEuler( &                                       ! translate back to Euler angles
+                                             math_mul33x33( &                                       ! pre-multiply
+                                               math_EulertoR(material_EulerAngles(1:3,c,i,e)), &    ! face-value orientation
+                                               texture_transformation(1:3,1:3,material_texture(c,i,e)) &          ! and transformation matrix
+                                             ) &
+                                             )
+      enddo
+    enddo
+  enddo
+
+ return
 !--------------------------------------------------------------------------------------------------
 ! precounting of elements for each homog/micro pair
  do e = 1_pInt, theMesh%Nelems
@@ -1075,47 +1064,17 @@ subroutine material_populateGrains
    elemsOfHomogMicro(homog,micro)%p(Nelems(homog,micro)) = e                                        ! remember elements active in this homog/micro pair
  enddo elementLooping
 
- allocate(volumeOfGrain(maxval(Ngrains)),       source=0.0_pReal)                                   ! reserve memory for maximum case
  allocate(phaseOfGrain(maxval(Ngrains)),        source=0_pInt)                                      ! reserve memory for maximum case
  allocate(textureOfGrain(maxval(Ngrains)),      source=0_pInt)                                      ! reserve memory for maximum case
  allocate(orientationOfGrain(3,maxval(Ngrains)),source=0.0_pReal)                                   ! reserve memory for maximum case
 
- if (iand(myDebug,debug_levelBasic) /= 0_pInt) then
+
-     write(6,'(/,a/)') ' MATERIAL grain population'
-     write(6,'(a32,1x,a32,1x,a6)') 'homogenization_name','microstructure_name','grain#'
- endif
  homogenizationLoop: do homog = 1_pInt,size(config_homogenization)
    dGrains = homogenization_Ngrains(homog)                                                          ! grain number per material point
-   microstructureLoop: do micro = 1_pInt,size(config_microstructure)                                                       ! all pairs of homog and micro
+   microstructureLoop: do micro = 1_pInt,size(config_microstructure)                                ! all pairs of homog and micro
      activePair: if (Ngrains(homog,micro) > 0_pInt) then
        myNgrains = Ngrains(homog,micro)                                                             ! assign short name for total number of grains to populate
        myNconstituents = microstructure_Nconstituents(micro)                                        ! assign short name for number of constituents
-       if (iand(myDebug,debug_levelBasic) /= 0_pInt) &
-         write(6,'(/,a32,1x,a32,1x,i6)') homogenization_name(homog),microstructure_name(micro),myNgrains
-
-
-!--------------------------------------------------------------------------------------------------
-! calculate volume of each grain
-
-       volumeOfGrain = 0.0_pReal
-       grain = 0_pInt
-
-       do hme = 1_pInt, Nelems(homog,micro)
-         e = elemsOfHomogMicro(homog,micro)%p(hme)                                                  ! my combination of homog and micro, only perform calculations for elements with homog, micro combinations which is indexed in cpElemsindex
-         if (microstructure_elemhomo(micro)) then                                                   ! homogeneous distribution of grains over each element's IPs
-           volumeOfGrain(grain+1_pInt:grain+dGrains) = sum(mesh_ipVolume(1:theMesh%elem%nIPs,e))/&
-                                                                         real(dGrains,pReal)        ! each grain combines size of all IPs in that element
-           grain = grain + dGrains                                                                  ! wind forward by Ngrains@IP
-         else
-           forall (i = 1_pInt:theMesh%elem%nIPs) &                                                   ! loop over IPs
-             volumeOfGrain(grain+(i-1)*dGrains+1_pInt:grain+i*dGrains) = &
-               mesh_ipVolume(i,e)/real(dGrains,pReal)                                               ! assign IPvolume/Ngrains@IP to all grains of IP
-           grain = grain + theMesh%elem%nIPs * dGrains                                               ! wind forward by Nips*Ngrains@IP
-         endif
-       enddo
-
-       if (grain /= myNgrains) &
-         call IO_error(0,el = homog,ip = micro,ext_msg = 'inconsistent grain count after volume calc')
 
 !--------------------------------------------------------------------------------------------------
 ! divide myNgrains as best over constituents
@@ -1165,8 +1124,7 @@ subroutine material_populateGrains
          phaseOfGrain  (grain+1_pInt:grain+NgrainsOfConstituent(i)) = phaseID                       ! assign resp. phase
          textureOfGrain(grain+1_pInt:grain+NgrainsOfConstituent(i)) = textureID                     ! assign resp. texture
 
-         myNorientations = ceiling(real(NgrainsOfConstituent(i),pReal)/&
+         myNorientations = ceiling(real(NgrainsOfConstituent(i),pReal)/1.0,pInt)                    ! max number of unique orientations (excl. symmetry)
-                                   real(texture_symmetry(textureID),pReal),pInt)                    ! max number of unique orientations (excl. symmetry)
 
 !--------------------------------------------------------------------------------------------------
 ! has texture components
@@ -1196,7 +1154,7 @@ subroutine material_populateGrains
 
          do j = 1_pInt,NgrainsOfConstituent(i)-1_pInt                                               ! walk thru grains of current constituent
            call random_number(rnd)
-           t = nint(rnd*real(NgrainsOfConstituent(i)-j,pReal)+real(j,pReal)+0.5_pReal,pInt)       ! select a grain in remaining list
+           t = nint(rnd*real(NgrainsOfConstituent(i)-j,pReal)+real(j,pReal)+0.5_pReal,pInt)         ! select a grain in remaining list
            m                               = phaseOfGrain(grain+t)                                  ! exchange current with random
            phaseOfGrain(grain+t)           = phaseOfGrain(grain+j)
            phaseOfGrain(grain+j)           = m
@@ -1244,7 +1202,6 @@ subroutine material_populateGrains
                                          currentGrainOfConstituent(c)))
                ipGrain = ipGrain + 1_pInt                                                           ! advance IP grain counter
                currentGrainOfConstituent(c)  = currentGrainOfConstituent(c) + 1_pInt                ! advance index of grain population for constituent c
-               material_volume(ipGrain,i,e)  = volumeOfGrain(grain+currentGrainOfConstituent(c))    ! assign properties
                material_phase(ipGrain,i,e)   = phaseOfGrain(grain+currentGrainOfConstituent(c))
                material_texture(ipGrain,i,e) = textureOfGrain(grain+currentGrainOfConstituent(c))
                material_EulerAngles(1:3,ipGrain,i,e) = orientationOfGrain(1:3,grain+currentGrainOfConstituent(c))
@@ -1254,7 +1211,6 @@ subroutine material_populateGrains
            grain = sum(NgrainsOfConstituent(1:c-1_pInt))                                            ! figure out actual starting index in overall/consecutive grain population
            do ipGrain = ipGrain + 1_pInt, dGrains                                                   ! ensure last constituent fills up to dGrains
              currentGrainOfConstituent(c)  = currentGrainOfConstituent(c) + 1_pInt
-             material_volume(ipGrain,i,e)  = volumeOfGrain(grain+currentGrainOfConstituent(c))
              material_phase(ipGrain,i,e)   = phaseOfGrain(grain+currentGrainOfConstituent(c))
              material_texture(ipGrain,i,e) = textureOfGrain(grain+currentGrainOfConstituent(c))
              material_EulerAngles(1:3,ipGrain,i,e) = orientationOfGrain(1:3,grain+currentGrainOfConstituent(c))
@@ -1263,7 +1219,6 @@ subroutine material_populateGrains
          enddo
 
          do i = i, theMesh%elem%nIPs                                                                 ! loop over IPs to (possibly) distribute copies from first IP
-           material_volume (1_pInt:dGrains,i,e) = material_volume (1_pInt:dGrains,1,e)
            material_phase  (1_pInt:dGrains,i,e) = material_phase  (1_pInt:dGrains,1,e)
            material_texture(1_pInt:dGrains,i,e) = material_texture(1_pInt:dGrains,1,e)
            material_EulerAngles(1:3,1_pInt:dGrains,i,e) = material_EulerAngles(1:3,1_pInt:dGrains,1,e)
@@ -1275,7 +1230,7 @@ subroutine material_populateGrains
  enddo homogenizationLoop
 
  deallocate(texture_transformation)
- deallocate(elemsOfHomogMicro)
+
  call config_deallocate('material.config/microstructure')
 
 end subroutine material_populateGrains

src/plastic_disloUCLA.f90

@@ -121,15 +121,7 @@ subroutine plastic_disloUCLA_init()
     math_expand
   use IO, only: &
     IO_error
-  use material, only: &
+  use material
-    phase_plasticity, &
-    phase_plasticityInstance, &
-    phase_Noutput, &
-    material_allocatePlasticState, &
-    PLASTICITY_DISLOUCLA_label, &
-    PLASTICITY_DISLOUCLA_ID, &
-    material_phase, &
-    plasticState
   use config, only: &
     config_phase
   use lattice

src/plastic_dislotwin.f90

@@ -185,15 +185,7 @@ subroutine plastic_dislotwin_init
    PI
  use IO, only: &
    IO_error
- use material, only: &
+ use material
-   phase_plasticity, &
-   phase_plasticityInstance, &
-   phase_Noutput, &
-   material_allocatePlasticState, &
-   PLASTICITY_DISLOTWIN_label, &
-   PLASTICITY_DISLOTWIN_ID, &
-   material_phase, &  
-   plasticState
  use config, only: &
    config_phase
  use lattice

src/plastic_isotropic.f90

@@ -91,19 +91,7 @@ subroutine plastic_isotropic_init
     debug_levelBasic
   use IO, only: &
     IO_error
-  use material, only: &
+  use material
-    phase_plasticity, &
-    phase_plasticityInstance, &
-    phase_Noutput, &
-    material_allocatePlasticState, &
-    PLASTICITY_ISOTROPIC_label, &
-    PLASTICITY_ISOTROPIC_ID, &
-    material_phase, &
-    plasticState
-#ifdef DEBUG
-  use material, only: &
-    phasememberAt
-#endif
   use config, only: &
     config_phase
   use lattice
@@ -464,7 +452,7 @@ function plastic_isotropic_postResults(Mp,instance,of) result(postResults)
         c = c + 1
       case (dot_gamma_ID)
         postResults(c+1) = prm%dot_gamma_0 &
-                         * (sqrt(1.5_pReal) * norm_Mp /(prm%M * stt%xi(of)))**prm%n
+                              * (sqrt(1.5_pReal) * norm_Mp /(prm%M * stt%xi(of)))**prm%n
         c = c + 1
  
     end select
@@ -479,11 +467,10 @@ end function plastic_isotropic_postResults
 !> @brief writes results to HDF5 output file
 !--------------------------------------------------------------------------------------------------
 subroutine plastic_isotropic_results(instance,group)
-#if defined(PETSc) || defined(DAMASK_HDF5)
+#if defined(PETSc) || defined(DAMASKHDF5)
-  use results, only: &
+  use results
-    results_writeDataset
 
-  integer,          intent(in) :: instance
+  integer, intent(in) :: instance
   character(len=*), intent(in) :: group
   
   integer :: o
@@ -496,10 +483,9 @@ subroutine plastic_isotropic_results(instance,group)
     end select
   enddo outputsLoop
   end associate
-  
 #else
-  integer,          intent(in) :: instance
+  integer, intent(in) :: instance
-  character(len=*), intent(in) :: group
+  character(len=*) :: group
 #endif
 
 end subroutine plastic_isotropic_results

src/plastic_kinematichardening.f90

@@ -112,19 +112,7 @@ subroutine plastic_kinehardening_init
    math_expand
  use IO, only: &
    IO_error
- use material, only: &
+ use material
-   phase_plasticity, &
-   phase_plasticityInstance, &
-   phase_Noutput, &
-   material_allocatePlasticState, &
-   PLASTICITY_kinehardening_label, &
-   PLASTICITY_kinehardening_ID, &
-   material_phase, &
-   plasticState
-#ifdef DEBUG
- use material, only: &
-   phasememberAt
-#endif
  use config, only: &
    config_phase
  use lattice

src/plastic_none.f90

@@ -23,13 +23,7 @@ subroutine plastic_none_init
    debug_level, &
    debug_constitutive, &
    debug_levelBasic
- use material, only: &
+ use material
-   phase_plasticity, &
-   material_allocatePlasticState, &
-   PLASTICITY_NONE_label, &
-   PLASTICITY_NONE_ID, &
-   material_phase, &
-   plasticState
 
  integer :: &
    Ninstance, &

src/plastic_nonlocal.f90

@@ -248,15 +248,7 @@ subroutine plastic_nonlocal_init
     debug_levelBasic
   use mesh, only: &
     theMesh
-  use material, only: &
+  use material
-    phase_plasticity, &
-    phase_plasticityInstance, &
-    phase_Noutput, &
-    PLASTICITY_NONLOCAL_label, &
-    PLASTICITY_NONLOCAL_ID, &
-    plasticState, &
-    material_phase, &
-    material_allocatePlasticState
   use config
   use lattice
 

src/plastic_phenopowerlaw.f90

@@ -116,15 +116,7 @@ subroutine plastic_phenopowerlaw_init
    math_expand
  use IO, only: &
    IO_error
- use material, only: &
+ use material
-   phase_plasticity, &
-   phase_plasticityInstance, &
-   phase_Noutput, &
-   material_allocatePlasticState, &
-   PLASTICITY_PHENOPOWERLAW_LABEL, &
-   PLASTICITY_PHENOPOWERLAW_ID, &
-   material_phase, &
-   plasticState
  use config, only: &
    config_phase
  use lattice

src/quit.f90

@@ -34,13 +34,13 @@ subroutine quit(stop_id)
 #endif
  
  call date_and_time(values = dateAndTime)
- write(6,'(/,a)') 'DAMASK terminated on:'
+ write(6,'(/,a)') ' DAMASK terminated on:'
- write(6,'(a,2(i2.2,a),i4.4)') 'Date:               ',dateAndTime(3),'/',&
+ write(6,'(a,2(i2.2,a),i4.4)') ' Date:               ',dateAndTime(3),'/',&
-                                                      dateAndTime(2),'/',&
+                                                       dateAndTime(2),'/',&
-                                                      dateAndTime(1)
+                                                       dateAndTime(1)
- write(6,'(a,2(i2.2,a),i2.2)') 'Time:               ',dateAndTime(5),':',&
+ write(6,'(a,2(i2.2,a),i2.2)') ' Time:               ',dateAndTime(5),':',&
-                                                      dateAndTime(6),':',&
+                                                       dateAndTime(6),':',&
-                                                      dateAndTime(7)
+                                                       dateAndTime(7)
 
  if (stop_id == 0 .and. ierr == 0 .and. error == 0) stop 0                                          ! normal termination
  if (stop_id == 2 .and. ierr == 0 .and. error == 0) stop 2                                          ! not all incs converged

src/results.f90

@@ -66,9 +66,9 @@ subroutine results_init
   write(6,'(a)')   ' https://doi.org/10.1007/s40192-018-0118-7'
 
   resultsFile = HDF5_openFile(trim(getSolverJobName())//'.hdf5','w',.true.)
-  call HDF5_addAttribute(resultsFile,'DADF5-version',0.1)
+  call HDF5_addAttribute(resultsFile,'DADF5-version',0.2)
   call HDF5_addAttribute(resultsFile,'DADF5-major',0)
-  call HDF5_addAttribute(resultsFile,'DADF5-minor',1)
+  call HDF5_addAttribute(resultsFile,'DADF5-minor',2)
   call HDF5_addAttribute(resultsFile,'DAMASK',DAMASKVERSION)
   call get_command(commandLine)
   call HDF5_addAttribute(resultsFile,'call',trim(commandLine))
@@ -103,7 +103,7 @@ end subroutine results_closeJobFile
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief closes the results file
+!> @brief creates the group of increment and adds time as attribute to the file
 !--------------------------------------------------------------------------------------------------
 subroutine results_addIncrement(inc,time)
  
@@ -250,6 +250,8 @@ subroutine results_writeScalarDataset_real(group,dataset,label,description,SIuni
   
 #ifdef PETSc
   call HDF5_write(groupHandle,dataset,label,.true.)
+#else
+  call HDF5_write(groupHandle,dataset,label,.false.)
 #endif
   
   if (HDF5_objectExists(groupHandle,label)) &
@@ -277,6 +279,8 @@ subroutine results_writeVectorDataset_real(group,dataset,label,description,SIuni
   
 #ifdef PETSc
   call HDF5_write(groupHandle,dataset,label,.true.)
+#else
+  call HDF5_write(groupHandle,dataset,label,.false.)
 #endif
   
   if (HDF5_objectExists(groupHandle,label)) &
@@ -305,6 +309,8 @@ subroutine results_writeTensorDataset_real(group,dataset,label,description,SIuni
   
 #ifdef PETSc
   call HDF5_write(groupHandle,dataset,label,.true.)
+#else
+  call HDF5_write(groupHandle,dataset,label,.false.)
 #endif
   
   if (HDF5_objectExists(groupHandle,label)) &
@@ -333,6 +339,8 @@ subroutine results_writeVectorDataset_int(group,dataset,label,description,SIunit
   
 #ifdef PETSc
   call HDF5_write(groupHandle,dataset,label,.true.)
+#else
+  call HDF5_write(groupHandle,dataset,label,.false.)
 #endif
   
   if (HDF5_objectExists(groupHandle,label)) &
@@ -347,7 +355,7 @@ end subroutine results_writeVectorDataset_int
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief stores a vector dataset in a group
+!> @brief stores a tensor dataset in a group
 !--------------------------------------------------------------------------------------------------
 subroutine results_writeTensorDataset_int(group,dataset,label,description,SIunit)
 
@@ -361,6 +369,8 @@ subroutine results_writeTensorDataset_int(group,dataset,label,description,SIunit
   
 #ifdef PETSc
   call HDF5_write(groupHandle,dataset,label,.true.)
+#else
+  call HDF5_write(groupHandle,dataset,label,.false.)
 #endif
   
   if (HDF5_objectExists(groupHandle,label)) &
@@ -375,7 +385,7 @@ end subroutine results_writeTensorDataset_int
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief stores a vector dataset in a group
+!> @brief stores a scalar dataset in a group
 !--------------------------------------------------------------------------------------------------
 subroutine results_writeScalarDataset_rotation(group,dataset,label,description,lattice_structure)
   use rotations, only: &
@@ -391,6 +401,8 @@ subroutine results_writeScalarDataset_rotation(group,dataset,label,description,l
   
 #ifdef PETSc
   call HDF5_write(groupHandle,dataset,label,.true.)
+#else
+  call HDF5_write(groupHandle,dataset,label,.false.)
 #endif
   
   if (HDF5_objectExists(groupHandle,label)) &