223031012
/
DAMASK_EICMD
mirror of https://github.com/achalhp/DAMASK_EICMD.git
1
0
Fork 0
Code Issues Projects Releases Wiki Activity

Merge branch 'improved-presicion-handling-2' into 'development' 

Improved presicion handling 2

See merge request damask/DAMASK!62
This commit is contained in:
Franz Roters 2019-03-07 16:58:15 +01:00
parent 507963d558 7a083c4098
commit 29812e320d
8 changed files with 522 additions and 488 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
src/CPFEM2.f90
View File

@ -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

91
src/DAMASK_interface.f90
View File

@ -10,13 +10,13 @@
!> 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
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
@ -42,8 +42,15 @@ contains
subroutine DAMASK_interface_init()
use, intrinsic :: &
iso_fortran_env
use :: &
use, intrinsic :: &
iso_c_binding
use PETScSys
use system_routines, only: &
signalusr1_C, &
signalusr2_C, &
getHostName, &
getCWD
#include <petsc/finclude/petscsys.h>
#if defined(__GFORTRAN__) && __GNUC__ < 5
===================================================================================================
@ -81,13 +88,6 @@ 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
@ -105,7 +105,7 @@ subroutine DAMASK_interface_init()
integer, allocatable, dimension(:) :: &
chunkPos
integer, dimension(8) :: &
dateAndTime ! type default integer
dateAndTime
PetscErrorCode :: ierr
external :: &
quit
@ -120,7 +120,7 @@ subroutine DAMASK_interface_init()
call MPI_Init_Thread(MPI_THREAD_FUNNELED,threadLevel,ierr);CHKERRQ(ierr)
if (threadLevel<MPI_THREAD_FUNNELED) then
write(6,'(a)') ' MPI library does not support OpenMP'
call quit(1_pInt)
call quit(1)
endif
#endif
call PETScInitialize(PETSC_NULL_CHARACTER,ierr) ! according to PETSc manual, that should be the first line in the code
@ -130,11 +130,11 @@ subroutine DAMASK_interface_init()
mainProcess: if (worldrank == 0) then
if (output_unit /= 6) then
write(output_unit,'(a)') ' STDOUT != 6'
call quit(1_pInt)
call quit(1)
endif
if (error_unit /= 0) then
write(output_unit,'(a)') ' STDERR != 0'
call quit(1_pInt)
call quit(1)
endif
else mainProcess
close(6) ! disable output for non-master processes (open 6 to rank specific file for debug)
@ -167,7 +167,7 @@ subroutine DAMASK_interface_init()
call get_command(commandLine)
chunkPos = IIO_stringPos(commandLine)
do i = 2_pInt, chunkPos(1)
do i = 2, chunkPos(1)
select case(IIO_stringValue(commandLine,chunkPos,i)) ! extract key
case ('-h','--help')
write(6,'(a)') ' #######################################################################'
@ -205,23 +205,23 @@ subroutine DAMASK_interface_init()
write(6,'(a)') ' Help:'
write(6,'(/,a)')' --help'
write(6,'(a,/)')' Prints this message and exits'
call quit(0_pInt) ! normal Termination
call quit(0) ! normal Termination
case ('-l', '--load', '--loadcase')
if ( i < chunkPos(1)) loadcaseArg = trim(IIO_stringValue(commandLine,chunkPos,i+1_pInt))
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_pInt))
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_pInt))
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_pInt)
interface_restartInc = IIO_IntValue(commandLine,chunkPos,i+1)
endif
end select
enddo
if (len_trim(loadcaseArg) == 0 .or. len_trim(geometryArg) == 0) then
write(6,'(a)') ' Please specify geometry AND load case (-h for help)'
call quit(1_pInt)
call quit(1)
endif
if (len_trim(workingDirArg) > 0) call setWorkingDirectory(trim(workingDirArg))
@ -243,7 +243,7 @@ subroutine DAMASK_interface_init()
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) &
if (interface_restartInc > 0) &
write(6,'(a,i6.6)') ' Restart from increment: ', interface_restartInc
call signalusr1_c(c_funloc(setSIGUSR1))
@ -280,8 +280,8 @@ subroutine setWorkingDirectory(workingDirectoryArg)
workingDirectory = trim(rectifyPath(workingDirectory))
error = setCWD(trim(workingDirectory))
if(error) then
write(6,'(a20,a,a16)') ' working directory "',trim(workingDirectory),'" does not exist'
call quit(1_pInt)
write(6,'(a20,a,a16)') ' Working directory "',trim(workingDirectory),'" does not exist'
call quit(1)
endif
end subroutine setWorkingDirectory
@ -331,7 +331,7 @@ character(len=1024) function getGeometryFile(geometryParameter)
inquire(file=trim(getGeometryFile), exist=file_exists)
if (.not. file_exists) then
write(6,'(a)') ' Geometry file does not exists ('//trim(getGeometryFile)//')'
call quit(1_pInt)
call quit(1)
endif
end function getGeometryFile
@ -355,8 +355,8 @@ character(len=1024) function getLoadCaseFile(loadCaseParameter)
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)
write(6,'(a)') ' Load case file does not exists ('//trim(getLoadCaseFile)//')'
call quit(1)
endif
end function getLoadCaseFile
@ -371,7 +371,7 @@ function rectifyPath(path)
implicit none
character(len=*) :: path
character(len=1024) :: rectifyPath
integer :: i,j,k,l ! no pInt
integer :: i,j,k,l
!--------------------------------------------------------------------------------------------------
! remove /./ from path
@ -416,7 +416,7 @@ 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
integer :: i,posLastCommonSlash,remainingSlashes
posLastCommonSlash = 0
remainingSlashes = 0
@ -435,6 +435,7 @@ character(len=1024) function makeRelativePath(a,b)
end function makeRelativePath
!--------------------------------------------------------------------------------------------------
!> @brief sets global variable SIGUSR1 to .true. if program receives SIGUSR1
!--------------------------------------------------------------------------------------------------
@ -471,8 +472,8 @@ end subroutine setSIGUSR2
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
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
@ -484,21 +485,21 @@ end function IIO_stringValue
!--------------------------------------------------------------------------------------------------
!> @brief taken from IO, check IO_intValue for documentation
!--------------------------------------------------------------------------------------------------
integer(pInt) pure function IIO_intValue(string,chunkPos,myChunk)
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(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, 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
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_pInt)
100 IIO_intValue = huge(1)
end function IIO_intValue
@ -509,20 +510,20 @@ end function IIO_intValue
pure function IIO_stringPos(string)
implicit none
integer(pInt), dimension(:), allocatable :: IIO_stringPos
integer, 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)
integer :: left, right
allocate(IIO_stringPos(1), source=0_pInt)
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,int(left, pInt), int(right, pInt)]
IIO_stringPos(1) = IIO_stringPos(1)+1_pInt
IIO_stringPos = [IIO_stringPos,left, right]
IIO_stringPos(1) = IIO_stringPos(1)+1
enddo
end function IIO_stringPos

