Skip to content

Commit 9f6b81f

Browse files
ChrisZYJChrisZYJ
authored
Robust floating point comparisons (#495)
1 parent 9630d35 commit 9f6b81f

23 files changed

+216
-168
lines changed

‎src/common/m_checker_common.fpp‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -10,6 +10,8 @@ module m_checker_common
1010

1111
use m_mpi_proxy !< Message passing interface (MPI) module proxy
1212

13+
use m_helper_basic !< Functions to compare floating point numbers
14+
1315
use m_helper
1416

1517
implicit none

‎src/common/m_helper.fpp‎

Lines changed: 1 addition & 56 deletions
Original file line numberDiff line numberDiff line change
@@ -36,10 +36,7 @@ module m_helper
3636
f_create_bbox, &
3737
s_print_2D_array, &
3838
f_xor, &
39-
f_logical_to_int, &
40-
f_approx_equal, &
41-
f_is_default, &
42-
f_all_default
39+
f_logical_to_int
4340

4441
contains
4542

@@ -535,56 +532,4 @@ contains
535532
end if
536533
end function f_logical_to_int
537534

538-
!> This procedure checks if two floating point numbers of kind(0d0) are within tolerance.
539-
!! @param a First number.
540-
!! @param b Second number.
541-
!! @param tol_input Relative error (default = 1d-6).
542-
!! @return Result of the comparison.
543-
logical function f_approx_equal(a, b, tol_input) result(res)
544-
! Reference: https://floating-point-gui.de/errors/comparison/
545-
546-
real(kind(0d0)), intent(in) :: a, b
547-
real(kind(0d0)), optional, intent(in) :: tol_input
548-
real(kind(0d0)) :: tol
549-
550-
if (present(tol_input)) then
551-
if (tol_input <= 0d0) then
552-
call s_mpi_abort('tol_input must be positive. Exiting ...')
553-
end if
554-
tol = tol_input
555-
else
556-
tol = 1d-6
557-
end if
558-
559-
if (a == b) then
560-
res = .true.
561-
else if (a == 0d0 .or. b == 0d0 .or. (abs(a) + abs(b) < tiny(a))) then
562-
res = (abs(a - b) < (tol*tiny(a)))
563-
else
564-
res = (abs(a - b)/min(abs(a) + abs(b), huge(a)) < tol)
565-
end if
566-
end function f_approx_equal
567-
568-
!> Checks if a real(kind(0d0)) variable is of default value.
569-
!! @param var Variable to check.
570-
logical function f_is_default(var) result(res)
571-
real(kind(0d0)), intent(in) :: var
572-
573-
res = f_approx_equal(var, dflt_real)
574-
end function f_is_default
575-
576-
!> Checks if ALL elements of a real(kind(0d0)) array are of default value.
577-
!! @param var_array Array to check.
578-
logical function f_all_default(var_array) result(res)
579-
real(kind(0d0)), intent(in) :: var_array(:)
580-
logical :: res_array(size(var_array))
581-
integer :: i
582-
583-
do i = 1, size(var_array)
584-
res_array(i) = f_is_default(var_array(i))
585-
end do
586-
587-
res = all(res_array)
588-
end function f_all_default
589-
590535
end module m_helper

‎src/common/m_helper_basic.f90‎

Lines changed: 71 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,71 @@
1+
!>
2+
!! @file m_helper_basic.f90
3+
!! @brief Contains module m_helper_basic
4+
5+
module m_helper_basic
6+
7+
! Dependencies =============================================================
8+
9+
use m_derived_types !< Definitions of the derived types
10+
! ==========================================================================
11+
12+
implicit none
13+
14+
private;
15+
public :: f_approx_equal, &
16+
f_is_default, &
17+
f_all_default
18+
19+
contains
20+
21+
!> This procedure checks if two floating point numbers of kind(0d0) are within tolerance.
22+
!! @param a First number.
23+
!! @param b Second number.
24+
!! @param tol_input Relative error (default = 1d-6).
25+
!! @return Result of the comparison.
26+
logical function f_approx_equal(a, b, tol_input) result(res)
27+
! Reference: https://floating-point-gui.de/errors/comparison/
28+
29+
real(kind(0d0)), intent(in) :: a, b
30+
real(kind(0d0)), optional, intent(in) :: tol_input
31+
real(kind(0d0)) :: tol
32+
33+
if (present(tol_input)) then
34+
tol = tol_input
35+
else
36+
tol = 1d-6
37+
end if
38+
39+
if (a == b) then
40+
res = .true.
41+
else if (a == 0d0 .or. b == 0d0 .or. (abs(a) + abs(b) < tiny(a))) then
42+
res = (abs(a - b) < (tol*tiny(a)))
43+
else
44+
res = (abs(a - b)/min(abs(a) + abs(b), huge(a)) < tol)
45+
end if
46+
end function f_approx_equal
47+
48+
!> Checks if a real(kind(0d0)) variable is of default value.
49+
!! @param var Variable to check.
50+
logical function f_is_default(var) result(res)
51+
!$acc routine seq
52+
real(kind(0d0)), intent(in) :: var
53+
54+
res = f_approx_equal(var, dflt_real)
55+
end function f_is_default
56+
57+
!> Checks if ALL elements of a real(kind(0d0)) array are of default value.
58+
!! @param var_array Array to check.
59+
logical function f_all_default(var_array) result(res)
60+
real(kind(0d0)), intent(in) :: var_array(:)
61+
logical :: res_array(size(var_array))
62+
integer :: i
63+
64+
do i = 1, size(var_array)
65+
res_array(i) = f_is_default(var_array(i))
66+
end do
67+
68+
res = all(res_array)
69+
end function f_all_default
70+
71+
end module m_helper_basic

‎src/common/m_variables_conversion.fpp‎

Lines changed: 4 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -19,6 +19,8 @@ module m_variables_conversion
1919

2020
use m_mpi_proxy !< Message passing interface (MPI) module proxy
2121

22+
use m_helper_basic !< Functions to compare floating point numbers
23+
2224
use m_helper
2325
! ==========================================================================
2426

@@ -1005,7 +1007,7 @@ contains
10051007
qK_prim_vf(i)%sf(j, k, l) = qK_cons_vf(i)%sf(j, k, l)
10061008
end do
10071009

1008-
if (sigma /= dflt_real) then
1010+
if (.not. f_is_default(sigma)) then
10091011
qK_prim_vf(c_idx)%sf(j, k, l) = qK_cons_vf(c_idx)%sf(j, k, l)
10101012
end if
10111013

@@ -1159,7 +1161,7 @@ contains
11591161
end do
11601162
end if
11611163

1162-
if (sigma /= dflt_real) then
1164+
if (.not. f_is_default(sigma)) then
11631165
q_cons_vf(c_idx)%sf(j, k, l) = q_prim_vf(c_idx)%sf(j, k, l)
11641166
end if
11651167

‎src/post_process/m_checker.f90‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -9,6 +9,8 @@ module m_checker
99

1010
use m_mpi_proxy !< Message passing interface (MPI) module proxy
1111

12+
use m_helper_basic !< Functions to compare floating point numbers
13+
1214
use m_helper
1315

1416
implicit none

‎src/post_process/m_global_parameters.fpp‎

Lines changed: 4 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -13,6 +13,8 @@ module m_global_parameters
1313
#endif
1414

1515
use m_derived_types !< Definitions of the derived types
16+
17+
use m_helper_basic !< Functions to compare floating point numbers
1618
! ==========================================================================
1719

1820
implicit none
@@ -492,7 +494,7 @@ contains
492494
sys_size = stress_idx%end
493495
end if
494496

495-
if (sigma /= dflt_real) then
497+
if (.not. f_is_default(sigma)) then
496498
c_idx = sys_size + 1
497499
sys_size = c_idx
498500
end if
@@ -518,7 +520,7 @@ contains
518520
sys_size = internalEnergies_idx%end
519521
alf_idx = 1 ! dummy, cannot actually have a void fraction
520522

521-
if (sigma /= dflt_real) then
523+
if (.not. f_is_default(sigma)) then
522524
c_idx = sys_size + 1
523525
sys_size = c_idx
524526
end if

‎src/pre_process/m_assign_variables.f90‎

Lines changed: 6 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -9,6 +9,9 @@ module m_assign_variables
99
use m_global_parameters ! Global parameters for the code
1010

1111
use m_variables_conversion ! Subroutines to change the state variables from
12+
13+
use m_helper_basic !< Functions to compare floating point numbers
14+
1215
! one form to another
1316
! ==========================================================================
1417

@@ -577,17 +580,17 @@ subroutine s_assign_patch_species_primitive_variables(patch_id, j, k, l, &
577580

578581
if (bubbles .and. (.not. polytropic) .and. (.not. qbmm)) then
579582
do i = 1, nb
580-
if (q_prim_vf(bub_idx%ps(i))%sf(j, k, l) == dflt_real) then
583+
if (f_is_default(q_prim_vf(bub_idx%ps(i))%sf(j, k, l))) then
581584
q_prim_vf(bub_idx%ps(i))%sf(j, k, l) = pb0(i)
582585
! print *, 'setting to pb0'
583586
end if
584-
if (q_prim_vf(bub_idx%ms(i))%sf(j, k, l) == dflt_real) then
587+
if (f_is_default(q_prim_vf(bub_idx%ms(i))%sf(j, k, l))) then
585588
q_prim_vf(bub_idx%ms(i))%sf(j, k, l) = mass_v0(i)
586589
end if
587590
end do
588591
end if
589592

590-
if (sigma /= dflt_real) then
593+
if (.not. f_is_default(sigma)) then
591594
q_prim_vf(c_idx)%sf(j, k, l) = eta*patch_icpp(patch_id)%cf_val + &
592595
(1d0 - eta)*patch_icpp(smooth_patch_id)%cf_val
593596
end if

‎src/pre_process/m_check_ib_patches.fpp‎

Lines changed: 30 additions & 28 deletions
Original file line numberDiff line numberDiff line change
@@ -18,6 +18,8 @@ module m_check_ib_patches
1818

1919
use m_compile_specific
2020

21+
use m_helper_basic !< Functions to compare floating point numbers
22+
2123
use m_helper
2224
! ==========================================================================
2325

@@ -82,9 +84,9 @@ contains
8284
! Constraints on the geometric parameters of the circle patch
8385
if (n == 0 .or. p > 0 .or. patch_ib(patch_id)%radius <= 0d0 &
8486
.or. &
85-
patch_ib(patch_id)%x_centroid == dflt_real &
87+
f_is_default(patch_ib(patch_id)%x_centroid) &
8688
.or. &
87-
patch_ib(patch_id)%y_centroid == dflt_real) then
89+
f_is_default(patch_ib(patch_id)%y_centroid)) then
8890

8991
call s_mpi_abort('Inconsistency(ies) detected in '// &
9092
'geometric parameters of circle '// &
@@ -107,8 +109,8 @@ contains
107109
! Constraints on the geometric parameters of the airfoil patch
108110
if (n == 0 .or. p > 0 .or. patch_ib(patch_id)%c <= 0d0 &
109111
.or. patch_ib(patch_id)%p <= 0d0 .or. patch_ib(patch_id)%t <= 0d0 &
110-
.or. patch_ib(patch_id)%m <= 0d0 .or. patch_ib(patch_id)%x_centroid == dflt_real &
111-
.or. patch_ib(patch_id)%y_centroid == dflt_real) then
112+
.or. patch_ib(patch_id)%m <= 0d0 .or. f_is_default(patch_ib(patch_id)%x_centroid) &
113+
.or. f_is_default(patch_ib(patch_id)%y_centroid)) then
112114

113115
call s_mpi_abort('Inconsistency(ies) detected in '// &
114116
'geometric parameters of airfoil '// &
@@ -131,9 +133,9 @@ contains
131133
! Constraints on the geometric parameters of the 3d airfoil patch
132134
if (n == 0 .or. p == 0 .or. patch_ib(patch_id)%c <= 0d0 &
133135
.or. patch_ib(patch_id)%p <= 0d0 .or. patch_ib(patch_id)%t <= 0d0 &
134-
.or. patch_ib(patch_id)%m <= 0d0 .or. patch_ib(patch_id)%x_centroid == dflt_real &
135-
.or. patch_ib(patch_id)%y_centroid == dflt_real .or. patch_ib(patch_id)%z_centroid == dflt_real &
136-
.or. patch_ib(patch_id)%length_z == dflt_real) then
136+
.or. patch_ib(patch_id)%m <= 0d0 .or. f_is_default(patch_ib(patch_id)%x_centroid) &
137+
.or. f_is_default(patch_ib(patch_id)%y_centroid) .or. f_is_default(patch_ib(patch_id)%z_centroid) &
138+
.or. f_is_default(patch_ib(patch_id)%length_z)) then
137139

138140
call s_mpi_abort('Inconsistency(ies) detected in '// &
139141
'geometric parameters of airfoil '// &
@@ -156,9 +158,9 @@ contains
156158
! Constraints on the geometric parameters of the rectangle patch
157159
if (n == 0 .or. p > 0 &
158160
.or. &
159-
patch_ib(patch_id)%x_centroid == dflt_real &
161+
f_is_default(patch_ib(patch_id)%x_centroid) &
160162
.or. &
161-
patch_ib(patch_id)%y_centroid == dflt_real &
163+
f_is_default(patch_ib(patch_id)%y_centroid) &
162164
.or. &
163165
patch_ib(patch_id)%length_x <= 0d0 &
164166
.or. &
@@ -185,11 +187,11 @@ contains
185187
! Constraints on the geometric parameters of the sphere patch
186188
if (n == 0 .or. p == 0 &
187189
.or. &
188-
patch_ib(patch_id)%x_centroid == dflt_real &
190+
f_is_default(patch_ib(patch_id)%x_centroid) &
189191
.or. &
190-
patch_ib(patch_id)%y_centroid == dflt_real &
192+
f_is_default(patch_ib(patch_id)%y_centroid) &
191193
.or. &
192-
patch_ib(patch_id)%z_centroid == dflt_real &
194+
f_is_default(patch_ib(patch_id)%z_centroid) &
193195
.or. &
194196
patch_ib(patch_id)%radius <= 0d0) then
195197

@@ -214,27 +216,27 @@ contains
214216
! Constraints on the geometric parameters of the cylinder patch
215217
if (p == 0 &
216218
.or. &
217-
patch_ib(patch_id)%x_centroid == dflt_real &
219+
f_is_default(patch_ib(patch_id)%x_centroid) &
218220
.or. &
219-
patch_ib(patch_id)%y_centroid == dflt_real &
221+
f_is_default(patch_ib(patch_id)%y_centroid) &
220222
.or. &
221-
patch_ib(patch_id)%z_centroid == dflt_real &
223+
f_is_default(patch_ib(patch_id)%z_centroid) &
222224
.or. &
223225
(patch_ib(patch_id)%length_x <= 0d0 .and. &
224226
patch_ib(patch_id)%length_y <= 0d0 .and. &
225227
patch_ib(patch_id)%length_z <= 0d0) &
226228
.or. &
227229
(patch_ib(patch_id)%length_x > 0d0 .and. &
228-
(patch_ib(patch_id)%length_y /= dflt_real .or. &
229-
patch_ib(patch_id)%length_z /= dflt_real)) &
230+
((.not. f_is_default(patch_ib(patch_id)%length_y)) .or. &
231+
(.not. f_is_default(patch_ib(patch_id)%length_z)))) &
230232
.or. &
231233
(patch_ib(patch_id)%length_y > 0d0 .and. &
232-
(patch_ib(patch_id)%length_x /= dflt_real .or. &
233-
patch_ib(patch_id)%length_z /= dflt_real)) &
234+
((.not. f_is_default(patch_ib(patch_id)%length_x)) .or. &
235+
(.not. f_is_default(patch_ib(patch_id)%length_z)))) &
234236
.or. &
235237
(patch_ib(patch_id)%length_z > 0d0 .and. &
236-
(patch_ib(patch_id)%length_x /= dflt_real .or. &
237-
patch_ib(patch_id)%length_y /= dflt_real)) &
238+
((.not. f_is_default(patch_ib(patch_id)%length_x)) .or. &
239+
(.not. f_is_default(patch_ib(patch_id)%length_y)))) &
238240
.or. &
239241
patch_ib(patch_id)%radius <= 0d0) then
240242

