Robust floating point comparisons

benchmarks/viscous_weno5_sgb_mono/case.py

@@ -28,7 +28,7 @@
 procs  = DICT["nodes"] * DICT["tasks_per_node"]
 ncells = math.floor(ppg * procs * ARGS["gbpp"])
 s      = math.floor((ncells / 2.0) ** (1/3))
-Nx, Ny, Nz = 2*s, s, s
+Nx, Ny, Nz = s, s, 2*s
 
 x0      = 10.E-04
 y0      = 10.E-04
@@ -86,8 +86,8 @@
 cfl     = 0.01
 Ldomain = 20.E-03
 L       = Ldomain/x0
-dx      = L/float(Nx)
-dt      = cfl*dx*c0/cact
+dz      = L/float(Nz)
+dt      = cfl*dz*c0/cact
 Lpulse  = 0.3*Ldomain
 Tpulse  = Lpulse/cact
 
@@ -98,12 +98,12 @@
     # ==========================================================
 
     # Computational Domain Parameters ==========================
-    'x_domain%beg'                 : -10.E-03/x0,
-    'x_domain%end'                 :  10.E-03/x0,
+    'x_domain%beg'                 : -5.E-03/x0,
+    'x_domain%end'                 :  5.E-03/x0,
     'y_domain%beg'                 : -5.E-03/y0,
-    'y_domain%end'                 : 5.E-03/y0,
-    'z_domain%beg'                 : -5.E-03/z0,
-    'z_domain%end'                 : 5.E-03/z0,
+    'y_domain%end'                 :  5.E-03/y0,
+    'z_domain%beg'                 : -10.E-03/z0,
+    'z_domain%end'                 :  10.E-03/z0,
     'stretch_x'                    : 'F',
     'cyl_coord'                    : 'F',
     'm'                            : Nx,
@@ -154,9 +154,9 @@
     'patch_icpp(1)%x_centroid'     : 0.,
     'patch_icpp(1)%y_centroid'     : 0.,
     'patch_icpp(1)%z_centroid'     : 0.,
-    'patch_icpp(1)%length_x'       : 20.E-03/x0,
+    'patch_icpp(1)%length_x'       : 10.E-03/x0,
     'patch_icpp(1)%length_y'       : 10.E-03/y0,
-    'patch_icpp(1)%length_z'       : 10.E-03/z0,
+    'patch_icpp(1)%length_z'       : 20.E-03/z0,
     'patch_icpp(1)%vel(1)'         : 0.0,
     'patch_icpp(1)%vel(2)'         : 0.0,
     'patch_icpp(1)%vel(3)'         : 0.0,
@@ -172,9 +172,9 @@
     'patch_icpp(2)%x_centroid'     : 0.,
     'patch_icpp(2)%y_centroid'     : 0.,
     'patch_icpp(2)%z_centroid'     : 0.,
-    'patch_icpp(2)%length_x'       : 5.E-03/x0,
+    'patch_icpp(2)%length_x'       : 10.E-03/x0,
     'patch_icpp(2)%length_y'       : 10.E-03/y0,
-    'patch_icpp(2)%length_z'       : 10.E-03/z0,
+    'patch_icpp(2)%length_z'       : 5.E-03/z0,
     'patch_icpp(2)%alter_patch(1)' : 'T',
     'patch_icpp(2)%vel(1)'         : 0.0,
     'patch_icpp(2)%vel(2)'         : 0.0,
@@ -230,7 +230,8 @@
     # Acoustic source ==========================================
     'Monopole'                     : 'T',
     'num_mono'                     : 1,
-    'Mono(1)%loc(1)'               : -5.E-03/x0,
+    'Mono(1)%support'              : 4,
+    'Mono(1)%loc(3)'               : -5.E-03/x0,
     'Mono(1)%npulse'               : 1,
     'Mono(1)%dir'                  : 1.,
     'Mono(1)%pulse'                : 1,

src/common/m_checker_common.fpp

@@ -10,6 +10,8 @@ module m_checker_common
 
     use m_mpi_proxy            !< Message passing interface (MPI) module proxy
 
+    use m_helper_basic         !< Functions to compare floating point numbers
+
     use m_helper
 
     implicit none

src/common/m_helper.fpp

