From 7d14a3746cd24b9cb47ec257238d40e0197cf29a Mon Sep 17 00:00:00 2001 From: wpbonelli Date: Wed, 7 Feb 2024 09:15:48 -0500 Subject: [PATCH] feat(timeselect): add module/type for time selections * support specification of times some event(s) should occur * like timeseries but no time-indexed variable, just times --- autotest/TestTimeSelect.f90 | 111 +++++++++++++++ autotest/meson.build | 3 +- autotest/tester.f90 | 4 +- make/makefile | 1 + msvs/mf6core.vfproj | 1 + src/Model/ModelUtilities/TimeSelect.f90 | 176 ++++++++++++++++++++++++ src/meson.build | 2 + 7 files changed, 296 insertions(+), 2 deletions(-) create mode 100644 autotest/TestTimeSelect.f90 create mode 100644 src/Model/ModelUtilities/TimeSelect.f90 diff --git a/autotest/TestTimeSelect.f90 b/autotest/TestTimeSelect.f90 new file mode 100644 index 00000000000..bf00c988d52 --- /dev/null +++ b/autotest/TestTimeSelect.f90 @@ -0,0 +1,111 @@ +module TestTimeSelect + use KindModule, only: I4B, DP, LGP + use testdrive, only: check, error_type, new_unittest, test_failed, & + to_string, unittest_type + use TimeSelectModule, only: TimeSelectType + use ConstantsModule, only: LINELENGTH + implicit none + private + public :: collect_timeselect + +contains + subroutine collect_timeselect(testsuite) + type(unittest_type), allocatable, intent(out) :: testsuite(:) + testsuite = [ & + new_unittest("is_increasing", test_is_increasing), & + new_unittest("slice", test_slice) & + ] + end subroutine collect_timeselect + + subroutine test_is_increasing(error) + type(error_type), allocatable, intent(out) :: error + type(TimeSelectType) :: ts + + call ts%expand(3) + + ! increasing + ts%times = (/0.0_DP, 1.0_DP, 2.0_DP/) + call check(error, ts%is_increasing()) + + ! not decreasing + ts%times = (/0.0_DP, 0.0_DP, 2.0_DP/) + call check(error,.not. ts%is_increasing()) + + ! decreasing + ts%times = (/2.0_DP, 1.0_DP, 0.0_DP/) + call check(error,.not. ts%is_increasing()) + end subroutine + + subroutine test_slice(error) + type(error_type), allocatable, intent(out) :: error + type(TimeSelectType) :: ts + logical(LGP) :: changed + + call ts%expand(3) + ts%times = (/0.0_DP, 1.0_DP, 2.0_DP/) + call check( & + error, & + size(ts%times) == 3, & + "expected size 3, got"//to_string(size(ts%times))) + + ! empty slice + call ts%set_slice(1.1_DP, 1.9_DP) + call check( & + error, & + ts%slice(1) == -1 .and. ts%slice(2) == -1, & + "empty slice failed, got ["// & + to_string(ts%slice(1))//","//to_string(ts%slice(2))//"]") + + ! single-item slice + call ts%set_slice(0.5_DP, 1.5_DP) + call check( & + error, & + ts%slice(1) == 2 .and. ts%slice(2) == 2, & + "1-item slice failed, got ["// & + to_string(ts%slice(1))//","//to_string(ts%slice(2))//"]") + + ! multi-item slice + changed = .false. + call ts%set_slice(0.5_DP, 2.5_DP, changed=changed) + call check(error, changed) + call check( & + error, & + ts%slice(1) == 2 .and. ts%slice(2) == 3, & + "2-item slice failed, got ["// & + to_string(ts%slice(1))//","//to_string(ts%slice(2))//"]") + + ! no-change + call ts%set_slice(0.1_DP, 2.5_DP, changed=changed) + call check(error,.not. changed) + call check( & + error, & + ts%slice(1) == 2 .and. ts%slice(2) == 3, & + "2-item slice failed, got ["// & + to_string(ts%slice(1))//","//to_string(ts%slice(2))//"]") + + ! lower bound equal to a time value + call ts%set_slice(0.0_DP, 2.5_DP) + call check( & + error, & + ts%slice(1) == 1 .and. ts%slice(2) == 3, & + "lb eq slice failed, got [" & + //to_string(ts%slice(1))//","//to_string(ts%slice(2))//"]") + + ! upper bound equal to a time value + call ts%set_slice(-0.5_DP, 2.0_DP) + call check( & + error, & + ts%slice(1) == 1 .and. ts%slice(2) == 3, & + "ub eq slice failed, got [" & + //to_string(ts%slice(1))//","//to_string(ts%slice(2))//"]") + + ! both bounds equal to a time value + call ts%set_slice(0.0_DP, 2.0_DP) + call check( & + error, & + ts%slice(1) == 1 .and. ts%slice(2) == 3, & + "lb ub eq slice failed, got [" & + //to_string(ts%slice(1))//","//to_string(ts%slice(2))//"]") + + end subroutine test_slice +end module TestTimeSelect diff --git a/autotest/meson.build b/autotest/meson.build index a36d6e1e82e..dc9ee68a221 100644 --- a/autotest/meson.build +++ b/autotest/meson.build @@ -9,7 +9,8 @@ if test_drive.found() and not fc_id.contains('intel') 'List', 'MathUtil', 'Message', - 'Sim' + 'Sim', + 'TimeSelect' ] test_srcs = files( diff --git a/autotest/tester.f90 b/autotest/tester.f90 index ed297ef643c..1c9a92d5466 100644 --- a/autotest/tester.f90 +++ b/autotest/tester.f90 @@ -11,6 +11,7 @@ program tester use TestMathUtil, only: collect_mathutil use TestMessage, only: collect_message use TestSim, only: collect_sim + use TestTimeSelect, only: collect_timeselect implicit none integer :: stat, is character(len=:), allocatable :: suite_name, test_name @@ -27,7 +28,8 @@ program tester new_testsuite("List", collect_list), & new_testsuite("MathUtil", collect_mathutil), & new_testsuite("Message", collect_message), & - new_testsuite("Sim", collect_sim) & + new_testsuite("Sim", collect_sim), & + new_testsuite("TimeSelect", collect_timeselect) & ] call get_argument(1, suite_name) diff --git a/make/makefile b/make/makefile index a70403d508b..4b685431f7f 100644 --- a/make/makefile +++ b/make/makefile @@ -387,6 +387,7 @@ $(OBJDIR)/sparsekit.o \ $(OBJDIR)/rcm.o \ $(OBJDIR)/blas1_d.o \ $(OBJDIR)/Iunit.o \ +$(OBJDIR)/TimeSelect.o \ $(OBJDIR)/RectangularGeometry.o \ $(OBJDIR)/CircularGeometry.o diff --git a/msvs/mf6core.vfproj b/msvs/mf6core.vfproj index 5da12df2fa7..3f4ddb891c0 100644 --- a/msvs/mf6core.vfproj +++ b/msvs/mf6core.vfproj @@ -264,6 +264,7 @@ + diff --git a/src/Model/ModelUtilities/TimeSelect.f90 b/src/Model/ModelUtilities/TimeSelect.f90 new file mode 100644 index 00000000000..600eafa7816 --- /dev/null +++ b/src/Model/ModelUtilities/TimeSelect.f90 @@ -0,0 +1,176 @@ +!> @brief Specify times for some event(s) to occur. +module TimeSelectModule + + use KindModule, only: DP, I4B, LGP + use ConstantsModule, only: DZERO, DONE + use ArrayHandlersModule, only: ExpandArray + use ErrorUtilModule, only: pstop + implicit none + public :: TimeSelectType + + !> @brief Represents a series of instants at which some event should occur. + !! + !! Supports slicing e.g. to filter times within a given period & time step. + !! Array storage can be expanded as needed. Note: array expansion must take + !! place before slicing; whenever expand() is invoked, the slice is cleared. + !! The series is assumed to strictly increase, is_increasing() checks this. + !< + type :: TimeSelectType + real(DP), allocatable :: times(:) + integer(I4B) :: slice(2) + contains + procedure :: destroy + procedure :: expand + procedure :: init + procedure :: is_increasing + procedure :: set_slice + end type TimeSelectType + +contains + + !> @brief Destroy the time selection object. + subroutine destroy(this) + ! -- dummy + class(TimeSelectType) :: this + deallocate (this%times) + end subroutine destroy + + !> @brief Expand capacity by the given amount. Resets the current slice. + subroutine expand(this, increment) + ! -- dummy + class(TimeSelectType) :: this + integer(I4B), optional, intent(in) :: increment + ! -- local + integer(I4B) :: inclocal, isize, newsize + real(DP), allocatable, dimension(:) :: temp + + if (present(increment)) then + inclocal = increment + else + inclocal = 1 + end if + + ! -- expand/repopulate array if already allocated, + ! otherwise allocate it to the given size. the + ! current slice is also reset here. + if (allocated(this%times)) then + isize = size(this%times) + newsize = isize + inclocal + allocate (temp(newsize)) + temp(1:isize) = this%times + deallocate (this%times) + call move_alloc(temp, this%times) + this%slice = (/1, newsize/) + else + allocate (this%times(inclocal)) + this%slice = (/1, inclocal/) + end if + end subroutine expand + + !> @brief Initialize or clear the time selection object. + subroutine init(this) + class(TimeSelectType) :: this + if (.not. allocated(this%times)) allocate (this%times(0)) + this%slice = (/0, 0/) + end subroutine + + !> @brief Determine if times strictly increase. + !! Returns true if empty or not yet allocated. + function is_increasing(this) result(inc) + class(TimeSelectType) :: this + logical(LGP) :: inc + integer(I4B) :: i + real(DP) :: l, t + + inc = .true. + if (.not. allocated(this%times)) return + do i = 1, size(this%times) + t = this%times(i) + if (i /= 1) then + if (l >= t) then + inc = .false. + return + end if + end if + l = t + end do + end function is_increasing + + !> @brief Slice the time selection between t0 and t1 (inclusive). + !! + !! Finds and stores the index of the first time at the same instant + !! as or following the start time, and of the last time at the same + !! instant as or preceding the end time. Allows filtering the times + !! for e.g. a particular stress period and time step. Array indices + !! are assumed to start at 1. If no times are found to fall within + !! the slice (i.e. the slice falls entirely between two consecutive + !! times or beyond the min/max range), indices are set to [-1, -1]. + !! + !! The given start and end times are first checked against currently + !! stored indices to avoid recalculating them if possible, allowing + !! multiple consuming components (e.g., subdomain particle tracking + !! solutions) to share the object efficiently, provided all proceed + !! through stress periods and time steps in lockstep, i.e. they all + !! solve any given period/step before any will proceed to the next. + !< + subroutine set_slice(this, t0, t1, changed) + ! -- dummy + class(TimeSelectType) :: this + real(DP), intent(in) :: t0, t1 + logical(LGP), intent(inout), optional :: changed + ! -- local + integer(I4B) :: i, i0, i1 + integer(I4B) :: l, u, lp, up + real(DP) :: t + + ! -- by default, need to iterate over all times + i0 = 1 + i1 = size(this%times) + + ! -- if no times fall within the slice, set to [-1, -1] + l = -1 + u = -1 + + ! -- previous bounding indices + lp = this%slice(1) + up = this%slice(2) + + ! -- Check if we can reuse either the lower or upper bound. + ! The lower doesn't need to change if it indexes the 1st + ! time simultaneous with or later than the slice's start. + ! The upper doesn't need to change if it indexes the last + ! time before or simultaneous with the slice's end. + if (lp > 0 .and. up > 0) then + if (lp > 1) then + if (this%times(lp - 1) < t0 .and. & + this%times(lp) >= t0) then + l = lp + i0 = l + end if + end if + if (up > 1 .and. up < i1) then + if (this%times(up + 1) > t1 .and. & + this%times(up) <= t1) then + u = up + i1 = u + end if + end if + if (l == lp .and. u == up) then + this%slice = (/l, u/) + if (present(changed)) changed = .false. + return + end if + end if + + ! -- recompute bounding indices if needed + do i = i0, i1 + t = this%times(i) + if (l < 0 .and. t >= t0 .and. t <= t1) l = i + if (l > 0 .and. t <= t1) u = i + end do + this%slice = (/l, u/) + if (present(changed)) changed = l /= lp .or. u /= up + + end subroutine + +end module TimeSelectModule diff --git a/src/meson.build b/src/meson.build index 92a1fef28bb..a9a193201e2 100644 --- a/src/meson.build +++ b/src/meson.build @@ -176,6 +176,8 @@ modflow_sources = files( 'Model' / 'ModelUtilities' / 'PackageMover.f90', 'Model' / 'ModelUtilities' / 'SfrCrossSectionManager.f90', 'Model' / 'ModelUtilities' / 'SfrCrossSectionUtils.f90', + 'Model' / 'ModelUtilities' / 'TimeSelect.f90', + 'Model' / 'ModelUtilities' / 'TspAdvOptions.f90', 'Model' / 'ModelUtilities' / 'SwfCxsUtils.f90', 'Model' / 'ModelUtilities' / 'TspAdvOptions.f90', 'Model' / 'ModelUtilities' / 'UzfCellGroup.f90',