Skip to content

Commit

Permalink
feat(timeselect): add module/type for time selections
Browse files Browse the repository at this point in the history
* support specification of times some event(s) should occur
* like timeseries but no time-indexed variable, just times
  • Loading branch information
wpbonelli committed Feb 19, 2024
1 parent 45744b5 commit 7d14a37
Show file tree
Hide file tree
Showing 7 changed files with 296 additions and 2 deletions.
111 changes: 111 additions & 0 deletions autotest/TestTimeSelect.f90
Original file line number Diff line number Diff line change
@@ -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
3 changes: 2 additions & 1 deletion autotest/meson.build
Original file line number Diff line number Diff line change
Expand Up @@ -9,7 +9,8 @@ if test_drive.found() and not fc_id.contains('intel')
'List',
'MathUtil',
'Message',
'Sim'
'Sim',
'TimeSelect'
]

test_srcs = files(
Expand Down
4 changes: 3 additions & 1 deletion autotest/tester.f90
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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)
Expand Down
1 change: 1 addition & 0 deletions make/makefile
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down
1 change: 1 addition & 0 deletions msvs/mf6core.vfproj
Original file line number Diff line number Diff line change
Expand Up @@ -264,6 +264,7 @@
<File RelativePath="..\src\Model\ModelUtilities\SfrCrossSectionManager.f90"/>
<File RelativePath="..\src\Model\ModelUtilities\SfrCrossSectionUtils.f90"/>
<File RelativePath="..\src\Model\ModelUtilities\SwfCxsUtils.f90"/>
<File RelativePath="..\src\Model\ModelUtilities\TimeSelect.f90"/>
<File RelativePath="..\src\Model\ModelUtilities\TspAdvOptions.f90"/>
<File RelativePath="..\src\Model\ModelUtilities\UzfCellGroup.f90"/>
<File RelativePath="..\src\Model\ModelUtilities\Xt3dAlgorithm.f90"/>
Expand Down
176 changes: 176 additions & 0 deletions src/Model/ModelUtilities/TimeSelect.f90
Original file line number Diff line number Diff line change
@@ -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
2 changes: 2 additions & 0 deletions src/meson.build
Original file line number Diff line number Diff line change
Expand Up @@ -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',
Expand Down

0 comments on commit 7d14a37

Please sign in to comment.