using dEq, dNeq, cEq, dEq

code/crystallite.f90

@@ -843,7 +843,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
            !$OMP DO PRIVATE(neighboring_e,neighboring_i)
            do e = FEsolving_execElem(1),FEsolving_execElem(2)
              do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e)
-               if (.not. crystallite_localPlasticity(1,i,e) .and. abs(crystallite_subFrac(1,i,e)) > tiny(0.0_pReal)) then
+               if (.not. crystallite_localPlasticity(1,i,e) .and. dNeq(crystallite_subFrac(1,i,e),0.0_pReal)) then
                  do n = 1_pInt,FE_NipNeighbors(FE_celltype(FE_geomtype(mesh_element(2,e))))
                    neighboring_e = mesh_ipNeighborhood(1,n,i,e)
                    neighboring_i = mesh_ipNeighborhood(2,n,i,e)
@@ -3623,6 +3623,7 @@ logical function crystallite_integrateStress(&
 
  !* inversion of Fi_current...
 
+ invFi_current = math_inv33(Fi_current)
  failedInversionFi: if (all(dEq(invFi_current,0.0_pReal))) then
 #ifndef _OPENMP
    if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) then
@@ -3883,7 +3884,7 @@ logical function crystallite_integrateStress(&
  invFp_new = math_mul33x33(invFp_current,B)
  invFp_new = invFp_new / math_det33(invFp_new)**(1.0_pReal/3.0_pReal)                               ! regularize by det
  Fp_new = math_inv33(invFp_new)
- failedInversionFp2: if (all(dEq(invFp_new,0.0_pReal))) then
+ failedInversionInvFp: if (all(dEq(Fp_new,0.0_pReal))) then
 #ifndef _OPENMP
    if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) then
      write(6,'(a,i8,1x,a,i8,a,1x,i2,1x,i3,a,i3)') '<< CRYST >> integrateStress failed on invFp_new inversion at el ip ipc ',&
@@ -3895,7 +3896,7 @@ logical function crystallite_integrateStress(&
    endif
 #endif
    return
- endif failedInversionFp2
+ endif failedInversionInvFp
  Fe_new = math_mul33x33(math_mul33x33(Fg_new,invFp_new),invFi_new)    ! calc resulting Fe
 
  !* calculate 1st Piola-Kirchhoff stress

code/homogenization_RGC.f90

@@ -386,6 +386,8 @@ end subroutine homogenization_RGC_partitionDeformation
 ! "happy" with result
 !--------------------------------------------------------------------------------------------------
 function homogenization_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el)
+ use prec, only: &
+   dEq
  use debug, only: &
    debug_level, &
    debug_homogenization,&
@@ -441,10 +443,10 @@ function homogenization_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el)
  real(pReal), dimension(:,:), allocatable :: tract,jmatrix,jnverse,smatrix,pmatrix,rmatrix
  real(pReal), dimension(:), allocatable   :: resid,relax,p_relax,p_resid,drelax
  
- if(abs(dt) < tiny(0.0_pReal)) then                                                                  ! zero time step
+ zeroTimeStep: if(dEq(dt,0.0_pReal)) then
-   homogenization_RGC_updateState = .true.                                                           ! pretend everything is fine and return
+   homogenization_RGC_updateState = .true.                                                          ! pretend everything is fine and return
    return                                                                    
- endif
+ endif zeroTimeStep
 
 !--------------------------------------------------------------------------------------------------
 ! get the dimension of the cluster (grains and interfaces)

code/math.f90

@@ -761,6 +761,8 @@ end function math_inv33
 !   returns error if not possible, i.e. if det close to zero
 !--------------------------------------------------------------------------------------------------
 pure subroutine math_invert33(A, InvA, DetA, error)
+ use prec, only: &
+   dEq
 
  implicit none
  logical, intent(out) :: error
@@ -774,7 +776,7 @@ pure subroutine math_invert33(A, InvA, DetA, error)
 
  DetA = A(1,1) * InvA(1,1) + A(1,2) * InvA(2,1) + A(1,3) * InvA(3,1)
 
- if (abs(DetA) <= tiny(DetA)) then
+ if (dEq(DetA,0.0_pReal)) then
    InvA = 0.0_pReal
    error = .true.
  else
@@ -1279,6 +1281,8 @@ end function math_qNorm
 !> @brief quaternion inversion
 !--------------------------------------------------------------------------------------------------
 pure function math_qInv(Q)
+ use prec, only: &
+   dNeq
 
  implicit none
  real(pReal), dimension(4), intent(in) ::  Q
@@ -1288,8 +1292,7 @@ pure function math_qInv(Q)
  math_qInv = 0.0_pReal
 
  squareNorm = math_qDot(Q,Q)
- if (abs(squareNorm) > tiny(squareNorm)) &
+ if (dNeq(squareNorm,0.0_pReal)) math_qInv = math_qConj(Q) / squareNorm
-   math_qInv = math_qConj(Q) / squareNorm
 
 end function math_qInv
 
@@ -2093,6 +2096,8 @@ end function math_eigenvectorBasisSym33
 !> @brief rotational part from polar decomposition of 33 tensor m
 !--------------------------------------------------------------------------------------------------
 function math_rotationalPart33(m)
+ use prec, only: &
+   dEq
  use IO, only: &
    IO_warning
 
@@ -2104,12 +2109,12 @@ function math_rotationalPart33(m)
  U = math_eigenvectorBasisSym33(math_mul33x33(transpose(m),m))
  Uinv = math_inv33(U)
 