17
src/HDF5_utilities.f90
View File

@ -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
write(6,'(/,a)') ' <<<+- HDF5_Utilities init -+>>>'
#include "compilation_info.f90"
!--------------------------------------------------------------------------------------------------
!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')
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_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')
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_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')
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

6
src/constitutive.f90
View File

@ -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, &

3
src/crystallite.f90
View File

@ -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, &

3
src/debug.f90
View File

@ -8,8 +8,7 @@
module debug
use prec, only: &
pInt, &
pReal, &
pLongInt
pReal
implicit none
private

1
src/material.f90
View File

@ -15,7 +15,6 @@ module material
tPlasticState, &
tSourceState, &
tHomogMapping, &
tPhaseMapping, &
group_float, &
group_int

113
src/prec.f90
View File

@ -7,28 +7,24 @@
!> @brief setting precision for real and int type
!--------------------------------------------------------------------------------------------------
module prec
! ToDo: use, intrinsic :: iso_fortran_env, only : I8 => int64, WP => real64
use, intrinsic :: IEEE_arithmetic, only:&
IEEE_selected_real_kind
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
! https://software.intel.com/en-us/blogs/2017/03/27/doctor-fortran-in-it-takes-all-kinds
#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)
integer, parameter, public :: pReal = IEEE_selected_real_kind(15,307) !< number with 15 significant digits, up to 1e+-300 (typically 64 bit)
#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
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
@ -83,9 +79,8 @@ module prec
integer(pInt), pointer, dimension(:,:) :: p
end type
type, public :: tPhaseMapping
integer(pInt), pointer, dimension(:,:,:) :: p
end type
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, &
@ -103,21 +98,21 @@ 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
integer(pInt), allocatable, dimension(:) :: realloc_lhs_test
external :: &
quit
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
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)
@ -136,9 +131,16 @@ 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
real(pReal) :: eps
if (present(tol)) then
eps = tol
else
eps = PREAL_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
@ -153,9 +155,16 @@ 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
real(pReal) :: eps
if (present(tol)) then
eps = tol
else
eps = PREAL_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
@ -170,9 +179,16 @@ 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
real(pReal) :: eps
if (present(tol)) then
eps = tol
else
eps = PREAL_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
@ -187,9 +203,16 @@ 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
real(pReal) :: eps
if (present(tol)) then
eps = tol
else
eps = PREAL_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
@ -205,9 +228,16 @@ 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
real(pReal) :: eps
if (present(tol)) then
eps = tol
else
eps = PREAL_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
@ -223,9 +253,16 @@ 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
real(pReal) :: eps
if (present(tol)) then
eps = tol
else
eps = PREAL_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