improved NaN checks

This commit is contained in:
Martin Diehl 2015-04-13 10:02:52 +00:00
parent 50998bd6a4
commit 5c1804e77d
4 changed files with 82 additions and 63 deletions

View File

@ -530,9 +530,14 @@ DAMASK_interface.o: DAMASK_spectral_interface.f90 \
# and no user-defined procedure with the same name as any intrinsic will be called except when it is explicitly declared external
# --> allows the use of 'getcwd'
prec.o: prec.f90
$(PREFIX) $(COMPILERNAME) $(COMPILE) -c prec.f90 -fno-range-check $(SUFFIX)
$(PREFIX) $(COMPILERNAME) $(COMPILE) -c prec.f90 -fno-range-check -fall-intrinsics -fno-fast-math $(SUFFIX)
# fno-range-check: Disable range checking on results of simplification of constant expressions during compilation
# --> allows the definition of DAMASK_NaN
#-fall-intrinsics: all intrinsic procedures (including the GNU-specific extensions) are accepted. -Wintrinsics-std will be ignored
# and no user-defined procedure with the same name as any intrinsic will be called except when it is explicitly declared external
# --> allows the use of 'isnan'
#-fno-fast-math:
# --> otherwise, when setting -ffast-math, isnan always evaluates to false (I would call it a bug)
else
DAMASK_interface.o: DAMASK_spectral_interface.f90 \
$(wildcard DAMASK_FEM_interface.f90) \

View File

@ -1,12 +1,12 @@
!##############################################################
!$Id$
#ifdef __GFORTRAN__
write(6,*) 'Compiled with ', compiler_version() !not supported by GFORTRAN 4.5 and ifort 12
write(6,*) 'With options ', compiler_options()
write(6,*) 'Compiled with ', compiler_version() !not supported by and ifort <= 15 (and old gfortran)
write(6,*) 'With options ', compiler_options()
#endif
#ifdef __INTEL_COMPILER
write(6,'(a,i4.4,a,i8.8)') ' Compiled with Intel fortran version ', __INTEL_COMPILER,&
', build date ', __INTEL_COMPILER_BUILD_DATE
', build date ', __INTEL_COMPILER_BUILD_DATE
#endif
write(6,*) 'Compiled on ', __DATE__,' at ',__TIME__
write(6,*)

View File

