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

using IEEE_is_NAN and IEEE_quiet_NaN instead of hand-written solution, will not work for gfortran < 5

This commit is contained in:
Martin Diehl 2017-05-04 00:32:44 +02:00
parent 923732ab4b
commit b35ff67f99
7 changed files with 92 additions and 139 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

70
src/crystallite.f90
View File

@ -1232,8 +1232,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, intrinsic :: &
IEEE_arithmetic
use numerics, only: &
numerics_integrationMode
use debug, only: &
@ -1331,9 +1331,9 @@ subroutine crystallite_integrateStateRK4()
if (crystallite_todo(g,i,e)) then
c = phasememberAt(g,i,e)
p = phaseAt(g,i,e)
NaN = any(prec_isNaN(plasticState(p)%dotState(:,c)))
NaN = any(IEEE_is_NaN(plasticState(p)%dotState(:,c)))
do mySource = 1_pInt, phase_Nsources(p)
NaN = NaN .or. any(prec_isNaN(sourceState(p)%p(mySource)%dotState(:,c)))
NaN = NaN .or. any(IEEE_is_NaN(sourceState(p)%p(mySource)%dotState(:,c)))
enddo
if (NaN) then ! NaN occured in any dotState
if (.not. crystallite_localPlasticity(g,i,e)) then ! if broken non-local...
@ -1475,9 +1475,9 @@ subroutine crystallite_integrateStateRK4()
p = phaseAt(g,i,e)
c = phasememberAt(g,i,e)
NaN = any(prec_isNaN(plasticState(p)%dotState(:,c)))
NaN = any(IEEE_is_NaN(plasticState(p)%dotState(:,c)))
do mySource = 1_pInt, phase_Nsources(p)
NaN = NaN .or. any(prec_isNaN(sourceState(p)%p(mySource)%dotState(:,c)))
NaN = NaN .or. any(IEEE_is_NaN(sourceState(p)%p(mySource)%dotState(:,c)))
enddo
if (NaN) then ! NaN occured in any dotState
if (.not. crystallite_localPlasticity(g,i,e)) then ! if broken non-local...
@ -1528,8 +1528,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, intrinsic :: &
IEEE_arithmetic
use debug, only: &
debug_level, &
debug_crystallite, &
@ -1647,9 +1647,9 @@ subroutine crystallite_integrateStateRKCK45()
if (crystallite_todo(g,i,e)) then
cc = phasememberAt(g,i,e)
p = phaseAt(g,i,e)
NaN = any(prec_isNaN(plasticState(p)%dotState(:,cc)))
NaN = any(IEEE_is_NaN(plasticState(p)%dotState(:,cc)))
do mySource = 1_pInt, phase_Nsources(p)
NaN = NaN .or. any(prec_isNaN(sourceState(p)%p(mySource)%dotState(:,cc)))
NaN = NaN .or. any(IEEE_is_NaN(sourceState(p)%p(mySource)%dotState(:,cc)))
enddo
if (NaN) then ! NaN occured in any dotState
if (.not. crystallite_localPlasticity(g,i,e)) then ! if broken non-local...
@ -1801,9 +1801,9 @@ subroutine crystallite_integrateStateRKCK45()
p = phaseAt(g,i,e)
cc = phasememberAt(g,i,e)
NaN = any(prec_isNaN(plasticState(p)%dotState(:,cc)))
NaN = any(IEEE_is_NaN(plasticState(p)%dotState(:,cc)))
do mySource = 1_pInt, phase_Nsources(p)
NaN = NaN .or. any(prec_isNaN(sourceState(p)%p(mySource)%dotState(:,cc)))
NaN = NaN .or. any(IEEE_is_NaN(sourceState(p)%p(mySource)%dotState(:,cc)))
enddo
if (NaN) then ! NaN occured in any dotState
if (.not. crystallite_localPlasticity(g,i,e)) then ! if broken non-local...
@ -2031,8 +2031,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, intrinsic :: &
IEEE_arithmetic
use debug, only: &
debug_level, &
debug_crystallite, &
@ -2133,9 +2133,9 @@ subroutine crystallite_integrateStateAdaptiveEuler()
if (crystallite_todo(g,i,e)) then
p = phaseAt(g,i,e)
c = phasememberAt(g,i,e)
NaN = any(prec_isNaN(plasticState(p)%dotState(:,c)))
NaN = any(IEEE_is_NaN(plasticState(p)%dotState(:,c)))
do mySource = 1_pInt, phase_Nsources(p)
NaN = NaN .or. any(prec_isNaN(sourceState(p)%p(mySource)%dotState(:,c)))
NaN = NaN .or. any(IEEE_is_NaN(sourceState(p)%p(mySource)%dotState(:,c)))
enddo
if (NaN) then ! NaN occured in any dotState
if (.not. crystallite_localPlasticity(g,i,e)) then ! if broken non-local...
@ -2254,9 +2254,9 @@ subroutine crystallite_integrateStateAdaptiveEuler()
if (crystallite_todo(g,i,e)) then
p = phaseAt(g,i,e)
c = phasememberAt(g,i,e)
NaN = any(prec_isNaN(plasticState(p)%dotState(:,c)))
NaN = any(IEEE_is_NaN(plasticState(p)%dotState(:,c)))
do mySource = 1_pInt, phase_Nsources(p)
NaN = NaN .or. any(prec_isNaN(sourceState(p)%p(mySource)%dotState(:,c)))
NaN = NaN .or. any(IEEE_is_NaN(sourceState(p)%p(mySource)%dotState(:,c)))
enddo
if (NaN) then ! NaN occured in any dotState
if (.not. crystallite_localPlasticity(g,i,e)) then ! if broken non-local...
@ -2391,8 +2391,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, intrinsic :: &
IEEE_arithmetic
use debug, only: &
debug_level, &
debug_crystallite, &
@ -2471,9 +2471,9 @@ eIter = FEsolving_execElem(1:2)
if (crystallite_todo(g,i,e) .and. .not. crystallite_converged(g,i,e)) then
c = phasememberAt(g,i,e)
p = phaseAt(g,i,e)
NaN = any(prec_isNaN(plasticState(p)%dotState(:,c)))
NaN = any(IEEE_is_NaN(plasticState(p)%dotState(:,c)))
do mySource = 1_pInt, phase_Nsources(p)
NaN = NaN .or. any(prec_isNaN(sourceState(p)%p(mySource)%dotState(:,c)))
NaN = NaN .or. any(IEEE_is_NaN(sourceState(p)%p(mySource)%dotState(:,c)))
enddo
if (NaN) then ! NaN occured in any dotState
if (.not. crystallite_localPlasticity(g,i,e) .and. .not. numerics_timeSyncing) then ! if broken non-local...
@ -2614,8 +2614,8 @@ end subroutine crystallite_integrateStateEuler
!> using Fixed Point Iteration to adapt the stepsize
!--------------------------------------------------------------------------------------------------
subroutine crystallite_integrateStateFPI()
use prec, only: &
prec_isNaN
use, intrinsic :: &
IEEE_arithmetic
use debug, only: &
debug_e, &
debug_i, &
@ -2737,9 +2737,9 @@ subroutine crystallite_integrateStateFPI()
if (crystallite_todo(g,i,e)) then
p = phaseAt(g,i,e)
c = phasememberAt(g,i,e)
NaN = any(prec_isNaN(plasticState(p)%dotState(:,c)))
NaN = any(IEEE_is_NaN(plasticState(p)%dotState(:,c)))
do mySource = 1_pInt, phase_Nsources(p)
NaN = NaN .or. any(prec_isNaN(sourceState(p)%p(mySource)%dotState(:,c)))
NaN = NaN .or. any(IEEE_is_NaN(sourceState(p)%p(mySource)%dotState(:,c)))
enddo
if (NaN) then ! NaN occured in any dotState
if (.not. crystallite_localPlasticity(g,i,e)) then ! if broken is a non-local...
@ -2851,9 +2851,9 @@ subroutine crystallite_integrateStateFPI()
if (crystallite_todo(g,i,e) .and. .not. crystallite_converged(g,i,e)) then
p = phaseAt(g,i,e)
c = phasememberAt(g,i,e)
NaN = any(prec_isNaN(plasticState(p)%dotState(:,c)))
NaN = any(IEEE_is_NaN(plasticState(p)%dotState(:,c)))
do mySource = 1_pInt, phase_Nsources(p)
NaN = NaN .or. any(prec_isNaN(sourceState(p)%p(mySource)%dotState(:,c)))
NaN = NaN .or. any(IEEE_is_NaN(sourceState(p)%p(mySource)%dotState(:,c)))
enddo
if (NaN) then ! NaN occured in any dotState
crystallite_todo(g,i,e) = .false. ! ... skip me next time
@ -3062,8 +3062,9 @@ end subroutine crystallite_integrateStateFPI
!> returns true, if state jump was successfull or not needed. false indicates NaN in delta state
!--------------------------------------------------------------------------------------------------
logical function crystallite_stateJump(ipc,ip,el)
use, intrinsic :: &
IEEE_arithmetic
use prec, only: &
prec_isNaN, &
dNeq0
use debug, only: &
debug_level, &
@ -3098,7 +3099,7 @@ logical function crystallite_stateJump(ipc,ip,el)
p = phaseAt(ipc,ip,el)
call constitutive_collectDeltaState(crystallite_Tstar_v(1:6,ipc,ip,el), crystallite_Fe(1:3,1:3,ipc,ip,el), ipc,ip,el)
mySizePlasticDeltaState = plasticState(p)%sizeDeltaState
if( any(prec_isNaN(plasticState(p)%deltaState(:,c)))) then ! NaN occured in deltaState
if( any(IEEE_is_NaN(plasticState(p)%deltaState(:,c)))) then ! NaN occured in deltaState
crystallite_stateJump = .false.
return
endif
@ -3106,7 +3107,7 @@ logical function crystallite_stateJump(ipc,ip,el)
plasticState(p)%deltaState(1:mySizePlasticDeltaState,c)
do mySource = 1_pInt, phase_Nsources(p)
mySizeSourceDeltaState = sourceState(p)%p(mySource)%sizeDeltaState
if( any(prec_isNaN(sourceState(p)%p(mySource)%deltaState(:,c)))) then ! NaN occured in deltaState
if( any(IEEE_is_NaN(sourceState(p)%p(mySource)%deltaState(:,c)))) then ! NaN occured in deltaState
crystallite_stateJump = .false.
return
endif
@ -3169,9 +3170,10 @@ logical function crystallite_integrateStress(&
el,& ! element number
timeFraction &
)
use, intrinsic :: &
IEEE_arithmetic
use prec, only: pLongInt, &
tol_math_check, &
prec_isNaN, &
dEq0
use numerics, only: nStress, &
aTol_crystalliteStress, &
@ -3430,7 +3432,7 @@ logical function crystallite_integrateStress(&
aTol_crystalliteStress) ! minimum lower cutoff
residuumLp = Lpguess - Lp_constitutive
if (any(prec_isNaN(residuumLp))) then ! NaN in residuum...
if (any(IEEE_is_NaN(residuumLp))) then ! NaN in residuum...
#ifndef _OPENMP
if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) &
write(6,'(a,i8,1x,a,i8,a,1x,i2,1x,i3,a,i3,a)') '<< CRYST >> integrateStress encountered NaN at el (elFE) ip ipc ', &
@ -3520,7 +3522,7 @@ logical function crystallite_integrateStress(&
aTolLi = max(rTol_crystalliteStress * max(norm2(Liguess),norm2(Li_constitutive)), & ! absolute tolerance from largest acceptable relative error
aTol_crystalliteStress) ! minimum lower cutoff
residuumLi = Liguess - Li_constitutive
if (any(prec_isNaN(residuumLi))) then ! NaN in residuum...
if (any(IEEE_is_NaN(residuumLi))) then ! NaN in residuum...
return ! ...me = .false. to inform integrator about problem
elseif (norm2(residuumLi) < aTolLi) then ! converged if below absolute tolerance
exit LiLoop ! ...leave iteration loop

6
src/lattice.f90
View File

@ -2170,8 +2170,8 @@ pure function lattice_symmetrize33(struct,T33)
!> @brief figures whether unit quat falls into stereographic standard triangle
!--------------------------------------------------------------------------------------------------
logical pure function lattice_qInSST(Q, struct)
use prec, only: &
prec_isNaN
use, intrinsic :: &
IEEE_arithmetic
use math, only: &
math_qToRodrig
@ -2181,7 +2181,7 @@ logical pure function lattice_qInSST(Q, struct)
real(pReal), dimension(3) :: Rodrig ! Rodrigues vector of Q
Rodrig = math_qToRodrig(Q)
if (any(prec_isNaN(Rodrig))) then
if (any(IEEE_is_NaN(Rodrig))) then
lattice_qInSST = .false.
else
select case (struct)

34
src/math.f90
View File

@ -1666,15 +1666,16 @@ end function math_qToEulerAxisAngle
!> @brief Rodrigues vector (x, y, z) from unit quaternion (w+ix+jy+kz)
!--------------------------------------------------------------------------------------------------
pure function math_qToRodrig(Q)
use, intrinsic :: &
IEEE_arithmetic
use prec, only: &
DAMASK_NaN, &
tol_math_check
implicit none
real(pReal), dimension(4), intent(in) :: Q
real(pReal), dimension(3) :: math_qToRodrig
math_qToRodrig = merge(Q(2:4)/Q(1),DAMASK_NaN,abs(Q(1)) > tol_math_check) ! NaN for 180 deg since Rodrig is unbound
math_qToRodrig = merge(Q(2:4)/Q(1),IEEE_value(1.0_pReal,IEEE_quiet_NaN),abs(Q(1)) > tol_math_check)! NaN for 180 deg since Rodrig is unbound
end function math_qToRodrig
@ -2095,8 +2096,8 @@ end function math_rotationalPart33
! will return NaN on error
!--------------------------------------------------------------------------------------------------
function math_eigenvaluesSym(m)
use prec, only: &
DAMASK_NaN
use, intrinsic :: &
IEEE_arithmetic
implicit none
real(pReal), dimension(:,:), intent(in) :: m
@ -2109,7 +2110,7 @@ function math_eigenvaluesSym(m)
vectors = m ! copy matrix to input (doubles as output) array
call dsyev('N','U',size(m,1),vectors,size(m,1),math_eigenvaluesSym,work,(64+2)*size(m,1),info)
if (info /= 0_pInt) math_eigenvaluesSym = DAMASK_NaN
if (info /= 0_pInt) math_eigenvaluesSym = IEEE_value(1.0_pReal,IEEE_quiet_NaN)
end function math_eigenvaluesSym
@ -2701,29 +2702,13 @@ pure function math_rotate_forward3333(tensor,rot_tensor)
end function math_rotate_forward3333
!--------------------------------------------------------------------------------------------------
!> @brief calculate average of tensor field
!--------------------------------------------------------------------------------------------------
function math_tensorAvg(field)
implicit none
real(pReal), dimension(3,3) :: math_tensorAvg
real(pReal), intent(in), dimension(:,:,:,:,:) :: field
real(pReal) :: wgt
wgt = 1.0_pReal/real(size(field,3)*size(field,4)*size(field,5), pReal)
math_tensorAvg = sum(sum(sum(field,dim=5),dim=4),dim=3)*wgt
end function math_tensorAvg
!--------------------------------------------------------------------------------------------------
!> @brief limits a scalar value to a certain range (either one or two sided)
! Will return NaN if left > right
!--------------------------------------------------------------------------------------------------
real(pReal) pure function math_limit(a, left, right)
use prec, only: &
DAMASK_NaN
use, intrinsic :: &
IEEE_arithmetic
implicit none
real(pReal), intent(in) :: a
@ -2735,7 +2720,8 @@ real(pReal) pure function math_limit(a, left, right)
merge(right, huge(a), present(right)) &
)
if (present(left) .and. present(right)) math_limit = merge (DAMASK_NaN,math_limit, left>right)
if (present(left) .and. present(right)) &
math_limit = merge (IEEE_value(1.0_pReal,IEEE_quiet_NaN),math_limit, left>right)
end function math_limit