- if (all(abs(Uinv) <= tiny(Uinv))) then                                                             ! math_inv33 returns zero when failed, avoid floating point equality comparison
+ inversionFailed: if (all(dEq(Uinv,0.0_pReal))) then
    math_rotationalPart33 = math_I3
    call IO_warning(650_pInt)
- else
+ else inversionFailed
    math_rotationalPart33 = math_mul33x33(m,Uinv)
- endif
+ endif inversionFailed
 
 end function math_rotationalPart33
 

code/plastic_disloUCLA.f90

@@ -981,7 +981,8 @@ end subroutine plastic_disloUCLA_LpAndItsTangent
 !--------------------------------------------------------------------------------------------------
 subroutine plastic_disloUCLA_dotState(Tstar_v,Temperature,ipc,ip,el)
  use prec, only: &
-   tol_math_check
+   tol_math_check, &
+   dEq
  use math, only: &
    pi
  use material, only: &
@@ -1119,7 +1120,7 @@ subroutine plastic_disloUCLA_dotState(Tstar_v,Temperature,ipc,ip,el)
      !* Dipole formation
      EdgeDipMinDistance = &
        plastic_disloUCLA_CEdgeDipMinDistance(instance)*plastic_disloUCLA_burgersPerSlipSystem(j,instance)
-     if (abs(tau_slip_pos) <= tiny(0.0_pReal)) then
+     if (dEq(tau_slip_pos,0.0_pReal)) then
        DotRhoDipFormation = 0.0_pReal
      else
        EdgeDipDistance = &
@@ -1147,7 +1148,7 @@ subroutine plastic_disloUCLA_dotState(Tstar_v,Temperature,ipc,ip,el)
         plastic_disloUCLA_CAtomicVolume(instance)*plastic_disloUCLA_burgersPerSlipSystem(j,instance)**(3.0_pReal)
      VacancyDiffusion = &
         plastic_disloUCLA_D0(instance)*exp(-plastic_disloUCLA_Qsd(instance)/(kB*Temperature))
-     if (abs(tau_slip_pos) <= tiny(0.0_pReal)) then
+     if (dEq(tau_slip_pos,0.0_pReal)) then
        DotRhoEdgeDipClimb = 0.0_pReal
      else
        ClimbVelocity = &
@@ -1180,7 +1181,8 @@ end subroutine plastic_disloUCLA_dotState
 !--------------------------------------------------------------------------------------------------
 function plastic_disloUCLA_postResults(Tstar_v,Temperature,ipc,ip,el)
  use prec, only: &
-   tol_math_check
+   tol_math_check, &
+   dEq
  use math, only: &
    pi
  use material, only: &
@@ -1408,7 +1410,7 @@ function plastic_disloUCLA_postResults(Tstar_v,Temperature,ipc,ip,el)
           c = c + ns
         elseif(plastic_disloUCLA_outputID(o,instance) == stress_exponent_ID) then
           do j = 1_pInt, ns
-            if (abs(gdot_slip_pos(j)+gdot_slip_neg(j))<=tiny(0.0_pReal)) then
+            if (dEq(gdot_slip_pos(j)+gdot_slip_neg(j),0.0_pReal)) then
               plastic_disloUCLA_postResults(c+j) = 0.0_pReal
             else
               plastic_disloUCLA_postResults(c+j) = (tau_slip_pos(j)+tau_slip_neg(j))/&

code/plastic_dislotwin.f90

@@ -200,6 +200,7 @@ contains
 subroutine plastic_dislotwin_init(fileUnit)
  use, intrinsic :: iso_fortran_env                                                                  ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment)
  use prec, only: &
+   dEq, &
    dNeq
  use debug, only: &
    debug_level,&
