diff --git a/.github/workflows/Fortran.yml b/.github/workflows/Fortran.yml index 5cc241e00..d93cdd27f 100644 --- a/.github/workflows/Fortran.yml +++ b/.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' diff --git a/.github/workflows/Python.yml b/.github/workflows/Python.yml index ca74c0016..3b6bb0df4 100644 --- a/.github/workflows/Python.yml +++ b/.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 diff --git a/CMakeLists.txt b/CMakeLists.txt index b4c405319..6b341dbbe 100644 --- a/CMakeLists.txt +++ b/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") -pkg_check_modules(PETSC REQUIRED PETSc>=3.12.0 PETSc<3.17.0) +set(ENV{PKG_CONFIG_PATH} "$ENV{PETSC_DIR}/$ENV{PETSC_ARCH}/lib/pkgconfig:$ENV{PKG_CONFIG_PATH}") +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_LINK_EXECUTABLE "${CMAKE_Fortran_LINK_EXECUTABLE} -o -L${PETSC_LIBRARY_DIRS} -lpetsc ${PETSC_EXTERNAL_LIB} -lz ${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} -o -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") diff --git a/PRIVATE b/PRIVATE index 459326e98..e3d9c0cc3 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit 459326e9840c843ade72b04cf28e50889c9779f1 +Subproject commit e3d9c0cc3ffa1435003f934a0f0e6c7d969a022a diff --git a/README b/README deleted file mode 100644 index 543de378f..000000000 --- a/README +++ /dev/null @@ -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 diff --git a/README.md b/README.md new file mode 100644 index 000000000..5b5edf677 --- /dev/null +++ b/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 + diff --git a/python/damask/VERSION b/python/damask/VERSION index 65787866a..f8dcddc79 100644 --- a/python/damask/VERSION +++ b/python/damask/VERSION @@ -1 +1 @@ -v3.0.0-alpha6-184-g1f98b04d4 +v3.0.0-alpha6-228-g758ad6072 diff --git a/python/damask/_grid.py b/python/damask/_grid.py index 6dc9879ce..d4a7ea02f 100644 --- a/python/damask/_grid.py +++ b/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³ diff --git a/python/damask/_result.py b/python/damask/_result.py index b8e025350..5ad9454e9 100644 --- a/python/damask/_result.py +++ b/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}) diff --git a/src/C_routines.c b/src/C_routines.c index 3d62a87c2..a25fde688 100644 --- a/src/C_routines.c +++ b/src/C_routines.c @@ -7,7 +7,11 @@ #include #include #include -#include "zlib.h" +#include + +#ifdef FYAML +#include +#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 diff --git a/src/DAMASK_interface.f90 b/src/DAMASK_interface.f90 index e958c8a29..1a6bb3ee6 100644 --- a/src/DAMASK_interface.f90 +++ b/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 diff --git a/src/HDF5_utilities.f90 b/src/HDF5_utilities.f90 index 93de4053b..86807507b 100644 --- a/src/HDF5_utilities.f90 +++ b/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 diff --git a/src/IO.f90 b/src/IO.f90 index 493970824..240c9e992 100644 --- a/src/IO.f90 +++ b/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 : pair' case (705) diff --git a/src/VTI.f90 b/src/VTI.f90 index d8b8f7778..a5918137e 100644 --- a/src/VTI.f90 +++ b/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, & - s + startPos, endPos cells = -1 diff --git a/src/YAML_parse.f90 b/src/YAML_parse.f90 index 8a3264cff..bfce06bc9 100644 --- a/src/YAML_parse.f90 +++ b/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 diff --git a/src/YAML_types.f90 b/src/YAML_types.f90 index d7ca71cbb..857b24bd0 100644 --- a/src/YAML_types.f90 +++ b/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) - item => self%first + allocate(item) + self%first => item + self%last => item else - item => self%first - do while (associated(item%next)) - item => item%next - enddo - allocate(item%next) - item => item%next + allocate(self%last%next) + item => self%last%next + self%last => item end if item%node => node diff --git a/src/grid/spectral_utilities.f90 b/src/grid/spectral_utilities.f90 index 525f68d7e..4375ec860 100644 --- a/src/grid/spectral_utilities.f90 +++ b/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 diff --git a/src/material.f90 b/src/material.f90 index 3ff280b0d..2258c40a1 100644 --- a/src/material.f90 +++ b/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, constituents%length - constituent => constituents%get(co) + do co = 1, size(ph_of(ma,:)>0) - 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 diff --git a/src/mesh/discretization_mesh.f90 b/src/mesh/discretization_mesh.f90 index 9baff52fb..70ee28343 100644 --- a/src/mesh/discretization_mesh.f90 +++ b/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) diff --git a/src/prec.f90 b/src/prec.f90 index d0753790e..73563d4f2 100644 --- a/src/prec.f90 +++ b/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 diff --git a/src/rotations.f90 b/src/rotations.f90 index 53f1fe976..7aa81fcab 100644 --- a/src/rotations.f90 +++ b/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 diff --git a/src/system_routines.f90 b/src/system_routines.f90 index 2eb0b7958..e0adf9dc0 100644 --- a/src/system_routines.f90 +++ b/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) - use, intrinsic :: ISO_C_Binding, only: C_INT, C_CHAR + function setCWD_C(cwd) bind(C) + use, intrinsic :: ISO_C_Binding, only: C_INT, C_CHAR - integer(C_INT) :: setCWD_C - character(kind=C_CHAR), dimension(*), intent(in) :: cwd - end function setCWD_C + integer(C_INT) :: setCWD_C + character(kind=C_CHAR), dimension(*), intent(in) :: cwd + end function setCWD_C - subroutine getCWD_C(cwd, stat) bind(C) - use, intrinsic :: ISO_C_Binding, only: C_INT, C_CHAR - use prec + subroutine getCWD_C(cwd, stat) bind(C) + use, intrinsic :: ISO_C_Binding, only: C_INT, C_CHAR + use prec - character(kind=C_CHAR), dimension(pPathLen+1), intent(out) :: cwd ! NULL-terminated array - integer(C_INT), intent(out) :: stat - end subroutine getCWD_C + character(kind=C_CHAR), dimension(pPathLen+1), intent(out) :: cwd ! NULL-terminated array + integer(C_INT), intent(out) :: stat + end subroutine getCWD_C - subroutine getHostName_C(hostname, stat) bind(C) - use, intrinsic :: ISO_C_Binding, only: C_INT, C_CHAR - use prec + subroutine getHostName_C(hostname, stat) bind(C) + use, intrinsic :: ISO_C_Binding, only: C_INT, C_CHAR + use prec - character(kind=C_CHAR), dimension(pStringLen+1), intent(out) :: hostname ! NULL-terminated array - integer(C_INT), intent(out) :: stat - end subroutine getHostName_C + character(kind=C_CHAR), dimension(pStringLen+1), intent(out) :: hostname ! NULL-terminated array + integer(C_INT), intent(out) :: stat + end subroutine getHostName_C - subroutine getUserName_C(username, stat) bind(C) - use, intrinsic :: ISO_C_Binding, only: C_INT, C_CHAR - use prec + subroutine getUserName_C(username, stat) bind(C) + use, intrinsic :: ISO_C_Binding, only: C_INT, C_CHAR + use prec - character(kind=C_CHAR), dimension(pStringLen+1), intent(out) :: username ! NULL-terminated array - integer(C_INT), intent(out) :: stat - end subroutine getUserName_C + character(kind=C_CHAR), dimension(pStringLen+1), intent(out) :: username ! NULL-terminated array + integer(C_INT), intent(out) :: stat + end subroutine getUserName_C - subroutine signalterm_C(handler) bind(C) - use, intrinsic :: ISO_C_Binding, only: C_FUNPTR + subroutine signalterm_C(handler) bind(C) + use, intrinsic :: ISO_C_Binding, only: C_FUNPTR - type(C_FUNPTR), intent(in), value :: handler - end subroutine signalterm_C + type(C_FUNPTR), intent(in), value :: handler + end subroutine signalterm_C - subroutine signalusr1_C(handler) bind(C) - use, intrinsic :: ISO_C_Binding, only: C_FUNPTR + subroutine signalusr1_C(handler) bind(C) + use, intrinsic :: ISO_C_Binding, only: C_FUNPTR - type(C_FUNPTR), intent(in), value :: handler - end subroutine signalusr1_C + type(C_FUNPTR), intent(in), value :: handler + end subroutine signalusr1_C - subroutine signalusr2_C(handler) bind(C) - use, intrinsic :: ISO_C_Binding, only: C_FUNPTR + subroutine signalusr2_C(handler) bind(C) + 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