@ -1524,6 +1524,8 @@ end subroutine crystallite_stressAndItsTangent
!> @brief integrate stress, state with 4th order explicit Runge Kutta method
!--------------------------------------------------------------------------------------------------
subroutine crystallite_integrateStateRK4()
use prec, only: &
prec_isNaN
use numerics, only: &
numerics_integrationMode
use debug, only: &
@ -1621,16 +1623,14 @@ subroutine crystallite_integrateStateRK4()
if (crystallite_todo(g,i,e)) then
c = mappingConstitutive(1,g,i,e)
p = mappingConstitutive(2,g,i,e)
if ( any(plasticState(p)%dotState(:,c) /= plasticState(p)%dotState(:,c)) .or.&
any(damageState(p)%dotState(:,c) /= damageState(p)%dotState(:,c)) .or.&
any(thermalState(p)%dotState(:,c) /= thermalState(p)%dotState(:,c)) .or.&
any(vacancyState(p)%dotState(:,c) /= vacancyState(p)%dotState(:,c))) then ! NaN occured in dotState
if (.not. crystallite_localPlasticity(g,i,e)) then ! if broken non-local...
if ( any(prec_isNaN([plasticState(p)%dotState(:,c), damageState(p)%dotState(:,c), &
thermalState(p)%dotState(:,c), vacancyState(p)%dotState(:,c)]))) then ! NaN occured in any dotState
if (.not. crystallite_localPlasticity(g,i,e)) then ! if broken non-local...
!$OMP CRITICAL (checkTodo)
crystallite_todo = crystallite_todo .and. crystallite_localPlasticity ! ...all non-locals skipped
crystallite_todo = crystallite_todo .and. crystallite_localPlasticity ! ...all non-locals skipped
!$OMP END CRITICAL (checkTodo)
else ! if broken local...
crystallite_todo(g,i,e) = .false. ! ... skip this one next time
else ! if broken local...
crystallite_todo(g,i,e) = .false. ! ... skip this one next time
endif
endif
endif
@ -1772,10 +1772,8 @@ subroutine crystallite_integrateStateRK4()
p = mappingConstitutive(2,g,i,e)
c = mappingConstitutive(1,g,i,e)
if ( any(plasticState(p)%dotState(:,c) /= plasticState(p)%dotState(:,c)) .or.&
any(damageState(p)%dotState(:,c) /= damageState(p)%dotState(:,c)) .or.&
any(thermalState(p)%dotState(:,c) /= thermalState(p)%dotState(:,c)) .or.&
any(vacancyState(p)%dotState(:,c) /= vacancyState(p)%dotState(:,c))) then ! NaN occured in dotState
if (any(prec_isNaN([plasticState(p)%dotState(:,c), damageState(p)%dotState(:,c), &
thermalState(p)%dotState(:,c), vacancyState(p)%dotState(:,c)]))) then ! NaN occured in any dotState
if (.not. crystallite_localPlasticity(g,i,e)) then ! if broken non-local...
!$OMP CRITICAL (checkTodo)
crystallite_todo = crystallite_todo .and. crystallite_localPlasticity ! ...all non-locals skipped
@ -1824,6 +1822,8 @@ end subroutine crystallite_integrateStateRK4
!> adaptive step size (use 5th order solution to advance = "local extrapolation")
!--------------------------------------------------------------------------------------------------
subroutine crystallite_integrateStateRKCK45()
use prec, only: &
prec_isNaN
use debug, only: &
debug_level, &
debug_crystallite, &
@ -1948,10 +1948,8 @@ subroutine crystallite_integrateStateRKCK45()
if (crystallite_todo(g,i,e)) then
cc = mappingConstitutive(1,g,i,e)
p = mappingConstitutive(2,g,i,e)
if ( any(plasticState(p)%dotState(:,cc) /= plasticState(p)%dotState(:,cc)) .or.&
any(damageState(p)%dotState(:,cc) /= damageState(p)%dotState(:,cc)) .or.&
any(thermalState(p)%dotState(:,cc) /= thermalState(p)%dotState(:,cc)) .or.&
any(vacancyState(p)%dotState(:,cc) /= vacancyState(p)%dotState(:,cc))) then ! NaN occured in dotState
if (any(prec_isNaN([plasticState(p)%dotState(:,cc), damageState(p)%dotState(:,cc), &
thermalState(p)%dotState(:,cc), vacancyState(p)%dotState(:,cc)]))) then ! NaN occured in any dotState
if (.not. crystallite_localPlasticity(g,i,e)) then ! if broken non-local...
!$OMP CRITICAL (checkTodo)
crystallite_todo = crystallite_todo .and. crystallite_localPlasticity ! ...all non-locals skipped
@ -2108,10 +2106,8 @@ subroutine crystallite_integrateStateRKCK45()
p = mappingConstitutive(2,g,i,e)
cc = mappingConstitutive(1,g,i,e)
if ( any(plasticState(p)%dotState(:,cc) /= plasticState(p)%dotState(:,cc)) .or.&
any(damageState(p)%dotState(:,cc) /= damageState(p)%dotState(:,cc)) .or.&
any(thermalState(p)%dotState(:,cc) /= thermalState(p)%dotState(:,cc)) .or.&
any(vacancyState(p)%dotState(:,cc) /= vacancyState(p)%dotState(:,cc))) then ! NaN occured in dotState
if (any(prec_isNaN([plasticState(p)%dotState(:,cc), damageState(p)%dotState(:,cc), &
thermalState(p)%dotState(:,cc), vacancyState(p)%dotState(:,cc)]))) then ! NaN occured in any dotState
if (.not. crystallite_localPlasticity(g,i,e)) then ! if broken non-local...
!$OMP CRITICAL (checkTodo)
crystallite_todo = crystallite_todo .and. crystallite_localPlasticity ! ...all non-locals skipped
@ -2359,7 +2355,8 @@ end subroutine crystallite_integrateStateRKCK45
!> @brief integrate stress, state with 1st order Euler method with adaptive step size
!--------------------------------------------------------------------------------------------------
subroutine crystallite_integrateStateAdaptiveEuler()
use prec, only: &
prec_isNaN
use debug, only: &
debug_level, &
debug_crystallite, &
@ -2470,10 +2467,8 @@ subroutine crystallite_integrateStateAdaptiveEuler()
if (crystallite_todo(g,i,e)) then
p = mappingConstitutive(2,g,i,e)
c = mappingConstitutive(1,g,i,e)
if ( any(plasticState(p)%dotState(:,c) /= plasticState(p)%dotState(:,c)) .or. &
any(damageState( p)%dotState(:,c) /= damageState( p)%dotState(:,c)) .or. &
any(thermalState(p)%dotState(:,c) /= thermalState(p)%dotState(:,c)) .or.&
any(vacancyState(p)%dotState(:,c) /= vacancyState(p)%dotState(:,c))) then ! NaN occured in dotState
if (any(prec_isNaN([plasticState(p)%dotState(:,c), damageState(p)%dotState(:,c), &
thermalState(p)%dotState(:,c), vacancyState(p)%dotState(:,c)]))) then ! NaN occured in any dotState
if (.not. crystallite_localPlasticity(g,i,e)) then ! if broken non-local...
!$OMP CRITICAL (checkTodo)
crystallite_todo = crystallite_todo .and. crystallite_localPlasticity ! ...all non-locals skipped
@ -2599,10 +2594,8 @@ subroutine crystallite_integrateStateAdaptiveEuler()
if (crystallite_todo(g,i,e)) then
p = mappingConstitutive(2,g,i,e)
c = mappingConstitutive(1,g,i,e)
if ( any(plasticState(p)%dotState(:,c) /= plasticState(p)%dotState(:,c)) .or.&
any(damageState( p)%dotState(:,c) /= damageState( p)%dotState(:,c)) .or.&
any(thermalState(p)%dotState(:,c) /= thermalState(p)%dotState(:,c)) .or.&
any(vacancyState(p)%dotState(:,c) /= vacancyState(p)%dotState(:,c))) then ! NaN occured in dotState
if (any(prec_isNaN([plasticState(p)%dotState(:,c), damageState(p)%dotState(:,c), &
thermalState(p)%dotState(:,c), vacancyState(p)%dotState(:,c)]))) then ! NaN occured in any dotState
if (.not. crystallite_localPlasticity(g,i,e)) then ! if broken non-local...
!$OMP CRITICAL (checkTodo)
crystallite_todo = crystallite_todo .and. crystallite_localPlasticity ! ...all non-locals skipped
@ -2749,6 +2742,8 @@ end subroutine crystallite_integrateStateAdaptiveEuler
!> @brief integrate stress, and state with 1st order explicit Euler method
!--------------------------------------------------------------------------------------------------
subroutine crystallite_integrateStateEuler()
use prec, only: &
prec_isNaN
use debug, only: &
debug_level, &
debug_crystallite, &
@ -2828,10 +2823,8 @@ eIter = FEsolving_execElem(1:2)
if (crystallite_todo(g,i,e) .and. .not. crystallite_converged(g,i,e)) then
c = mappingConstitutive(1,g,i,e)
p = mappingConstitutive(2,g,i,e)
if ( any(plasticState(p)%dotState(:,c) /= plasticState(p)%dotState(:,c)) .or. &
any(damageState( p)%dotState(:,c) /= damageState( p)%dotState(:,c)) .or. &
any(thermalState(p)%dotState(:,c) /= thermalState(p)%dotState(:,c)) .or. &
any(vacancyState(p)%dotState(:,c) /= vacancyState(p)%dotState(:,c))) then ! NaN occured in dotState
if (any(prec_isNaN([plasticState(p)%dotState(:,c), damageState(p)%dotState(:,c), &
thermalState(p)%dotState(:,c), vacancyState(p)%dotState(:,c)]))) then ! NaN occured in any dotState
if (.not. crystallite_localPlasticity(g,i,e) .and. .not. numerics_timeSyncing) then ! if broken non-local...
!$OMP CRITICAL (checkTodo)
crystallite_todo = crystallite_todo .and. crystallite_localPlasticity ! ...all non-locals skipped
@ -2975,6 +2968,8 @@ end subroutine crystallite_integrateStateEuler
!> using Fixed Point Iteration to adapt the stepsize
!--------------------------------------------------------------------------------------------------
subroutine crystallite_integrateStateFPI()
use prec, only: &
prec_isNaN
use debug, only: &
debug_e, &
debug_i, &
@ -3116,11 +3111,8 @@ subroutine crystallite_integrateStateFPI()
if (crystallite_todo(g,i,e)) then
p = mappingConstitutive(2,g,i,e)
c = mappingConstitutive(1,g,i,e)
if ( any(plasticState(p)%dotState(:,c) /= plasticState(p)%dotState(:,c)) .or. &
any(damageState( p)%dotState(:,c) /= damageState( p)%dotState(:,c)) .or. &
any(thermalState(p)%dotState(:,c) /= thermalState(p)%dotState(:,c)) .or. &
any(vacancyState(p)%dotState(:,c) /= vacancyState(p)%dotState(:,c))) then !NaN occured in dotState
if (any(prec_isNaN([plasticState(p)%dotState(:,c), damageState(p)%dotState(:,c), &
thermalState(p)%dotState(:,c), vacancyState(p)%dotState(:,c)]))) then ! NaN occured in any dotState
if (.not. crystallite_localPlasticity(g,i,e)) then ! if broken is a non-local...
!$OMP CRITICAL (checkTodo)
crystallite_todo = crystallite_todo .and. crystallite_localPlasticity ! ...all non-locals done (and broken)
@ -3237,10 +3229,8 @@ subroutine crystallite_integrateStateFPI()
if (crystallite_todo(g,i,e) .and. .not. crystallite_converged(g,i,e)) then
p = mappingConstitutive(2,g,i,e)
c = mappingConstitutive(1,g,i,e)
if ( any(plasticState(p)%dotState(:,c) /= plasticState(p)%dotState(:,c)) .or. &
any(damageState( p)%dotState(:,c) /= damageState( p)%dotState(:,c)) .or. &
any(thermalState(p)%dotState(:,c) /= thermalState(p)%dotState(:,c)) .or. &
any(vacancyState(p)%dotState(:,c) /= vacancyState(p)%dotState(:,c))) then ! NaN occured in dotState
if (any(prec_isNaN([plasticState(p)%dotState(:,c), damageState(p)%dotState(:,c), &
thermalState(p)%dotState(:,c), vacancyState(p)%dotState(:,c)]))) then ! NaN occured in any dotState
crystallite_todo(g,i,e) = .false. ! ... skip me next time
if (.not. crystallite_localPlasticity(g,i,e)) then ! if me is non-local...
!$OMP CRITICAL (checkTodo)

