Wrong integer numbers (#887)

Co-authored-by: Xuzheng Tian <xtian64@login-phoenix-rh9-3.pace.gatech.edu>
Co-authored-by: Spencer Bryngelson <sbryngelson@gmail.com>

Author: 3 people

src/common/m_checker_common.fpp

@@ -322,7 +322,7 @@ contains
         @:PROHIBIT(surface_tension .and. sigma < 0._wp, &
             "sigma must be greater than or equal to zero")
 
-        @:PROHIBIT(surface_tension .and. sigma == dflt_real, &
+        @:PROHIBIT(surface_tension .and. f_approx_equal(sigma, dflt_real), &
             "sigma must be set if surface_tension is enabled")
 
         @:PROHIBIT(.not. f_is_default(sigma) .and. .not. surface_tension, &
@@ -347,9 +347,12 @@ contains
         !! Called by s_check_inputs_common for all three stages
     impure subroutine s_check_inputs_moving_bc
         #:for X, VB2, VB3 in [('x', 'vb2', 'vb3'), ('y', 'vb3', 'vb1'), ('z', 'vb1', 'vb2')]
-            if (any((/bc_${X}$%vb1, bc_${X}$%vb2, bc_${X}$%vb3/) /= 0._wp)) then
+            if (.not. (f_approx_equal(bc_${X}$%vb1, 0._wp) .and. &
+                       f_approx_equal(bc_${X}$%vb2, 0._wp) .and. &
+                       f_approx_equal(bc_${X}$%vb3, 0._wp))) then
                 if (bc_${X}$%beg == BC_SLIP_WALL) then
-                    if (any((/bc_${X}$%${VB2}$, bc_${X}$%${VB3}$/) /= 0._wp)) then
+                    if (.not. (f_approx_equal(bc_${X}$%${VB2}$, 0._wp) .and. &
+                               f_approx_equal(bc_${X}$%${VB3}$, 0._wp))) then
                         call s_mpi_abort("bc_${X}$%beg must be -15 if "// &
                                          "bc_${X}$%${VB2}$ or bc_${X}$%${VB3}$ "// &
                                          "is set. Exiting.", CASE_FILE_ERROR_CODE)
@@ -362,9 +365,12 @@ contains
         #:endfor
 
         #:for X, VE2, VE3 in [('x', 've2', 've3'), ('y', 've3', 've1'), ('z', 've1', 've2')]
-            if (any((/bc_${X}$%ve1, bc_${X}$%ve2, bc_${X}$%ve3/) /= 0._wp)) then
+            if (.not. (f_approx_equal(bc_${X}$%ve1, 0._wp) .and. &
+                       f_approx_equal(bc_${X}$%ve2, 0._wp) .and. &
+                       f_approx_equal(bc_${X}$%ve3, 0._wp))) then
                 if (bc_${X}$%end == BC_SLIP_WALL) then
-                    if (any((/bc_${X}$%${VE2}$, bc_${X}$%${VE3}$/) /= 0._wp)) then
+                    if (.not. (f_approx_equal(bc_${X}$%${VE2}$, 0._wp) .and. &
+                               f_approx_equal(bc_${X}$%${VE3}$, 0._wp))) then
                         call s_mpi_abort("bc_${X}$%end must be -15 if "// &
                                          "bc_${X}$%${VE2}$ or bc_${X}$%${VE3}$ "// &
                                          "is set. Exiting.", CASE_FILE_ERROR_CODE)

src/common/m_eigen_solver.f90

@@ -10,6 +10,8 @@ module m_eigen_solver
 
     use m_precision_select
 
+    use m_helper_basic         !< Functions to compare floating point numbers
+
     implicit none
 
     private; 
@@ -124,7 +126,7 @@ pure subroutine cbal(nm, nl, ar, ai, low, igh, scale)
 
             do 110 i = 1, l
                 if (i == j) go to 110
-                if (ar(j, i) /= 0.0_wp .or. ai(j, i) /= 0.0_wp) go to 120
+                if (.not. f_approx_equal(ar(j, i), 0.0_wp) .or. .not. f_approx_equal(ai(j, i), 0.0_wp)) go to 120
 110         end do
 
             ml = l
@@ -140,7 +142,7 @@ pure subroutine cbal(nm, nl, ar, ai, low, igh, scale)
 
             do 150 i = k, l
                 if (i == j) go to 150
-                if (ar(i, j) /= 0.0_wp .or. ai(i, j) /= 0.0_wp) go to 170
+                if (.not. f_approx_equal(ar(i, j), 0.0_wp) .or. .not. f_approx_equal(ai(i, j), 0.0_wp)) go to 170
 150         end do
 
             ml = k
@@ -164,7 +166,7 @@ pure subroutine cbal(nm, nl, ar, ai, low, igh, scale)
                 r = r + abs(ar(i, j)) + abs(ai(i, j))
 200         end do
 !     guard against zero c or r due to underflow
-            if (c == 0.0_wp .or. r == 0.0_wp) go to 270
+            if (f_approx_equal(c, 0.0_wp) .or. f_approx_equal(r, 0.0_wp)) go to 270
             g = r/radix
             f = 1.0_wp
             s = c + r
@@ -242,7 +244,7 @@ pure subroutine corth(nm, nl, low, igh, ar, ai, ortr, orti)
             do 90 i = ml, igh
                 scale = scale + abs(ar(i, ml - 1)) + abs(ai(i, ml - 1))
 90          end do
-            if (scale == 0._wp) go to 180
+            if (f_approx_equal(scale, 0._wp)) go to 180
             mp = ml + igh
 !     for i=igh step -1 until ml do
             do 100 ii = ml, igh
@@ -254,7 +256,7 @@ pure subroutine corth(nm, nl, low, igh, ar, ai, ortr, orti)
 
             g = sqrt(h)
             call pythag(ortr(ml), orti(ml), f)
-            if (f == 0._wp) go to 103
+            if (f_approx_equal(f, 0._wp)) go to 103
             h = h + f*g
             g = g/f
             ortr(ml) = (1.0_wp + g)*ortr(ml)
@@ -374,8 +376,8 @@ pure subroutine comqr2(nm, nl, low, igh, ortr, orti, hr, hi, wr, wi, zr, zi, ier
 !     for i=igh-1 step -1 until low+1 do
 105     do 140 ii = 1, iend
             i = igh - ii
-            if (abs(ortr(i)) == 0._wp .and. abs(orti(i)) == 0._wp) go to 140
-            if (abs(hr(i, i - 1)) == 0._wp .and. abs(hi(i, i - 1)) == 0._wp) go to 140
+            if (f_approx_equal(abs(ortr(i)), 0._wp) .and. f_approx_equal(abs(orti(i)), 0._wp)) go to 140
+            if (f_approx_equal(abs(hr(i, i - 1)), 0._wp) .and. f_approx_equal(abs(hi(i, i - 1)), 0._wp)) go to 140
 !     norm below is negative of h formed in corth
             norm = hr(i, i - 1)*ortr(i) + hi(i, i - 1)*orti(i)
             ip1 = i + 1
@@ -410,7 +412,7 @@ pure subroutine comqr2(nm, nl, low, igh, ortr, orti, hr, hi, wr, wi, zr, zi, ier
 
         do 170 i = l, igh
             ll = min0(i + 1, igh)
-            if (abs(hi(i, i - 1)) == 0._wp) go to 170
+            if (f_approx_equal(abs(hi(i, i - 1)), 0._wp)) go to 170
             call pythag(hr(i, i - 1), hi(i, i - 1), norm)
             yr = hr(i, i - 1)/norm
             yi = hi(i, i - 1)/norm
@@ -458,7 +460,7 @@ pure subroutine comqr2(nm, nl, low, igh, ortr, orti, hr, hi, wr, wi, zr, zi, ier
             tst1 = abs(hr(l - 1, l - 1)) + abs(hi(l - 1, l - 1)) &
                    + abs(hr(l, l)) + abs(hi(l, l))
             tst2 = tst1 + abs(hr(l, l - 1))
-            if (tst2 == tst1) go to 300
+            if (f_approx_equal(tst2, tst1)) go to 300
 260     end do
 !     form shift
 300     if (l == en) go to 660
@@ -468,7 +470,7 @@ pure subroutine comqr2(nm, nl, low, igh, ortr, orti, hr, hi, wr, wi, zr, zi, ier
         si = hi(en, en)
         xr = hr(enm1, en)*hr(en, enm1)
         xi = hi(enm1, en)*hr(en, enm1)
-        if (xr == 0.0_wp .and. xi == 0.0_wp) go to 340
+        if (f_approx_equal(xr, 0.0_wp) .and. f_approx_equal(xi, 0.0_wp)) go to 340
         yr = (hr(enm1, enm1) - sr)/2.0_wp
         yi = (hi(enm1, enm1) - si)/2.0_wp
         call csroot(yr**2 - yi**2 + xr, 2.0_wp*yr*yi + xi, zzr, zzi)
@@ -522,7 +524,7 @@ pure subroutine comqr2(nm, nl, low, igh, ortr, orti, hr, hi, wr, wi, zr, zi, ier
 500     end do
 
         si = hi(en, en)
-        if (abs(si) == 0._wp) go to 540
+        if (f_approx_equal(abs(si), 0._wp)) go to 540
         call pythag(hr(en, en), si, norm)
         sr = hr(en, en)/norm
         si = si/norm
@@ -567,7 +569,7 @@ pure subroutine comqr2(nm, nl, low, igh, ortr, orti, hr, hi, wr, wi, zr, zi, ier
 590         end do
 600     end do
 
-        if (abs(si) == 0._wp) go to 240
+        if (f_approx_equal(abs(si), 0._wp)) go to 240
 
         do 630 i = 1, en
             yr = hr(i, en)
@@ -602,7 +604,7 @@ pure subroutine comqr2(nm, nl, low, igh, ortr, orti, hr, hi, wr, wi, zr, zi, ier
             end do
         end do
 
-        if (nl == 1 .or. norm == 0._wp) go to 1001
+        if (nl == 1 .or. f_approx_equal(norm, 0._wp)) go to 1001
 !     for en=nl step -1 until 2 do
         do 800 nn = 2, nl
             en = nl + 2 - nn
@@ -625,7 +627,7 @@ pure subroutine comqr2(nm, nl, low, igh, ortr, orti, hr, hi, wr, wi, zr, zi, ier
 
                 yr = xr - wr(i)
                 yi = xi - wi(i)
-                if (yr /= 0.0_wp .or. yi /= 0.0_wp) go to 765
+                if (.not. f_approx_equal(yr, 0.0_wp) .or. .not. f_approx_equal(yi, 0.0_wp)) go to 765
                 tst1 = norm
                 yr = tst1
 760             yr = 0.01_wp*yr
@@ -635,7 +637,7 @@ pure subroutine comqr2(nm, nl, low, igh, ortr, orti, hr, hi, wr, wi, zr, zi, ier
                 call cdiv(zzr, zzi, yr, yi, hr(i, en), hi(i, en))
 !     overflow control
                 tr = abs(hr(i, en)) + abs(hi(i, en))
-                if (tr == 0.0_wp) go to 780
+                if (f_approx_equal(tr, 0.0_wp)) go to 780
                 tst1 = tr
                 tst2 = tst1 + 1.0_wp/tst1
                 if (tst2 > tst1) go to 780
@@ -796,11 +798,11 @@ pure elemental subroutine pythag(a, b, c)
 
         real(wp) :: p, r, s, t, u
         p = max(abs(a), abs(b))
-        if (p == 0.0_wp) go to 20
+        if (f_approx_equal(p, 0.0_wp)) go to 20
         r = (min(abs(a), abs(b))/p)**2
 10      continue
         t = 4.0_wp + r
-        if (t == 4.0_wp) go to 20
+        if (f_approx_equal(t, 4.0_wp)) go to 20
         s = r/t
         u = 1.0_wp + 2.0_wp*s
         p = u*p

src/common/m_helper.fpp

@@ -163,7 +163,7 @@ contains
         phi_nv = (1._wp + sqrt(mu_n/mu_v)*(M_v/M_n)**(0.25_wp))**2 &
                  /(sqrt(8._wp)*sqrt(1._wp + M_n/M_v))
         ! internal bubble pressure
-        pb0 = pl0 + 2._wp*ss/(R0ref*R0)
+        pb0(:) = pl0 + 2._wp*ss/(R0ref*R0(:))
 
         ! mass fraction of vapor
         chi_vw0 = 1._wp/(1._wp + R_v/R_n*(pb0/pv - 1._wp))
@@ -179,10 +179,10 @@ contains
         rho_m0 = pv/(chi_vw0*R_v*temp)
 
         ! mass of gas/vapor computed using dimensional quantities
-        mass_n0 = 4._wp*(pb0 - pv)*pi/(3._wp*R_n*temp*rhol0)*R0**3
-        mass_v0 = 4._wp*pv*pi/(3._wp*R_v*temp*rhol0)*R0**3
+        mass_n0(:) = 4._wp*(pb0(:) - pv)*pi/(3._wp*R_n*temp*rhol0)*R0(:)**3
+        mass_v0(:) = 4._wp*pv*pi/(3._wp*R_v*temp*rhol0)*R0(:)**3
         ! Peclet numbers
-        Pe_T = rho_m0*cp_m0*uu*R0ref/k_m0
+        Pe_T(:) = rho_m0*cp_m0(:)*uu*R0ref/k_m0(:)
         Pe_c = uu*R0ref/D_m
 
         Tw = temp
@@ -191,8 +191,8 @@ contains
         !if(.not. qbmm) then
         R_n = rhol0*R_n*temp/pl0
         R_v = rhol0*R_v*temp/pl0
-        k_n = k_n/k_m0
-        k_v = k_v/k_m0
+        k_n(:) = k_n(:)/k_m0(:)
+        k_v(:) = k_v(:)/k_m0(:)
         pb0 = pb0/pl0
         pv = pv/pl0
         Tw = 1._wp
@@ -203,7 +203,7 @@ contains
         !end if
 
         ! natural frequencies
-        omegaN = sqrt(3._wp*k_poly*Ca + 2._wp*(3._wp*k_poly - 1._wp)/(Web*R0))/R0
+        omegaN(:) = sqrt(3._wp*k_poly*Ca + 2._wp*(3._wp*k_poly - 1._wp)/(Web*R0))/R0
         do ir = 1, Nb
             call s_transcoeff(omegaN(ir)*R0(ir), Pe_T(ir)*R0(ir), &
                               Re_trans_T(ir), Im_trans_T(ir))
@@ -488,7 +488,7 @@ contains
         real(wp) :: P
 
         if (l == 0) then
-            P = 1_wp
+            P = 1._wp
         else if (l == 1) then
             P = x
         else

src/common/m_helper_basic.f90

@@ -10,6 +10,7 @@ module m_helper_basic
 
     private; 
     public :: f_approx_equal, &
+              f_approx_in_array, &
               f_is_default, &
               f_all_default, &
               f_is_integer, &
@@ -20,7 +21,7 @@ module m_helper_basic
     !> This procedure checks if two floating point numbers of wp are within tolerance.
     !! @param a First number.
     !! @param b Second number.
-    !! @param tol_input Relative error (default = 1e-6_wp).
+    !! @param tol_input Relative error (default = 1e-10_wp).
     !! @return Result of the comparison.
     logical pure elemental function f_approx_equal(a, b, tol_input) result(res)
         !$acc routine seq
@@ -31,7 +32,7 @@ logical pure elemental function f_approx_equal(a, b, tol_input) result(res)
         if (present(tol_input)) then
             tol = tol_input
         else
-            tol = 1e-6_wp
+            tol = 1e-10_wp
         end if
 
         if (a == b) then
@@ -43,6 +44,35 @@ logical pure elemental function f_approx_equal(a, b, tol_input) result(res)
         end if
     end function f_approx_equal
 
+    !> This procedure checks if the point numbers of wp belongs to another array are within tolerance.
+    !! @param a First number.
+    !! @param b Array that contains several point numbers.
+    !! @param tol_input Relative error (default = 1e-10_wp).
+    !! @return Result of the comparison.
+    logical pure function f_approx_in_array(a, b, tol_input) result(res)
+        !$acc routine seq
+        real(wp), intent(in) :: a
+        real(wp), intent(in) :: b(:)
+        real(wp), optional, intent(in) :: tol_input
+        real(wp) :: tol
+        integer :: i
+
+        res = .false.
+
+        if (present(tol_input)) then
+            tol = tol_input
+        else
+            tol = 1e-10_wp
+        end if
+
+        do i = 1, size(b)
+            if (f_approx_equal(a, b(i), tol)) then
+                res = .true.
+                exit
+            end if
+        end do
+    end function f_approx_in_array
+
     !> Checks if a real(wp) variable is of default value.
     !! @param var Variable to check.
     logical pure elemental function f_is_default(var) result(res)

src/common/m_mpi_common.fpp

@@ -608,8 +608,6 @@ contains
 
         integer :: pack_offset, unpack_offset
 
-        real(wp), pointer :: p_send, p_recv
-
 #ifdef MFC_MPI
 
         call nvtxStartRange("RHS-COMM-PACKBUF")

src/common/m_phase_change.fpp

@@ -17,6 +17,8 @@ module m_phase_change
 
     use ieee_arithmetic
 
+    use m_helper_basic         !< Functions to compare floating point numbers
+
     implicit none
 
     private; 
@@ -748,7 +750,7 @@ contains
         ! Generic loop iterators
         integer :: ns
 
-        if ((pSat == 0.0_wp) .and. (TSIn == 0.0_wp)) then
+        if ((f_approx_equal(pSat, 0.0_wp)) .and. (f_approx_equal(TSIn, 0.0_wp))) then
 
             ! assigning Saturation temperature
             TSat = 0.0_wp