17
src/plastic_isotropic.f90
View File

@ -9,8 +9,7 @@
module plastic_isotropic
use prec, only: &
pReal,&
pInt, &
DAMASK_NaN
pInt
implicit none
private
@ -36,14 +35,14 @@ module plastic_isotropic
integer(kind(undefined_ID)), allocatable, dimension(:) :: &
outputID
real(pReal) :: &
fTaylor = DAMASK_NaN, &
tau0 = DAMASK_NaN, &
gdot0 = DAMASK_NaN, &
n = DAMASK_NaN, &
h0 = DAMASK_NaN, &
fTaylor, &
tau0, &
gdot0, &
n, &
h0, &
h0_slopeLnRate = 0.0_pReal, &
tausat = DAMASK_NaN, &
a = DAMASK_NaN, &
tausat, &
a, &
aTolFlowstress = 1.0_pReal, &
aTolShear = 1.0e-6_pReal, &
tausat_SinhFitA= 0.0_pReal, &

10
src/plastic_nonlocal.f90
View File

@ -2377,9 +2377,9 @@ end subroutine plastic_nonlocal_deltaState
!---------------------------------------------------------------------------------------------------
subroutine plastic_nonlocal_dotState(Tstar_v, Fe, Fp, Temperature, &
timestep,subfrac, ip,el)
use prec, only: DAMASK_NaN, &
dNeq0, &
use, intrinsic :: &
IEEE_arithmetic
use prec, only: dNeq0, &
dNeq, &
dEq0
use numerics, only: numerics_integrationMode, &
@ -2701,7 +2701,7 @@ if (.not. phase_localPlasticity(material_phase(1_pInt,ip,el))) then
write(6,'(a)') '<< CONST >> enforcing cutback !!!'
endif
#endif
plasticState(p)%dotState = DAMASK_NaN ! -> return NaN and, hence, enforce cutback
plasticState(p)%dotState = IEEE_value(1.0_pReal,IEEE_quiet_NaN) ! -> return NaN and, hence, enforce cutback
return
endif
@ -2984,7 +2984,7 @@ if ( any(rhoSglOriginal(1:ns,1:4) + rhoDot(1:ns,1:4) * timestep < -aTolRho(in
write(6,'(a)') '<< CONST >> enforcing cutback !!!'
endif
#endif
plasticState(p)%dotState = DAMASK_NaN
plasticState(p)%dotState = IEEE_value(1.0_pReal,IEEE_quiet_NaN)
return
else
forall (s = 1:ns, t = 1_pInt:4_pInt)

36
src/prec.f90
View File

@ -5,22 +5,12 @@
!> @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
!> @details setting precision for real and int type and for DAMASK_NaN. Definition is made
!! depending on makro "INT" defined during compilation
!! for details on NaN see https://software.intel.com/en-us/forums/topic/294680
!--------------------------------------------------------------------------------------------------
module prec
#if !(defined(__GFORTRAN__) && __GNUC__ < 5)
use, intrinsic :: & ! unfortunately not avialable in gfortran <= 5
IEEE_arithmetic
#endif
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)
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)
#else
NO SUITABLE PRECISION FOR REAL SELECTED, STOPPING COMPILATION
#endif
@ -106,7 +96,6 @@ module prec
public :: &
prec_init, &
prec_isNaN, &
dEq, &
dEq0, &
cEq, &
@ -118,7 +107,7 @@ contains
!--------------------------------------------------------------------------------------------------
!> @brief reporting precision and checking if DAMASK_NaN is set correctly
!> @brief reporting precision
!--------------------------------------------------------------------------------------------------
subroutine prec_init
use, intrinsic :: &
@ -133,36 +122,13 @@ subroutine prec_init
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,e10.3)') ' NaN: ', DAMASK_NaN
write(6,'(a,l3)') ' NaN != NaN: ',DAMASK_NaN /= DAMASK_NaN
write(6,'(a,l3,/)') ' NaN check passed ',prec_isNAN(DAMASK_NaN)
if ((.not. prec_isNaN(DAMASK_NaN)) .or. (DAMASK_NaN == DAMASK_NaN)) call quit(9000)
realloc_lhs_test = [1_pInt,2_pInt]
if (realloc_lhs_test(2)/=2_pInt) call quit(9000)
end subroutine prec_init
!--------------------------------------------------------------------------------------------------
!> @brief figures out if a floating point number is NaN
! basically just a small wrapper, because gfortran < 5.0 does not have the IEEE module
!--------------------------------------------------------------------------------------------------
logical elemental pure function prec_isNaN(a)
implicit none
real(pReal), intent(in) :: a
#if (defined(__GFORTRAN__) && __GNUC__ < 5)
intrinsic :: isNaN
prec_isNaN = isNaN(a)
#else
prec_isNaN = IEEE_is_NaN(a)
#endif
end function prec_isNaN
!--------------------------------------------------------------------------------------------------
!> @brief equality comparison for float with double precision
! replaces "==" but for certain (relative) tolerance. Counterpart to dNeq

6
src/spectral_utilities.f90
View File

@ -741,8 +741,8 @@ end function utilities_curlRMS
!> @brief calculates mask compliance tensor used to adjust F to fullfill stress BC
!--------------------------------------------------------------------------------------------------
function utilities_maskedCompliance(rot_BC,mask_stress,C)
use prec, only: &
prec_isNaN
use, intrinsic :: &
IEEE_arithmetic
use IO, only: &
IO_error
use math, only: &
@ -794,7 +794,7 @@ function utilities_maskedCompliance(rot_BC,mask_stress,C)
endif; enddo; endif; enddo
call math_invert(size_reduced, c_reduced, s_reduced, errmatinv) ! invert reduced stiffness
if (any(prec_isNaN(s_reduced))) errmatinv = .true.
if (any(IEEE_is_NaN(s_reduced))) errmatinv = .true.
if(errmatinv) call IO_error(error_ID=400_pInt,ext_msg='utilities_maskedCompliance')
temp99_Real = 0.0_pReal ! fill up compliance with zeros
k = 0_pInt