feat(MathUtil): add modular arithmetic utilities

Author: wpbonelli

autotest/TestMathUtil.f90

@@ -0,0 +1,17 @@
+module TestMathUtil
+  use KindModule, only: I4B, DP
+  use testdrive, only: check, error_type, new_unittest, test_failed, &
+                       to_string, unittest_type
+  use MathUtilModule, only: is_same, is_close
+  implicit none
+  private
+  public :: collect_mathutil
+
+contains
+
+  subroutine collect_mathutil(testsuite)
+    type(unittest_type), allocatable, intent(out) :: testsuite(:)
+    allocate(testsuite(0))
+  end subroutine collect_mathutil
+
+end module TestMathUtil

autotest/meson.build

@@ -5,6 +5,7 @@ if test_drive.found() and not fc_id.contains('intel')
     'DevFeature',
     'GeomUtil',
     'InputOutput',
+    'MathUtil',
     'Sim'
   ]
 

autotest/tester.f90

@@ -6,6 +6,7 @@ program tester
   use TestDevFeature, only: collect_dev_feature
   use TestGeomUtil, only: collect_geomutil
   use TestInputOutput, only: collect_inputoutput
+  use TestMathUtil, only: collect_mathutil
   use TestSim, only: collect_sim
   implicit none
   integer :: stat, is
@@ -19,6 +20,7 @@ program tester
                new_testsuite("DevFeature", collect_dev_feature), &
                new_testsuite("GeomUtil", collect_geomutil), &
                new_testsuite("InputOutput", collect_inputoutput), &
+               new_testsuite("MathUtil", collect_mathutil), &
                new_testsuite("Sim", collect_sim) &
                ]
 

make/makefile

@@ -86,6 +86,7 @@ $(OBJDIR)/version.o \
 $(OBJDIR)/Message.o \
 $(OBJDIR)/Sim.o \
 $(OBJDIR)/OpenSpec.o \
+$(OBJDIR)/MathUtil.o \
 $(OBJDIR)/InputOutput.o \
 $(OBJDIR)/TableTerm.o \
 $(OBJDIR)/Table.o \

msvs/mf6core.vfproj

@@ -395,6 +395,7 @@
 		<File RelativePath="..\src\Utilities\List.f90"/>
 		<File RelativePath="..\src\Utilities\ListReader.f90"/>
 		<File RelativePath="..\src\Utilities\LongLineReader.f90"/>
+		<File RelativePath="..\src\Utilities\MathUtil.f90"/>
 		<File RelativePath="..\src\Utilities\Message.f90"/>
 		<File RelativePath="..\src\Utilities\OpenSpec.f90"/>
 		<File RelativePath="..\src\Utilities\PackageBudget.f90"/>

src/Solution/NumericalSolution.f90

@@ -13,7 +13,8 @@ module NumericalSolutionModule
                              LENMEMPATH
   use MemoryHelperModule, only: create_mem_path
   use TableModule, only: TableType, table_cr
-  use GenericUtilitiesModule, only: is_same, sim_message
+  use GenericUtilitiesModule, only: sim_message
+  use MathUtilModule, only: is_same
   use VersionModule, only: IDEVELOPMODE
   use BaseModelModule, only: BaseModelType
   use BaseExchangeModule, only: BaseExchangeType

src/Utilities/MathUtil.f90

@@ -0,0 +1,87 @@
+module MathUtilModule
+  use KindModule, only: DP, I4B, LGP
+  use ErrorUtilModule, only: pstop
+  use ConstantsModule, only: MAXCHARLEN, LENHUGELINE, &
+                             DZERO, DPREC, DSAME, &
+                             LINELENGTH, LENHUGELINE, VSUMMARY
+
+  implicit none
+  private
+
+  public :: add_wrap, subtr_wrap, incr_wrap, decr_wrap
+
+contains
+
+!> @brief Add a given integer to another, "wrapping around" if needed.
+!!
+!! This function implements modular addition on the non-negative open interval
+!! from 0 to ilimit. If (istart - iadd) > ilimit, the result wraps around to 1.
+!<
+function add_wrap(istart, iadd, ilimit) result(iwrapped)
+  integer(I4B) :: istart, iadd, ilimit, iwrapped
+
+  ! -- A few conditions must hold to ensure only one wraparound is needed
+  if (istart > ilimit) &
+    call pstop(1, 'istart must not exceed ilimit')
+  if (iadd < 0 .or. iadd > ilimit) &
+    call pstop(1, 'iadd most neither be negative nor exceed ilimit')
+
+  ! -- Add iadd to istart.  If the result exceeds ilimit, wrap around
+  ! -- starting with 1.
+  iwrapped = istart + iadd
+  if (iwrapped .gt. ilimit) iwrapped = iwrapped - ilimit
+end function add_wrap
+
+!> @brief Subtract a given integer from another, "wrapping around" if needed.
+!!
+!! This function implements modular subtraction on the non-negative open interval
+!! from 0 to ilimit. If (istart - isubtr) < 0, the result wraps around to ilimit.
+!<
+function subtr_wrap(istart, isubtr, ilimit) result(iwrapped)
+  integer(I4B) :: istart, isubtr, ilimit, iwrapped
+
+  ! -- A few conditions must hold to ensure only one wraparound is needed:
+  if (istart < 1) &
+    call pstop(1, 'istart must be positive')
+  if (isubtr < 0 .or. isubtr > ilimit) &
+    call pstop(1, 'isubtr must neither be negative nor exceed ilimit')
+
+  ! -- Subtract isubtr from istart.  If the result is less than 1, wrap
+  ! -- around starting with ilimit.
+  iwrapped = istart - isubtr
+  if (iwrapped .lt. 1) iwrapped = iwrapped + ilimit
+end function subtr_wrap
+
+!> @brief Increments an index, "wrapping around" if needed
+function incr_wrap(istart, ilimit) result(iwrapped)
+  integer(I4B) :: istart, ilimit, iwrapped
+
+  ! -- Require istart <= ilimit
+  if (istart > ilimit) &
+    call pstop(1, 'istart must be less than or equal to ilimit')
+
+  ! -- Increment istart. If istart is at ilimit, wrap around to 1.
+  if (istart .ne. ilimit) then
+    iwrapped = istart + 1
+  else
+    iwrapped = 1
+  end if
+end function incr_wrap
+
+!> @brief Decrements an index, "wrapping around" if needed
+function decr_wrap(istart, ilimit) result(iwrapped)
+  integer(I4B) :: istart, ilimit, iwrapped
+
+  ! -- Require istart >= 1
+  if (istart < 1) &
+    call pstop(1, 'istart must be greater than or equal to 1')
+
+  ! -- Decrement istart. If istart is 1, wrap around to ilimit.
+  if (istart .ne. 1) then
+    iwrapped = istart - 1
+  else
+    iwrapped = ilimit
+  end if
+end function decr_wrap
+
+end module MathUtilModule

src/meson.build

@@ -237,6 +237,7 @@ modflow_sources = files(
     'Utilities' / 'List.f90',
     'Utilities' / 'ListReader.f90',
     'Utilities' / 'LongLineReader.f90',
+    'Utilities' / 'MathUtil.f90',
     'Utilities' / 'Message.f90',
     'Utilities' / 'OpenSpec.f90',
     'Utilities' / 'PackageBudget.f90',