@@ -751,8 +752,8 @@ subroutine plastic_dislotwin_init(fileUnit)
       if (plastic_dislotwin_Qsd(instance) <= 0.0_pReal) &
         call IO_error(211_pInt,el=instance,ext_msg='Qsd ('//PLASTICITY_DISLOTWIN_label//')')
       if (sum(plastic_dislotwin_Ntwin(:,instance)) > 0_pInt) then
-        if (abs(plastic_dislotwin_SFE_0K(instance))  <= tiny(0.0_pReal) .and. &
+        if (dEq(plastic_dislotwin_SFE_0K(instance), 0.0_pReal) .and. &
-            abs(plastic_dislotwin_dSFE_dT(instance)) <= tiny(0.0_pReal) .and. &
+            dEq(plastic_dislotwin_dSFE_dT(instance),0.0_pReal) .and. &
                             lattice_structure(phase) == LATTICE_fcc_ID) &
           call IO_error(211_pInt,el=instance,ext_msg='SFE0K ('//PLASTICITY_DISLOTWIN_label//')')
         if (plastic_dislotwin_aTolRho(instance) <= 0.0_pReal) &
@@ -761,8 +762,8 @@ subroutine plastic_dislotwin_init(fileUnit)
           call IO_error(211_pInt,el=instance,ext_msg='aTolTwinFrac ('//PLASTICITY_DISLOTWIN_label//')')
       endif
       if (sum(plastic_dislotwin_Ntrans(:,instance)) > 0_pInt) then
-        if (abs(plastic_dislotwin_SFE_0K(instance))  <= tiny(0.0_pReal) .and. &
+        if (dEq(plastic_dislotwin_SFE_0K(instance), 0.0_pReal) .and. &
-            abs(plastic_dislotwin_dSFE_dT(instance)) <= tiny(0.0_pReal) .and. &
+            dEq(plastic_dislotwin_dSFE_dT(instance),0.0_pReal) .and. &
                             lattice_structure(phase) == LATTICE_fcc_ID) &
           call IO_error(211_pInt,el=instance,ext_msg='SFE0K ('//PLASTICITY_DISLOTWIN_label//')')
         if (plastic_dislotwin_aTolTransFrac(instance) <= 0.0_pReal) &
@@ -1630,7 +1631,8 @@ end subroutine plastic_dislotwin_microstructure
 !--------------------------------------------------------------------------------------------------
 subroutine plastic_dislotwin_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,Temperature,ipc,ip,el)
  use prec, only: &
-   tol_math_check
+   tol_math_check, &
+   dNeq
  use math, only: &
    math_Plain3333to99, &
    math_Mandel6to33, &
@@ -1777,8 +1779,8 @@ subroutine plastic_dislotwin_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,Temperature
  
 !--------------------------------------------------------------------------------------------------
 ! Shear banding (shearband) part
- if(abs(plastic_dislotwin_sbVelocity(instance))   > tiny(0.0_pReal) .and. &
+ if(dNeq(plastic_dislotwin_sbVelocity(instance),  0.0_pReal) .and. &
-    abs(plastic_dislotwin_sbResistance(instance)) > tiny(0.0_pReal)) then
+    dNeq(plastic_dislotwin_sbResistance(instance),0.0_pReal)) then
    gdot_sb = 0.0_pReal
    dgdot_dtausb = 0.0_pReal
    call math_eigenValuesVectorsSym(math_Mandel6to33(Tstar_v),eigValues,eigVectors,error)
@@ -1944,7 +1946,8 @@ end subroutine plastic_dislotwin_LpAndItsTangent
 !--------------------------------------------------------------------------------------------------
 subroutine plastic_dislotwin_dotState(Tstar_v,Temperature,ipc,ip,el)
  use prec, only: &
-   tol_math_check
+   tol_math_check, &
+   dEq
  use math, only: &
    pi
  use material, only: &
@@ -2045,7 +2048,7 @@ subroutine plastic_dislotwin_dotState(Tstar_v,Temperature,ipc,ip,el)
       !* Dipole formation
       EdgeDipMinDistance = &
         plastic_dislotwin_CEdgeDipMinDistance(instance)*plastic_dislotwin_burgersPerSlipSystem(j,instance)
-      if (abs(tau_slip(j)) <= tiny(0.0_pReal)) then
+      if (dEq(tau_slip(j),0.0_pReal)) then
         DotRhoDipFormation = 0.0_pReal
       else
         EdgeDipDistance = &
@@ -2073,10 +2076,10 @@ subroutine plastic_dislotwin_dotState(Tstar_v,Temperature,ipc,ip,el)
         plastic_dislotwin_CAtomicVolume(instance)*plastic_dislotwin_burgersPerSlipSystem(j,instance)**(3.0_pReal)
       VacancyDiffusion = &
         plastic_dislotwin_D0(instance)*exp(-plastic_dislotwin_Qsd(instance)/(kB*Temperature))
-      if (abs(tau_slip(j)) <= tiny(0.0_pReal)) then
+      if (dEq(tau_slip(j),0.0_pReal)) then
         DotRhoEdgeDipClimb = 0.0_pReal
       else
-        if (EdgeDipDistance-EdgeDipMinDistance <= tiny(0.0_pReal)) then
+        if (dEq(EdgeDipDistance-EdgeDipMinDistance,0.0_pReal)) then
           DotRhoEdgeDipClimb = 0.0_pReal
         else
           ClimbVelocity = 3.0_pReal*lattice_mu(ph)*VacancyDiffusion*AtomicVolume/ &
@@ -2191,7 +2194,8 @@ end subroutine plastic_dislotwin_dotState
 !--------------------------------------------------------------------------------------------------
 function plastic_dislotwin_postResults(Tstar_v,Temperature,ipc,ip,el)
  use prec, only: &
-   tol_math_check
+   tol_math_check, &
+   dEq
  use math, only: &
    pi, &
    math_Mandel6to33, &
@@ -2506,11 +2510,8 @@ function plastic_dislotwin_postResults(Tstar_v,Temperature,ipc,ip,el)
              endif
 
              !* Stress exponent
-             if (abs(gdot_slip(j))<=tiny(0.0_pReal)) then
+             plastic_dislotwin_postResults(c+j) = &
-               plastic_dislotwin_postResults(c+j) = 0.0_pReal
+               merge(0.0_pReal,(tau/gdot_slip(j))*dgdot_dtauslip,dEq(gdot_slip(j),0.0_pReal))
-             else
-               plastic_dislotwin_postResults(c+j) = (tau/gdot_slip(j))*dgdot_dtauslip
-             endif
          enddo ; enddo
          c = c + ns
       case (sb_eigenvalues_ID)

code/plastic_isotropic.f90

@@ -524,6 +524,8 @@ end subroutine plastic_isotropic_LiAndItsTangent
 !> @brief calculates the rate of change of microstructure
 !--------------------------------------------------------------------------------------------------
 subroutine plastic_isotropic_dotState(Tstar_v,ipc,ip,el)
+ use prec, only: &
+   dEq
  use math, only: &
    math_mul6x6
  use material, only: &
@@ -570,7 +572,7 @@ subroutine plastic_isotropic_dotState(Tstar_v,ipc,ip,el)
 !--------------------------------------------------------------------------------------------------
 ! hardening coefficient
  if (abs(gamma_dot) > 1e-12_pReal) then
-   if (abs(param(instance)%tausat_SinhFitA) <= tiny(0.0_pReal)) then
+   if (dEq(param(instance)%tausat_SinhFitA,0.0_pReal)) then
      saturation = param(instance)%tausat
    else
      saturation = (  param(instance)%tausat &

code/plastic_nonlocal.f90

@@ -1549,6 +1549,8 @@ end subroutine plastic_nonlocal_aTolState
 !> @brief calculates quantities characterizing the microstructure
 !--------------------------------------------------------------------------------------------------
 subroutine plastic_nonlocal_microstructure(Fe, Fp, ip, el)
+use prec, only: &
+  dEq
 use IO, only: &
   IO_error
 use math, only: &
@@ -1792,7 +1794,7 @@ if (.not. phase_localPlasticity(ph) .and. shortRangeStressCorrection(instance))
                                   - neighbor_rhoExcess(c,s,neighbors(2))
       enddo
       invConnections = math_inv33(connections)
-      if (all(abs(invConnections) <= tiny(0.0_pReal))) &                                            ! check for failed in version (math_inv33 returns 0) and avoid floating point equality comparison
+      if (all(dEq(invConnections,0.0_pReal))) &
         call IO_error(-1_pInt,ext_msg='back stress calculation: inversion error')
       rhoExcessGradient(c) = math_mul3x3(m(1:3,s,c), &
                                          math_mul33x3(invConnections,rhoExcessDifferences))
@@ -2200,6 +2202,8 @@ end subroutine plastic_nonlocal_LpAndItsTangent
 !> @brief (instantaneous) incremental change of microstructure
 !--------------------------------------------------------------------------------------------------
 subroutine plastic_nonlocal_deltaState(Tstar_v,ip,el)
+use prec, only: &
+  dNeq
 use debug,    only: debug_level, &
                     debug_constitutive, &
                     debug_levelBasic, &
@@ -2322,8 +2326,8 @@ dUpper(1:ns,2) = lattice_mu(ph) * burgers(1:ns,instance) / (4.0_pReal * pi * abs
 
 
 forall (c = 1_pInt:2_pInt)
-  where(sqrt(rhoSgl(1:ns,2*c-1)+rhoSgl(1:ns,2*c)+&
+  where(dNeq(sqrt(rhoSgl(1:ns,2*c-1)+rhoSgl(1:ns,2*c)+&
-        abs(rhoSgl(1:ns,2*c+3))+abs(rhoSgl(1:ns,2*c+4))+rhoDip(1:ns,c)) >= tiny(0.0_pReal)) &
+             abs(rhoSgl(1:ns,2*c+3))+abs(rhoSgl(1:ns,2*c+4))+rhoDip(1:ns,c)),0.0_pReal)) &
     dUpper(1:ns,c) = min(1.0_pReal / sqrt(rhoSgl(1:ns,2*c-1) + rhoSgl(1:ns,2*c) & 
                        + abs(rhoSgl(1:ns,2*c+3)) + abs(rhoSgl(1:ns,2*c+4)) + rhoDip(1:ns,c)), &
                        dUpper(1:ns,c))
@@ -2335,7 +2339,7 @@ deltaDUpper = dUpper - dUpperOld
 !*** dissociation by stress increase
 deltaRhoDipole2SingleStress = 0.0_pReal
 forall (c=1_pInt:2_pInt, s=1_pInt:ns, deltaDUpper(s,c) < 0.0_pReal .and. &
-                                        abs(dUpperOld(s,c) - dLower(s,c)) > tiny(0.0_pReal)) &
+                                        dNeq(dUpperOld(s,c) - dLower(s,c),0.0_pReal)) &
   deltaRhoDipole2SingleStress(s,8_pInt+c) = rhoDip(s,c) * deltaDUpper(s,c) &
                                            / (dUpperOld(s,c) - dLower(s,c))
 
@@ -2383,7 +2387,8 @@ subroutine plastic_nonlocal_dotState(Tstar_v, Fe, Fp, Temperature, &
                                      timestep,subfrac, ip,el)
 
 use prec,     only: DAMASK_NaN, &
-                    dNeq
+                    dNeq, &
+                    dEq
 use numerics, only: numerics_integrationMode, &
                     numerics_timeSyncing
 use IO,       only: IO_error
@@ -2617,8 +2622,8 @@ dUpper(1:ns,1) = lattice_mu(ph) * burgers(1:ns,instance) &
 dUpper(1:ns,2) = lattice_mu(ph) * burgers(1:ns,instance) &
                / (4.0_pReal * pi * abs(tau))
 forall (c = 1_pInt:2_pInt)
-  where(sqrt(rhoSgl(1:ns,2*c-1)+rhoSgl(1:ns,2*c)+&
+  where(dNeq(sqrt(rhoSgl(1:ns,2*c-1)+rhoSgl(1:ns,2*c)+&
-        abs(rhoSgl(1:ns,2*c+3))+abs(rhoSgl(1:ns,2*c+4))+rhoDip(1:ns,c)) >= tiny(0.0_pReal)) &
+             abs(rhoSgl(1:ns,2*c+3))+abs(rhoSgl(1:ns,2*c+4))+rhoDip(1:ns,c)),0.0_pReal)) &
     dUpper(1:ns,c) = min(1.0_pReal / sqrt(rhoSgl(1:ns,2*c-1) + rhoSgl(1:ns,2*c) & 
                        + abs(rhoSgl(1:ns,2*c+3)) + abs(rhoSgl(1:ns,2*c+4)) + rhoDip(1:ns,c)), &
                        dUpper(1:ns,c))
@@ -2830,11 +2835,11 @@ if (.not. phase_localPlasticity(material_phase(1_pInt,ip,el))) then
       my_rhoSgl = rhoSgl
       my_v = v
       if(numerics_timeSyncing) then
-        if (abs(subfrac(1_pInt,ip,el))<= tiny(0.0_pReal)) then
+        if (dEq(subfrac(1_pInt,ip,el),0.0_pReal)) then
           my_rhoSgl = rhoSgl0
           my_v = v0
         elseif (neighbor_n > 0_pInt) then
-          if (abs(subfrac(1_pInt,neighbor_ip,neighbor_el))<= tiny(0.0_pReal)) then
+          if (dEq(subfrac(1_pInt,neighbor_ip,neighbor_el),0.0_pReal)) then
             my_rhoSgl = rhoSgl0
             my_v = v0
           endif
@@ -3172,6 +3177,8 @@ end subroutine plastic_nonlocal_updateCompatibility
 !* calculates quantities characterizing the microstructure           *
 !*********************************************************************
 function plastic_nonlocal_dislocationstress(Fe, ip, el)
+use prec, only: &
+  dEq
 use math,     only: math_mul33x33, &
                     math_mul33x3, &
                     math_inv33, &
@@ -3384,7 +3391,7 @@ if (.not. phase_localPlasticity(ph)) then
                     Rsquare = R * R
                     Rcube = Rsquare * R 
                     denominator = R * (R + flipSign * lambda)
-                    if (abs(denominator)<= tiny(0.0_pReal)) exit ipLoop
+                    if (dEq(denominator,0.0_pReal)) exit ipLoop
                       
                     sigma(1,1) = sigma(1,1) - real(side,pReal) &
                                             * flipSign * z / denominator &
@@ -3429,7 +3436,7 @@ if (.not. phase_localPlasticity(ph)) then
                     Rsquare = R * R
                     Rcube = Rsquare * R 
                     denominator = R * (R + flipSign * lambda)
-                    if (abs(denominator)<= tiny(0.0_pReal)) exit ipLoop
+                    if (dEq(denominator,0.0_pReal)) exit ipLoop
                     
                     sigma(1,2) = sigma(1,2) - real(side,pReal) * flipSign * z &
                                                                * (1.0_pReal - lattice_nu(ph)) / denominator &
@@ -3518,6 +3525,8 @@ end function plastic_nonlocal_dislocationstress
 !> @brief return array of constitutive results
 !--------------------------------------------------------------------------------------------------
 function plastic_nonlocal_postResults(Tstar_v,Fe,ip,el)
+ use prec, only: &
+   dNeq
  use math, only: &
    math_mul6x6, &
    math_mul33x3, &
@@ -3634,8 +3643,8 @@ dUpper(1:ns,1) = lattice_mu(ph) * burgers(1:ns,instance) &
 dUpper(1:ns,2) = lattice_mu(ph) * burgers(1:ns,instance) &
                / (4.0_pReal * pi * abs(tau))
 forall (c = 1_pInt:2_pInt)
-  where(sqrt(rhoSgl(1:ns,2*c-1)+rhoSgl(1:ns,2*c)+&
+  where(dNeq(sqrt(rhoSgl(1:ns,2*c-1)+rhoSgl(1:ns,2*c)+&
-        abs(rhoSgl(1:ns,2*c+3))+abs(rhoSgl(1:ns,2*c+4))+rhoDip(1:ns,c)) >= tiny(0.0_pReal)) &
+             abs(rhoSgl(1:ns,2*c+3))+abs(rhoSgl(1:ns,2*c+4))+rhoDip(1:ns,c)),0.0_pReal)) &
     dUpper(1:ns,c) = min(1.0_pReal / sqrt(rhoSgl(1:ns,2*c-1) + rhoSgl(1:ns,2*c) & 
                        + abs(rhoSgl(1:ns,2*c+3)) + abs(rhoSgl(1:ns,2*c+4)) + rhoDip(1:ns,c)), &
                        dUpper(1:ns,c))

code/plastic_phenoplus.f90

@@ -112,6 +112,8 @@ contains
 !--------------------------------------------------------------------------------------------------
 subroutine plastic_phenoplus_init(fileUnit)
  use, intrinsic :: iso_fortran_env                                                                  ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment)
+ use prec, only: &
+   dEq
  use debug, only: &
    debug_level, &
    debug_constitutive,&
@@ -481,7 +483,7 @@ subroutine plastic_phenoplus_init(fileUnit)
      if (any(plastic_phenoplus_tausat_slip(:,instance) <= 0.0_pReal .and. &
              plastic_phenoplus_Nslip(:,instance) > 0)) &
        call IO_error(211_pInt,el=instance,ext_msg='tausat_slip ('//PLASTICITY_PHENOPLUS_label//')')
-     if (any(abs(plastic_phenoplus_a_slip(instance)) <= tiny(0.0_pReal) .and. &
+     if (any(dEq(plastic_phenoplus_a_slip(instance),0.0_pReal) .and. &
              plastic_phenoplus_Nslip(:,instance) > 0)) &
        call IO_error(211_pInt,el=instance,ext_msg='a_slip ('//PLASTICITY_PHENOPLUS_label//')')
      if (any(plastic_phenoplus_tau0_twin(:,instance) < 0.0_pReal .and. &
@@ -923,6 +925,8 @@ end subroutine plastic_phenoplus_microstructure
 !> @brief calculates plastic velocity gradient and its tangent
 !--------------------------------------------------------------------------------------------------
 subroutine plastic_phenoplus_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,ipc,ip,el)
+ use prec, only: &
+   dNeq
  use math, only: &
    math_Plain3333to99, &
    math_Mandel6to33
@@ -1038,7 +1042,7 @@ subroutine plastic_phenoplus_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,ipc,ip,el)
                (gdot_slip_pos+gdot_slip_neg)*lattice_Sslip(1:3,1:3,1,index_myFamily+i,ph)
 
      ! Calculation of the tangent of Lp
-     if (abs(gdot_slip_pos) > tiny(0.0_pReal)) then
+     if (dNeq(gdot_slip_pos,0.0_pReal)) then
        dgdot_dtauslip_pos = gdot_slip_pos*plastic_phenoplus_n_slip(instance)/tau_slip_pos
        forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) &
          dLp_dTstar3333(k,l,m,n) = dLp_dTstar3333(k,l,m,n) + &
@@ -1046,7 +1050,7 @@ subroutine plastic_phenoplus_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,ipc,ip,el)
                                                      nonSchmid_tensor(m,n,1)
      endif
 
-     if (abs(gdot_slip_neg) > tiny(0.0_pReal)) then
+     if (dNeq(gdot_slip_neg,0.0_pReal)) then
        dgdot_dtauslip_neg = gdot_slip_neg*plastic_phenoplus_n_slip(instance)/tau_slip_neg
        forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) &
          dLp_dTstar3333(k,l,m,n) = dLp_dTstar3333(k,l,m,n) + &
@@ -1073,7 +1077,7 @@ subroutine plastic_phenoplus_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,ipc,ip,el)
      Lp = Lp + gdot_twin*lattice_Stwin(1:3,1:3,index_myFamily+i,ph)
 
      ! Calculation of the tangent of Lp
-     if (abs(gdot_twin) > tiny(0.0_pReal)) then
+     if (dNeq(gdot_twin,0.0_pReal)) then
        dgdot_dtautwin = gdot_twin*plastic_phenoplus_n_twin(instance)/tau_twin
        forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) &
          dLp_dTstar3333(k,l,m,n) = dLp_dTstar3333(k,l,m,n) + &

code/plastic_phenopowerlaw.f90

@@ -124,6 +124,8 @@ contains
 !--------------------------------------------------------------------------------------------------
 subroutine plastic_phenopowerlaw_init(fileUnit)
  use, intrinsic :: iso_fortran_env                                                                  ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment)
+ use prec, only: &
+   dEq
  use debug, only: &
    debug_level, &
    debug_constitutive,&
@@ -487,7 +489,7 @@ subroutine plastic_phenopowerlaw_init(fileUnit)
      if (any(plastic_phenopowerlaw_tausat_slip(:,instance) <= 0.0_pReal .and. &
              plastic_phenopowerlaw_Nslip(:,instance) > 0)) &
        call IO_error(211_pInt,el=instance,ext_msg='tausat_slip ('//PLASTICITY_PHENOPOWERLAW_label//')')
-     if (any(abs(plastic_phenopowerlaw_a_slip(instance)) <= tiny(0.0_pReal) .and. &
+     if (any(dEq(plastic_phenopowerlaw_a_slip(instance),0.0_pReal) .and. &
              plastic_phenopowerlaw_Nslip(:,instance) > 0)) &
        call IO_error(211_pInt,el=instance,ext_msg='a_slip ('//PLASTICITY_PHENOPOWERLAW_label//')')
      if (any(plastic_phenopowerlaw_tau0_twin(:,instance) < 0.0_pReal .and. &
@@ -774,6 +776,8 @@ end subroutine plastic_phenopowerlaw_aTolState
 !> @brief calculates plastic velocity gradient and its tangent
 !--------------------------------------------------------------------------------------------------
 subroutine plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,ipc,ip,el)
+ use prec, only: &
+   dNeq
  use math, only: &
    math_Plain3333to99, &
    math_Mandel6to33
@@ -863,7 +867,7 @@ subroutine plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,ipc,ip,
                (gdot_slip_pos+gdot_slip_neg)*lattice_Sslip(1:3,1:3,1,index_myFamily+i,ph)
 
      ! Calculation of the tangent of Lp
-     if (abs(gdot_slip_pos) > tiny(0.0_pReal)) then
+     if (dNeq(gdot_slip_pos,0.0_pReal)) then
        dgdot_dtauslip_pos = gdot_slip_pos*plastic_phenopowerlaw_n_slip(instance)/tau_slip_pos
        forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) &
          dLp_dTstar3333(k,l,m,n) = dLp_dTstar3333(k,l,m,n) + &
@@ -871,7 +875,7 @@ subroutine plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,ipc,ip,
                                                      nonSchmid_tensor(m,n,1)
      endif
 
-     if (abs(gdot_slip_neg) > tiny(0.0_pReal)) then
+     if (dNeq(gdot_slip_neg,0.0_pReal)) then
        dgdot_dtauslip_neg = gdot_slip_neg*plastic_phenopowerlaw_n_slip(instance)/tau_slip_neg
        forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) &
          dLp_dTstar3333(k,l,m,n) = dLp_dTstar3333(k,l,m,n) + &
@@ -898,7 +902,7 @@ subroutine plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,ipc,ip,
      Lp = Lp + gdot_twin*lattice_Stwin(1:3,1:3,index_myFamily+i,ph)
 
      ! Calculation of the tangent of Lp
-     if (abs(gdot_twin) > tiny(0.0_pReal)) then
+     if (dNeq(gdot_twin,0.0_pReal)) then
        dgdot_dtautwin = gdot_twin*plastic_phenopowerlaw_n_twin(instance)/tau_twin
        forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) &
          dLp_dTstar3333(k,l,m,n) = dLp_dTstar3333(k,l,m,n) + &

code/prec.f90

@@ -115,7 +115,9 @@ module prec
    prec_init, &
    prec_isNaN, &
    dEq, &
-   dNeq
+   cEq, &
+   dNeq, &
+   cNeq
  
 contains
 
@@ -180,23 +182,23 @@ end function prec_isNaN
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief equality comparison for double precision
+!> @brief equality comparison for float with double precision
 ! replaces "==" but for certain (relative) tolerance. Counterpart to dNeq
 ! http://www.cygnus-software.com/papers/comparingfloats/comparingfloats.htm
 !--------------------------------------------------------------------------------------------------
 logical elemental pure function dEq(a,b,tol)
 
  implicit none
- real(pReal), intent(in) :: a,b
+ real(pReal), intent(in)           :: a,b
  real(pReal), intent(in), optional :: tol
- real(pReal), parameter  :: eps = 2.220446049250313E-16                                             ! DBL_EPSILON in C
+ real(pReal), parameter            :: eps = 2.220446049250313E-16                                   ! DBL_EPSILON in C
 
  dEq = merge(.True., .False.,abs(a-b) <= merge(tol,eps,present(tol))*maxval(abs([a,b])))
 end function dEq
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief inequality comparison for double precision
+!> @brief inequality comparison for float with double precision
 ! replaces "!=" but for certain (relative) tolerance. Counterpart to dEq
 ! http://www.cygnus-software.com/papers/comparingfloats/comparingfloats.htm
 !--------------------------------------------------------------------------------------------------
@@ -205,9 +207,44 @@ logical elemental pure function dNeq(a,b,tol)
  implicit none
  real(pReal), intent(in)           :: a,b
  real(pReal), intent(in), optional :: tol
- real(pReal), parameter  :: eps = 2.220446049250313E-16                                             ! DBL_EPSILON in C
+ real(pReal), parameter            :: eps = 2.220446049250313E-16                                   ! DBL_EPSILON in C
 
  dNeq = merge(.False., .True.,abs(a-b) <= merge(tol,eps,present(tol))*maxval(abs([a,b])))
 end function dNeq
 
+!--------------------------------------------------------------------------------------------------
+!> @brief equality comparison for complex with double precision
+! replaces "==" but for certain (relative) tolerance. Counterpart to cNeq
+! http://www.cygnus-software.com/papers/comparingfloats/comparingfloats.htm
+! probably a component wise comparison would be more accurate than the comparsion of the absolute
+! value
+!--------------------------------------------------------------------------------------------------
+logical elemental pure function cEq(a,b,tol)
+
+ implicit none
+ complex(pReal), intent(in)           :: a,b
+ real(pReal),    intent(in), optional :: tol
+ real(pReal),    parameter            :: eps = 2.220446049250313E-16                                ! DBL_EPSILON in C
+
+ cEq = merge(.True., .False.,abs(a-b) <= merge(tol,eps,present(tol))*maxval(abs([a,b])))
+end function cEq
+
+
+!--------------------------------------------------------------------------------------------------
+!> @brief inequality comparison for complex with double precision
+! replaces "!=" but for certain (relative) tolerance. Counterpart to cEq
+! http://www.cygnus-software.com/papers/comparingfloats/comparingfloats.htm
+! probably a component wise comparison would be more accurate than the comparsion of the absolute
+! value
+!--------------------------------------------------------------------------------------------------
+logical elemental pure function cNeq(a,b,tol)
+
+ implicit none
+ complex(pReal), intent(in)           :: a,b
+ real(pReal),    intent(in), optional :: tol
+ real(pReal),    parameter            :: eps = 2.220446049250313E-16                                ! DBL_EPSILON in C
+
+ cNeq = merge(.False., .True.,abs(a-b) <= merge(tol,eps,present(tol))*maxval(abs([a,b])))
+end function cNeq
+
 end module prec

code/spectral_utilities.f90

@@ -402,7 +402,7 @@ subroutine utilities_updateGamma(C,saveReference)
  integer(pInt) :: &
    i, j, k, &
    l, m, n, o
- logical :: ierr
+ logical :: err
 
  C_ref = C
  if (saveReference) then
@@ -426,7 +426,7 @@ subroutine utilities_updateGamma(C,saveReference)
        matA(1:3,1:3) = real(temp33_complex); matA(4:6,4:6) = real(temp33_complex)
        matA(1:3,4:6) = aimag(temp33_complex); matA(4:6,1:3) = -aimag(temp33_complex)
        if (abs(math_det33(matA(1:3,1:3))) > 1e-16) then
-         call math_invert(6_pInt, matA, matInvA, ierr)
+         call math_invert(6_pInt, matA, matInvA, err)
          temp33_complex = cmplx(matInvA(1:3,1:3),matInvA(1:3,4:6),pReal)
          forall(l=1_pInt:3_pInt, m=1_pInt:3_pInt, n=1_pInt:3_pInt, o=1_pInt:3_pInt) &
            gamma_hat(l,m,n,o,i,j,k-grid3Offset) =  temp33_complex(l,n)* &
@@ -938,6 +938,10 @@ end subroutine utilities_fourierTensorDivergence
 !--------------------------------------------------------------------------------------------------
 subroutine utilities_constitutiveResponse(F_lastInc,F,timeinc, &
                                           P,C_volAvg,C_minmaxAvg,P_av,forwardData,rotation_BC)
+ use prec, only: &
+   dNeq
+ use IO, only: &
+   IO_error
  use debug, only: &
    debug_reset, &
    debug_info
@@ -976,10 +980,9 @@ subroutine utilities_constitutiveResponse(F_lastInc,F,timeinc, &
    age
 
  integer(pInt) :: &
-   j,k
+   j,k,ierr
  real(pReal), dimension(3,3,3,3) :: max_dPdF, min_dPdF
  real(pReal)   :: max_dPdF_norm, min_dPdF_norm, defgradDetMin, defgradDetMax, defgradDet
- PetscErrorCode :: ierr
 
  external :: &
    MPI_Reduce, &
@@ -1013,7 +1016,9 @@ subroutine utilities_constitutiveResponse(F_lastInc,F,timeinc, &
      defgradDetMin = min(defgradDetMin,defgradDet)
    end do
    call MPI_reduce(MPI_IN_PLACE,defgradDetMax,1,MPI_DOUBLE,MPI_MAX,0,PETSC_COMM_WORLD,ierr)
+   if(ierr /=0_pInt) call IO_error(894_pInt, ext_msg='MPI_Allreduce max')
    call MPI_reduce(MPI_IN_PLACE,defgradDetMin,1,MPI_DOUBLE,MPI_MIN,0,PETSC_COMM_WORLD,ierr)
+   if(ierr /=0_pInt) call IO_error(894_pInt, ext_msg='MPI_Allreduce min')
    if (worldrank == 0_pInt) then
      write(6,'(a,1x,es11.4)') ' max determinant of deformation =', defgradDetMax
      write(6,'(a,1x,es11.4)') ' min determinant of deformation =', defgradDetMin
@@ -1039,7 +1044,9 @@ subroutine utilities_constitutiveResponse(F_lastInc,F,timeinc, &
  end do
 
  call MPI_Allreduce(MPI_IN_PLACE,max_dPdF,81,MPI_DOUBLE,MPI_MAX,PETSC_COMM_WORLD,ierr)
+ if(ierr /=0_pInt) call IO_error(894_pInt, ext_msg='MPI_Allreduce max')
  call MPI_Allreduce(MPI_IN_PLACE,min_dPdF,81,MPI_DOUBLE,MPI_MIN,PETSC_COMM_WORLD,ierr)
+ if(ierr /=0_pInt) call IO_error(894_pInt, ext_msg='MPI_Allreduce min')
 
  C_minmaxAvg = 0.5_pReal*(max_dPdF + min_dPdF)
  C_volAvg = sum(sum(materialpoint_dPdF,dim=6),dim=5) * wgt
@@ -1194,6 +1201,10 @@ end function utilities_getFreqDerivative
 ! convolution
 !--------------------------------------------------------------------------------------------------
 subroutine utilities_updateIPcoords(F)
+ use prec, only: &
+   cNeq
+ use IO, only: &
+   IO_error
  use math, only: &
    math_mul33x3
  use mesh, only: &
@@ -1205,10 +1216,9 @@ subroutine utilities_updateIPcoords(F)
  implicit none
 
  real(pReal),   dimension(3,3,grid(1),grid(2),grid3), intent(in) :: F
- integer(pInt) :: i, j, k, m
+ integer(pInt) :: i, j, k, m, ierr
  real(pReal),   dimension(3) :: step, offset_coords
  real(pReal),   dimension(3,3) :: Favg
- PetscErrorCode :: ierr
  external &
    MPI_Bcast
 
@@ -1219,8 +1229,8 @@ subroutine utilities_updateIPcoords(F)
  call utilities_FFTtensorForward()
  call utilities_fourierTensorDivergence()
 
- do k = 1_pInt, grid3; do j = 1_pInt, grid(2) ;do i = 1_pInt, grid1Red
+ do k = 1_pInt, grid3; do j = 1_pInt, grid(2); do i = 1_pInt, grid1Red
-   if (any(abs(xi1st(1:3,i,j,k)) > tiny(0.0_pReal))) &
+   if (any(cNeq(xi1st(1:3,i,j,k),cmplx(0.0_pReal,0.0_pReal)))) &
      vectorField_fourier(1:3,i,j,k) = vectorField_fourier(1:3,i,j,k)/ &
                                       sum(conjg(-xi1st(1:3,i,j,k))*xi1st(1:3,i,j,k))
  enddo; enddo; enddo
@@ -1230,12 +1240,14 @@ subroutine utilities_updateIPcoords(F)
 ! average F
  if (grid3Offset == 0_pInt) Favg = real(tensorField_fourier(1:3,1:3,1,1,1),pReal)*wgt
  call MPI_Bcast(Favg,9,MPI_DOUBLE,0,PETSC_COMM_WORLD,ierr)
+ if(ierr /=0_pInt) call IO_error(894_pInt, ext_msg='update_IPcoords')
 
 !--------------------------------------------------------------------------------------------------
 ! add average to fluctuation and put (0,0,0) on (0,0,0)
  step = geomSize/real(grid, pReal)
  if (grid3Offset == 0_pInt) offset_coords = vectorField_real(1:3,1,1,1)
  call MPI_Bcast(offset_coords,3,MPI_DOUBLE,0,PETSC_COMM_WORLD,ierr)
+ if(ierr /=0_pInt) call IO_error(894_pInt, ext_msg='update_IPcoords')
  offset_coords = math_mul33x3(Favg,step/2.0_pReal) - offset_coords
  m = 1_pInt
  do k = 1_pInt,grid3; do j = 1_pInt,grid(2); do i = 1_pInt,grid(1)