From 1a2382e36ec15d8ca31e552c2a31a690bf1aa318 Mon Sep 17 00:00:00 2001 From: wpbonelli Date: Thu, 8 Feb 2024 22:16:05 -0500 Subject: [PATCH] add try_advance() routine to timeselect, simplify subcel methods --- autotest/TestTimeSelect.f90 | 48 ++++++++--------- src/Model/ModelUtilities/TimeSelect.f90 | 52 ++++++++++++------- src/Model/ParticleTracking/prt1oc1.f90 | 4 +- src/Model/ParticleTracking/prt1prp1.f90 | 4 +- .../ParticleTracker/MethodSubcellPollock.f90 | 26 ++++------ .../ParticleTracker/MethodSubcellTernary.f90 | 26 ++++------ 6 files changed, 84 insertions(+), 76 deletions(-) diff --git a/autotest/TestTimeSelect.f90 b/autotest/TestTimeSelect.f90 index bf00c988d52..a1b80a9de10 100644 --- a/autotest/TestTimeSelect.f90 +++ b/autotest/TestTimeSelect.f90 @@ -25,15 +25,15 @@ subroutine test_is_increasing(error) ! increasing ts%times = (/0.0_DP, 1.0_DP, 2.0_DP/) - call check(error, ts%is_increasing()) + call check(error, ts%increasing()) ! not decreasing ts%times = (/0.0_DP, 0.0_DP, 2.0_DP/) - call check(error,.not. ts%is_increasing()) + call check(error,.not. ts%increasing()) ! decreasing ts%times = (/2.0_DP, 1.0_DP, 0.0_DP/) - call check(error,.not. ts%is_increasing()) + call check(error,.not. ts%increasing()) end subroutine subroutine test_slice(error) @@ -49,63 +49,63 @@ subroutine test_slice(error) "expected size 3, got"//to_string(size(ts%times))) ! empty slice - call ts%set_slice(1.1_DP, 1.9_DP) + call ts%select(1.1_DP, 1.9_DP) call check( & error, & - ts%slice(1) == -1 .and. ts%slice(2) == -1, & + ts%selection(1) == -1 .and. ts%selection(2) == -1, & "empty slice failed, got ["// & - to_string(ts%slice(1))//","//to_string(ts%slice(2))//"]") + to_string(ts%selection(1))//","//to_string(ts%selection(2))//"]") ! single-item slice - call ts%set_slice(0.5_DP, 1.5_DP) + call ts%select(0.5_DP, 1.5_DP) call check( & error, & - ts%slice(1) == 2 .and. ts%slice(2) == 2, & + ts%selection(1) == 2 .and. ts%selection(2) == 2, & "1-item slice failed, got ["// & - to_string(ts%slice(1))//","//to_string(ts%slice(2))//"]") + to_string(ts%selection(1))//","//to_string(ts%selection(2))//"]") ! multi-item slice changed = .false. - call ts%set_slice(0.5_DP, 2.5_DP, changed=changed) + call ts%select(0.5_DP, 2.5_DP, changed=changed) call check(error, changed) call check( & error, & - ts%slice(1) == 2 .and. ts%slice(2) == 3, & + ts%selection(1) == 2 .and. ts%selection(2) == 3, & "2-item slice failed, got ["// & - to_string(ts%slice(1))//","//to_string(ts%slice(2))//"]") + to_string(ts%selection(1))//","//to_string(ts%selection(2))//"]") ! no-change - call ts%set_slice(0.1_DP, 2.5_DP, changed=changed) + call ts%select(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, & + ts%selection(1) == 2 .and. ts%selection(2) == 3, & "2-item slice failed, got ["// & - to_string(ts%slice(1))//","//to_string(ts%slice(2))//"]") + to_string(ts%selection(1))//","//to_string(ts%selection(2))//"]") ! lower bound equal to a time value - call ts%set_slice(0.0_DP, 2.5_DP) + call ts%select(0.0_DP, 2.5_DP) call check( & error, & - ts%slice(1) == 1 .and. ts%slice(2) == 3, & + ts%selection(1) == 1 .and. ts%selection(2) == 3, & "lb eq slice failed, got [" & - //to_string(ts%slice(1))//","//to_string(ts%slice(2))//"]") + //to_string(ts%selection(1))//","//to_string(ts%selection(2))//"]") ! upper bound equal to a time value - call ts%set_slice(-0.5_DP, 2.0_DP) + call ts%select(-0.5_DP, 2.0_DP) call check( & error, & - ts%slice(1) == 1 .and. ts%slice(2) == 3, & + ts%selection(1) == 1 .and. ts%selection(2) == 3, & "ub eq slice failed, got [" & - //to_string(ts%slice(1))//","//to_string(ts%slice(2))//"]") + //to_string(ts%selection(1))//","//to_string(ts%selection(2))//"]") ! both bounds equal to a time value - call ts%set_slice(0.0_DP, 2.0_DP) + call ts%select(0.0_DP, 2.0_DP) call check( & error, & - ts%slice(1) == 1 .and. ts%slice(2) == 3, & + ts%selection(1) == 1 .and. ts%selection(2) == 3, & "lb ub eq slice failed, got [" & - //to_string(ts%slice(1))//","//to_string(ts%slice(2))//"]") + //to_string(ts%selection(1))//","//to_string(ts%selection(2))//"]") end subroutine test_slice end module TestTimeSelect diff --git a/src/Model/ModelUtilities/TimeSelect.f90 b/src/Model/ModelUtilities/TimeSelect.f90 index 1939933a221..3d8057b971b 100644 --- a/src/Model/ModelUtilities/TimeSelect.f90 +++ b/src/Model/ModelUtilities/TimeSelect.f90 @@ -10,20 +10,21 @@ module TimeSelectModule !> @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. + !! Supports selection e.g. to filter times in a selected 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. + !! place before selection; when expand() is invoked the selection is cleared. + !! The time series is assumed to strictly increase, increasing() checks this. !< type :: TimeSelectType real(DP), allocatable :: times(:) - integer(I4B) :: slice(2) + integer(I4B) :: selection(2) contains procedure :: destroy procedure :: expand procedure :: init - procedure :: is_increasing - procedure :: set_slice + procedure :: increasing + procedure :: select + procedure :: try_advance end type TimeSelectType contains @@ -39,19 +40,19 @@ subroutine expand(this, increment) class(TimeSelectType) :: this integer(I4B), optional, intent(in) :: increment call ExpandArray(this%times, increment=increment) - this%slice = (/1, size(this%times)/) + this%selection = (/1, size(this%times)/) 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/) + this%selection = (/0, 0/) end subroutine !> @brief Determine if times strictly increase. !! Returns true if empty or not yet allocated. - function is_increasing(this) result(inc) + function increasing(this) result(inc) class(TimeSelectType) :: this logical(LGP) :: inc integer(I4B) :: i @@ -69,17 +70,17 @@ function is_increasing(this) result(inc) end if l = t end do - end function is_increasing + end function increasing - !> @brief Slice the time selection between t0 and t1 (inclusive). + !> @brief Select times 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 selection (i.e. it falls entirely between two consecutive + !! times or beyond the time 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 @@ -88,7 +89,7 @@ end function is_increasing !! 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) + subroutine select(this, t0, t1, changed) ! -- dummy class(TimeSelectType) :: this real(DP), intent(in) :: t0, t1 @@ -107,8 +108,8 @@ subroutine set_slice(this, t0, t1, changed) u = -1 ! -- previous bounding indices - lp = this%slice(1) - up = this%slice(2) + lp = this%selection(1) + up = this%selection(2) ! -- Check if we can reuse either the lower or upper bound. ! The lower doesn't need to change if it indexes the 1st @@ -131,7 +132,7 @@ subroutine set_slice(this, t0, t1, changed) end if end if if (l == lp .and. u == up) then - this%slice = (/l, u/) + this%selection = (/l, u/) if (present(changed)) changed = .false. return end if @@ -143,9 +144,24 @@ subroutine set_slice(this, t0, t1, changed) 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/) + this%selection = (/l, u/) if (present(changed)) changed = l /= lp .or. u /= up end subroutine + !> @brief Update the selection to match the current time step. + subroutine try_advance(this) + ! -- modules + use TdisModule, only: kper, kstp, nper, nstp, totimc, delt + ! -- dummy + class(TimeSelectType) :: this + ! -- local + real(DP) :: l, u + l = minval(this%times) + u = maxval(this%times) + if (.not. (kper == 1 .and. kstp == 1)) l = totimc + if (.not. (kper == nper .and. kstp == nstp(kper))) u = totimc + delt + call this%select(l, u) + end subroutine try_advance + end module TimeSelectModule diff --git a/src/Model/ParticleTracking/prt1oc1.f90 b/src/Model/ParticleTracking/prt1oc1.f90 index 104a995c7e8..5eabcc00c28 100644 --- a/src/Model/ParticleTracking/prt1oc1.f90 +++ b/src/Model/ParticleTracking/prt1oc1.f90 @@ -304,7 +304,7 @@ subroutine prt_oc_read_options(this) call this%tracktimes%expand() this%tracktimes%times(size(this%tracktimes%times)) = dval end do ttloop - if (.not. this%tracktimes%is_increasing()) then + if (.not. this%tracktimes%increasing()) then errmsg = "TRACK TIMES MUST STRICTLY INCREASE" call store_error(errmsg) call this%parser%StoreErrorUnit(terminate=.true.) @@ -328,7 +328,7 @@ subroutine prt_oc_read_options(this) read (line, '(f30.0)') dval this%tracktimes%times(i) = dval end do - if (.not. this%tracktimes%is_increasing()) then + if (.not. this%tracktimes%increasing()) then errmsg = "TRACK TIMES MUST STRICTLY INCREASE" call store_error(errmsg) call this%parser%StoreErrorUnit(terminate=.true.) diff --git a/src/Model/ParticleTracking/prt1prp1.f90 b/src/Model/ParticleTracking/prt1prp1.f90 index 2309f668109..eb7c60e14a7 100644 --- a/src/Model/ParticleTracking/prt1prp1.f90 +++ b/src/Model/ParticleTracking/prt1prp1.f90 @@ -725,7 +725,7 @@ subroutine prp_options(this, option, found) call this%releasetimes%expand() this%releasetimes%times(size(this%releasetimes%times)) = dval end do rtloop - if (.not. this%releasetimes%is_increasing()) then + if (.not. this%releasetimes%increasing()) then errmsg = "RELEASE TIMES MUST STRICTLY INCREASE" call store_error(errmsg) call this%parser%StoreErrorUnit(terminate=.true.) @@ -749,7 +749,7 @@ subroutine prp_options(this, option, found) read (line, '(f30.0)') dval this%releasetimes%times(i) = dval end do - if (.not. this%releasetimes%is_increasing()) then + if (.not. this%releasetimes%increasing()) then errmsg = "RELEASE TIMES MUST STRICTLY INCREASE" call store_error(errmsg) call this%parser%StoreErrorUnit(terminate=.true.) diff --git a/src/Solution/ParticleTracker/MethodSubcellPollock.f90 b/src/Solution/ParticleTracker/MethodSubcellPollock.f90 index 85823e3cee5..fd9013a86df 100644 --- a/src/Solution/ParticleTracker/MethodSubcellPollock.f90 +++ b/src/Solution/ParticleTracker/MethodSubcellPollock.f90 @@ -83,7 +83,7 @@ end subroutine apply_msp subroutine track_subcell(this, subcell, particle, tmax) ! modules use ParticleModule, only: get_particle_id - use TdisModule, only: kper, kstp, nper, nstp, totimc, delt + use TdisModule, only: kper, kstp ! dummy class(MethodSubcellPollockType), intent(inout) :: this class(SubcellRectType), intent(in) :: subcell @@ -92,11 +92,12 @@ subroutine track_subcell(this, subcell, particle, tmax) ! local doubleprecision :: vx, dvxdx, vy, dvydy, vz, dvzdz doubleprecision :: dtexitx, dtexity, dtexitz, dtexit, texit, dt, t, t0 - doubleprecision :: x, y, z, l, u + doubleprecision :: x, y, z integer :: statusVX, statusVY, statusVZ, i doubleprecision :: initialX, initialY, initialZ integer :: exitFace integer :: reason + integer(I4B) :: tslice(2) !< user-time slice for the current time step reason = -1 @@ -157,19 +158,14 @@ subroutine track_subcell(this, subcell, particle, tmax) texit = particle%ttrack + dtexit t0 = particle%ttrack - if (all(this%tracktimes%slice > 0)) then - ! -- Select user tracking times to solve. If this is the first time step - ! of the simulation, include all times before it begins; if it is the - ! last time step, include all times after it ends. Otherwise take the - ! times within the current period and time step only. - l = minval(this%tracktimes%times) - u = maxval(this%tracktimes%times) - if (.not. (kper == 1 .and. kstp == 1)) l = totimc - if (.not. (kper == nper .and. kstp == nstp(kper))) u = totimc + delt - call this%tracktimes%set_slice(l, u) - - ! -- Solve for tracking times in this slice - do i = this%tracktimes%slice(1), this%tracktimes%slice(2) + ! -- Select user tracking times to solve. If this is the first time step + ! of the simulation, include all times before it begins; if it is the + ! last time step include all times after it ends, otherwise the times + ! within the current period and time step only. + call this%tracktimes%try_advance() + tslice = this%tracktimes%selection + if (all(tslice > 0)) then + do i = tslice(1), tslice(2) t = this%tracktimes%times(i) if (t >= particle%ttrack .and. t < texit .and. t < tmax) then dt = t - t0 diff --git a/src/Solution/ParticleTracker/MethodSubcellTernary.f90 b/src/Solution/ParticleTracker/MethodSubcellTernary.f90 index e600b6bb42c..f0edb1229a9 100644 --- a/src/Solution/ParticleTracker/MethodSubcellTernary.f90 +++ b/src/Solution/ParticleTracker/MethodSubcellTernary.f90 @@ -66,7 +66,7 @@ end subroutine apply_mst !> @brief Track a particle across a triangular subcell using the ternary method subroutine track_subcell(this, subcell, particle, tmax) ! modules - use TdisModule, only: kper, kstp, nper, nstp, totimc, delt + use TdisModule, only: kper, kstp ! dummy class(MethodSubcellTernaryType), intent(inout) :: this class(SubcellTriType), intent(in) :: subcell @@ -82,12 +82,13 @@ subroutine track_subcell(this, subcell, particle, tmax) double precision :: rot(2, 2), res(2), loc(2) double precision :: alp, bet, alp0, bet0, alp1, bet1, alp2, bet2, alpi, beti double precision :: vzbot, vztop, vzi, vziodz, az, dtexitz - double precision :: dt, t, t0, dtexitxy, texit, x, y, z, l, u + double precision :: dt, t, t0, dtexitxy, texit, x, y, z integer :: izstatus, itopbotexit integer :: ntmax, nsave, isolv, itrifaceenter, itrifaceexit double precision :: tol, step, dtexit, alpexit, betexit integer :: ntdebug ! kluge integer :: reason, i + integer(I4B) :: tslice(2) lbary = .true. ! kluge ntmax = 10000 @@ -196,19 +197,14 @@ subroutine track_subcell(this, subcell, particle, tmax) texit = particle%ttrack + dtexit t0 = particle%ttrack - if (all(this%tracktimes%slice > 0)) then - ! -- Select user tracking times to solve. If this is the first time step - ! of the simulation, include all times before it begins; if it is the - ! last time step, include all times after it ends. Otherwise take the - ! times within the current period and time step only. - l = minval(this%tracktimes%times) - u = maxval(this%tracktimes%times) - if (.not. (kper == 1 .and. kstp == 1)) l = totimc - if (.not. (kper == nper .and. kstp == nstp(kper))) u = totimc + delt - call this%tracktimes%set_slice(l, u) - - ! -- Solve for tracking times in this slice - do i = this%tracktimes%slice(1), this%tracktimes%slice(2) + ! -- Select user tracking times to solve. If this is the first time step + ! of the simulation, include all times before it begins; if it is the + ! last time step, include all times after it ends. Otherwise take the + ! times within the current period and time step only. + call this%tracktimes%try_advance() + tslice = this%tracktimes%selection + if (all(tslice > 0)) then + do i = tslice(1), tslice(2) t = this%tracktimes%times(i) if (t >= particle%ttrack .and. t < texit .and. t < tmax) then dt = t - t0