@@ -36,10 +36,7 @@ module m_helper
               f_create_bbox, &
               s_print_2D_array, &
               f_xor, &
-              f_logical_to_int, &
-              f_approx_equal, &
-              f_is_default, &
-              f_all_default
+              f_logical_to_int
 
 contains
 
@@ -535,56 +532,4 @@ contains
         end if
     end function f_logical_to_int
 
-    !> This procedure checks if two floating point numbers of kind(0d0) are within tolerance.
-    !! @param a First number.
-    !! @param b Second number.
-    !! @param tol_input Relative error (default = 1d-6).
-    !! @return Result of the comparison.
-    logical function f_approx_equal(a, b, tol_input) result(res)
-        ! Reference: https://floating-point-gui.de/errors/comparison/
-
-        real(kind(0d0)), intent(in) :: a, b
-        real(kind(0d0)), optional, intent(in) :: tol_input
-        real(kind(0d0)) :: tol
-
-        if (present(tol_input)) then
-            if (tol_input <= 0d0) then
-                call s_mpi_abort('tol_input must be positive. Exiting ...')
-            end if
-            tol = tol_input
-        else
-            tol = 1d-6
-        end if
-
-        if (a == b) then
-            res = .true.
-        else if (a == 0d0 .or. b == 0d0 .or. (abs(a) + abs(b) < tiny(a))) then
-            res = (abs(a - b) < (tol*tiny(a)))
-        else
-            res = (abs(a - b)/min(abs(a) + abs(b), huge(a)) < tol)
-        end if
-    end function f_approx_equal
-
-    !> Checks if a real(kind(0d0)) variable is of default value.
-    !! @param var Variable to check.
-    logical function f_is_default(var) result(res)
-        real(kind(0d0)), intent(in) :: var
-
-        res = f_approx_equal(var, dflt_real)
-    end function f_is_default
-
-    !> Checks if ALL elements of a real(kind(0d0)) array are of default value.
-    !! @param var_array Array to check.
-    logical function f_all_default(var_array) result(res)
-        real(kind(0d0)), intent(in) :: var_array(:)
-        logical :: res_array(size(var_array))
-        integer :: i
-
-        do i = 1, size(var_array)
-            res_array(i) = f_is_default(var_array(i))
-        end do
-
-        res = all(res_array)
-    end function f_all_default
-
 end module m_helper

src/common/m_helper_basic.f90

@@ -0,0 +1,71 @@
+!>
+!! @file m_helper_basic.f90
+!! @brief Contains module m_helper_basic
+
+module m_helper_basic
+
+    ! Dependencies =============================================================
+
+    use m_derived_types        !< Definitions of the derived types
+    ! ==========================================================================
+
+    implicit none
+
+    private; 
+    public :: f_approx_equal, &
+              f_is_default, &
+              f_all_default
+
+contains
+
+    !> This procedure checks if two floating point numbers of kind(0d0) are within tolerance.
+    !! @param a First number.
+    !! @param b Second number.
+    !! @param tol_input Relative error (default = 1d-6).
+    !! @return Result of the comparison.
+    logical function f_approx_equal(a, b, tol_input) result(res)
+        ! Reference: https://floating-point-gui.de/errors/comparison/
+
+        real(kind(0d0)), intent(in) :: a, b
+        real(kind(0d0)), optional, intent(in) :: tol_input
+        real(kind(0d0)) :: tol
+
+        if (present(tol_input)) then
+            tol = tol_input
+        else
+            tol = 1d-6
+        end if
+
+        if (a == b) then
+            res = .true.
+        else if (a == 0d0 .or. b == 0d0 .or. (abs(a) + abs(b) < tiny(a))) then
+            res = (abs(a - b) < (tol*tiny(a)))
+        else
+            res = (abs(a - b)/min(abs(a) + abs(b), huge(a)) < tol)
+        end if
+    end function f_approx_equal
+
+    !> Checks if a real(kind(0d0)) variable is of default value.
+    !! @param var Variable to check.
+    logical function f_is_default(var) result(res)
+        !$acc routine seq
+        real(kind(0d0)), intent(in) :: var
+
+        res = f_approx_equal(var, dflt_real)
+    end function f_is_default
+
+    !> Checks if ALL elements of a real(kind(0d0)) array are of default value.
+    !! @param var_array Array to check.
+    logical function f_all_default(var_array) result(res)
+        real(kind(0d0)), intent(in) :: var_array(:)
+        logical :: res_array(size(var_array))
+        integer :: i
+
+        do i = 1, size(var_array)
+            res_array(i) = f_is_default(var_array(i))
+        end do
+
+        res = all(res_array)
+    end function f_all_default
+
+end module m_helper_basic

