Merge remote-tracking branch 'origin/development' into result.incs

.github/workflows/Fortran.yml

@@ -2,7 +2,7 @@ name: Grid and Mesh Solver
 on: [push]
 
 env:
-  PETSC_VERSION: '3.16.2'
+  PETSC_VERSION: '3.17.0'
   HOMEBREW_NO_ANALYTICS: 'ON'     # Make Homebrew installation a little quicker
   HOMEBREW_NO_AUTO_UPDATE: 'ON'
   HOMEBREW_NO_BOTTLE_SOURCE_FALLBACK: 'ON'

.github/workflows/Python.yml

@@ -11,6 +11,7 @@ jobs:
       matrix:
         python-version: ['3.8', '3.9'] #, '3.10']
         os: [ubuntu-latest, macos-latest, windows-latest]
+      fail-fast: false
 
     steps:
       - uses: actions/checkout@v2

CMakeLists.txt

@@ -9,8 +9,10 @@ endif()
 
 # Dummy project to determine compiler names and version
 project(Prerequisites LANGUAGES)
-set(ENV{PKG_CONFIG_PATH} "$ENV{PETSC_DIR}/$ENV{PETSC_ARCH}/lib/pkgconfig")
+set(ENV{PKG_CONFIG_PATH} "$ENV{PETSC_DIR}/$ENV{PETSC_ARCH}/lib/pkgconfig:$ENV{PKG_CONFIG_PATH}")
-pkg_check_modules(PETSC REQUIRED PETSc>=3.12.0 PETSc<3.17.0)
+pkg_check_modules(PETSC_MIN REQUIRED PETSc>=3.12.0 QUIET) #CMake does not support version range
+pkg_check_modules(PETSC REQUIRED PETSc<3.18.0)
+
 pkg_get_variable(CMAKE_Fortran_COMPILER PETSc fcompiler)
 pkg_get_variable(CMAKE_C_COMPILER PETSc ccompiler)
 
@@ -25,6 +27,13 @@ else()
 endif()
 add_definitions("-D${DAMASK_SOLVER}")
 
+# EXPERIMENTAL: This might help to detect HDF5 and FFTW3 in the future if PETSc is not aware of them
+set(ENV{PKG_CONFIG_PATH} "$ENV{PETSC_DIR}/$ENV{PETSC_ARCH}/externalpackages:$ENV{PKG_CONFIG_PATH}")
+pkg_check_modules(HDF5 hdf5)
+pkg_check_modules(FFTW3 fftw3)
+pkg_check_modules(fYAML libfyaml)
+pkg_check_modules(zlib zlib)
+
 file(STRINGS ${PROJECT_SOURCE_DIR}/VERSION DAMASK_VERSION)
 
 message("\nBuilding ${CMAKE_PROJECT_NAME} ${DAMASK_VERSION}\n")
@@ -32,6 +41,9 @@ message("\nBuilding ${CMAKE_PROJECT_NAME} ${DAMASK_VERSION}\n")
 add_definitions(-DPETSC)
 add_definitions(-DDAMASKVERSION="${DAMASK_VERSION}")
 add_definitions(-DCMAKE_SYSTEM="${CMAKE_SYSTEM}")
+if(PETSC_VERSION VERSION_GREATER_EQUAL 3.17.0)
+   add_definitions("-DCHKERRQ=PetscCall")
+endif()
 
 if(CMAKE_BUILD_TYPE STREQUAL "")
     set(CMAKE_BUILD_TYPE "RELEASE")
@@ -104,8 +116,18 @@ if(CMAKE_BUILD_TYPE STREQUAL "DEBUG")
     set(CMAKE_Fortran_LINK_EXECUTABLE           "${CMAKE_Fortran_LINK_EXECUTABLE} ${DEBUG_FLAGS}")
 endif()
 
-set(CMAKE_Fortran_FLAGS_${CMAKE_BUILD_TYPE}   "${CMAKE_Fortran_FLAGS_${CMAKE_BUILD_TYPE}} ${PETSC_INCLUDES} ${BUILDCMD_POST}")
+set(CMAKE_Fortran_FLAGS_${CMAKE_BUILD_TYPE} "${CMAKE_Fortran_FLAGS_${CMAKE_BUILD_TYPE}} ${PETSC_INCLUDES} ${BUILDCMD_POST}")
-set(CMAKE_Fortran_LINK_EXECUTABLE             "${CMAKE_Fortran_LINK_EXECUTABLE} <OBJECTS> -o <TARGET> <LINK_LIBRARIES> -L${PETSC_LIBRARY_DIRS} -lpetsc ${PETSC_EXTERNAL_LIB} -lz ${BUILDCMD_POST}")
+
+set(CMAKE_Fortran_LINK_EXECUTABLE "${CMAKE_Fortran_LINK_EXECUTABLE} <OBJECTS> -o <TARGET> <LINK_LIBRARIES> -L${PETSC_LIBRARY_DIRS} -lpetsc ${PETSC_EXTERNAL_LIB} -lz")
+
+if(fYAML_FOUND STREQUAL "1")
+    set(CMAKE_Fortran_LINK_EXECUTABLE "${CMAKE_Fortran_LINK_EXECUTABLE} -L${fYAML_LIBRARY_DIRS} -l${fYAML_LIBRARIES}")
+    add_definitions(-DFYAML)
+    pkg_get_variable(fYAML_INCLUDE_DIR libfyaml includedir) # fYAML_INCLUDE_DIRS and fYAML_CFLAGS are not working
+    set(CMAKE_C_FLAGS_${CMAKE_BUILD_TYPE} "${CMAKE_C_FLAGS_${CMAKE_BUILD_TYPE}} -I${fYAML_INCLUDE_DIR}")
+endif()
+
+set(CMAKE_Fortran_LINK_EXECUTABLE "${CMAKE_Fortran_LINK_EXECUTABLE} ${BUILDCMD_POST}")
 
 message("Fortran Compiler Flags:\n${CMAKE_Fortran_FLAGS_${CMAKE_BUILD_TYPE}}\n")
 message("C Compiler Flags:\n${CMAKE_C_FLAGS_${CMAKE_BUILD_TYPE}}\n")

