diff --git a/src/crystallite.f90 b/src/crystallite.f90 index efeab9dcd..8d2ae5f23 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -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: & @@ -1316,7 +1316,7 @@ subroutine crystallite_integrateStateRK4() ! first Runge-Kutta step !$OMP PARALLEL !$OMP DO - do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains + do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains if (crystallite_todo(g,i,e)) & call constitutive_collectDotState(crystallite_Tstar_v(1:6,g,i,e), & crystallite_Fe, & @@ -1326,22 +1326,22 @@ subroutine crystallite_integrateStateRK4() !$OMP ENDDO !$OMP DO PRIVATE(p,c,NaN) - do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains + do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains !$OMP FLUSH(crystallite_todo) 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... + if (NaN) 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 @@ -1356,7 +1356,7 @@ subroutine crystallite_integrateStateRK4() !$OMP PARALLEL !$OMP DO PRIVATE(p,c) - do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains + do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains if (crystallite_todo(g,i,e)) then p = phaseAt(g,i,e) c = phasememberAt(g,i,e) @@ -1371,7 +1371,7 @@ subroutine crystallite_integrateStateRK4() !$OMP ENDDO !$OMP DO PRIVATE(mySizePlasticDotState,mySizeSourceDotState,p,c) - do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains + do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains if (crystallite_todo(g,i,e)) then p = phaseAt(g,i,e) @@ -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, & @@ -2732,22 +2732,22 @@ subroutine crystallite_integrateStateFPI() !$OMP ENDDO !$OMP DO PRIVATE(p,c,NaN) - do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains + do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains !$OMP FLUSH(crystallite_todo) 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... + if (NaN) 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) + crystallite_todo = crystallite_todo .and. crystallite_localPlasticity ! ...all non-locals done (and broken) !$OMP END CRITICAL (checkTodo) - else ! broken one was local... - crystallite_todo(g,i,e) = .false. ! ... done (and broken) + else ! broken one was local... + crystallite_todo(g,i,e) = .false. ! ... done (and broken) endif endif endif @@ -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, & @@ -3426,11 +3428,11 @@ logical function crystallite_integrateStress(& !* update current residuum and check for convergence of loop - aTolLp = max(rTol_crystalliteStress * max(norm2(Lpguess),norm2(Lp_constitutive)), & ! absolute tolerance from largest acceptable relative error - aTol_crystalliteStress) ! minimum lower cutoff + aTolLp = max(rTol_crystalliteStress * max(norm2(Lpguess),norm2(Lp_constitutive)), & ! absolute tolerance from largest acceptable relative error + 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 ', & @@ -3438,16 +3440,16 @@ logical function crystallite_integrateStress(& ' ; iteration ', NiterationStressLp,& ' >> returning..!' #endif - return ! ...me = .false. to inform integrator about problem - elseif (norm2(residuumLp) < aTolLp) then ! converged if below absolute tolerance - exit LpLoop ! ...leave iteration loop + return ! ...me = .false. to inform integrator about problem + elseif (norm2(residuumLp) < aTolLp) then ! converged if below absolute tolerance + exit LpLoop ! ...leave iteration loop elseif ( NiterationStressLp == 1_pInt & - .or. norm2(residuumLp) < norm2(residuumLp_old)) then ! not converged, but improved norm of residuum (always proceed in first iteration)... - residuumLp_old = residuumLp ! ...remember old values and... + .or. norm2(residuumLp) < norm2(residuumLp_old)) then ! not converged, but improved norm of residuum (always proceed in first iteration)... + residuumLp_old = residuumLp ! ...remember old values and... Lpguess_old = Lpguess - steplengthLp = steplengthLp0 ! ...proceed with normal step length (calculate new search direction) - else ! not converged and residuum not improved... - steplengthLp = 0.5_pReal * steplengthLp ! ...try with smaller step length in same direction + steplengthLp = steplengthLp0 ! ...proceed with normal step length (calculate new search direction) + else ! not converged and residuum not improved... + steplengthLp = 0.5_pReal * steplengthLp ! ...try with smaller step length in same direction Lpguess = Lpguess_old + steplengthLp * deltaLp cycle LpLoop endif @@ -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 diff --git a/src/lattice.f90 b/src/lattice.f90 index 7726c772d..a970ed85a 100644 --- a/src/lattice.f90 +++ b/src/lattice.f90 @@ -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) diff --git a/src/math.f90 b/src/math.f90 index 9629401ca..3a0533ddf 100644 --- a/src/math.f90 +++ b/src/math.f90 @@ -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 @@ -2092,11 +2093,11 @@ end function math_rotationalPart33 !-------------------------------------------------------------------------------------------------- !> @brief Eigenvalues of symmetric matrix m -! will return NaN on error +! 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 diff --git a/src/plastic_isotropic.f90 b/src/plastic_isotropic.f90 index a82b9c820..bea9c616e 100644 --- a/src/plastic_isotropic.f90 +++ b/src/plastic_isotropic.f90 @@ -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, & diff --git a/src/plastic_nonlocal.f90 b/src/plastic_nonlocal.f90 index 0c0a48523..0a8c4c3f9 100644 --- a/src/plastic_nonlocal.f90 +++ b/src/plastic_nonlocal.f90 @@ -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) diff --git a/src/prec.f90 b/src/prec.f90 index c6ae656b8..671e15990 100644 --- a/src/prec.f90 +++ b/src/prec.f90 @@ -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 diff --git a/src/spectral_utilities.f90 b/src/spectral_utilities.f90 index a33c7a9f5..0773a9065 100644 --- a/src/spectral_utilities.f90 +++ b/src/spectral_utilities.f90 @@ -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