@@ -256,19 +258,19 @@ contains
256258
call s_int_to_str(patch_id, iStr)
257259

258260
! Constraints on the geometric parameters of the inactive patch
259-
if (patch_ib(patch_id)%x_centroid /= dflt_real &
261+
if ((.not. f_is_default(patch_ib(patch_id)%x_centroid)) &
260262
.or. &
261-
patch_ib(patch_id)%y_centroid /= dflt_real &
263+
(.not. f_is_default(patch_ib(patch_id)%y_centroid)) &
262264
.or. &
263-
patch_ib(patch_id)%z_centroid /= dflt_real &
265+
(.not. f_is_default(patch_ib(patch_id)%z_centroid)) &
264266
.or. &
265-
patch_ib(patch_id)%length_x /= dflt_real &
267+
(.not. f_is_default(patch_ib(patch_id)%length_x)) &
266268
.or. &
267-
patch_ib(patch_id)%length_y /= dflt_real &
269+
(.not. f_is_default(patch_ib(patch_id)%length_y)) &
268270
.or. &
269-
patch_ib(patch_id)%length_z /= dflt_real &
271+
(.not. f_is_default(patch_ib(patch_id)%length_z)) &
270272
.or. &
271-
patch_ib(patch_id)%radius /= dflt_real) then
273+
(.not. f_is_default(patch_ib(patch_id)%radius))) then
272274

273275
call s_mpi_abort('Inconsistency(ies) detected in '// &
274276
'geometric parameters of inactive '// &

‎src/pre_process/m_check_patches.fpp‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -18,6 +18,8 @@ module m_check_patches
1818

1919
use m_compile_specific
2020

21+
use m_helper_basic !< Functions to compare floating point numbers
22+
2123
use m_helper
2224
! ==========================================================================
2325

‎src/pre_process/m_checker.fpp‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -9,6 +9,8 @@ module m_checker
99

1010
use m_mpi_proxy !< Message passing interface (MPI) module proxy
1111

12+
use m_helper_basic !< Functions to compare floating point numbers
13+
1214
use m_helper
1315

1416
implicit none

0 commit comments

Comments
 (0)