PRIVATE

@@ -1 +1 @@
-Subproject commit 459326e9840c843ade72b04cf28e50889c9779f1
+Subproject commit e3d9c0cc3ffa1435003f934a0f0e6c7d969a022a

README

@@ -1,13 +0,0 @@
-DAMASK - The Düsseldorf Advanced Material Simulation Kit
-Visit damask.mpie.de for installation and usage instructions
-
-CONTACT INFORMATION
-
-Max-Planck-Institut für Eisenforschung GmbH
-Max-Planck-Str. 1
-40237 Düsseldorf
-Germany
-
-damask@mpie.de
-https://damask.mpie.de
-https://git.damask.mpie.de

README.md

@@ -0,0 +1,15 @@
+# DAMASK - The Düsseldorf Advanced Material Simulation Kit
+
+Visit [damask.mpie.de](https://damask.mpie.de) for installation and usage instructions
+
+## Contact Information
+
+Max-Planck-Institut für Eisenforschung GmbH  
+Max-Planck-Str. 1  
+40237 Düsseldorf  
+Germany  
+
+damask@mpie.de  
+https://damask.mpie.de  
+https://git.damask.mpie.de  
+

python/damask/VERSION

@@ -1 +1 @@
-v3.0.0-alpha6-184-g1f98b04d4
+v3.0.0-alpha6-228-g758ad6072

python/damask/_grid.py

@@ -292,7 +292,7 @@ class Grid:
         >>> import damask
         >>> N_grains = 20
         >>> cells = (32,32,32)
-        >>> damask.util.run(f'neper -T -n {N_grains} -tesrsize {cells[0]}:{cells[1]}:{cells[2]} -periodicity "all" -format "vtk"')
+        >>> damask.util.run(f'neper -T -n {N_grains} -tesrsize {cells[0]}:{cells[1]}:{cells[2]} -periodicity all -format vtk')
         >>> damask.Grid.load_Neper(f'n{N_grains}-id1.vtk')
         cells:  32 × 32 × 32
         size:   1.0 × 1.0 × 1.0 m³

python/damask/_result.py

@@ -1154,14 +1154,14 @@ class Result:
 
         >>> import damask
         >>> r = damask.Result('my_file.hdf5')
-        >>> r.strain(t='U',m=0.5)
+        >>> r.add_strain(t='U',m=0.5)
 
         Add the plastic Euler-Almansi strain based on the
         plastic deformation gradient 'F_p':
 
         >>> import damask
         >>> r = damask.Result('my_file.hdf5')
-        >>> r.strain('F_p','V',-1)
+        >>> r.add_strain('F_p','V',-1)
 
         """
         self._add_generic_pointwise(self._add_strain,{'F':F},{'t':t,'m':m})

src/C_routines.c

@@ -7,7 +7,11 @@
 #include <pwd.h>
 #include <sys/types.h>
 #include <sys/stat.h>
-#include "zlib.h"
+#include <zlib.h>
+
+#ifdef FYAML
+#include <libfyaml.h>
+#endif
 
 #define PATHLEN 4096
 #define STRLEN   256
@@ -80,3 +84,26 @@ void inflate_c(const uLong *s_deflated, const uLong *s_inflated, const Byte defl
     }
   }
 }
+
+#ifdef FYAML
+void to_flow_c(char **flow, int* length_flow, const char *mixed){
+  struct fy_document *fyd = NULL;
+  enum fy_emitter_cfg_flags emit_flags = FYECF_MODE_FLOW_ONELINE | FYECF_STRIP_LABELS | FYECF_STRIP_TAGS |FYECF_STRIP_DOC;
+
+  fyd = fy_document_build_from_string(NULL, mixed, -1);
+  if (!fyd) {
+    *length_flow = -1;
+    return;
+  }
+  int err = fy_document_resolve(fyd);
+  if (err) {
+    *length_flow = -1;
+    return;
+  }
+
+  *flow = fy_emit_document_to_string(fyd,emit_flags);
+  *length_flow = strlen(*flow);
+
+  fy_document_destroy(fyd);
+}
+#endif

src/DAMASK_interface.f90

@@ -11,7 +11,7 @@
 !--------------------------------------------------------------------------------------------------
 #define PETSC_MAJOR 3
 #define PETSC_MINOR_MIN 12
-#define PETSC_MINOR_MAX 16
+#define PETSC_MINOR_MAX 17
 
 module DAMASK_interface
   use, intrinsic :: ISO_fortran_env

src/HDF5_utilities.f90

@@ -48,6 +48,7 @@ module HDF5_utilities
 !> @details for parallel IO, all dimension except for the last need to match
 !--------------------------------------------------------------------------------------------------
   interface HDF5_write
+#if defined(__GFORTRAN__) && __GNUC__<11
     module procedure HDF5_write_real1
     module procedure HDF5_write_real2
     module procedure HDF5_write_real3
@@ -55,7 +56,6 @@ module HDF5_utilities
     module procedure HDF5_write_real5
     module procedure HDF5_write_real6
     module procedure HDF5_write_real7
-
     module procedure HDF5_write_int1
     module procedure HDF5_write_int2
     module procedure HDF5_write_int3
@@ -63,6 +63,10 @@ module HDF5_utilities
     module procedure HDF5_write_int5
     module procedure HDF5_write_int6
     module procedure HDF5_write_int7
+#else
+    module procedure HDF5_write_real
+    module procedure HDF5_write_int
+#endif
   end interface HDF5_write
 
 !--------------------------------------------------------------------------------------------------
@@ -1210,6 +1214,7 @@ subroutine HDF5_read_int7(dataset,loc_id,datasetName,parallel)
 
 end subroutine HDF5_read_int7
 
+#if defined(__GFORTRAN__) && __GNUC__<11
 
 !--------------------------------------------------------------------------------------------------
 !> @brief write dataset of type real with 1 dimension
@@ -1499,6 +1504,71 @@ subroutine HDF5_write_real7(dataset,loc_id,datasetName,parallel)
 
 end subroutine HDF5_write_real7
 
+#else
+
+!--------------------------------------------------------------------------------------------------
+!> @brief write dataset of type real with 1-7 dimension
+!--------------------------------------------------------------------------------------------------
+subroutine HDF5_write_real(dataset,loc_id,datasetName,parallel)
+
+  real(pReal),      intent(in), dimension(..) :: dataset                                            !< data written to file
+  integer(HID_T),   intent(in)  :: loc_id                                                           !< file or group handle
+  character(len=*), intent(in)  :: datasetName                                                      !< name of the dataset in the file
+  logical, intent(in), optional :: parallel                                                         !< dataset is distributed over multiple processes
+
+
+  integer :: hdferr
+  integer(HID_T)   :: dset_id, filespace_id, memspace_id, plist_id
+  integer(HSIZE_T), dimension(rank(dataset)) :: &
+    myStart, &
+    myShape, &                                                                                      !< shape of the dataset (this process)
+    totalShape                                                                                      !< shape of the dataset (all processes)
+
+!---------------------------------------------------------------------------------------------------
+! determine shape of dataset
+  myShape = int(shape(dataset),HSIZE_T)
+  if (any(myShape(1:size(myShape)-1) == 0)) return                                                  !< empty dataset (last dimension can be empty)
+
+  if (present(parallel)) then
+    call initialize_write(dset_id, filespace_id, memspace_id, plist_id, &
+                          myStart, totalShape,loc_id,myShape,datasetName,H5T_NATIVE_DOUBLE,parallel)
+  else
+    call initialize_write(dset_id, filespace_id, memspace_id, plist_id, &
+                          myStart, totalShape,loc_id,myShape,datasetName,H5T_NATIVE_DOUBLE,parallel_default)
+  end if
+
+  if (product(totalShape) /= 0) then
+    select rank(dataset)
+      rank (1)
+        call H5Dwrite_f(dset_id, H5T_NATIVE_DOUBLE,dataset,int(totalShape,HSIZE_T), hdferr,&
+                        file_space_id = filespace_id, mem_space_id = memspace_id, xfer_prp = plist_id)
+      rank (2)
+        call H5Dwrite_f(dset_id, H5T_NATIVE_DOUBLE,dataset,int(totalShape,HSIZE_T), hdferr,&
+                        file_space_id = filespace_id, mem_space_id = memspace_id, xfer_prp = plist_id)
+      rank (3)
+        call H5Dwrite_f(dset_id, H5T_NATIVE_DOUBLE,dataset,int(totalShape,HSIZE_T), hdferr,&
+                        file_space_id = filespace_id, mem_space_id = memspace_id, xfer_prp = plist_id)
+      rank (4)
+        call H5Dwrite_f(dset_id, H5T_NATIVE_DOUBLE,dataset,int(totalShape,HSIZE_T), hdferr,&
+                        file_space_id = filespace_id, mem_space_id = memspace_id, xfer_prp = plist_id)
+      rank (5)
+        call H5Dwrite_f(dset_id, H5T_NATIVE_DOUBLE,dataset,int(totalShape,HSIZE_T), hdferr,&
+                        file_space_id = filespace_id, mem_space_id = memspace_id, xfer_prp = plist_id)
+      rank (6)
+        call H5Dwrite_f(dset_id, H5T_NATIVE_DOUBLE,dataset,int(totalShape,HSIZE_T), hdferr,&
+                        file_space_id = filespace_id, mem_space_id = memspace_id, xfer_prp = plist_id)
+      rank (7)
+        call H5Dwrite_f(dset_id, H5T_NATIVE_DOUBLE,dataset,int(totalShape,HSIZE_T), hdferr,&
+                        file_space_id = filespace_id, mem_space_id = memspace_id, xfer_prp = plist_id)
+    end select
+    if(hdferr < 0) error stop 'HDF5 error'
+  end if
+
+  call finalize_write(plist_id, dset_id, filespace_id, memspace_id)
+
+end subroutine HDF5_write_real
+#endif
+
 
 !--------------------------------------------------------------------------------------------------
 !> @brief Write dataset of type string (scalar).
@@ -1561,6 +1631,7 @@ subroutine HDF5_write_str(dataset,loc_id,datasetName)
 
 end subroutine HDF5_write_str
 
+#if defined(__GFORTRAN__) && __GNUC__<11
 
 !--------------------------------------------------------------------------------------------------
 !> @brief write dataset of type integer with 1 dimension
@@ -1849,6 +1920,70 @@ subroutine HDF5_write_int7(dataset,loc_id,datasetName,parallel)
 
 end subroutine HDF5_write_int7
 
+#else
+
+!--------------------------------------------------------------------------------------------------
+!> @brief write dataset of type integer with 1-7 dimension
+!--------------------------------------------------------------------------------------------------
+subroutine HDF5_write_int(dataset,loc_id,datasetName,parallel)
+
+  integer,          intent(in), dimension(..) :: dataset                                            !< data written to file
+  integer(HID_T),   intent(in)  :: loc_id                                                           !< file or group handle
+  character(len=*), intent(in)  :: datasetName                                                      !< name of the dataset in the file
+  logical, intent(in), optional :: parallel                                                         !< dataset is distributed over multiple processes
+
+
+  integer :: hdferr
+  integer(HID_T)   :: dset_id, filespace_id, memspace_id, plist_id
+  integer(HSIZE_T), dimension(rank(dataset)) :: &
+    myStart, &
+    myShape, &                                                                                      !< shape of the dataset (this process)
+    totalShape                                                                                      !< shape of the dataset (all processes)
+
+!---------------------------------------------------------------------------------------------------
+! determine shape of dataset
+  myShape = int(shape(dataset),HSIZE_T)
+  if (any(myShape(1:size(myShape)-1) == 0)) return                                                  !< empty dataset (last dimension can be empty)
+
+  if (present(parallel)) then
+    call initialize_write(dset_id, filespace_id, memspace_id, plist_id, &
+                          myStart, totalShape, loc_id,myShape,datasetName,H5T_NATIVE_INTEGER,parallel)
+  else
+    call initialize_write(dset_id, filespace_id, memspace_id, plist_id, &
+                          myStart, totalShape, loc_id,myShape,datasetName,H5T_NATIVE_INTEGER,parallel_default)
+  end if
+
+  if (product(totalShape) /= 0) then
+    select rank(dataset)
+      rank(1)
+        call H5Dwrite_f(dset_id, H5T_NATIVE_INTEGER,dataset,int(totalShape,HSIZE_T), hdferr,&
+                       file_space_id = filespace_id, mem_space_id = memspace_id, xfer_prp = plist_id)
+      rank(2)
+        call H5Dwrite_f(dset_id, H5T_NATIVE_INTEGER,dataset,int(totalShape,HSIZE_T), hdferr,&
+                       file_space_id = filespace_id, mem_space_id = memspace_id, xfer_prp = plist_id)
+      rank(3)
+        call H5Dwrite_f(dset_id, H5T_NATIVE_INTEGER,dataset,int(totalShape,HSIZE_T), hdferr,&
+                       file_space_id = filespace_id, mem_space_id = memspace_id, xfer_prp = plist_id)
+      rank(4)
+        call H5Dwrite_f(dset_id, H5T_NATIVE_INTEGER,dataset,int(totalShape,HSIZE_T), hdferr,&
+                       file_space_id = filespace_id, mem_space_id = memspace_id, xfer_prp = plist_id)
+      rank(5)
+        call H5Dwrite_f(dset_id, H5T_NATIVE_INTEGER,dataset,int(totalShape,HSIZE_T), hdferr,&
+                       file_space_id = filespace_id, mem_space_id = memspace_id, xfer_prp = plist_id)
+      rank(6)
+        call H5Dwrite_f(dset_id, H5T_NATIVE_INTEGER,dataset,int(totalShape,HSIZE_T), hdferr,&
+                       file_space_id = filespace_id, mem_space_id = memspace_id, xfer_prp = plist_id)
+      rank(7)
+        call H5Dwrite_f(dset_id, H5T_NATIVE_INTEGER,dataset,int(totalShape,HSIZE_T), hdferr,&
+                       file_space_id = filespace_id, mem_space_id = memspace_id, xfer_prp = plist_id)
+    end select
+    if(hdferr < 0) error stop 'HDF5 error'
+  end if
+
+  call finalize_write(plist_id, dset_id, filespace_id, memspace_id)
+
+end subroutine HDF5_write_int
+#endif
 
 !--------------------------------------------------------------------------------------------------
 !> @brief initialize HDF5 handles, determines global shape and start for parallel read

src/IO.f90

@@ -483,7 +483,9 @@ subroutine IO_error(error_ID,el,ip,g,instance,ext_msg)
     case (701)
       msg = 'Incorrect indent/Null value not allowed'
     case (702)
-      msg = 'Invalid use of flow yaml'
+      msg = 'Invalid use of flow YAML'
+    case (703)
+      msg = 'Invalid YAML'
     case (704)
       msg = 'Space expected after a colon for <key>: <value> pair'
     case (705)

src/VTI.f90

@@ -81,7 +81,7 @@ subroutine VTI_readDataset_raw(base64_str,dataType,headerType,compressed, &
   character(len=:), allocatable, intent(out) :: dataType, headerType, base64_str
   logical, intent(out) :: compressed
 
-  logical :: inFile,inImage,gotCellData
+  logical :: inFile, inImage
   integer(pI64) :: &
     startPos, endPos, &
     s
@@ -152,10 +152,9 @@ subroutine VTI_readCellsSizeOrigin(cells,geomSize,origin, &
     fileContent
 
   character(len=:), allocatable :: dataType, headerType
-  logical :: inFile,inImage,gotCellData,compressed
+  logical :: inFile, inImage, compressed
   integer(pI64) :: &
-    startPos, endPos, &
+    startPos, endPos
-    s
 
 
   cells = -1

src/YAML_parse.f90

@@ -8,6 +8,9 @@ module YAML_parse
   use prec
   use IO
   use YAML_types
+#ifdef FYAML
+  use system_routines
+#endif
 
   implicit none
   private
@@ -16,14 +19,34 @@ module YAML_parse
     YAML_parse_init, &
     YAML_parse_str
 
+#ifdef FYAML
+  interface
+
+    subroutine to_flow_C(flow,length_flow,mixed) bind(C)
+      use, intrinsic :: ISO_C_Binding, only: C_INT, C_CHAR, C_PTR
+
+      type(C_PTR), intent(out) :: flow
+      integer(C_INT), intent(out) :: length_flow
+      character(kind=C_CHAR), dimension(*), intent(in) :: mixed
+    end subroutine to_flow_C
+
+  end interface
+#endif
+
+
 contains
 
 !--------------------------------------------------------------------------------------------------
 !> @brief Do sanity checks.
 !--------------------------------------------------------------------------------------------------
-subroutine YAML_parse_init
+subroutine YAML_parse_init()
 
-  call selfTest
+  print'(/,1x,a)', '<<<+-  YAML_parse init  -+>>>'
+#ifdef FYAML
+  print'(/,1x,a)', 'libfyaml powered'
+#else
+  call selfTest()
+#endif
 
 end subroutine YAML_parse_init
 
@@ -60,7 +83,7 @@ recursive function parse_flow(YAML_flow) result(node)
     s, &                                                                                            ! start position of dictionary or list
     d                                                                                               ! position of key: value separator (':')
 
-  flow_string = trim(adjustl(YAML_flow(:)))
+  flow_string = trim(adjustl(YAML_flow))
   if (len_trim(flow_string) == 0) then
     node => emptyDict
     return
@@ -145,8 +168,11 @@ logical function quotedString(line)
 
   character(len=*), intent(in) :: line
 
+
   quotedString = .false.
 
+  if (len(line) == 0) return
+
   if (scan(line(:1),IO_QUOTES) == 1) then
     quotedString = .true.
     if(line(len(line):len(line)) /= line(:1)) call IO_error(710,ext_msg=line)
@@ -155,8 +181,37 @@ logical function quotedString(line)
 end function quotedString
 
 
+#ifdef FYAML
 !--------------------------------------------------------------------------------------------------
-! @brief Returns Indentation.
+! @brief Convert all block style YAML parts to flow style.
+!--------------------------------------------------------------------------------------------------
+function to_flow(mixed) result(flow)
+
+  character(len=*), intent(in) :: mixed
+  character(:,C_CHAR), allocatable :: flow
+
+  type(C_PTR) :: str_ptr
+  integer(C_INT) :: strlen
+
+
+  call to_flow_C(str_ptr,strlen,f_c_string(mixed))
+  if (strlen < 1) call IO_error(703,ext_msg='libyfaml')
+  allocate(character(len=strlen,kind=c_char) :: flow)
+
+  block
+    character(len=strlen,kind=c_char), pointer :: s
+    call c_f_pointer(str_ptr,s)
+    flow = s(:len(s)-1)
+  end block
+
+  call free_C(str_ptr)
+
+end function to_flow
+
+
+#else
+!--------------------------------------------------------------------------------------------------
+! @brief Determine Indentation.
 ! @details It determines the indentation level for a given block/line.
 ! In cases for nested lists, an offset is added to determine the indent of the item block (skip
 ! leading dashes)
@@ -280,7 +335,7 @@ subroutine skip_empty_lines(blck,s_blck)
   enddo
 
 end subroutine skip_empty_lines
- 
+
 
 !--------------------------------------------------------------------------------------------------
 ! @brief skip file header
@@ -303,7 +358,7 @@ subroutine skip_file_header(blck,s_blck)
       call IO_error(708,ext_msg = line)
     end if
   end if
- 
+
 end subroutine skip_file_header
 
 
@@ -371,7 +426,7 @@ subroutine list_item_inline(blck,s_blck,inline,offset)
 
   character(len=:), allocatable :: line
   integer :: indent,indent_next
- 
+
   indent = indentDepth(blck(s_blck:),offset)
   line   = IO_rmComment(blck(s_blck:s_blck + index(blck(s_blck:),IO_EOL) - 2))
   inline = line(indent-offset+3:)
@@ -385,7 +440,7 @@ subroutine list_item_inline(blck,s_blck,inline,offset)
     indent_next = indentDepth(blck(s_blck:))
   enddo
 
-  if(scan(inline,",") > 0) inline = '"'//inline//'"' 
+  if(scan(inline,",") > 0) inline = '"'//inline//'"'
 
 end subroutine list_item_inline
 
@@ -737,7 +792,7 @@ end subroutine
 
 
 !--------------------------------------------------------------------------------------------------
-! @brief convert all block style YAML parts to flow style
+! @brief Convert all block style YAML parts to flow style.
 !--------------------------------------------------------------------------------------------------
 function to_flow(blck)
 
@@ -749,7 +804,7 @@ function to_flow(blck)
                                    s_flow, &                                                        !< start position in flow
                                    offset, &                                                        !< counts leading '- ' in nested lists
                                    end_line
- 
+
   allocate(character(len=len(blck)*2)::to_flow)
   s_flow = 1
   s_blck = 1
@@ -876,7 +931,7 @@ subroutine selfTest
 
   character(len=*), parameter :: flow = &
     '{a: ["b", {c: "d"}, "e"]}'
-  
+
   if( .not. to_flow(flow_multi)        == flow) error stop 'to_flow'
   end block multi_line_flow1
 
@@ -889,7 +944,7 @@ subroutine selfTest
     "[c,"//IO_EOL//&
     "d"//IO_EOL//&
     "e, f]}"//IO_EOL
- 
+
   character(len=*), parameter :: flow = &
     "[{a: {b: [c, d e, f]}}]"
 
@@ -921,5 +976,6 @@ subroutine selfTest
   end block basic_mixed
 
 end subroutine selfTest
+#endif
 
 end module YAML_parse

src/YAML_types.f90

@@ -119,7 +119,8 @@ module YAML_types
 
   type, extends(tNode), public :: tList
 
-    class(tItem), pointer  :: first => NULL()
+    class(tItem), pointer :: first => NULL(), &
+                             last => NULL()
 
     contains
     procedure :: asFormattedString => tList_asFormattedString
@@ -144,7 +145,7 @@ module YAML_types
   end type tDict
 
 
-  type :: tItem
+  type, public :: tItem
     character(len=:), allocatable :: key
     class(tNode),     pointer     :: node => NULL()
     class(tItem),     pointer     :: next => NULL()
@@ -1348,15 +1349,13 @@ subroutine tList_append(self,node)
   type(tItem), pointer :: item
 
   if (.not. associated(self%first)) then
-    allocate(self%first)
+    allocate(item)
-    item => self%first
+    self%first => item
+    self%last => item
   else
-    item => self%first
+    allocate(self%last%next)
-    do while (associated(item%next))
+    item => self%last%next
-      item => item%next
+    self%last => item
-    enddo
-    allocate(item%next)
-    item => item%next
   end if
 
   item%node => node

src/grid/spectral_utilities.f90

@@ -142,7 +142,7 @@ contains
 !> level chosen.
 !> Initializes FFTW.
 !--------------------------------------------------------------------------------------------------
-subroutine spectral_utilities_init
+subroutine spectral_utilities_init()
 
   PetscErrorCode :: err_PETSc
   integer        :: i, j, k, &
@@ -350,6 +350,8 @@ subroutine spectral_utilities_init
     allocate (gamma_hat(3,3,3,3,cells1Red,cells(2),cells3), source = cmplx(0.0_pReal,0.0_pReal,pReal))
   endif
 
+  call selfTest()
+
 end subroutine spectral_utilities_init
 
 
@@ -1146,4 +1148,41 @@ subroutine utilities_saveReferenceStiffness
 
 end subroutine utilities_saveReferenceStiffness
 
+
+!--------------------------------------------------------------------------------------------------
+!> @brief Check correctness of forward-backward transform.
+!--------------------------------------------------------------------------------------------------
+subroutine selfTest()
+
+  real(pReal), allocatable, dimension(:,:,:,:,:) :: tensorField_real_
+  real(pReal), allocatable, dimension(:,:,:,:) :: vectorField_real_
+  real(pReal), allocatable, dimension(:,:,:) :: scalarField_real_
+
+
+  call random_number(tensorField_real)
+  tensorField_real(1:3,1:3,cells(1)+1:cells1Red*2,:,:) = 0.0_pReal
+  tensorField_real_ = tensorField_real
+  call utilities_FFTtensorForward()
+  call utilities_FFTtensorBackward()
+  tensorField_real(1:3,1:3,cells(1)+1:cells1Red*2,:,:) = 0.0_pReal
+  if (maxval(abs(tensorField_real_ - tensorField_real))>5.0e-15_pReal) error stop 'tensorField'
+
+  call random_number(vectorField_real)
+  vectorField_real(1:3,cells(1)+1:cells1Red*2,:,:) = 0.0_pReal
+  vectorField_real_ = vectorField_real
+  call utilities_FFTvectorForward()
+  call utilities_FFTvectorBackward()
+  vectorField_real(1:3,cells(1)+1:cells1Red*2,:,:) = 0.0_pReal
+  if (maxval(abs(vectorField_real_ - vectorField_real))>5.0e-15_pReal) error stop 'vectorField'
+
+  call random_number(scalarField_real)
+  scalarField_real(cells(1)+1:cells1Red*2,:,:) = 0.0_pReal
+  scalarField_real_ = scalarField_real
+  call utilities_FFTscalarForward()
+  call utilities_FFTscalarBackward()
+  scalarField_real(cells(1)+1:cells1Red*2,:,:) = 0.0_pReal
+  if (maxval(abs(scalarField_real_ - scalarField_real))>5.0e-15_pReal) error stop 'scalarField'
+
+end subroutine selfTest
+
 end module spectral_utilities

src/material.f90

@@ -91,9 +91,13 @@ subroutine parse()
                            homogenizations, &
                            homogenization
 
+  class(tItem), pointer :: item
   integer, dimension(:), allocatable :: &
     counterPhase, &
-    counterHomogenization
+    counterHomogenization, &
+    ho_of
+  integer, dimension(:,:), allocatable :: ph_of
+  real(pReal), dimension(:,:), allocatable :: v_of
 
   real(pReal) :: v
   integer :: &
@@ -102,11 +106,14 @@ subroutine parse()
     co, ce, &
     ma
 
+
   materials       => config_material%get('material')
   phases          => config_material%get('phase')
   homogenizations => config_material%get('homogenization')
 
-  call sanityCheck(materials, homogenizations)
+
+  if (maxval(discretization_materialAt) > materials%length) &
+    call IO_error(155,ext_msg='More materials requested than found in material.yaml')
 
 #if defined (__GFORTRAN__)
   material_name_phase          = getKeys(phases)
@@ -123,6 +130,49 @@ subroutine parse()
   end do
   homogenization_maxNconstituents = maxval(homogenization_Nconstituents)
 
+  allocate(material_v(homogenization_maxNconstituents,discretization_Ncells),source=0.0_pReal)
+
+  allocate(material_O_0(materials%length))
+  allocate(material_F_i_0(materials%length))
+
+  allocate(ho_of(materials%length))
+  allocate(ph_of(materials%length,homogenization_maxNconstituents),source=-1)
+  allocate( v_of(materials%length,homogenization_maxNconstituents),source=0.0_pReal)
+
+  ! parse YAML structure
+  select type(materials)
+
+    class is(tList)
+
+      item => materials%first
+      do ma = 1, materials%length
+        material => item%node
+        ho_of(ma) = homogenizations%getIndex(material%get_asString('homogenization'))
+        constituents => material%get('constituents')
+
+        homogenization => homogenizations%get(ho_of(ma))
+        if (constituents%length /= homogenization%get_asInt('N_constituents')) call IO_error(148)
+
+        allocate(material_O_0(ma)%data(constituents%length))
+        allocate(material_F_i_0(ma)%data(1:3,1:3,constituents%length))
+
+        do co = 1, constituents%length
+          constituent => constituents%get(co)
+           v_of(ma,co) = constituent%get_asFloat('v')
+          ph_of(ma,co) = phases%getIndex(constituent%get_asString('phase'))
+
+          call material_O_0(ma)%data(co)%fromQuaternion(constituent%get_as1dFloat('O',requiredSize=4))
+          material_F_i_0(ma)%data(1:3,1:3,co) = constituent%get_as2dFloat('F_i',defaultVal=math_I3,requiredShape=[3,3])
+
+        end do
+        if (dNeq(sum(v_of(ma,:)),1.0_pReal,1.e-9_pReal)) call IO_error(153,ext_msg='constituent')
+
+        item => item%next
+      end do
+
+  end select
+
+
   allocate(counterPhase(phases%length),source=0)
   allocate(counterHomogenization(homogenizations%length),source=0)
 
@@ -132,12 +182,13 @@ subroutine parse()
   allocate(material_phaseID(homogenization_maxNconstituents,discretization_Ncells),source=0)
   allocate(material_phaseEntry(homogenization_maxNconstituents,discretization_Ncells),source=0)
 
-  allocate(material_v(homogenization_maxNconstituents,discretization_Ncells),source=0.0_pReal)
 
+  ! build mappings
   do el = 1, discretization_Nelems
-    material => materials%get(discretization_materialAt(el))
 
-    ho = homogenizations%getIndex(material%get_asString('homogenization'))
+    ma = discretization_materialAt(el)
+    ho = ho_of(ma)
+
     do ip = 1, discretization_nIPs
       ce = (el-1)*discretization_nIPs + ip
       material_homogenizationID(ce) = ho
@@ -145,13 +196,11 @@ subroutine parse()
       material_homogenizationEntry(ce) = counterHomogenization(ho)
     end do
 
-    constituents => material%get('constituents')
+    do co = 1, size(ph_of(ma,:)>0)
-    do co = 1, constituents%length
-      constituent => constituents%get(co)
 
-      v = constituent%get_asFloat('v')
+      v  =  v_of(ma,co)
+      ph = ph_of(ma,co)
 
-      ph = phases%getIndex(constituent%get_asString('phase'))
       do ip = 1, discretization_nIPs
         ce = (el-1)*discretization_nIPs + ip
         material_phaseID(co,ce) = ph
@@ -161,54 +210,11 @@ subroutine parse()
       end do
 
     end do
-    if (dNeq(sum(material_v(1:constituents%length,ce)),1.0_pReal,1.e-9_pReal)) &
-      call IO_error(153,ext_msg='constituent')
-
   end do
 
-  allocate(material_O_0(materials%length))
-  allocate(material_F_i_0(materials%length))
-
-  do ma = 1, materials%length
-    material     => materials%get(ma)
-    constituents => material%get('constituents')
-    allocate(material_O_0(ma)%data(constituents%length))
-    allocate(material_F_i_0(ma)%data(1:3,1:3,constituents%length))
-    do co = 1, constituents%length
-      constituent => constituents%get(co)
-      call material_O_0(ma)%data(co)%fromQuaternion(constituent%get_as1dFloat('O',requiredSize=4))
-      material_F_i_0(ma)%data(1:3,1:3,co) = constituent%get_as2dFloat('F_i',defaultVal=math_I3,requiredShape=[3,3])
-    enddo
- enddo
-
 end subroutine parse
 
 
-!--------------------------------------------------------------------------------------------------
-!> @brief Check if material.yaml is consistent and contains sufficient # of materials
-!--------------------------------------------------------------------------------------------------
-subroutine sanityCheck(materials,homogenizations)
-
-  class(tNode), intent(in) :: materials, &
-                              homogenizations
-
-  class(tNode), pointer :: material, &
-                           homogenization, &
-                           constituents
-  integer :: m
-
-  if (maxval(discretization_materialAt) > materials%length) &
-    call IO_error(155,ext_msg='More materials requested than found in material.yaml')
-
-  do m = 1, materials%length
-    material => materials%get(m)
-    constituents   => material%get('constituents')
-    homogenization => homogenizations%get(material%get_asString('homogenization'))
-    if (constituents%length /= homogenization%get_asInt('N_constituents')) call IO_error(148)
-  end do
-
-end subroutine sanityCheck
-
 #if defined (__GFORTRAN__)
 !--------------------------------------------------------------------------------------------------
 !> @brief %keys() is broken on gfortran

src/mesh/discretization_mesh.f90

@@ -100,7 +100,11 @@ subroutine discretization_mesh_init(restart)
   debug_element = config_debug%get_asInt('element',defaultVal=1)
   debug_ip      = config_debug%get_asInt('integrationpoint',defaultVal=1)
 
+#if (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR>16)
+  call DMPlexCreateFromFile(PETSC_COMM_WORLD,interface_geomFile,'n/a',PETSC_TRUE,globalMesh,err_PETSc)
+#else
   call DMPlexCreateFromFile(PETSC_COMM_WORLD,interface_geomFile,PETSC_TRUE,globalMesh,err_PETSc)
+#endif
   CHKERRQ(err_PETSc)
   call DMGetDimension(globalMesh,dimPlex,err_PETSc)
   CHKERRQ(err_PETSc)

src/prec.f90

@@ -48,7 +48,7 @@ contains
 !--------------------------------------------------------------------------------------------------
 !> @brief Report precision and do self test.
 !--------------------------------------------------------------------------------------------------
-subroutine prec_init
+subroutine prec_init()
 
   print'(/,1x,a)', '<<<+-  prec init  -+>>>'
 
@@ -60,7 +60,7 @@ subroutine prec_init
   print'(  a,e10.3)', '   epsilon value:      ',PREAL_EPSILON
   print'(  a,i3)',    '   decimal precision:  ',precision(0.0_pReal)
 
-  call selfTest
+  call selfTest()
 
 end subroutine prec_init
 
@@ -245,7 +245,7 @@ end function prec_bytesToC_INT64_T
 !--------------------------------------------------------------------------------------------------
 !> @brief Check correctness of some prec functions.
 !--------------------------------------------------------------------------------------------------
-subroutine selfTest
+subroutine selfTest()
 
   integer, allocatable, dimension(:) :: realloc_lhs_test
   real(pReal),   dimension(1) :: f

src/rotations.f90

@@ -658,11 +658,7 @@ function om2ax(om) result(ax)
   else
     call dgeev('N','V',3,om_,3,Wr,Wi,devNull,3,VR,3,work,size(work,1),ierr)
     if (ierr /= 0) error stop 'LAPACK error'
-#if defined(__GFORTRAN__) &&  __GNUC__<9
-    i = maxloc(merge(1,0,cEq(cmplx(Wr,Wi,pReal),cmplx(1.0_pReal,0.0_pReal,pReal),tol=1.0e-14_pReal)),dim=1)
-#else
     i = findloc(cEq(cmplx(Wr,Wi,pReal),cmplx(1.0_pReal,0.0_pReal,pReal),tol=1.0e-14_pReal),.true.,dim=1) !find eigenvalue (1,0)
-#endif
     if (i == 0) error stop 'om2ax conversion failed'
     ax(1:3) = VR(1:3,i)
     where (                dNeq0([om(2,3)-om(3,2), om(3,1)-om(1,3), om(1,2)-om(2,1)])) &
@@ -1427,10 +1423,6 @@ subroutine selfTest()
 
   do i = 1, 20
 
-#if defined(__GFORTRAN__) && __GNUC__<9
-    if(i<7) cycle
-#endif
-
     if(i==1) then
       qu = om2qu(math_I3)
     elseif(i==2) then

src/system_routines.f90

@@ -17,59 +17,67 @@ module system_routines
     getUserName, &
     signalterm_C, &
     signalusr1_C, &
-    signalusr2_C
+    signalusr2_C, &
+    f_c_string, &
+    free_C
 
 
   interface
 
-  function setCWD_C(cwd) bind(C)
+    function setCWD_C(cwd) bind(C)
-    use, intrinsic :: ISO_C_Binding, only: C_INT, C_CHAR
+      use, intrinsic :: ISO_C_Binding, only: C_INT, C_CHAR
 
-    integer(C_INT) :: setCWD_C
+      integer(C_INT) :: setCWD_C
-    character(kind=C_CHAR), dimension(*), intent(in) :: cwd
+      character(kind=C_CHAR), dimension(*), intent(in) :: cwd
-  end function setCWD_C
+    end function setCWD_C
 
-  subroutine getCWD_C(cwd, stat) bind(C)
+    subroutine getCWD_C(cwd, stat) bind(C)
-    use, intrinsic :: ISO_C_Binding, only: C_INT, C_CHAR
+      use, intrinsic :: ISO_C_Binding, only: C_INT, C_CHAR
-    use prec
+      use prec
 
-    character(kind=C_CHAR), dimension(pPathLen+1), intent(out) :: cwd                               ! NULL-terminated array
+      character(kind=C_CHAR), dimension(pPathLen+1), intent(out) :: cwd                             ! NULL-terminated array
-    integer(C_INT),                                intent(out) :: stat
+      integer(C_INT),                                intent(out) :: stat
-  end subroutine getCWD_C
+    end subroutine getCWD_C
 
-  subroutine getHostName_C(hostname, stat) bind(C)
+    subroutine getHostName_C(hostname, stat) bind(C)
-    use, intrinsic :: ISO_C_Binding, only: C_INT, C_CHAR
+      use, intrinsic :: ISO_C_Binding, only: C_INT, C_CHAR
-    use prec
+      use prec
 
-    character(kind=C_CHAR), dimension(pStringLen+1), intent(out) :: hostname                        ! NULL-terminated array
+      character(kind=C_CHAR), dimension(pStringLen+1), intent(out) :: hostname                        ! NULL-terminated array
-    integer(C_INT),                                  intent(out) :: stat
+      integer(C_INT),                                  intent(out) :: stat
-  end subroutine getHostName_C
+    end subroutine getHostName_C
 
-  subroutine getUserName_C(username, stat) bind(C)
+    subroutine getUserName_C(username, stat) bind(C)
-    use, intrinsic :: ISO_C_Binding, only: C_INT, C_CHAR
+      use, intrinsic :: ISO_C_Binding, only: C_INT, C_CHAR
-    use prec
+      use prec
 
-    character(kind=C_CHAR), dimension(pStringLen+1), intent(out) :: username                        ! NULL-terminated array
+      character(kind=C_CHAR), dimension(pStringLen+1), intent(out) :: username                        ! NULL-terminated array
-    integer(C_INT),                                  intent(out) :: stat
+      integer(C_INT),                                  intent(out) :: stat
-  end subroutine getUserName_C
+    end subroutine getUserName_C
 
-  subroutine signalterm_C(handler) bind(C)
+    subroutine signalterm_C(handler) bind(C)
-    use, intrinsic :: ISO_C_Binding, only: C_FUNPTR
+      use, intrinsic :: ISO_C_Binding, only: C_FUNPTR
 
-    type(C_FUNPTR), intent(in), value :: handler
+      type(C_FUNPTR), intent(in), value :: handler
-  end subroutine signalterm_C
+    end subroutine signalterm_C
 
-  subroutine signalusr1_C(handler) bind(C)
+    subroutine signalusr1_C(handler) bind(C)
-    use, intrinsic :: ISO_C_Binding, only: C_FUNPTR
+      use, intrinsic :: ISO_C_Binding, only: C_FUNPTR
 
-    type(C_FUNPTR), intent(in), value :: handler
+      type(C_FUNPTR), intent(in), value :: handler
-  end subroutine signalusr1_C
+    end subroutine signalusr1_C
 
-  subroutine signalusr2_C(handler) bind(C)
+    subroutine signalusr2_C(handler) bind(C)
-    use, intrinsic :: ISO_C_Binding, only: C_FUNPTR
+      use, intrinsic :: ISO_C_Binding, only: C_FUNPTR
+
+      type(C_FUNPTR), intent(in), value :: handler
+    end subroutine signalusr2_C
+
+    subroutine free_C(ptr) bind(C,name='free')
+      import c_ptr
+      type(c_ptr), value :: ptr
+    end subroutine free_C
 
-    type(C_FUNPTR), intent(in), value :: handler
-  end subroutine signalusr2_C
 
   end interface