refactor: introduce MathUtilModule, move is_same, add tests

autotest/TestMathUtil.f90

@@ -1,8 +1,9 @@
 module TestMathUtil
   use KindModule, only: I4B, DP
+  use ConstantsModule, only: DNODATA, DZERO
   use testdrive, only: check, error_type, new_unittest, test_failed, &
                        to_string, unittest_type
-  use MathUtilModule, only: mod_offset
+  use MathUtilModule, only: is_close, mod_offset
   implicit none
   private
   public :: collect_mathutil
@@ -12,6 +13,9 @@ module TestMathUtil
   subroutine collect_mathutil(testsuite)
     type(unittest_type), allocatable, intent(out) :: testsuite(:)
     testsuite = [ &
+                new_unittest("is_close_symmetric", test_is_close_symmetric), &
+                new_unittest("is_close_symmetric_near_0", &
+                             test_is_close_symmetric_near_0), &
                 new_unittest("mod_offset", &
                              test_mod_offset) &
                 ]
@@ -41,4 +45,123 @@ subroutine test_mod_offset(error)
     call check(error, mod_offset(2.0_DP, 3.0_DP, 3.0_DP) == 5.0_DP)
   end subroutine test_mod_offset
 
+  subroutine test_is_close_symmetric(error)
+    type(error_type), allocatable, intent(out) :: error
+    real(DP) :: a, b, rtol
+
+    ! exact match
+    a = 1.0_DP
+    b = 1.0_DP
+    call check(error, is_close(a, b), &
+               "exp eq: a="//to_string(a)// &
+               ", b="//to_string(b)// &
+               ", eps=default")
+    if (allocated(error)) return
+
+    ! mismatch with default rtol
+    b = 1.0001_DP
+    call check(error, (.not. (is_close(a, b))), &
+               "exp ne: a="//to_string(a)// &
+               ", b="//to_string(b)// &
+               ", eps=default")
+    if (allocated(error)) return
+
+    ! inexact match with large rtol
+    rtol = 1d-2
+    call check(error, is_close(a, b, rtol=rtol), &
+               "exp eq: a="//to_string(a)// &
+               ", b="//to_string(b)// &
+               ", rtol="//to_string(rtol))
+    if (allocated(error)) return
+
+    ! mismatch when we reduce rtol
+    rtol = 0.5d-5
+    call check(error, (.not. is_close(a, b, rtol=rtol)), &
+               "exp ne: a="//to_string(a)// &
+               ", b="//to_string(b)// &
+               ", rtol="//to_string(rtol))
+    if (allocated(error)) return
+
+    ! +/-0
+    call check(error, is_close(0.0_DP, -0.0_DP), &
+               "exp ne: a="//to_string(a)// &
+               ", b="//to_string(b)// &
+               ", eps=default")
+
+    ! DNODATA
+    call check(error, (.not. is_close(0.0_DP, DNODATA)), &
+               "exp ne: a="//to_string(a)// &
+               ", b="//to_string(b)// &
+               ", eps=default")
+    call check(error, is_close(DNODATA, DNODATA), &
+               "exp ne: a="//to_string(a)// &
+               ", b="//to_string(b)// &
+               ", eps=default")
+    call check(error, (.not. is_close(DNODATA, DNODATA / 10)), &
+               "exp ne: a="//to_string(a)// &
+               ", b="//to_string(b)// &
+               ", eps=default")
+    call check(error, (.not. is_close(DNODATA, DNODATA * 10)), &
+               "exp ne: a="//to_string(a)// &
+               ", b="//to_string(b)// &
+               ", eps=default")
+
+  end subroutine test_is_close_symmetric
+
+  subroutine test_is_close_symmetric_near_0(error)
+    type(error_type), allocatable, intent(out) :: error
+    real(DP) :: a, b, rtol, atol
+
+    a = 0.0_DP
+    b = 0.0_DP
+    call check(error, is_close(a, b), &
+               "exp eq: a="//to_string(a)// &
+               ", b="//to_string(b)// &
+               ", rtol=default")
+    if (allocated(error)) return
+
+    a = DZERO
+    b = DZERO
+    call check(error, is_close(a, b), &
+               "exp eq: a="//to_string(a)// &
+               ", b="//to_string(b)// &
+               ", rtol=default")
+    if (allocated(error)) return
+
+    b = 1d-4
+    call check(error, (.not. is_close(a, b)), &
+               "exp eq: a="//to_string(a)// &
+               ", b="//to_string(b)// &
+               ", rtol=default")
+    if (allocated(error)) return
+
+    rtol = 0.999_DP
+    call check(error, &
+               ! expect failure, see above
+               (.not. is_close(a, b, rtol=rtol)), &
+               "exp eq: a="//to_string(a)// &
+               ", b="//to_string(b)// &
+               ", rtol="//to_string(rtol))
+    if (allocated(error)) return
+
+    ! absolute comparison is appropriate when a and/or b are near or equal to 0
+    b = 1d-4
+    atol = 1d-3
+    call check(error, is_close(a, b, atol=atol), &
+               "exp eq: a="//to_string(a)// &
+               ", b="//to_string(b)// &
+               ", atol="//to_string(atol))
+    if (allocated(error)) return
+
+    ! make sure the absolute tolerance is applied
+    b = 1d-4
+    atol = 1d-5
+    call check(error, (.not. is_close(a, b, atol=atol)), &
+               "exp eq: a="//to_string(a)// &
+               ", b="//to_string(b)// &
+               ", atol="//to_string(atol))
+    if (allocated(error)) return
+
+  end subroutine test_is_close_symmetric_near_0
+
 end module TestMathUtil

make/makefile

@@ -85,7 +85,7 @@ $(OBJDIR)/compilerversion.o \
 $(OBJDIR)/version.o \
 $(OBJDIR)/Sim.o \
 $(OBJDIR)/OpenSpec.o \
-$(OBJDIR)/genericutils.o \
+$(OBJDIR)/MathUtil.o \
 $(OBJDIR)/InputOutput.o \
 $(OBJDIR)/TableTerm.o \
 $(OBJDIR)/Table.o \
@@ -324,7 +324,6 @@ $(OBJDIR)/BaseGeometry.o \
 $(OBJDIR)/mf6.o \
 $(OBJDIR)/StringList.o \
 $(OBJDIR)/MemorySetHandler.o \
-$(OBJDIR)/MathUtil.o \
 $(OBJDIR)/ilut.o \
 $(OBJDIR)/sparsekit.o \
 $(OBJDIR)/rcm.o \

msvs/mf6core.vfproj

@@ -384,7 +384,6 @@
 			<FileConfiguration Name="Release|Win32">
 				<Tool Name="VFFortranCompilerTool" Preprocess="preprocessYes"/></FileConfiguration></File>
 		<File RelativePath="..\src\Utilities\DevFeature.f90"/>
-		<File RelativePath="..\src\Utilities\genericutils.f90"/>
 		<File RelativePath="..\src\Utilities\ErrorUtil.f90"/>
 		<File RelativePath="..\src\Utilities\GeomUtil.f90"/>
 		<File RelativePath="..\src\Utilities\HashTable.f90"/>

msvs/mf6lib.vfproj

@@ -143,7 +143,6 @@
 		<File RelativePath="..\src\Utilities\comarg.f90"/>
 		<File RelativePath="..\src\Utilities\compilerversion.fpp"/>
 		<File RelativePath="..\src\Utilities\Constants.f90"/>
-		<File RelativePath="..\src\Utilities\genericutils.f90"/>
 		<File RelativePath="..\src\Utilities\ErrorUtil.f90"/>
 		<File RelativePath="..\src\Utilities\GeomUtil.f90"/>
 		<File RelativePath="..\src\Utilities\HashTable.f90"/>

pymake/makefile

@@ -114,7 +114,6 @@ $(OBJDIR)/GwfNpfOptions.o \
 $(OBJDIR)/GwfBuyInputData.o \
 $(OBJDIR)/ims8misc.o \
 $(OBJDIR)/LinearSolverBase.o \
-$(OBJDIR)/genericutils.o \
 $(OBJDIR)/ArrayHandlers.o \
 $(OBJDIR)/IndexMap.o \
 $(OBJDIR)/version.o \

src/Model/Connection/GridSorting.f90

@@ -2,7 +2,7 @@ module GridSorting
   use KindModule, only: I4B, DP, LGP
   use ConstantsModule, only: DHALF
   use CellWithNbrsModule, only: GlobalCellType
-  use GenericUtilitiesModule, only: is_same
+  use MathUtilModule, only: is_close
   use BaseDisModule, only: dis_transform_xy
   implicit none
   private
@@ -75,11 +75,11 @@ function lessThan(n, m) result(isLess)
                     dis_bot_m(gcm%index))
 
       ! compare
-      if (.not. is_same(zn, zm, 10 * epsilon(zn))) then
+      if (.not. is_close(zn, zm, 10 * epsilon(zn))) then
         isLess = zn > zm
-      else if (.not. is_same(yn, ym, 10 * epsilon(yn))) then
+      else if (.not. is_close(yn, ym, 10 * epsilon(yn))) then
         isLess = yn > ym
-      else if (.not. is_same(xn, xm, 10 * epsilon(xn))) then
+      else if (.not. is_close(xn, xm, 10 * epsilon(xn))) then
         isLess = xn < xm
       else
         isLess = .false.

src/Model/GroundWaterFlow/gwf3csub8.f90

@@ -18,7 +18,7 @@ module GwfCsubModule
                              TABLEFT, TABCENTER, TABRIGHT, &
                              TABSTRING, TABUCSTRING, TABINTEGER, TABREAL
   use MemoryHelperModule, only: create_mem_path
-  use GenericUtilitiesModule, only: is_same
+  use MathUtilModule, only: is_close
   use MessageModule, only: write_message
   use SmoothingModule, only: sQuadraticSaturation, &
                              sQuadraticSaturationDerivative, &
@@ -2875,7 +2875,7 @@ subroutine csub_fc(this, kiter, hold, hnew, matrix_sln, idxglo, rhs)
           rhs(node) = rhs(node) + rhsterm
           !
           ! -- calculate interbed water compressibility terms
-          if (this%brg /= DZERO .and. idelay == 0) then
+          if (.not. is_close(this%brg, DZERO) .and. idelay == 0) then
             call this%csub_nodelay_wcomp_fc(ib, node, tled, area, &
                                             hnew(node), hold(node), &
                                             hcof, rhsterm)

src/Model/GroundWaterFlow/gwf3lak8.f90

@@ -28,7 +28,7 @@ module LakModule
   use BaseDisModule, only: DisBaseType
   use SimModule, only: count_errors, store_error, store_error_unit, &
                        deprecation_warning
-  use GenericUtilitiesModule, only: is_same
+  use MathUtilModule, only: is_close
   use BlockParserModule, only: BlockParserType
   use BaseDisModule, only: DisBaseType
   use SimVariablesModule, only: errmsg, warnmsg
@@ -833,7 +833,7 @@ subroutine lak_read_lake_connections(this)
                                    warnmsg, this%parser%GetUnit())
         case default
           read (keyword, *) rval
-          if (is_same(rval, DNODATA)) then
+          if (is_close(rval, DNODATA)) then
             is_lake_bed = .FALSE.
           else
             is_lake_bed = .TRUE.
@@ -1832,7 +1832,7 @@ subroutine lak_read_initial_attr(this)
           end if
           length = this%connlength(j)
         end if
-        if (is_same(this%bedleak(j), DNODATA)) then
+        if (is_close(this%bedleak(j), DNODATA)) then
           clb(j) = DNODATA
         else if (this%bedleak(j) > DZERO) then
           clb(j) = DONE / this%bedleak(j)
@@ -1844,7 +1844,7 @@ subroutine lak_read_initial_attr(this)
         else
           caq(j) = DZERO
         end if
-        if (is_same(this%bedleak(j), DNODATA)) then
+        if (is_close(this%bedleak(j), DNODATA)) then
           this%satcond(j) = area / caq(j)
         else if (clb(j) * caq(j) > DZERO) then
           this%satcond(j) = area / (clb(j) + caq(j))
@@ -1879,7 +1879,7 @@ subroutine lak_read_initial_attr(this)
           nn = this%cellid(j)
           area = this%warea(j)
           c1 = DZERO
-          if (is_same(clb(j), DNODATA)) then
+          if (is_close(clb(j), DNODATA)) then
             cbedleak = '     NONE     '
             cbedcond = '     NONE     '
           else if (clb(j) > DZERO) then
@@ -2604,7 +2604,7 @@ subroutine lak_calculate_evaporation(this, ilak, stage, avail, ev)
     call this%lak_calculate_sarea(ilak, stage, sa)
     ev = sa * this%evaporation(ilak)
     if (ev > avail) then
-      if (is_same(avail, DPREC)) then
+      if (is_close(avail, DPREC)) then
         ev = DZERO
       else
         ev = -avail

src/Solution/LinearMethods/ims8base.f90

@@ -9,7 +9,7 @@ MODULE IMSLinearBaseModule
   use KindModule, only: DP, I4B
   use ConstantsModule, only: LINELENGTH, IZERO, &
                              DZERO, DPREC, DEM6, DEM3, DHALF, DONE
-  use GenericUtilitiesModule, only: is_same
+  use MathUtilModule, only: is_close
   use BlockParserModule, only: BlockParserType
   use IMSReorderingModule, only: ims_odrv
   use ConvergenceSummaryModule
@@ -211,7 +211,7 @@ SUBROUTINE ims_base_cg(ICNVG, ITMAX, INNERIT, &
       IF (ICNVG .NE. 0) EXIT INNER
       !
       ! -- CHECK THAT CURRENT AND PREVIOUS rho ARE DIFFERENT
-      lsame = is_same(rho, rho0)
+      lsame = is_close(rho, rho0)
       IF (lsame) THEN
         EXIT INNER
       END IF
@@ -513,15 +513,15 @@ SUBROUTINE ims_base_bcgs(ICNVG, ITMAX, INNERIT, &
       !
       ! -- CHECK THAT CURRENT AND PREVIOUS rho, alpha, AND omega ARE
       !    DIFFERENT
-      lsame = is_same(rho, rho0)
+      lsame = is_close(rho, rho0)
       IF (lsame) THEN
         EXIT INNER
       END IF
-      lsame = is_same(alpha, alpha0)
+      lsame = is_close(alpha, alpha0)
       IF (lsame) THEN
         EXIT INNER
       END IF
-      lsame = is_same(omega, omega0)
+      lsame = is_close(omega, omega0)
       IF (lsame) THEN
         EXIT INNER
       END IF

src/Solution/NumericalSolution.f90

@@ -13,8 +13,8 @@ module NumericalSolutionModule
                              LENMEMPATH
   use MemoryHelperModule, only: create_mem_path
   use TableModule, only: TableType, table_cr
-  use GenericUtilitiesModule, only: is_same
   Use MessageModule, only: write_message
+  use MathUtilModule, only: is_close
   use VersionModule, only: IDEVELOPMODE
   use BaseModelModule, only: BaseModelType
   use BaseExchangeModule, only: BaseExchangeType
@@ -2476,7 +2476,7 @@ subroutine sln_ls(this, kiter, kstp, kper, in_iter, iptc, ptcf)
           end if
         end if
       else
-        lsame = is_same(l2norm, this%l2norm0)
+        lsame = is_close(l2norm, this%l2norm0)
         if (lsame) then
           iptc = 0
         end if

src/Utilities/InputOutput.f90

@@ -11,7 +11,6 @@ module InputOutputModule
                              TABLEFT, TABCENTER, TABRIGHT,                     &
                              TABSTRING, TABUCSTRING, TABINTEGER, TABREAL,      &
                              DZERO
-  use GenericUtilitiesModule, only: is_same
   use MessageModule, only: write_message
   private
   public :: GetUnit,                                                           &