View File

@ -9,9 +9,13 @@
!> @brief setting precision for real and int type depending on makros "FLOAT" and "INT"
!> @details setting precision for real and int type and for DAMASK_NaN. Definition is made
!! depending on makros "FLOAT" and "INT" defined during compilation
!! for details on NaN see https://software.intel.com/en-us/forums/topic/294680
!--------------------------------------------------------------------------------------------------
module prec
#ifndef __GFORTRAN__
use, intrinsic :: & ! unfortunately not in commonly used gfortran versions
IEEE_arithmetic
#endif
implicit none
private
@ -20,19 +24,17 @@ module prec
SPECTRAL SOLVER AND OWN FEM DO NOT SUPPORT SINGLE PRECISION, STOPPING COMPILATION
#endif
integer, parameter, public :: pReal = 4 !< floating point single precition (was selected_real_kind(6,37), number with 6 significant digits, up to 1e+-37)
#ifdef __INTEL_COMPILER
real(pReal), parameter, public :: DAMASK_NaN = Z'7F800001' !< quiet NaN for single precision (from http://www.hpc.unimelb.edu.au/doc/f90lrm/dfum_035.html, copy can be found in documentation/Code/Fortran)
#endif
#ifdef __GFORTRAN__
#ifndef __GFORTRAN__
real(pReal), parameter, public :: DAMASK_NaN = IEEE_value(IEEE_quiet_NaN) !< quiet NaN
#else
real(pReal), parameter, public :: DAMASK_NaN = real(Z'7F800001', pReal) !< quiet NaN for single precision (from http://www.hpc.unimelb.edu.au/doc/f90lrm/dfum_035.html, copy can be found in documentation/Code/Fortran)
#endif
#elif (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)
#ifdef __INTEL_COMPILER
real(pReal), parameter, public :: DAMASK_NaN = Z'7FF8000000000000' !< quiet NaN for double precision (from http://www.hpc.unimelb.edu.au/doc/f90lrm/dfum_035.html, copy can be found in documentation/Code/Fortran)
#endif
#ifdef __GFORTRAN__
real(pReal), parameter, public :: DAMASK_NaN = real(Z'7FF8000000000000', pReal) !< quiet NaN for double precision (from http://www.hpc.unimelb.edu.au/doc/f90lrm/dfum_035.html, copy can be found in documentation/Code/Fortran)
#ifndef __GFORTRAN__
real(pReal), parameter, public :: DAMASK_NaN = IEEE_value(IEEE_quiet_NaN) !< quiet NaN
#else
real(pReal), parameter, public :: DAMASK_NaN = real(Z'7FF8000000000000',pReal) !< quiet NaN for double precision (from http://www.hpc.unimelb.edu.au/doc/f90lrm/dfum_035.html, copy can be found in documentation/Code/Fortran)
#endif
#else
NO SUITABLE PRECISION FOR REAL SELECTED, STOPPING COMPILATION
@ -52,6 +54,7 @@ module prec
type, public :: p_vec !< variable length datatype used for storage of state
real(pReal), dimension(:), allocatable :: p
end type p_vec
type, public :: p_intvec
integer(pInt), dimension(:), allocatable :: p
end type p_intvec
@ -66,7 +69,7 @@ type, public :: p_intvec
nonlocal = .false. !< absolute tolerance for state integration
real(pReal), allocatable, dimension(:) :: &
atolState
real(pReal), pointer, dimension(:,:) :: &
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
state, & !< state
dotState !< state rate
real(pReal), allocatable, dimension(:,:) :: &
@ -88,7 +91,7 @@ type, public :: p_intvec
nSlip = 0_pInt , &
nTwin = 0_pInt, &
nTrans = 0_pInt
real(pReal), pointer, dimension(:,:) :: &
real(pReal), pointer, dimension(:,:), contiguous :: &
slipRate, & !< slip rate
accumulatedSlip !< accumulated plastic slip
end type
@ -112,16 +115,19 @@ type, public :: p_intvec
#endif
public :: &
prec_init
prec_init, &
prec_isNaN
contains
!--------------------------------------------------------------------------------------------------
!> @brief reporting precision and checking if DAMASK_NaN is set correctly
!--------------------------------------------------------------------------------------------------
subroutine prec_init
use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment)
use, intrinsic :: &
iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment)
implicit none
integer(pInt) :: worldrank = 0_pInt
#ifdef PETSc
@ -145,11 +151,29 @@ subroutine prec_init
write(6,'(a,i3)') ' Bytes for pInt: ',pInt
write(6,'(a,i3)') ' Bytes for pLongInt: ',pLongInt
write(6,'(a,e10.3)') ' NaN: ', DAMASK_NaN
write(6,'(a,l3,/)') ' NaN /= NaN: ',DAMASK_NaN/=DAMASK_NaN
write(6,'(a,l3)') ' NaN != NaN: ',DAMASK_NaN /= DAMASK_NaN
write(6,'(a,l3,/)') ' NaN check passed ',prec_isNAN(DAMASK_NaN)
endif mainProcess
if (DAMASK_NaN == DAMASK_NaN) call quit(9000)
if ((.not. prec_isNaN(DAMASK_NaN)) .or. (DAMASK_NaN == DAMASK_NaN)) call quit(9000)
end subroutine prec_init
!--------------------------------------------------------------------------------------------------
!> @brief figures out if a floating point number is NaN
! basically just a small wrapper, because gfortran < 4.9 does not have the IEEE module
!--------------------------------------------------------------------------------------------------
logical elemental function prec_isNaN(a)
implicit none
real(pReal), intent(in) :: a
#ifndef __GFORTRAN__
prec_isNaN = IEEE_is_NaN(a)
#else
prec_isNaN = isNaN(a)
#endif
end function prec_isNaN
end module prec