From 977f61452b12d2ecd7cca847a93693f984f50e4a Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 6 Mar 2019 15:24:42 +0100 Subject: [PATCH 1/7] compiler-independent defintion of real and integer kinds real(8) does not neccessarily mean a real with 8 byte (but for gfortran and ifort it does) --- src/material.f90 | 1 - src/prec.f90 | 195 +++++++++++++++++++++++------------------------ 2 files changed, 95 insertions(+), 101 deletions(-) diff --git a/src/material.f90 b/src/material.f90 index 3ef22a199..edee30d20 100644 --- a/src/material.f90 +++ b/src/material.f90 @@ -15,7 +15,6 @@ module material tPlasticState, & tSourceState, & tHomogMapping, & - tPhaseMapping, & group_float, & group_int diff --git a/src/prec.f90 b/src/prec.f90 index ea539011f..8a625b97e 100644 --- a/src/prec.f90 +++ b/src/prec.f90 @@ -7,94 +7,89 @@ !> @brief setting precision for real and int type !-------------------------------------------------------------------------------------------------- module prec -! ToDo: use, intrinsic :: iso_fortran_env, only : I8 => int64, WP => real64 - implicit none - private -#if (FLOAT==8) - integer, parameter, public :: pReal = 8 !< floating point double precision (was selected_real_kind(15,300), number with 15 significant digits, up to 1e+-300) -#else - NO SUITABLE PRECISION FOR REAL SELECTED, STOPPING COMPILATION -#endif - -#if (INT==4) - integer, parameter, public :: pInt = 4 !< integer representation 32 bit (was selected_int_kind(9), number with at least up to +- 1e9) -#elif (INT==8) - integer, parameter, public :: pInt = 8 !< integer representation 64 bit (was selected_int_kind(12), number with at least up to +- 1e12) -#else - NO SUITABLE PRECISION FOR INTEGER SELECTED, STOPPING COMPILATION -#endif - - integer, parameter, public :: pStringLen = 256 !< default string lenth - integer, parameter, public :: pLongInt = 8 !< integer representation 64 bit (was selected_int_kind(12), number with at least up to +- 1e12) - real(pReal), parameter, public :: tol_math_check = 1.0e-8_pReal !< tolerance for internal math self-checks (rotation) - - integer(pInt), allocatable, dimension(:) :: realloc_lhs_test - - type, public :: group_float !< variable length datatype used for storage of state - real(pReal), dimension(:), pointer :: p - end type group_float - - type, public :: group_int - integer(pInt), dimension(:), pointer :: p - end type group_int - -!http://stackoverflow.com/questions/3948210/can-i-have-a-pointer-to-an-item-in-an-allocatable-array - type, public :: tState - integer(pInt) :: & - sizeState = 0_pInt, & !< size of state - sizeDotState = 0_pInt, & !< size of dot state, i.e. state(1:sizeDot) follows time evolution by dotState rates - offsetDeltaState = 0_pInt, & !< index offset of delta state - sizeDeltaState = 0_pInt, & !< size of delta state, i.e. state(offset+1:offset+sizeDelta) follows time evolution by deltaState increments - sizePostResults = 0_pInt !< size of output data - real(pReal), pointer, dimension(:), contiguous :: & - atolState - real(pReal), pointer, dimension(:,:), contiguous :: & ! a pointer is needed here because we might point to state/doState. However, they will never point to something, but are rather allocated and, hence, contiguous - state0, & - state, & !< state - dotState, & !< rate of state change - deltaState !< increment of state change - real(pReal), allocatable, dimension(:,:) :: & - partionedState0, & - subState0, & - previousDotState, & !< state rate of previous xxxx - previousDotState2, & !< state rate two xxxx ago - RK4dotState - real(pReal), allocatable, dimension(:,:,:) :: & - RKCK45dotState - end type - - type, extends(tState), public :: tPlasticState - integer(pInt) :: & - nSlip = 0_pInt , & - nTwin = 0_pInt, & - nTrans = 0_pInt - logical :: & - nonlocal = .false. - real(pReal), pointer, dimension(:,:) :: & - slipRate, & !< slip rate - accumulatedSlip !< accumulated plastic slip - end type - - type, public :: tSourceState - type(tState), dimension(:), allocatable :: p !< tState for each active source mechanism in a phase - end type + use, intrinsic :: IEEE_arithmetic, only:& + IEEE_selected_real_kind - type, public :: tHomogMapping - integer(pInt), pointer, dimension(:,:) :: p - end type + implicit none + private + ! https://software.intel.com/en-us/blogs/2017/03/27/doctor-fortran-in-it-takes-all-kinds - type, public :: tPhaseMapping - integer(pInt), pointer, dimension(:,:,:) :: p - end type + integer, parameter, public :: pReal = IEEE_selected_real_kind(15,307) !< number with 15 significant digits, up to 1e+-300 (typically 64 bit) +#if (INT==4) + integer, parameter, public :: pInt = selected_int_kind(9) !< number with at least up to +-1e9 (typically 32 bit) +#elif (INT==8) + integer, parameter, public :: pInt = selected_int_kind(18) !< number with at least up to +-1e18 (typically 64 bit) +#else + NO SUITABLE PRECISION FOR INTEGER SELECTED, STOPPING COMPILATION +#endif + integer, parameter, public :: pLongInt = selected_int_kind(18) !< number with at least up to +-1e18 (typically 64 bit) + integer, parameter, public :: pStringLen = 256 !< default string length - public :: & - prec_init, & - dEq, & - dEq0, & - cEq, & - dNeq, & - dNeq0, & - cNeq + real(pReal), parameter, public :: tol_math_check = 1.0e-8_pReal !< tolerance for internal math self-checks (rotation) + + + type, public :: group_float !< variable length datatype used for storage of state + real(pReal), dimension(:), pointer :: p + end type group_float + + type, public :: group_int + integer(pInt), dimension(:), pointer :: p + end type group_int + + ! http://stackoverflow.com/questions/3948210/can-i-have-a-pointer-to-an-item-in-an-allocatable-array + type, public :: tState + integer(pInt) :: & + sizeState = 0_pInt, & !< size of state + sizeDotState = 0_pInt, & !< size of dot state, i.e. state(1:sizeDot) follows time evolution by dotState rates + offsetDeltaState = 0_pInt, & !< index offset of delta state + sizeDeltaState = 0_pInt, & !< size of delta state, i.e. state(offset+1:offset+sizeDelta) follows time evolution by deltaState increments + sizePostResults = 0_pInt !< size of output data + real(pReal), pointer, dimension(:), contiguous :: & + atolState + real(pReal), pointer, dimension(:,:), contiguous :: & ! a pointer is needed here because we might point to state/doState. However, they will never point to something, but are rather allocated and, hence, contiguous + state0, & + state, & !< state + dotState, & !< rate of state change + deltaState !< increment of state change + real(pReal), allocatable, dimension(:,:) :: & + partionedState0, & + subState0, & + previousDotState, & !< state rate of previous xxxx + previousDotState2, & !< state rate two xxxx ago + RK4dotState + real(pReal), allocatable, dimension(:,:,:) :: & + RKCK45dotState + end type + + type, extends(tState), public :: tPlasticState + integer(pInt) :: & + nSlip = 0_pInt , & + nTwin = 0_pInt, & + nTrans = 0_pInt + logical :: & + nonlocal = .false. + real(pReal), pointer, dimension(:,:) :: & + slipRate, & !< slip rate + accumulatedSlip !< accumulated plastic slip + end type + + type, public :: tSourceState + type(tState), dimension(:), allocatable :: p !< tState for each active source mechanism in a phase + end type + + type, public :: tHomogMapping + integer(pInt), pointer, dimension(:,:) :: p + end type + + + public :: & + prec_init, & + dEq, & + dEq0, & + cEq, & + dNeq, & + dNeq0, & + cNeq contains @@ -103,24 +98,24 @@ contains !> @brief reporting precision !-------------------------------------------------------------------------------------------------- subroutine prec_init -#if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800 - use, intrinsic :: iso_fortran_env, only: & - compiler_version, & - compiler_options -#endif - implicit none - external :: & - quit + implicit none + integer(pInt), allocatable, dimension(:) :: realloc_lhs_test - write(6,'(/,a)') ' <<<+- prec init -+>>>' -#include "compilation_info.f90" - write(6,'(a,i3)') ' Bytes for pReal: ',pReal - write(6,'(a,i3)') ' Bytes for pInt: ',pInt - write(6,'(a,i3)') ' Bytes for pLongInt: ',pLongInt + external :: & + quit - realloc_lhs_test = [1_pInt,2_pInt] - if (realloc_lhs_test(2)/=2_pInt) call quit(9000) + write(6,'(/,a)') ' <<<+- prec init -+>>>' + + write(6,'(a,i3)') ' Size of integer in bit: ',bit_size(0_pInt) + write(6,'(a,i19)') ' Maximum value: ',huge(0_pInt) + write(6,'(/,a,i3)') ' Size of float in bit: ',storage_size(0.0_pReal) + write(6,'(a,e10.3)') ' Maximum value: ',huge(0.0_pReal) + write(6,'(a,e10.3)') ' Minimum value: ',tiny(0.0_pReal) + write(6,'(a,i3)') ' Decimal precision: ',precision(0.0_pReal) + + realloc_lhs_test = [1_pInt,2_pInt] + if (realloc_lhs_test(2)/=2_pInt) call quit(9000) end subroutine prec_init From c9e7311b428dfd01488328be4b6bdbdc5650af02 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 6 Mar 2019 15:47:48 +0100 Subject: [PATCH 2/7] no need to use pInt here --- src/DAMASK_interface.f90 | 667 ++++++++++++++++++++------------------- 1 file changed, 334 insertions(+), 333 deletions(-) diff --git a/src/DAMASK_interface.f90 b/src/DAMASK_interface.f90 index 630b5b921..b2930510b 100644 --- a/src/DAMASK_interface.f90 +++ b/src/DAMASK_interface.f90 @@ -6,33 +6,33 @@ !> @brief Interfacing between the PETSc-based solvers and the material subroutines provided !! by DAMASK !> @details Interfacing between the PETSc-based solvers and the material subroutines provided -!> by DAMASK. Interpretating the command line arguments to get load case, geometry file, +!> by DAMASK. Interpretating the command line arguments to get load case, geometry file, !> and working directory. !-------------------------------------------------------------------------------------------------- module DAMASK_interface - use prec, only: & - pInt - implicit none - private - logical, public, protected :: SIGUSR1,SIGUSR2 - integer(pInt), public, protected :: & - interface_restartInc = 0_pInt !< Increment at which calculation starts - character(len=1024), public, protected :: & - geometryFile = '', & !< parameter given for geometry file - loadCaseFile = '' !< parameter given for load case file + implicit none + private + logical, public, protected :: & + SIGUSR1, & !< user-defined signal 1 + SIGUSR2 !< user-defined signal 2 + integer, public, protected :: & + interface_restartInc = 0 !< Increment at which calculation starts + character(len=1024), public, protected :: & + geometryFile = '', & !< parameter given for geometry file + loadCaseFile = '' !< parameter given for load case file - public :: & - getSolverJobName, & - DAMASK_interface_init - private :: & - setWorkingDirectory, & - getGeometryFile, & - getLoadCaseFile, & - rectifyPath, & - makeRelativePath, & - IIO_stringValue, & - IIO_intValue, & - IIO_stringPos + public :: & + getSolverJobName, & + DAMASK_interface_init + private :: & + setWorkingDirectory, & + getGeometryFile, & + getLoadCaseFile, & + rectifyPath, & + makeRelativePath, & + IIO_stringValue, & + IIO_intValue, & + IIO_stringPos contains !-------------------------------------------------------------------------------------------------- @@ -40,10 +40,17 @@ contains !! information on computation to screen !-------------------------------------------------------------------------------------------------- subroutine DAMASK_interface_init() - use, intrinsic :: & - iso_fortran_env - use :: & - iso_c_binding + use, intrinsic :: & + iso_fortran_env + use, intrinsic :: & + iso_c_binding + use PETScSys + use system_routines, only: & + signalusr1_C, & + signalusr2_C, & + getHostName, & + getCWD + #include #if defined(__GFORTRAN__) && __GNUC__ < 5 =================================================================================================== @@ -59,13 +66,13 @@ subroutine DAMASK_interface_init() #if defined(__INTEL_COMPILER) && __INTEL_COMPILER < 1600 =================================================================================================== - 16.0 16.0 16.0 16.0 16.0 16.0 16.0 16.0 16.0 16.0 16.0 16.0 16.0 16.0 16.0 16.0 16.0 16.0 16.0 + 16.0 16.0 16.0 16.0 16.0 16.0 16.0 16.0 16.0 16.0 16.0 16.0 16.0 16.0 16.0 16.0 16.0 16.0 16.0 =================================================================================================== ================== THIS VERSION OF DAMASK REQUIRES ifort > 16.0 ================================ ====================== THIS VERSION OF DAMASK REQUIRES ifort > 16.0 =========================== ========================= THIS VERSION OF DAMASK REQUIRES ifort > 16.0 ======================== =================================================================================================== - 16.0 16.0 16.0 16.0 16.0 16.0 16.0 16.0 16.0 16.0 16.0 16.0 16.0 16.0 16.0 16.0 16.0 16.0 16.0 + 16.0 16.0 16.0 16.0 16.0 16.0 16.0 16.0 16.0 16.0 16.0 16.0 16.0 16.0 16.0 16.0 16.0 16.0 16.0 =================================================================================================== #endif @@ -81,175 +88,168 @@ subroutine DAMASK_interface_init() =================================================================================================== #endif - use PETScSys - use system_routines, only: & - signalusr1_C, & - signalusr2_C, & - getHostName, & - getCWD - - implicit none - character(len=1024) :: & - commandLine, & !< command line call as string - loadcaseArg = '', & !< -l argument given to the executable - geometryArg = '', & !< -g argument given to the executable - workingDirArg = '', & !< -w argument given to the executable - userName !< name of user calling the executable - integer :: & - i, & + implicit none + character(len=1024) :: & + commandLine, & !< command line call as string + loadcaseArg = '', & !< -l argument given to the executable + geometryArg = '', & !< -g argument given to the executable + workingDirArg = '', & !< -w argument given to the executable + userName !< name of user calling the executable + integer :: & + i, & #ifdef _OPENMP - threadLevel, & + threadLevel, & #endif - worldrank = 0, & - worldsize = 0 - integer, allocatable, dimension(:) :: & - chunkPos - integer, dimension(8) :: & - dateAndTime ! type default integer - PetscErrorCode :: ierr - external :: & - quit + worldrank = 0, & + worldsize = 0 + integer, allocatable, dimension(:) :: & + chunkPos + integer, dimension(8) :: & + dateAndTime + PetscErrorCode :: ierr + external :: & + quit - open(6, encoding='UTF-8') ! for special characters in output + open(6, encoding='UTF-8') ! for special characters in output !-------------------------------------------------------------------------------------------------- ! PETSc Init #ifdef _OPENMP - ! If openMP is enabled, check if the MPI libary supports it and initialize accordingly. - ! Otherwise, the first call to PETSc will do the initialization. - call MPI_Init_Thread(MPI_THREAD_FUNNELED,threadLevel,ierr);CHKERRQ(ierr) - if (threadLevel>>' - write(6,'(/,a)') ' Roters et al., Computational Materials Science 158, 2018, 420-478' - write(6,'(a,/)') ' https://doi.org/10.1016/j.commatsci.2018.04.030' + call date_and_time(values = dateAndTime) + write(6,'(/,a)') ' <<<+- DAMASK_interface init -+>>>' + write(6,'(/,a)') ' Roters et al., Computational Materials Science 158, 2018, 420-478' + write(6,'(a,/)') ' https://doi.org/10.1016/j.commatsci.2018.04.030' - write(6,'(a,/)') ' Version: '//DAMASKVERSION + write(6,'(a,/)') ' Version: '//DAMASKVERSION -! https://github.com/jeffhammond/HPCInfo/blob/master/docs/Preprocessor-Macros.md + ! https://github.com/jeffhammond/HPCInfo/blob/master/docs/Preprocessor-Macros.md #if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800 - write(6,*) 'Compiled with: ', compiler_version() - write(6,*) 'Compiler options: ', compiler_options() + write(6,*) 'Compiled with: ', compiler_version() + write(6,*) 'Compiler options: ', compiler_options() #elif defined(__INTEL_COMPILER) - write(6,'(a,i4.4,a,i8.8)') ' Compiled with Intel fortran version :', __INTEL_COMPILER,& - ', build date :', __INTEL_COMPILER_BUILD_DATE + write(6,'(a,i4.4,a,i8.8)') ' Compiled with Intel fortran version :', __INTEL_COMPILER,& + ', build date :', __INTEL_COMPILER_BUILD_DATE #elif defined(__PGI) - write(6,'(a,i4.4,a,i8.8)') ' Compiled with PGI fortran version :', __PGIC__,& - '.', __PGIC_MINOR__ + write(6,'(a,i4.4,a,i8.8)') ' Compiled with PGI fortran version :', __PGIC__,& + '.', __PGIC_MINOR__ #endif - write(6,*) 'Compiled on ', __DATE__,' at ',__TIME__ + write(6,*) 'Compiled on ', __DATE__,' at ',__TIME__ - write(6,'(a,2(i2.2,a),i4.4)') ' Date: ',dateAndTime(3),'/',dateAndTime(2),'/', dateAndTime(1) - write(6,'(a,2(i2.2,a),i2.2)') ' Time: ',dateAndTime(5),':', dateAndTime(6),':', dateAndTime(7) - - call get_command(commandLine) - chunkPos = IIO_stringPos(commandLine) - do i = 2_pInt, chunkPos(1) - select case(IIO_stringValue(commandLine,chunkPos,i)) ! extract key - case ('-h','--help') - write(6,'(a)') ' #######################################################################' - write(6,'(a)') ' DAMASK Command Line Interface:' - write(6,'(a)') ' For PETSc-based solvers for the Düsseldorf Advanced Material Simulation Kit' - write(6,'(a,/)')' #######################################################################' - write(6,'(a,/)')' Valid command line switches:' - write(6,'(a)') ' --geom (-g, --geometry)' - write(6,'(a)') ' --load (-l, --loadcase)' - write(6,'(a)') ' --workingdir (-w, --wd, --workingdirectory, -d, --directory)' - write(6,'(a)') ' --restart (-r, --rs)' - write(6,'(a)') ' --help (-h)' - write(6,'(/,a)')' -----------------------------------------------------------------------' - write(6,'(a)') ' Mandatory arguments:' - write(6,'(/,a)')' --geom PathToGeomFile/NameOfGeom' - write(6,'(a)') ' Specifies the location of the geometry definition file.' - write(6,'(/,a)')' --load PathToLoadFile/NameOfLoadFile' - write(6,'(a)') ' Specifies the location of the load case definition file.' - write(6,'(/,a)')' -----------------------------------------------------------------------' - write(6,'(a)') ' Optional arguments:' - write(6,'(/,a)')' --workingdirectory PathToWorkingDirectory' - write(6,'(a)') ' Specifies the working directory and overwrites the default ./' - write(6,'(a)') ' Make sure the file "material.config" exists in the working' - write(6,'(a)') ' directory.' - write(6,'(a)') ' For further configuration place "numerics.config"' - write(6,'(a)')' and "debug.config" in that directory.' - write(6,'(/,a)')' --restart XX' - write(6,'(a)') ' Reads in increment XX and continues with calculating' - write(6,'(a)') ' increment XX+1 based on this.' - write(6,'(a)') ' Appends to existing results file' - write(6,'(a)') ' "NameOfGeom_NameOfLoadFile".' - write(6,'(a)') ' Works only if the restart information for increment XX' - write(6,'(a)') ' is available in the working directory.' - write(6,'(/,a)')' -----------------------------------------------------------------------' - write(6,'(a)') ' Help:' - write(6,'(/,a)')' --help' - write(6,'(a,/)')' Prints this message and exits' - call quit(0_pInt) ! normal Termination - case ('-l', '--load', '--loadcase') - if ( i < chunkPos(1)) loadcaseArg = trim(IIO_stringValue(commandLine,chunkPos,i+1_pInt)) - case ('-g', '--geom', '--geometry') - if (i < chunkPos(1)) geometryArg = trim(IIO_stringValue(commandLine,chunkPos,i+1_pInt)) - case ('-w', '-d', '--wd', '--directory', '--workingdir', '--workingdirectory') - if (i < chunkPos(1)) workingDirArg = trim(IIO_stringValue(commandLine,chunkPos,i+1_pInt)) - case ('-r', '--rs', '--restart') - if (i < chunkPos(1)) then - interface_restartInc = IIO_IntValue(commandLine,chunkPos,i+1_pInt) - endif - end select - enddo + write(6,'(a,2(i2.2,a),i4.4)') ' Date: ',dateAndTime(3),'/',dateAndTime(2),'/', dateAndTime(1) + write(6,'(a,2(i2.2,a),i2.2)') ' Time: ',dateAndTime(5),':', dateAndTime(6),':', dateAndTime(7) - if (len_trim(loadcaseArg) == 0 .or. len_trim(geometryArg) == 0) then - write(6,'(a)') ' Please specify geometry AND load case (-h for help)' - call quit(1_pInt) - endif + call get_command(commandLine) + chunkPos = IIO_stringPos(commandLine) + do i = 2, chunkPos(1) + select case(IIO_stringValue(commandLine,chunkPos,i)) ! extract key + case ('-h','--help') + write(6,'(a)') ' #######################################################################' + write(6,'(a)') ' DAMASK Command Line Interface:' + write(6,'(a)') ' For PETSc-based solvers for the Düsseldorf Advanced Material Simulation Kit' + write(6,'(a,/)')' #######################################################################' + write(6,'(a,/)')' Valid command line switches:' + write(6,'(a)') ' --geom (-g, --geometry)' + write(6,'(a)') ' --load (-l, --loadcase)' + write(6,'(a)') ' --workingdir (-w, --wd, --workingdirectory, -d, --directory)' + write(6,'(a)') ' --restart (-r, --rs)' + write(6,'(a)') ' --help (-h)' + write(6,'(/,a)')' -----------------------------------------------------------------------' + write(6,'(a)') ' Mandatory arguments:' + write(6,'(/,a)')' --geom PathToGeomFile/NameOfGeom' + write(6,'(a)') ' Specifies the location of the geometry definition file.' + write(6,'(/,a)')' --load PathToLoadFile/NameOfLoadFile' + write(6,'(a)') ' Specifies the location of the load case definition file.' + write(6,'(/,a)')' -----------------------------------------------------------------------' + write(6,'(a)') ' Optional arguments:' + write(6,'(/,a)')' --workingdirectory PathToWorkingDirectory' + write(6,'(a)') ' Specifies the working directory and overwrites the default ./' + write(6,'(a)') ' Make sure the file "material.config" exists in the working' + write(6,'(a)') ' directory.' + write(6,'(a)') ' For further configuration place "numerics.config"' + write(6,'(a)')' and "debug.config" in that directory.' + write(6,'(/,a)')' --restart XX' + write(6,'(a)') ' Reads in increment XX and continues with calculating' + write(6,'(a)') ' increment XX+1 based on this.' + write(6,'(a)') ' Appends to existing results file' + write(6,'(a)') ' "NameOfGeom_NameOfLoadFile".' + write(6,'(a)') ' Works only if the restart information for increment XX' + write(6,'(a)') ' is available in the working directory.' + write(6,'(/,a)')' -----------------------------------------------------------------------' + write(6,'(a)') ' Help:' + write(6,'(/,a)')' --help' + write(6,'(a,/)')' Prints this message and exits' + call quit(0) ! normal Termination + case ('-l', '--load', '--loadcase') + if ( i < chunkPos(1)) loadcaseArg = trim(IIO_stringValue(commandLine,chunkPos,i+1)) + case ('-g', '--geom', '--geometry') + if (i < chunkPos(1)) geometryArg = trim(IIO_stringValue(commandLine,chunkPos,i+1)) + case ('-w', '-d', '--wd', '--directory', '--workingdir', '--workingdirectory') + if (i < chunkPos(1)) workingDirArg = trim(IIO_stringValue(commandLine,chunkPos,i+1)) + case ('-r', '--rs', '--restart') + if (i < chunkPos(1)) then + interface_restartInc = IIO_IntValue(commandLine,chunkPos,i+1) + endif + end select + enddo - if (len_trim(workingDirArg) > 0) call setWorkingDirectory(trim(workingDirArg)) - geometryFile = getGeometryFile(geometryArg) - loadCaseFile = getLoadCaseFile(loadCaseArg) + if (len_trim(loadcaseArg) == 0 .or. len_trim(geometryArg) == 0) then + write(6,'(a)') ' Please specify geometry AND load case (-h for help)' + call quit(1) + endif - call get_environment_variable('USER',userName) - ! ToDo: https://stackoverflow.com/questions/8953424/how-to-get-the-username-in-c-c-in-linux - write(6,'(/,a,i4.1)') ' MPI processes: ',worldsize - write(6,'(a,a)') ' Host name: ', trim(getHostName()) - write(6,'(a,a)') ' User name: ', trim(userName) + if (len_trim(workingDirArg) > 0) call setWorkingDirectory(trim(workingDirArg)) + geometryFile = getGeometryFile(geometryArg) + loadCaseFile = getLoadCaseFile(loadCaseArg) - write(6,'(/a,a)') ' Command line call: ', trim(commandLine) - if (len(trim(workingDirArg)) > 0) & - write(6,'(a,a)') ' Working dir argument: ', trim(workingDirArg) - write(6,'(a,a)') ' Geometry argument: ', trim(geometryArg) - write(6,'(a,a)') ' Loadcase argument: ', trim(loadcaseArg) - write(6,'(a,a)') ' Working directory: ', trim(getCWD()) - write(6,'(a,a)') ' Geometry file: ', trim(geometryFile) - write(6,'(a,a)') ' Loadcase file: ', trim(loadCaseFile) - write(6,'(a,a)') ' Solver job name: ', trim(getSolverJobName()) - if (interface_restartInc > 0_pInt) & - write(6,'(a,i6.6)') ' Restart from increment: ', interface_restartInc + call get_environment_variable('USER',userName) + ! ToDo: https://stackoverflow.com/questions/8953424/how-to-get-the-username-in-c-c-in-linux + write(6,'(/,a,i4.1)') ' MPI processes: ',worldsize + write(6,'(a,a)') ' Host name: ', trim(getHostName()) + write(6,'(a,a)') ' User name: ', trim(userName) - call signalusr1_c(c_funloc(setSIGUSR1)) - call signalusr2_c(c_funloc(setSIGUSR2)) - SIGUSR1 = .false. - SIGUSR2 = .false. + write(6,'(/a,a)') ' Command line call: ', trim(commandLine) + if (len(trim(workingDirArg)) > 0) & + write(6,'(a,a)') ' Working dir argument: ', trim(workingDirArg) + write(6,'(a,a)') ' Geometry argument: ', trim(geometryArg) + write(6,'(a,a)') ' Load case argument: ', trim(loadcaseArg) + write(6,'(a,a)') ' Working directory: ', trim(getCWD()) + write(6,'(a,a)') ' Geometry file: ', trim(geometryFile) + write(6,'(a,a)') ' Loadcase file: ', trim(loadCaseFile) + write(6,'(a,a)') ' Solver job name: ', trim(getSolverJobName()) + if (interface_restartInc > 0) & + write(6,'(a,i6.6)') ' Restart from increment: ', interface_restartInc + + call signalusr1_c(c_funloc(setSIGUSR1)) + call signalusr2_c(c_funloc(setSIGUSR2)) + SIGUSR1 = .false. + SIGUSR2 = .false. end subroutine DAMASK_interface_init @@ -260,29 +260,29 @@ end subroutine DAMASK_interface_init !! possibly converting relative arguments to absolut path !-------------------------------------------------------------------------------------------------- subroutine setWorkingDirectory(workingDirectoryArg) - use system_routines, only: & - getCWD, & - setCWD + use system_routines, only: & + getCWD, & + setCWD - implicit none - character(len=*), intent(in) :: workingDirectoryArg !< working directory argument - character(len=1024) :: workingDirectory !< working directory argument - logical :: error - external :: quit + implicit none + character(len=*), intent(in) :: workingDirectoryArg !< working directory argument + character(len=1024) :: workingDirectory !< working directory argument + logical :: error + external :: quit - absolutePath: if (workingDirectoryArg(1:1) == '/') then - workingDirectory = workingDirectoryArg - else absolutePath - workingDirectory = getCWD() - workingDirectory = trim(workingDirectory)//'/'//workingDirectoryArg - endif absolutePath + absolutePath: if (workingDirectoryArg(1:1) == '/') then + workingDirectory = workingDirectoryArg + else absolutePath + workingDirectory = getCWD() + workingDirectory = trim(workingDirectory)//'/'//workingDirectoryArg + endif absolutePath - workingDirectory = trim(rectifyPath(workingDirectory)) - error = setCWD(trim(workingDirectory)) - if(error) then - write(6,'(a20,a,a16)') ' working directory "',trim(workingDirectory),'" does not exist' - call quit(1_pInt) - endif + workingDirectory = trim(rectifyPath(workingDirectory)) + error = setCWD(trim(workingDirectory)) + if(error) then + write(6,'(a20,a,a16)') ' Working directory "',trim(workingDirectory),'" does not exist' + call quit(1) + endif end subroutine setWorkingDirectory @@ -292,22 +292,22 @@ end subroutine setWorkingDirectory !-------------------------------------------------------------------------------------------------- character(len=1024) function getSolverJobName() - implicit none - integer :: posExt,posSep - character(len=1024) :: tempString + implicit none + integer :: posExt,posSep + character(len=1024) :: tempString - tempString = geometryFile - posExt = scan(tempString,'.',back=.true.) - posSep = scan(tempString,'/',back=.true.) + tempString = geometryFile + posExt = scan(tempString,'.',back=.true.) + posSep = scan(tempString,'/',back=.true.) - getSolverJobName = tempString(posSep+1:posExt-1) + getSolverJobName = tempString(posSep+1:posExt-1) - tempString = loadCaseFile - posExt = scan(tempString,'.',back=.true.) - posSep = scan(tempString,'/',back=.true.) + tempString = loadCaseFile + posExt = scan(tempString,'.',back=.true.) + posSep = scan(tempString,'/',back=.true.) - getSolverJobName = trim(getSolverJobName)//'_'//tempString(posSep+1:posExt-1) + getSolverJobName = trim(getSolverJobName)//'_'//tempString(posSep+1:posExt-1) end function getSolverJobName @@ -316,23 +316,23 @@ end function getSolverJobName !> @brief basename of geometry file with extension from command line arguments !-------------------------------------------------------------------------------------------------- character(len=1024) function getGeometryFile(geometryParameter) - use system_routines, only: & - getCWD + use system_routines, only: & + getCWD - implicit none - character(len=1024), intent(in) :: geometryParameter - logical :: file_exists - external :: quit + implicit none + character(len=1024), intent(in) :: geometryParameter + logical :: file_exists + external :: quit - getGeometryFile = trim(geometryParameter) - if (scan(getGeometryFile,'/') /= 1) getGeometryFile = trim(getCWD())//'/'//trim(getGeometryFile) - getGeometryFile = makeRelativePath(trim(getCWD()), getGeometryFile) + getGeometryFile = trim(geometryParameter) + if (scan(getGeometryFile,'/') /= 1) getGeometryFile = trim(getCWD())//'/'//trim(getGeometryFile) + getGeometryFile = makeRelativePath(trim(getCWD()), getGeometryFile) - inquire(file=trim(getGeometryFile), exist=file_exists) - if (.not. file_exists) then - write(6,'(a)') ' Geometry file does not exists ('//trim(getGeometryFile)//')' - call quit(1_pInt) - endif + inquire(file=trim(getGeometryFile), exist=file_exists) + if (.not. file_exists) then + write(6,'(a)') ' Geometry file does not exists ('//trim(getGeometryFile)//')' + call quit(1) + endif end function getGeometryFile @@ -341,23 +341,23 @@ end function getGeometryFile !> @brief relative path of loadcase from command line arguments !-------------------------------------------------------------------------------------------------- character(len=1024) function getLoadCaseFile(loadCaseParameter) - use system_routines, only: & - getCWD + use system_routines, only: & + getCWD - implicit none - character(len=1024), intent(in) :: loadCaseParameter - logical :: file_exists - external :: quit + implicit none + character(len=1024), intent(in) :: loadCaseParameter + logical :: file_exists + external :: quit - getLoadCaseFile = trim(loadCaseParameter) - if (scan(getLoadCaseFile,'/') /= 1) getLoadCaseFile = trim(getCWD())//'/'//trim(getLoadCaseFile) - getLoadCaseFile = makeRelativePath(trim(getCWD()), getLoadCaseFile) + getLoadCaseFile = trim(loadCaseParameter) + if (scan(getLoadCaseFile,'/') /= 1) getLoadCaseFile = trim(getCWD())//'/'//trim(getLoadCaseFile) + getLoadCaseFile = makeRelativePath(trim(getCWD()), getLoadCaseFile) - inquire(file=trim(getLoadCaseFile), exist=file_exists) - if (.not. file_exists) then - write(6,'(a)') ' Geometry file does not exists ('//trim(getLoadCaseFile)//')' - call quit(1_pInt) - endif + inquire(file=trim(getLoadCaseFile), exist=file_exists) + if (.not. file_exists) then + write(6,'(a)') ' Load case file does not exists ('//trim(getLoadCaseFile)//')' + call quit(1) + endif end function getLoadCaseFile @@ -368,162 +368,163 @@ end function getLoadCaseFile !-------------------------------------------------------------------------------------------------- function rectifyPath(path) - implicit none - character(len=*) :: path - character(len=1024) :: rectifyPath - integer :: i,j,k,l ! no pInt + implicit none + character(len=*) :: path + character(len=1024) :: rectifyPath + integer :: i,j,k,l !-------------------------------------------------------------------------------------------------- ! remove /./ from path - rectifyPath = trim(path) - l = len_trim(rectifyPath) - do i = l,3,-1 - if (rectifyPath(i-2:i) == '/./') rectifyPath(i-1:l) = rectifyPath(i+1:l)//' ' - enddo + rectifyPath = trim(path) + l = len_trim(rectifyPath) + do i = l,3,-1 + if (rectifyPath(i-2:i) == '/./') rectifyPath(i-1:l) = rectifyPath(i+1:l)//' ' + enddo !-------------------------------------------------------------------------------------------------- ! remove // from path - l = len_trim(rectifyPath) - do i = l,2,-1 - if (rectifyPath(i-1:i) == '//') rectifyPath(i-1:l) = rectifyPath(i:l)//' ' - enddo + l = len_trim(rectifyPath) + do i = l,2,-1 + if (rectifyPath(i-1:i) == '//') rectifyPath(i-1:l) = rectifyPath(i:l)//' ' + enddo !-------------------------------------------------------------------------------------------------- ! remove ../ and corresponding directory from rectifyPath - l = len_trim(rectifyPath) - i = index(rectifyPath(i:l),'../') - j = 0 - do while (i > j) - j = scan(rectifyPath(1:i-2),'/',back=.true.) - rectifyPath(j+1:l) = rectifyPath(i+3:l)//repeat(' ',2+i-j) - if (rectifyPath(j+1:j+1) == '/') then !search for '//' that appear in case of XXX/../../XXX - k = len_trim(rectifyPath) - rectifyPath(j+1:k-1) = rectifyPath(j+2:k) - rectifyPath(k:k) = ' ' - endif - i = j+index(rectifyPath(j+1:l),'../') - enddo - if(len_trim(rectifyPath) == 0) rectifyPath = '/' + l = len_trim(rectifyPath) + i = index(rectifyPath(i:l),'../') + j = 0 + do while (i > j) + j = scan(rectifyPath(1:i-2),'/',back=.true.) + rectifyPath(j+1:l) = rectifyPath(i+3:l)//repeat(' ',2+i-j) + if (rectifyPath(j+1:j+1) == '/') then !search for '//' that appear in case of XXX/../../XXX + k = len_trim(rectifyPath) + rectifyPath(j+1:k-1) = rectifyPath(j+2:k) + rectifyPath(k:k) = ' ' + endif + i = j+index(rectifyPath(j+1:l),'../') + enddo + if(len_trim(rectifyPath) == 0) rectifyPath = '/' end function rectifyPath - + !-------------------------------------------------------------------------------------------------- !> @brief relative path from absolute a to absolute b !-------------------------------------------------------------------------------------------------- character(len=1024) function makeRelativePath(a,b) - implicit none - character (len=*), intent(in) :: a,b - character (len=1024) :: a_cleaned,b_cleaned - integer :: i,posLastCommonSlash,remainingSlashes !no pInt + implicit none + character (len=*), intent(in) :: a,b + character (len=1024) :: a_cleaned,b_cleaned + integer :: i,posLastCommonSlash,remainingSlashes - posLastCommonSlash = 0 - remainingSlashes = 0 - a_cleaned = rectifyPath(trim(a)//'/') - b_cleaned = rectifyPath(b) + posLastCommonSlash = 0 + remainingSlashes = 0 + a_cleaned = rectifyPath(trim(a)//'/') + b_cleaned = rectifyPath(b) - do i = 1, min(1024,len_trim(a_cleaned),len_trim(rectifyPath(b_cleaned))) - if (a_cleaned(i:i) /= b_cleaned(i:i)) exit - if (a_cleaned(i:i) == '/') posLastCommonSlash = i - enddo - do i = posLastCommonSlash+1,len_trim(a_cleaned) - if (a_cleaned(i:i) == '/') remainingSlashes = remainingSlashes + 1 - enddo + do i = 1, min(1024,len_trim(a_cleaned),len_trim(rectifyPath(b_cleaned))) + if (a_cleaned(i:i) /= b_cleaned(i:i)) exit + if (a_cleaned(i:i) == '/') posLastCommonSlash = i + enddo + do i = posLastCommonSlash+1,len_trim(a_cleaned) + if (a_cleaned(i:i) == '/') remainingSlashes = remainingSlashes + 1 + enddo - makeRelativePath = repeat('..'//'/',remainingSlashes)//b_cleaned(posLastCommonSlash+1:len_trim(b_cleaned)) + makeRelativePath = repeat('..'//'/',remainingSlashes)//b_cleaned(posLastCommonSlash+1:len_trim(b_cleaned)) end function makeRelativePath + !-------------------------------------------------------------------------------------------------- !> @brief sets global variable SIGUSR1 to .true. if program receives SIGUSR1 !-------------------------------------------------------------------------------------------------- subroutine setSIGUSR1(signal) bind(C) - use :: iso_c_binding - - implicit none - integer(C_INT), value :: signal - SIGUSR1 = .true. - - write(6,*) 'received signal ',signal, 'set SIGUSR1' - + use :: iso_c_binding + + implicit none + integer(C_INT), value :: signal + SIGUSR1 = .true. + + write(6,*) 'received signal ',signal, 'set SIGUSR1' + end subroutine setSIGUSR1 - - + + !-------------------------------------------------------------------------------------------------- !> @brief sets global variable SIGUSR2 to .true. if program receives SIGUSR2 !-------------------------------------------------------------------------------------------------- subroutine setSIGUSR2(signal) bind(C) - use :: iso_c_binding + use :: iso_c_binding - implicit none - integer(C_INT), value :: signal - SIGUSR2 = .true. + implicit none + integer(C_INT), value :: signal + SIGUSR2 = .true. + + write(6,*) 'received signal ',signal, 'set SIGUSR2' - write(6,*) 'received signal ',signal, 'set SIGUSR2' - end subroutine setSIGUSR2 !-------------------------------------------------------------------------------------------------- -!> @brief taken from IO, check IO_stringValue for documentation +!> @brief taken from IO, check IO_stringValue for documentation !-------------------------------------------------------------------------------------------------- pure function IIO_stringValue(string,chunkPos,myChunk) - - implicit none - integer(pInt), dimension(:), intent(in) :: chunkPos !< positions of start and end of each tag/chunk in given string - integer(pInt), intent(in) :: myChunk !< position number of desired chunk - character(len=chunkPos(myChunk*2+1)-chunkPos(myChunk*2)+1) :: IIO_stringValue - character(len=*), intent(in) :: string !< raw input with known start and end of each chunk - IIO_stringValue = string(chunkPos(myChunk*2):chunkPos(myChunk*2+1)) + implicit none + integer, dimension(:), intent(in) :: chunkPos !< positions of start and end of each tag/chunk in given string + integer, intent(in) :: myChunk !< position number of desired chunk + character(len=chunkPos(myChunk*2+1)-chunkPos(myChunk*2)+1) :: IIO_stringValue + character(len=*), intent(in) :: string !< raw input with known start and end of each chunk + + IIO_stringValue = string(chunkPos(myChunk*2):chunkPos(myChunk*2+1)) end function IIO_stringValue !-------------------------------------------------------------------------------------------------- -!> @brief taken from IO, check IO_intValue for documentation +!> @brief taken from IO, check IO_intValue for documentation !-------------------------------------------------------------------------------------------------- -integer(pInt) pure function IIO_intValue(string,chunkPos,myChunk) - - implicit none - character(len=*), intent(in) :: string !< raw input with known start and end of each chunk - integer(pInt), intent(in) :: myChunk !< position number of desired sub string - integer(pInt), dimension(:), intent(in) :: chunkPos !< positions of start and end of each tag/chunk in given string +integer pure function IIO_intValue(string,chunkPos,myChunk) + + implicit none + character(len=*), intent(in) :: string !< raw input with known start and end of each chunk + integer, intent(in) :: myChunk !< position number of desired sub string + integer, dimension(:), intent(in) :: chunkPos !< positions of start and end of each tag/chunk in given string - valuePresent: if (myChunk > chunkPos(1) .or. myChunk < 1_pInt) then - IIO_intValue = 0_pInt - else valuePresent - read(UNIT=string(chunkPos(myChunk*2):chunkPos(myChunk*2+1)),ERR=100,FMT=*) IIO_intValue - endif valuePresent - return -100 IIO_intValue = huge(1_pInt) + valuePresent: if (myChunk > chunkPos(1) .or. myChunk < 1) then + IIO_intValue = 0 + else valuePresent + read(UNIT=string(chunkPos(myChunk*2):chunkPos(myChunk*2+1)),ERR=100,FMT=*) IIO_intValue + endif valuePresent + return +100 IIO_intValue = huge(1) end function IIO_intValue !-------------------------------------------------------------------------------------------------- -!> @brief taken from IO, check IO_stringPos for documentation +!> @brief taken from IO, check IO_stringPos for documentation !-------------------------------------------------------------------------------------------------- pure function IIO_stringPos(string) - implicit none - integer(pInt), dimension(:), allocatable :: IIO_stringPos - character(len=*), intent(in) :: string !< string in which chunks are searched for - - character(len=*), parameter :: SEP=achar(44)//achar(32)//achar(9)//achar(10)//achar(13) ! comma and whitespaces - integer :: left, right ! no pInt (verify and scan return default integer) + implicit none + integer, dimension(:), allocatable :: IIO_stringPos + character(len=*), intent(in) :: string !< string in which chunks are searched for - allocate(IIO_stringPos(1), source=0_pInt) - right = 0 - - do while (verify(string(right+1:),SEP)>0) - left = right + verify(string(right+1:),SEP) - right = left + scan(string(left:),SEP) - 2 - IIO_stringPos = [IIO_stringPos,int(left, pInt), int(right, pInt)] - IIO_stringPos(1) = IIO_stringPos(1)+1_pInt - enddo + character(len=*), parameter :: SEP=achar(44)//achar(32)//achar(9)//achar(10)//achar(13) ! comma and whitespaces + integer :: left, right + + allocate(IIO_stringPos(1), source=0) + right = 0 + + do while (verify(string(right+1:),SEP)>0) + left = right + verify(string(right+1:),SEP) + right = left + scan(string(left:),SEP) - 2 + IIO_stringPos = [IIO_stringPos,left, right] + IIO_stringPos(1) = IIO_stringPos(1)+1 + enddo end function IIO_stringPos From 05eb80d38c9cdce8b895956c9bbb998c2a86d7fb Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 6 Mar 2019 15:49:31 +0100 Subject: [PATCH 3/7] pLongInt was not used --- src/CPFEM2.f90 | 2 +- src/constitutive.f90 | 6 ++---- src/crystallite.f90 | 3 +-- src/debug.f90 | 3 +-- 4 files changed, 5 insertions(+), 9 deletions(-) diff --git a/src/CPFEM2.f90 b/src/CPFEM2.f90 index b2aa2f598..087328bf6 100644 --- a/src/CPFEM2.f90 +++ b/src/CPFEM2.f90 @@ -93,7 +93,7 @@ subroutine CPFEM_init compiler_options #endif use prec, only: & - pInt, pReal, pLongInt + pInt, pReal use IO, only: & IO_timeStamp, & IO_error diff --git a/src/constitutive.f90 b/src/constitutive.f90 index b3b5ae875..36886da18 100644 --- a/src/constitutive.f90 +++ b/src/constitutive.f90 @@ -755,8 +755,7 @@ end subroutine constitutive_hooke_SandItsTangents !-------------------------------------------------------------------------------------------------- subroutine constitutive_collectDotState(S, FeArray, Fi, FpArray, subdt, ipc, ip, el) use prec, only: & - pReal, & - pLongInt + pReal use debug, only: & debug_level, & debug_constitutive, & @@ -896,8 +895,7 @@ end subroutine constitutive_collectDotState !-------------------------------------------------------------------------------------------------- subroutine constitutive_collectDeltaState(S, Fe, Fi, ipc, ip, el) use prec, only: & - pReal, & - pLongInt + pReal use debug, only: & debug_level, & debug_constitutive, & diff --git a/src/crystallite.f90 b/src/crystallite.f90 index 29724ced1..dbb3484d8 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -1088,8 +1088,7 @@ logical function integrateStress(& ) use, intrinsic :: & IEEE_arithmetic - use prec, only: pLongInt, & - tol_math_check, & + use prec, only: tol_math_check, & dEq0 use numerics, only: nStress, & aTol_crystalliteStress, & diff --git a/src/debug.f90 b/src/debug.f90 index 6debf84c2..7dcc018d3 100644 --- a/src/debug.f90 +++ b/src/debug.f90 @@ -8,8 +8,7 @@ module debug use prec, only: & pInt, & - pReal, & - pLongInt + pReal implicit none private From 51f8b1961f566475baec747aab5220c239d94caa Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 6 Mar 2019 15:52:52 +0100 Subject: [PATCH 4/7] simplify integer precision handling essentially, it should be ok to always use the default integer (which is 32 bit unless using MSC.Marc) and use 64 bit integer only for special cases where an overflow could happen --- src/prec.f90 | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/src/prec.f90 b/src/prec.f90 index 8a625b97e..ab47f8e34 100644 --- a/src/prec.f90 +++ b/src/prec.f90 @@ -15,12 +15,10 @@ module prec ! https://software.intel.com/en-us/blogs/2017/03/27/doctor-fortran-in-it-takes-all-kinds integer, parameter, public :: pReal = IEEE_selected_real_kind(15,307) !< number with 15 significant digits, up to 1e+-300 (typically 64 bit) -#if (INT==4) - integer, parameter, public :: pInt = selected_int_kind(9) !< number with at least up to +-1e9 (typically 32 bit) -#elif (INT==8) +#if(INT==8) integer, parameter, public :: pInt = selected_int_kind(18) !< number with at least up to +-1e18 (typically 64 bit) #else - NO SUITABLE PRECISION FOR INTEGER SELECTED, STOPPING COMPILATION + integer, parameter, public :: pInt = selected_int_kind(9) !< number with at least up to +-1e9 (typically 32 bit) #endif integer, parameter, public :: pLongInt = selected_int_kind(18) !< number with at least up to +-1e18 (typically 64 bit) integer, parameter, public :: pStringLen = 256 !< default string length From a965c46025718a4e0a51c6e1279a856d21226467 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 6 Mar 2019 15:57:42 +0100 Subject: [PATCH 5/7] improved functions for floating point comparison - less stric tolerance for comparison to zero - better readable - avoid "merge" on optional arguments (not safe) - "equal" and "notEqual" are now symmetric (assignment of <= and < is arbitrary) --- src/prec.f90 | 104 ++++++++++++++++++++++++++++++++++++--------------- 1 file changed, 74 insertions(+), 30 deletions(-) diff --git a/src/prec.f90 b/src/prec.f90 index ab47f8e34..36f92b7bd 100644 --- a/src/prec.f90 +++ b/src/prec.f90 @@ -79,6 +79,8 @@ module prec integer(pInt), pointer, dimension(:,:) :: p end type + real(pReal), private, parameter :: DBL_EPSILON = 2.220446049250313E-16_pReal !< minimum positive number such that 1.0 + DBL_EPSILON /= 1.0. + real(pReal), private, parameter :: DBL_MIN = 2.2250738585072014E-308_pReal !< smallest normalized floating point number public :: & prec_init, & @@ -126,12 +128,19 @@ end subroutine prec_init !-------------------------------------------------------------------------------------------------- logical elemental pure function dEq(a,b,tol) - implicit none - real(pReal), intent(in) :: a,b - real(pReal), intent(in), optional :: tol - real(pReal), parameter :: eps = 2.220446049250313E-16 ! DBL_EPSILON in C + implicit none + real(pReal), intent(in) :: a,b + real(pReal), intent(in), optional :: tol + real(pReal) :: eps + + if (present(tol)) then + eps = tol + else + eps = DBL_EPSILON * maxval(abs([a,b])) + endif + + dEq = merge(.True.,.False.,abs(a-b) < eps) - dEq = merge(.True.,.False.,abs(a-b) <= merge(tol,eps,present(tol))*maxval(abs([a,b]))) end function dEq @@ -143,12 +152,19 @@ end function dEq !-------------------------------------------------------------------------------------------------- logical elemental pure function dNeq(a,b,tol) - implicit none - real(pReal), intent(in) :: a,b - real(pReal), intent(in), optional :: tol - real(pReal), parameter :: eps = 2.220446049250313E-16 ! DBL_EPSILON in C + implicit none + real(pReal), intent(in) :: a,b + real(pReal), intent(in), optional :: tol + real(pReal) :: eps + + if (present(tol)) then + eps = tol + else + eps = DBL_EPSILON * maxval(abs([a,b])) + endif + + dNeq = merge(.False.,.True.,abs(a-b) <= eps) - dNeq = merge(.False.,.True.,abs(a-b) <= merge(tol,eps,present(tol))*maxval(abs([a,b]))) end function dNeq @@ -160,12 +176,19 @@ end function dNeq !-------------------------------------------------------------------------------------------------- logical elemental pure function dEq0(a,tol) - implicit none - real(pReal), intent(in) :: a - real(pReal), intent(in), optional :: tol - real(pReal), parameter :: eps = 2.2250738585072014E-308 ! smallest non-denormalized number + implicit none + real(pReal), intent(in) :: a + real(pReal), intent(in), optional :: tol + real(pReal) :: eps + + if (present(tol)) then + eps = tol + else + eps = DBL_MIN * 10.0_pReal + endif + + dEq0 = merge(.True.,.False.,abs(a) < eps) - dEq0 = merge(.True.,.False.,abs(a) <= merge(tol,eps,present(tol))) end function dEq0 @@ -177,12 +200,19 @@ end function dEq0 !-------------------------------------------------------------------------------------------------- logical elemental pure function dNeq0(a,tol) - implicit none - real(pReal), intent(in) :: a - real(pReal), intent(in), optional :: tol - real(pReal), parameter :: eps = 2.2250738585072014E-308 ! smallest non-denormalized number + implicit none + real(pReal), intent(in) :: a + real(pReal), intent(in), optional :: tol + real(pReal) :: eps + + if (present(tol)) then + eps = tol + else + eps = DBL_MIN * 10.0_pReal + endif + + dNeq0 = merge(.False.,.True.,abs(a) <= eps) - dNeq0 = merge(.False.,.True.,abs(a) <= merge(tol,eps,present(tol))) end function dNeq0 @@ -195,12 +225,19 @@ end function dNeq0 !-------------------------------------------------------------------------------------------------- logical elemental pure function cEq(a,b,tol) - implicit none - complex(pReal), intent(in) :: a,b - real(pReal), intent(in), optional :: tol - real(pReal), parameter :: eps = 2.220446049250313E-16 ! DBL_EPSILON in C + implicit none + complex(pReal), intent(in) :: a,b + real(pReal), intent(in), optional :: tol + real(pReal) :: eps + + if (present(tol)) then + eps = tol + else + eps = DBL_EPSILON * maxval(abs([a,b])) + endif + + cEq = merge(.True.,.False.,abs(a-b) < eps) - cEq = merge(.True.,.False.,abs(a-b) <= merge(tol,eps,present(tol))*maxval(abs([a,b]))) end function cEq @@ -213,12 +250,19 @@ end function cEq !-------------------------------------------------------------------------------------------------- logical elemental pure function cNeq(a,b,tol) - implicit none - complex(pReal), intent(in) :: a,b - real(pReal), intent(in), optional :: tol - real(pReal), parameter :: eps = 2.220446049250313E-16 ! DBL_EPSILON in C + implicit none + complex(pReal), intent(in) :: a,b + real(pReal), intent(in), optional :: tol + real(pReal) :: eps + + if (present(tol)) then + eps = tol + else + eps = DBL_EPSILON * maxval(abs([a,b])) + endif + + cNeq = merge(.False.,.True.,abs(a-b) <= eps) - cNeq = merge(.False.,.True.,abs(a-b) <= merge(tol,eps,present(tol))*maxval(abs([a,b]))) end function cNeq end module prec From ccb62da24ae4a00a861e3abbd80e0f94a1f2752c Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 6 Mar 2019 16:33:39 +0100 Subject: [PATCH 6/7] kind-ID does not need to coincide with the number of bytes --- src/HDF5_utilities.f90 | 31 ++++++++++++++++--------------- 1 file changed, 16 insertions(+), 15 deletions(-) diff --git a/src/HDF5_utilities.f90 b/src/HDF5_utilities.f90 index 0582318ce..9878f1c76 100644 --- a/src/HDF5_utilities.f90 +++ b/src/HDF5_utilities.f90 @@ -89,26 +89,27 @@ module HDF5_utilities contains subroutine HDF5_utilities_init - use, intrinsic :: & - iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment) - implicit none - integer(HDF5_ERR_TYPE) :: hdferr - integer(SIZE_T) :: typeSize + implicit none + integer(HDF5_ERR_TYPE) :: hdferr + integer(SIZE_T) :: typeSize - write(6,'(/,a)') ' <<<+- HDF5_Utilities init -+>>>' -#include "compilation_info.f90" + write(6,'(/,a)') ' <<<+- HDF5_Utilities init -+>>>' !-------------------------------------------------------------------------------------------------- !initialize HDF5 library and check if integer and float type size match - call h5open_f(hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_Utilities_init: h5open_f') - call h5tget_size_f(H5T_NATIVE_INTEGER,typeSize, hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_Utilities_init: h5tget_size_f (int)') - if (int(pInt,SIZE_T)/=typeSize) call IO_error(0_pInt,ext_msg='pInt does not match H5T_NATIVE_INTEGER') - call h5tget_size_f(H5T_NATIVE_DOUBLE,typeSize, hdferr) - if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_Utilities_init: h5tget_size_f (double)') - if (int(pReal,SIZE_T)/=typeSize) call IO_error(0_pInt,ext_msg='pReal does not match H5T_NATIVE_DOUBLE') + call h5open_f(hdferr) + if (hdferr < 0) call IO_error(1,ext_msg='HDF5_Utilities_init: h5open_f') + + call h5tget_size_f(H5T_NATIVE_INTEGER,typeSize, hdferr) + if (hdferr < 0) call IO_error(1,ext_msg='HDF5_Utilities_init: h5tget_size_f (int)') + if (int(bit_size(0),SIZE_T)/=typeSize*8) & + call IO_error(0_pInt,ext_msg='Default integer size does not match H5T_NATIVE_INTEGER') + + call h5tget_size_f(H5T_NATIVE_DOUBLE,typeSize, hdferr) + if (hdferr < 0) call IO_error(1,ext_msg='HDF5_Utilities_init: h5tget_size_f (double)') + if (int(storage_size(0.0_pReal),SIZE_T)/=typeSize*8) & + call IO_error(0,ext_msg='pReal does not match H5T_NATIVE_DOUBLE') end subroutine HDF5_utilities_init From 7a083c4098ddaa43352e89c4c16aed1a46f16359 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 7 Mar 2019 11:02:27 +0100 Subject: [PATCH 7/7] [skip ci] more general names and procedure --- src/prec.f90 | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/src/prec.f90 b/src/prec.f90 index 36f92b7bd..b134e5b46 100644 --- a/src/prec.f90 +++ b/src/prec.f90 @@ -79,8 +79,8 @@ module prec integer(pInt), pointer, dimension(:,:) :: p end type - real(pReal), private, parameter :: DBL_EPSILON = 2.220446049250313E-16_pReal !< minimum positive number such that 1.0 + DBL_EPSILON /= 1.0. - real(pReal), private, parameter :: DBL_MIN = 2.2250738585072014E-308_pReal !< smallest normalized floating point number + real(pReal), private, parameter :: PREAL_EPSILON = epsilon(0.0_pReal) !< minimum positive number such that 1.0 + EPSILON /= 1.0. + real(pReal), private, parameter :: PREAL_MIN = tiny(0.0_pReal) !< smallest normalized floating point number public :: & prec_init, & @@ -136,7 +136,7 @@ logical elemental pure function dEq(a,b,tol) if (present(tol)) then eps = tol else - eps = DBL_EPSILON * maxval(abs([a,b])) + eps = PREAL_EPSILON * maxval(abs([a,b])) endif dEq = merge(.True.,.False.,abs(a-b) < eps) @@ -160,7 +160,7 @@ logical elemental pure function dNeq(a,b,tol) if (present(tol)) then eps = tol else - eps = DBL_EPSILON * maxval(abs([a,b])) + eps = PREAL_EPSILON * maxval(abs([a,b])) endif dNeq = merge(.False.,.True.,abs(a-b) <= eps) @@ -184,7 +184,7 @@ logical elemental pure function dEq0(a,tol) if (present(tol)) then eps = tol else - eps = DBL_MIN * 10.0_pReal + eps = PREAL_MIN * 10.0_pReal endif dEq0 = merge(.True.,.False.,abs(a) < eps) @@ -208,7 +208,7 @@ logical elemental pure function dNeq0(a,tol) if (present(tol)) then eps = tol else - eps = DBL_MIN * 10.0_pReal + eps = PREAL_MIN * 10.0_pReal endif dNeq0 = merge(.False.,.True.,abs(a) <= eps) @@ -233,7 +233,7 @@ logical elemental pure function cEq(a,b,tol) if (present(tol)) then eps = tol else - eps = DBL_EPSILON * maxval(abs([a,b])) + eps = PREAL_EPSILON * maxval(abs([a,b])) endif cEq = merge(.True.,.False.,abs(a-b) < eps) @@ -258,7 +258,7 @@ logical elemental pure function cNeq(a,b,tol) if (present(tol)) then eps = tol else - eps = DBL_EPSILON * maxval(abs([a,b])) + eps = PREAL_EPSILON * maxval(abs([a,b])) endif cNeq = merge(.False.,.True.,abs(a-b) <= eps)