diff --git a/autotest/TestMathUtil.f90 b/autotest/TestMathUtil.f90 index a1fc90a1233..b28e9050d8e 100644 --- a/autotest/TestMathUtil.f90 +++ b/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 diff --git a/make/makefile b/make/makefile index cd2e35d8612..e0c00d9cbf2 100644 --- a/make/makefile +++ b/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 \ diff --git a/msvs/mf6core.vfproj b/msvs/mf6core.vfproj index 6c78c5313af..7d547f60d79 100644 --- a/msvs/mf6core.vfproj +++ b/msvs/mf6core.vfproj @@ -384,7 +384,6 @@ - diff --git a/msvs/mf6lib.vfproj b/msvs/mf6lib.vfproj index 7186195f1e3..dee2f209199 100644 --- a/msvs/mf6lib.vfproj +++ b/msvs/mf6lib.vfproj @@ -143,7 +143,6 @@ - diff --git a/pymake/makefile b/pymake/makefile index 484e7798cf0..8e7065a0d6d 100644 --- a/pymake/makefile +++ b/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 \ diff --git a/src/Model/Connection/GridSorting.f90 b/src/Model/Connection/GridSorting.f90 index 3360b62809d..5d15fa2382c 100644 --- a/src/Model/Connection/GridSorting.f90 +++ b/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. diff --git a/src/Model/GroundWaterFlow/gwf3csub8.f90 b/src/Model/GroundWaterFlow/gwf3csub8.f90 index 73029f582db..db87b8ff411 100644 --- a/src/Model/GroundWaterFlow/gwf3csub8.f90 +++ b/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) diff --git a/src/Model/GroundWaterFlow/gwf3lak8.f90 b/src/Model/GroundWaterFlow/gwf3lak8.f90 index 52ffc9e2916..3834f640263 100644 --- a/src/Model/GroundWaterFlow/gwf3lak8.f90 +++ b/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 diff --git a/src/Solution/LinearMethods/ims8base.f90 b/src/Solution/LinearMethods/ims8base.f90 index 72bc98eb32b..4c82d9f87af 100644 --- a/src/Solution/LinearMethods/ims8base.f90 +++ b/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 diff --git a/src/Solution/NumericalSolution.f90 b/src/Solution/NumericalSolution.f90 index 456b58d54bd..41bc43e576b 100644 --- a/src/Solution/NumericalSolution.f90 +++ b/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 diff --git a/src/Utilities/InputOutput.f90 b/src/Utilities/InputOutput.f90 index 7253b996f3e..913ae044d5b 100644 --- a/src/Utilities/InputOutput.f90 +++ b/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, & diff --git a/src/Utilities/MathUtil.f90 b/src/Utilities/MathUtil.f90 index a33d3332d39..ad6af25ddd0 100644 --- a/src/Utilities/MathUtil.f90 +++ b/src/Utilities/MathUtil.f90 @@ -7,7 +7,7 @@ module MathUtilModule implicit none private - public :: mod_offset + public :: mod_offset, is_close interface mod_offset module procedure :: mod_offset_int, mod_offset_dbl @@ -15,6 +15,67 @@ module MathUtilModule contains + !> @brief Check if a real value is approximately equal to another. + !! + !! By default the determination is symmetric in a and b, as in + !! Python's math.isclose, with relative difference scaled by a + !! factor of the larger absolute value of a and b. The formula + !! is: abs(a - b) <= max(rtol * max(abs(a), abs(b)), atol). + !! + !! If symmetric is set to false the test is asymmetric in a and + !! b, with b taken to be the reference value, and the alternate + !! formula (abs(a - b) <= (atol + rtol * abs(b))) is used. This + !! is the approach taken by numpy.allclose. + !! + !! Defaults for rtol and atol are DSAME and DZERO, respectively. + !! If a or b are near 0 (especially if either is 0), an absolute + !! tolerance suitable for the particular case should be provided. + !! For a justification of a zero absolute tolerance default see: + !! https://peps.python.org/pep-0485/#absolute-tolerance-default + !< + pure logical function is_close(a, b, rtol, atol, symmetric) + ! dummy + real(DP), intent(in) :: a !< first real + real(DP), intent(in) :: b !< second real (reference value if asymmetric) + real(DP), intent(in), optional :: rtol !< relative tolerance (default=DSAME) + real(DP), intent(in), optional :: atol !< absolute tolerance (default=DZERO) + logical(LGP), intent(in), optional :: symmetric !< toggle (a)symmetric comparison + ! local + real(DP) :: lrtol, latol + logical(LGP) :: lsymmetric + + ! check for exact equality + if (a == b) then + is_close = .true. + return + end if + + ! process optional arguments + if (.not. present(rtol)) then + lrtol = DSAME + else + lrtol = rtol + end if + if (.not. present(atol)) then + latol = DZERO + else + latol = atol + end if + if (.not. present(symmetric)) then + lsymmetric = .true. + else + lsymmetric = symmetric + end if + + if (lsymmetric) then + ! "weak" symmetric test, https://peps.python.org/pep-0485/#which-symmetric-test + is_close = abs(a - b) <= max(lrtol * max(abs(a), abs(b)), latol) + else + ! asymmetric, https://numpy.org/doc/stable/reference/generated/numpy.isclose.html + is_close = (abs(a - b) <= (latol + lrtol * abs(b))) + end if + end function is_close + !> @brief Modulo with offset for integer values. pure function mod_offset_int(a, n, d) result(mo) ! -- dummy diff --git a/src/Utilities/TimeSeries/TimeArraySeries.f90 b/src/Utilities/TimeSeries/TimeArraySeries.f90 index 24e1125c31c..7c8c4607bbc 100644 --- a/src/Utilities/TimeSeries/TimeArraySeries.f90 +++ b/src/Utilities/TimeSeries/TimeArraySeries.f90 @@ -4,7 +4,7 @@ module TimeArraySeriesModule use BlockParserModule, only: BlockParserType use ConstantsModule, only: LINELENGTH, UNDEFINED, STEPWISE, LINEAR, & LENTIMESERIESNAME, LENMODELNAME, DZERO, DONE - use GenericUtilitiesModule, only: is_same + use MathUtilModule, only: is_close use InputOutputModule, only: GetUnit, openfile use KindModule, only: DP, I4B use ListModule, only: ListType, ListNodeType @@ -474,7 +474,7 @@ subroutine get_values_at_time(this, nvals, values, time) ierr = 1 end if else - if (is_same(taEarlier%taTime, time)) then + if (is_close(taEarlier%taTime, time)) then values = taEarlier%taArray else ! -- Only earlier time is available, and it is not time of interest; @@ -488,7 +488,7 @@ subroutine get_values_at_time(this, nvals, values, time) end if else if (associated(taLater)) then - if (is_same(taLater%taTime, time)) then + if (is_close(taLater%taTime, time)) then values = taLater%taArray else ! -- only later time is available, and it is not time of interest @@ -717,7 +717,7 @@ subroutine get_latest_preceding_node(this, time, tslNode) if (associated(currNode%nextNode)) then obj => currNode%nextNode%GetItem() ta => CastAsTimeArrayType(obj) - if (ta%taTime < time .or. is_same(ta%taTime, time)) then + if (ta%taTime < time .or. is_close(ta%taTime, time)) then currNode => currNode%nextNode else exit diff --git a/src/Utilities/TimeSeries/TimeSeries.f90 b/src/Utilities/TimeSeries/TimeSeries.f90 index 50904b0b9cd..b4a5d9d48f2 100644 --- a/src/Utilities/TimeSeries/TimeSeries.f90 +++ b/src/Utilities/TimeSeries/TimeSeries.f90 @@ -5,7 +5,7 @@ module TimeSeriesModule use ConstantsModule, only: LINELENGTH, UNDEFINED, STEPWISE, LINEAR, & LINEAREND, LENTIMESERIESNAME, LENHUGELINE, & DZERO, DONE, DNODATA - use GenericUtilitiesModule, only: is_same + use MathUtilModule, only: is_close use InputOutputModule, only: GetUnit, openfile, ParseLine, upcase use ListModule, only: ListType, ListNodeType use SimVariablesModule, only: errmsg @@ -319,7 +319,7 @@ subroutine get_surrounding_records(this, time, tsrecEarlier, tsrecLater) if (associated(currNode%nextNode)) then obj => currNode%nextNode%GetItem() tsr => CastAsTimeSeriesRecordType(obj) - if (tsr%tsrTime < time .and. .not. is_same(tsr%tsrTime, time)) then + if (tsr%tsrTime < time .and. .not. is_close(tsr%tsrTime, time)) then currNode => currNode%nextNode else exit @@ -356,7 +356,7 @@ subroutine get_surrounding_records(this, time, tsrecEarlier, tsrecLater) obj => tsNode1%GetItem() tsrec1 => CastAsTimeSeriesRecordType(obj) time1 = tsrec1%tsrTime - do while (time1 < time .and. .not. is_same(time1, time)) + do while (time1 < time .and. .not. is_close(time1, time)) if (associated(tsNode1%nextNode)) then tsNode1 => tsNode1%nextNode obj => tsNode1%GetItem() @@ -373,8 +373,8 @@ subroutine get_surrounding_records(this, time, tsrecEarlier, tsrecLater) ! end if ! - if (time0 < time .or. is_same(time0, time)) tsrecEarlier => tsrec0 - if (time1 > time .or. is_same(time1, time)) tsrecLater => tsrec1 + if (time0 < time .or. is_close(time0, time)) tsrecEarlier => tsrec0 + if (time1 > time .or. is_close(time1, time)) tsrecLater => tsrec1 ! ! -- Return return @@ -418,7 +418,7 @@ subroutine get_surrounding_nodes(this, time, nodeEarlier, nodeLater) if (associated(currNode%nextNode)) then obj => currNode%nextNode%GetItem() tsr => CastAsTimeSeriesRecordType(obj) - if (tsr%tsrTime < time .and. .not. is_same(tsr%tsrTime, time)) then + if (tsr%tsrTime < time .and. .not. is_close(tsr%tsrTime, time)) then currNode => currNode%nextNode else exit @@ -454,7 +454,7 @@ subroutine get_surrounding_nodes(this, time, nodeEarlier, nodeLater) obj => tsNode1%GetItem() tsrec1 => CastAsTimeSeriesRecordType(obj) time1 = tsrec1%tsrTime - do while (time1 < time .and. .not. is_same(time1, time)) + do while (time1 < time .and. .not. is_close(time1, time)) if (associated(tsNode1%nextNode)) then tsNode1 => tsNode1%nextNode obj => tsNode1%GetItem() @@ -467,11 +467,11 @@ subroutine get_surrounding_nodes(this, time, nodeEarlier, nodeLater) ! end if ! - if (time0 < time .or. is_same(time0, time)) then + if (time0 < time .or. is_close(time0, time)) then tsrecEarlier => tsrec0 nodeEarlier => tsNode0 end if - if (time1 > time .or. is_same(time1, time)) then + if (time1 > time .or. is_close(time1, time)) then tsrecLater => tsrec1 nodeLater => tsNode1 end if @@ -552,7 +552,7 @@ function get_value_at_time(this, time, extendToEndOfSimulation) ierr = 1 end if else - if (extendToEndOfSimulation .or. is_same(tsrEarlier%tsrTime, time)) then + if (extendToEndOfSimulation .or. is_close(tsrEarlier%tsrTime, time)) then get_value_at_time = tsrEarlier%tsrValue else ! -- Only earlier time is available, and it is not time of interest; @@ -566,7 +566,7 @@ function get_value_at_time(this, time, extendToEndOfSimulation) end if else if (associated(tsrLater)) then - if (is_same(tsrLater%tsrTime, time)) then + if (is_close(tsrLater%tsrTime, time)) then get_value_at_time = tsrLater%tsrValue else ! -- only later time is available, and it is not time of interest @@ -624,7 +624,7 @@ function get_integrated_value(this, time0, time1, extendToEndOfSimulation) currObj => currNode%GetItem() currRecord => CastAsTimeSeriesRecordType(currObj) currTime = currRecord%tsrTime - if (is_same(currTime, time1)) then + if (is_close(currTime, time1)) then ! Current node time = time1 so should be ldone ldone = .true. elseif (currTime < time1) then @@ -657,12 +657,12 @@ function get_integrated_value(this, time0, time1, extendToEndOfSimulation) if (lprocess) then ! -- determine lower and upper limits of time span of interest ! within current interval - if (currTime > time0 .or. is_same(currTime, time0)) then + if (currTime > time0 .or. is_close(currTime, time0)) then t0 = currTime else t0 = time0 end if - if (nextTime < time1 .or. is_same(nextTime, time1)) then + if (nextTime < time1 .or. is_close(nextTime, time1)) then t1 = nextTime else t1 = time1 @@ -697,7 +697,7 @@ function get_integrated_value(this, time0, time1, extendToEndOfSimulation) ! -- Are we done yet? if (t1 > time1) then ldone = .true. - elseif (is_same(t1, time1)) then + elseif (is_close(t1, time1)) then ldone = .true. else ! -- We are not done yet @@ -796,7 +796,7 @@ subroutine get_latest_preceding_node(this, time, tslNode) if (associated(currNode%nextNode)) then obj => currNode%nextNode%GetItem() tsr => CastAsTimeSeriesRecordType(obj) - if (tsr%tsrTime < time .or. is_same(tsr%tsrTime, time)) then + if (tsr%tsrTime < time .or. is_close(tsr%tsrTime, time)) then currNode => currNode%nextNode else exit @@ -829,7 +829,7 @@ subroutine get_latest_preceding_node(this, time, tslNode) end do end if ! - if (time0 < time .or. is_same(time0, time)) tslNode => tsNode0 + if (time0 < time .or. is_close(time0, time)) tslNode => tsNode0 ! ! -- Return return @@ -946,7 +946,7 @@ function GetTimeSeriesRecord(this, time, epsi) result(res) do tsr => this%GetNextTimeSeriesRecord() if (associated(tsr)) then - if (is_same(tsr%tsrTime, time)) then + if (is_close(tsr%tsrTime, time)) then res => tsr exit end if diff --git a/src/Utilities/genericutils.f90 b/src/Utilities/genericutils.f90 deleted file mode 100644 index e6921da9dd6..00000000000 --- a/src/Utilities/genericutils.f90 +++ /dev/null @@ -1,66 +0,0 @@ -!> @brief This module contains generic utilties -!! -!! This module contains generic utilities that have -!! limited dependencies. -!! -!< -module GenericUtilitiesModule - use KindModule, only: DP, I4B, LGP - use ConstantsModule, only: MAXCHARLEN, LENHUGELINE, & - DZERO, DPREC, DSAME, & - LINELENGTH, LENHUGELINE, VSUMMARY - use SimVariablesModule, only: istdout, isim_level - ! - implicit none - - private - - public :: is_same - -contains - - !> @brief Function to determine if two reals are the same - !! - !! Function to evaluate if the difference between a and b are less than eps - !! (i.e. a and b are the same). - !! - !< - function is_same(a, b, eps) result(lvalue) - ! -- return variable - logical(LGP) :: lvalue !< boolean indicating if a and b are the same - ! -- dummy variables - real(DP), intent(in) :: a !< first number to evaluate - real(DP), intent(in) :: b !< second number to evaluate - real(DP), intent(in), optional :: eps !< optional maximum difference between a abd b (default=DSAME) - ! -- local variables - real(DP) :: epsloc - real(DP) :: denom - real(DP) :: rdiff - ! - ! -- evaluate optioanl arguments - if (present(eps)) then - epsloc = eps - else - epsloc = DSAME - end if - lvalue = .FALSE. - if (a == b) then - lvalue = .TRUE. - else - if (abs(b) > abs(a)) then - denom = b - else - denom = a - if (abs(denom) == DZERO) then - denom = DPREC - end if - end if - rdiff = abs((a - b) / denom) - if (rdiff <= epsloc) then - lvalue = .TRUE. - end if - end if - - end function is_same - -end module GenericUtilitiesModule diff --git a/src/meson.build b/src/meson.build index 800accc6553..e77d9891412 100644 --- a/src/meson.build +++ b/src/meson.build @@ -227,7 +227,6 @@ modflow_sources = files( 'Utilities' / 'defmacro.F90', 'Utilities' / 'DevFeature.f90', 'Utilities' / 'ErrorUtil.f90', - 'Utilities' / 'genericutils.f90', 'Utilities' / 'GeomUtil.f90', 'Utilities' / 'HashTable.f90', 'Utilities' / 'HeadFileReader.f90', diff --git a/utils/mf5to6/make/makefile b/utils/mf5to6/make/makefile index 00f87156629..fd923f00ba0 100644 --- a/utils/mf5to6/make/makefile +++ b/utils/mf5to6/make/makefile @@ -37,14 +37,14 @@ $(OBJDIR)/compilerversion.o \ $(OBJDIR)/version.o \ $(OBJDIR)/OpenSpec.o \ $(OBJDIR)/GlobalVariables.o \ +$(OBJDIR)/ErrorUtil.o \ $(OBJDIR)/SimPHMF.o \ -$(OBJDIR)/genericutils.o \ +$(OBJDIR)/MathUtil.o \ $(OBJDIR)/InputOutput.o \ $(OBJDIR)/TableTerm.o \ $(OBJDIR)/Table.o \ $(OBJDIR)/MemoryHelper.o \ $(OBJDIR)/CharString.o \ -$(OBJDIR)/ErrorUtil.o \ $(OBJDIR)/Memory.o \ $(OBJDIR)/List.o \ $(OBJDIR)/MemoryList.o \ diff --git a/utils/mf5to6/msvs/mf5to6.vfproj b/utils/mf5to6/msvs/mf5to6.vfproj index f206d7d30aa..10ff7233751 100644 --- a/utils/mf5to6/msvs/mf5to6.vfproj +++ b/utils/mf5to6/msvs/mf5to6.vfproj @@ -103,7 +103,6 @@ - diff --git a/utils/mf5to6/pymake/extrafiles.txt b/utils/mf5to6/pymake/extrafiles.txt index 6b5a9b976c4..f8a3888f945 100644 --- a/utils/mf5to6/pymake/extrafiles.txt +++ b/utils/mf5to6/pymake/extrafiles.txt @@ -9,7 +9,6 @@ ../../../src/Utilities/Constants.f90 ../../../src/Utilities/SimVariables.f90 ../../../src/Utilities/compilerversion.F90 -../../../src/Utilities/genericutils.f90 ../../../src/Utilities/Message.f90 ../../../src/Utilities/ErrorUtil.f90 ../../../src/Utilities/GeomUtil.f90 @@ -17,6 +16,7 @@ ../../../src/Utilities/kind.f90 ../../../src/Utilities/List.f90 ../../../src/Utilities/LongLineReader.f90 +../../../src/Utilities/MathUtil.f90 ../../../src/Utilities/OpenSpec.f90 ../../../src/Utilities/defmacro.F90 ../../../src/Utilities/version.f90 diff --git a/utils/mf5to6/src/MultiLayerObsModule.f90 b/utils/mf5to6/src/MultiLayerObsModule.f90 index f8c77064f46..e2c87430a8e 100644 --- a/utils/mf5to6/src/MultiLayerObsModule.f90 +++ b/utils/mf5to6/src/MultiLayerObsModule.f90 @@ -1,8 +1,8 @@ module MultiLayerObs - - use ConstantsModule, only: DONE, MAXCHARLEN - use ConstantsPHMFModule, only: LENOBSNAMENEW - use GenericUtilitiesModule, only: is_same + + use ConstantsModule, only: DONE, MAXCHARLEN + use ConstantsPHMFModule, only: LENOBSNAMENEW + use MathUtilModule, only: is_close use ListModule, only: ListType use SimModule, only: store_error, ustop @@ -26,7 +26,7 @@ module MultiLayerObs procedure, public :: CheckWeightSum end type - contains +contains ! Non-type-bound procedures @@ -43,7 +43,7 @@ function CastAsLayerObsType(obj) result(res) end select return end function CastAsLayerObsType - + subroutine ConstructLayerObs(newLayerObs, layobname, layer, weight) ! dummy type(LayerObsType), pointer :: newLayerObs @@ -51,7 +51,7 @@ subroutine ConstructLayerObs(newLayerObs, layobname, layer, weight) integer, intent(in) :: layer double precision, intent(in) :: weight ! - allocate(newLayerObs) + allocate (newLayerObs) newLayerObs%lobsname = layobname newLayerObs%layer = layer newLayerObs%weight = weight @@ -71,7 +71,7 @@ subroutine AddLayerObsToList(list, layerobs) ! return end subroutine AddLayerObsToList - + function GetLayerObsFromList(list, indx) result(res) ! dummy type(ListType), intent(inout) :: list @@ -88,18 +88,18 @@ function GetLayerObsFromList(list, indx) result(res) ! return end function GetLayerObsFromList - + subroutine ConstructMLObs(newMLObs, obsname) ! dummy type(MLObsType), pointer, intent(inout) :: newMLObs character(len=LENOBSNAMENEW), intent(in) :: obsname ! - allocate(newMLObs) + allocate (newMLObs) newMLObs%mlobsname = obsname ! return end subroutine ConstructMLObs - + subroutine AddMLObsToList(list, mlobs) ! dummy type(ListType), intent(inout) :: list @@ -112,7 +112,7 @@ subroutine AddMLObsToList(list, mlobs) ! return end subroutine AddMLObsToList - + function GetMLObsFromList(list, indx) result(res) ! dummy type(ListType), intent(inout) :: list @@ -129,37 +129,37 @@ function GetMLObsFromList(list, indx) result(res) ! return end function GetMLObsFromList - + ! Type-bound procedures of MLObsType - + subroutine CheckWeightSum(this) ! dummy class(MLObsType) :: this ! local double precision :: weightsum integer :: i, nlayers - type(LayerObsType), pointer :: layobs => null() + type(LayerObsType), pointer :: layobs => null() character(len=MAXCHARLEN) :: ermsg ! formats - 10 format('Weights of layer observations do not sum to 1.0 for', & - ' multilayer observation: ',a) +10 format('Weights of layer observations do not sum to 1.0 for', & + ' multilayer observation: ', a) ! if (this%summ) return ! weightsum = 0.0d0 nlayers = this%LayerObsList%Count() - do i=1,nlayers + do i = 1, nlayers layobs => GetLayerObsFromList(this%LayerObsList, i) weightsum = weightsum + layobs%weight - enddo + end do ! - if (.not. is_same(weightsum, DONE)) then - write(ermsg,10)trim(this%mlobsname) + if (.not. is_close(weightsum, DONE)) then + write (ermsg, 10) trim(this%mlobsname) call store_error(ermsg) call ustop() - endif + end if ! return end subroutine - + end module MultiLayerObs diff --git a/utils/mf5to6/src/Preproc/ObsBlock.f90 b/utils/mf5to6/src/Preproc/ObsBlock.f90 index 61f3e5bad78..c4e51a3562e 100644 --- a/utils/mf5to6/src/Preproc/ObsBlock.f90 +++ b/utils/mf5to6/src/Preproc/ObsBlock.f90 @@ -3,7 +3,7 @@ module ObsBlockModule use BlockParserModule, only: BlockParserType use ConstantsModule, only: DONE, DZERO, & LINELENGTH, MAXCHARLEN, LENOBSNAME - use GenericUtilitiesModule, only: is_same + use MathUtilModule, only: is_close use ConstantsPHMFModule, only: CONTINUOUS, SINGLE, LENOBSNAMENEW use DnmDis3dModule, only: Dis3dType use GlobalVariablesPHMFModule, only: verbose @@ -184,7 +184,7 @@ subroutine process_block(this, insertLine, WriteBeginEnd, parser) iadjrow = 0 jadjcol = 0 ! - if (.not. is_same(xoff, DZERO)) then + if (.not. is_close(xoff, DZERO)) then if (xoff > DZERO) then if (jcol < ncol) then if (dis3d%idomain(jcol+1, irow, layer) == 1) then @@ -202,7 +202,7 @@ subroutine process_block(this, insertLine, WriteBeginEnd, parser) endif endif ! - if (.not. is_same(yoff, DZERO)) then + if (.not. is_close(yoff, DZERO)) then if (yoff > DZERO) then if (irow > 1) then if (dis3d%idomain(jcol, irow-1, layer) == 1) then diff --git a/utils/mf5to6/src/Preproc/ObservePHMF.f90 b/utils/mf5to6/src/Preproc/ObservePHMF.f90 index d35fd18079d..c50c7ad5f8d 100644 --- a/utils/mf5to6/src/Preproc/ObservePHMF.f90 +++ b/utils/mf5to6/src/Preproc/ObservePHMF.f90 @@ -15,7 +15,7 @@ module ObserveModule use ConstantsModule, only: DONE, DZERO, LENOBSNAME, & LENOBSTYPE, MAXCHARLEN use ConstantsPHMFModule, only: LENOBSNAMENEW, HUGEDBL, HDRYDEFAULT - use GenericUtilitiesModule, only: is_same + use MathUtilModule, only: is_close use ListModule, only: ListType use SimModule, only: store_warning, store_error, & store_error_unit, ustop @@ -206,7 +206,7 @@ subroutine CalcSimVal(this, itime) sumweights = DZERO k = 0 do i=1,nsrc - if (is_same(this%srcvals(itime, i), this%hdry)) then + if (is_close(this%srcvals(itime, i), this%hdry)) then k = k + 1 weights(i) = DZERO else diff --git a/utils/zonebudget/make/makefile b/utils/zonebudget/make/makefile index db96d3c3471..0a848b02b10 100644 --- a/utils/zonebudget/make/makefile +++ b/utils/zonebudget/make/makefile @@ -25,7 +25,7 @@ $(OBJDIR)/compilerversion.o \ $(OBJDIR)/version.o \ $(OBJDIR)/Sim.o \ $(OBJDIR)/OpenSpec.o \ -$(OBJDIR)/genericutils.o \ +$(OBJDIR)/MathUtil.o \ $(OBJDIR)/InputOutput.o \ $(OBJDIR)/LongLineReader.o \ $(OBJDIR)/DevFeature.o \ diff --git a/utils/zonebudget/msvs/zonebudget.vfproj b/utils/zonebudget/msvs/zonebudget.vfproj index f24872d4846..f6e8613efd2 100644 --- a/utils/zonebudget/msvs/zonebudget.vfproj +++ b/utils/zonebudget/msvs/zonebudget.vfproj @@ -44,7 +44,6 @@ - diff --git a/utils/zonebudget/pymake/extrafiles.txt b/utils/zonebudget/pymake/extrafiles.txt index 846b560f0d5..54d11147ff3 100644 --- a/utils/zonebudget/pymake/extrafiles.txt +++ b/utils/zonebudget/pymake/extrafiles.txt @@ -4,9 +4,9 @@ ../../../src/Utilities/Budget.f90 ../../../src/Utilities/Constants.f90 ../../../src/Utilities/compilerversion.F90 -../../../src/Utilities/genericutils.f90 ../../../src/Utilities/ErrorUtil.f90 ../../../src/Utilities/GeomUtil.f90 +../../../src/Utilities/MathUtil.f90 ../../../src/Utilities/InputOutput.f90 ../../../src/Utilities/kind.f90 ../../../src/Utilities/LongLineReader.f90