src/common/m_variables_conversion.fpp

@@ -19,6 +19,8 @@ module m_variables_conversion
 
     use m_mpi_proxy            !< Message passing interface (MPI) module proxy
 
+    use m_helper_basic         !< Functions to compare floating point numbers
+
     use m_helper
     ! ==========================================================================
 
@@ -1005,7 +1007,7 @@ contains
                         qK_prim_vf(i)%sf(j, k, l) = qK_cons_vf(i)%sf(j, k, l)
                     end do
 
-                    if (sigma /= dflt_real) then
+                    if (.not. f_is_default(sigma)) then
                         qK_prim_vf(c_idx)%sf(j, k, l) = qK_cons_vf(c_idx)%sf(j, k, l)
                     end if
 
@@ -1159,7 +1161,7 @@ contains
                         end do
                     end if
 
-                    if (sigma /= dflt_real) then
+                    if (.not. f_is_default(sigma)) then
                         q_cons_vf(c_idx)%sf(j, k, l) = q_prim_vf(c_idx)%sf(j, k, l)
                     end if
 

src/post_process/m_checker.f90

@@ -9,6 +9,8 @@ module m_checker
 
     use m_mpi_proxy            !< Message passing interface (MPI) module proxy
 
+    use m_helper_basic         !< Functions to compare floating point numbers
+
     use m_helper
 
     implicit none

src/post_process/m_global_parameters.fpp

@@ -13,6 +13,8 @@ module m_global_parameters
 #endif
 
     use m_derived_types         !< Definitions of the derived types
+
+    use m_helper_basic          !< Functions to compare floating point numbers
     ! ==========================================================================
 
     implicit none
@@ -489,7 +491,7 @@ contains
                 sys_size = stress_idx%end
             end if
 
-            if (sigma /= dflt_real) then
+            if (.not. f_is_default(sigma)) then
                 c_idx = sys_size + 1
                 sys_size = c_idx
             end if
@@ -515,7 +517,7 @@ contains
             sys_size = internalEnergies_idx%end
             alf_idx = 1 ! dummy, cannot actually have a void fraction
 
-            if (sigma /= dflt_real) then
+            if (.not. f_is_default(sigma)) then
                 c_idx = sys_size + 1
                 sys_size = c_idx
             end if

src/pre_process/m_assign_variables.f90

@@ -9,6 +9,9 @@ module m_assign_variables
     use m_global_parameters     ! Global parameters for the code
 
     use m_variables_conversion  ! Subroutines to change the state variables from
+
+    use m_helper_basic         !< Functions to compare floating point numbers
+
     ! one form to another
     ! ==========================================================================
 
@@ -577,17 +580,17 @@ subroutine s_assign_patch_species_primitive_variables(patch_id, j, k, l, &
 
         if (bubbles .and. (.not. polytropic) .and. (.not. qbmm)) then
             do i = 1, nb
-                if (q_prim_vf(bub_idx%ps(i))%sf(j, k, l) == dflt_real) then
+                if (f_is_default(q_prim_vf(bub_idx%ps(i))%sf(j, k, l))) then
                     q_prim_vf(bub_idx%ps(i))%sf(j, k, l) = pb0(i)
                     ! print *, 'setting to pb0'
                 end if
-                if (q_prim_vf(bub_idx%ms(i))%sf(j, k, l) == dflt_real) then
+                if (f_is_default(q_prim_vf(bub_idx%ms(i))%sf(j, k, l))) then
                     q_prim_vf(bub_idx%ms(i))%sf(j, k, l) = mass_v0(i)
                 end if
             end do
         end if
 
-        if (sigma /= dflt_real) then
+        if (.not. f_is_default(sigma)) then
             q_prim_vf(c_idx)%sf(j, k, l) = eta*patch_icpp(patch_id)%cf_val + &
                                            (1d0 - eta)*patch_icpp(smooth_patch_id)%cf_val
         end if