DAMASK_EICMD/src/prec.f90

264 lines
11 KiB
Fortran
Raw Normal View History

!--------------------------------------------------------------------------------------------------
!> @author Franz Roters, Max-Planck-Institut für Eisenforschung GmbH
!> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH
!> @author Christoph Kords, Max-Planck-Institut für Eisenforschung GmbH
!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
!> @author Luv Sharma, Max-Planck-Institut für Eisenforschung GmbH
!> @brief setting precision for real and int type
!--------------------------------------------------------------------------------------------------
module prec
2019-05-11 01:18:43 +05:30
use, intrinsic :: IEEE_arithmetic
2019-03-07 23:59:48 +05:30
implicit none
2019-09-23 12:23:56 +05:30
public
! https://software.intel.com/en-us/blogs/2017/03/27/doctor-fortran-in-it-takes-all-kinds
2019-03-07 23:07:22 +05:30
integer, parameter, public :: pReal = IEEE_selected_real_kind(15,307) !< number with 15 significant digits, up to 1e+-307 (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
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 :: pPathLen = 4096 !< maximum length of a path name on linux
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
2019-03-07 23:59:48 +05:30
integer, 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
2019-03-07 23:59:48 +05:30
integer :: &
sizeState = 0, & !< size of state
sizeDotState = 0, & !< size of dot state, i.e. state(1:sizeDot) follows time evolution by dotState rates
offsetDeltaState = 0, & !< index offset of delta state
2019-12-11 04:40:02 +05:30
sizeDeltaState = 0 !< size of delta state, i.e. state(offset+1:offset+sizeDelta) follows time evolution by deltaState increments
real(pReal), pointer, dimension(:), contiguous :: &
atolState
2019-03-07 23:59:48 +05:30
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
2019-03-07 23:59:48 +05:30
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
2019-03-07 23:59:48 +05:30
type, public :: tHomogMapping
2019-03-07 23:59:48 +05:30
integer, pointer, dimension(:,:) :: p
end type
2019-03-07 23:59:48 +05:30
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
2020-01-02 23:14:51 +05:30
integer, dimension(0), parameter, public :: &
emptyIntArray = [integer::]
real(pReal), dimension(0), parameter, public :: &
emptyRealArray = [real(pReal)::]
character(len=pStringLen), dimension(0), parameter, public :: &
emptyStringArray = [character(len=pStringLen)::]
2019-12-21 15:49:33 +05:30
2019-09-23 12:23:56 +05:30
private :: &
unitTest
2019-03-07 23:59:48 +05:30
contains
2015-04-13 15:32:52 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief reporting precision
!--------------------------------------------------------------------------------------------------
subroutine prec_init
2015-04-13 15:32:52 +05:30
write(6,'(/,a)') ' <<<+- prec init -+>>>'
2019-03-07 23:59:48 +05:30
write(6,'(a,i3)') ' Size of integer in bit: ',bit_size(0)
write(6,'(a,i19)') ' Maximum value: ',huge(0)
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)
2019-09-23 12:23:56 +05:30
call unitTest
2019-03-07 23:59:48 +05:30
end subroutine prec_init
2015-04-13 15:32:52 +05:30
!--------------------------------------------------------------------------------------------------
2016-05-29 14:15:03 +05:30
!> @brief equality comparison for float with double precision
! replaces "==" but for certain (relative) tolerance. Counterpart to dNeq
! https://randomascii.wordpress.com/2012/02/25/comparing-floating-point-numbers-2012-edition/
! AlmostEqualRelative
!--------------------------------------------------------------------------------------------------
logical elemental pure function dEq(a,b,tol)
real(pReal), intent(in) :: a,b
real(pReal), intent(in), optional :: tol
real(pReal) :: eps
2019-03-07 23:59:48 +05:30
if (present(tol)) then
eps = tol
else
eps = PREAL_EPSILON * maxval(abs([a,b]))
endif
dEq = merge(.True.,.False.,abs(a-b) <= eps)
end function dEq
!--------------------------------------------------------------------------------------------------
2016-05-29 14:15:03 +05:30
!> @brief inequality comparison for float with double precision
! replaces "!=" but for certain (relative) tolerance. Counterpart to dEq
! https://randomascii.wordpress.com/2012/02/25/comparing-floating-point-numbers-2012-edition/
! AlmostEqualRelative NOT
!--------------------------------------------------------------------------------------------------
logical elemental pure function dNeq(a,b,tol)
real(pReal), intent(in) :: a,b
real(pReal), intent(in), optional :: tol
2019-03-07 23:59:48 +05:30
if (present(tol)) then
2019-09-23 12:23:56 +05:30
dNeq = .not. dEq(a,b,tol)
else
2019-09-23 12:23:56 +05:30
dNeq = .not. dEq(a,b)
endif
end function dNeq
!--------------------------------------------------------------------------------------------------
!> @brief equality to 0 comparison for float with double precision
! replaces "==0" but everything not representable as a normal number is treated as 0. Counterpart to dNeq0
! https://de.mathworks.com/help/matlab/ref/realmin.html
! https://docs.oracle.com/cd/E19957-01/806-3568/ncg_math.html
!--------------------------------------------------------------------------------------------------
logical elemental pure function dEq0(a,tol)
real(pReal), intent(in) :: a
real(pReal), intent(in), optional :: tol
real(pReal) :: eps
2019-03-07 23:59:48 +05:30
if (present(tol)) then
eps = tol
else
eps = PREAL_MIN * 10.0_pReal
endif
dEq0 = merge(.True.,.False.,abs(a) <= eps)
end function dEq0
!--------------------------------------------------------------------------------------------------
!> @brief inequality to 0 comparison for float with double precision
! replaces "!=0" but everything not representable as a normal number is treated as 0. Counterpart to dEq0
! https://de.mathworks.com/help/matlab/ref/realmin.html
! https://docs.oracle.com/cd/E19957-01/806-3568/ncg_math.html
!--------------------------------------------------------------------------------------------------
logical elemental pure function dNeq0(a,tol)
real(pReal), intent(in) :: a
real(pReal), intent(in), optional :: tol
2019-03-07 23:59:48 +05:30
if (present(tol)) then
2019-09-23 12:23:56 +05:30
dNeq0 = .not. dEq0(a,tol)
else
2019-09-23 12:23:56 +05:30
dNeq0 = .not. dEq0(a)
endif
end function dNeq0
2016-05-29 14:15:03 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief equality comparison for complex with double precision
! replaces "==" but for certain (relative) tolerance. Counterpart to cNeq
! https://randomascii.wordpress.com/2012/02/25/comparing-floating-point-numbers-2012-edition/
2016-05-29 14:15:03 +05:30
! probably a component wise comparison would be more accurate than the comparsion of the absolute
! value
!--------------------------------------------------------------------------------------------------
logical elemental pure function cEq(a,b,tol)
complex(pReal), intent(in) :: a,b
real(pReal), intent(in), optional :: tol
real(pReal) :: eps
2019-03-07 23:59:48 +05:30
if (present(tol)) then
eps = tol
else
eps = PREAL_EPSILON * maxval(abs([a,b]))
endif
cEq = merge(.True.,.False.,abs(a-b) <= eps)
2016-05-29 14:15:03 +05:30
end function cEq
!--------------------------------------------------------------------------------------------------
!> @brief inequality comparison for complex with double precision
! replaces "!=" but for certain (relative) tolerance. Counterpart to cEq
! https://randomascii.wordpress.com/2012/02/25/comparing-floating-point-numbers-2012-edition/
2016-05-29 14:15:03 +05:30
! probably a component wise comparison would be more accurate than the comparsion of the absolute
! value
!--------------------------------------------------------------------------------------------------
logical elemental pure function cNeq(a,b,tol)
complex(pReal), intent(in) :: a,b
real(pReal), intent(in), optional :: tol
2019-03-07 23:59:48 +05:30
if (present(tol)) then
2019-09-23 12:23:56 +05:30
cNeq = .not. cEq(a,b,tol)
else
2019-09-23 12:23:56 +05:30
cNeq = .not. cEq(a,b)
endif
2016-05-29 14:15:03 +05:30
end function cNeq
2019-09-23 12:23:56 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief check correctness of (some) prec functions
!--------------------------------------------------------------------------------------------------
subroutine unitTest
integer, allocatable, dimension(:) :: realloc_lhs_test
real(pReal), dimension(2) :: r
external :: &
quit
call random_number(r)
r = r/minval(r)
if(.not. all(dEq(r,r+PREAL_EPSILON))) call quit(9000)
if(dEq(r(1),r(2)) .and. dNeq(r(1),r(2))) call quit(9000)
if(.not. all(dEq0(r-(r+PREAL_MIN)))) call quit(9000)
2019-09-23 12:23:56 +05:30
realloc_lhs_test = [1,2]
if (any(realloc_lhs_test/=[1,2])) call quit(9000)
end subroutine unitTest
end module prec