From 5c1804e77d40dd7684eda9ef0b5b2f07deb82f46 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Mon, 13 Apr 2015 10:02:52 +0000 Subject: [PATCH] improved NaN checks --- code/Makefile | 7 +++- code/compilation_info.f90 | 6 ++-- code/crystallite.f90 | 74 +++++++++++++++++---------------------- code/prec.f90 | 58 +++++++++++++++++++++--------- 4 files changed, 82 insertions(+), 63 deletions(-) diff --git a/code/Makefile b/code/Makefile index 51c5b4900..7384d60e4 100644 --- a/code/Makefile +++ b/code/Makefile @@ -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) \ diff --git a/code/compilation_info.f90 b/code/compilation_info.f90 index fc1dae1f2..64e6b136c 100644 --- a/code/compilation_info.f90 +++ b/code/compilation_info.f90 @@ -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,*) diff --git a/code/crystallite.f90 b/code/crystallite.f90 index 660599efa..2b7e0ced6 100644 --- a/code/crystallite.f90 +++ b/code/crystallite.f90 @@ -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) diff --git a/code/prec.f90 b/code/prec.f90 index e86b6f277..307855510 100644 --- a/code/prec.f90 +++ b/code/prec.f90 @@ -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