diff --git a/autotest/__snapshots__/test_prt_dry/test_mf6model[4-prtdry_stay].npy b/autotest/__snapshots__/test_prt_dry/test_mf6model[4-prtdry_stay].npy index 5a921e5f6f3..d8140e81665 100644 Binary files a/autotest/__snapshots__/test_prt_dry/test_mf6model[4-prtdry_stay].npy and b/autotest/__snapshots__/test_prt_dry/test_mf6model[4-prtdry_stay].npy differ diff --git a/autotest/test_prt_dry.py b/autotest/test_prt_dry.py index d68af937af1..72fb1b7ee0d 100644 --- a/autotest/test_prt_dry.py +++ b/autotest/test_prt_dry.py @@ -70,6 +70,8 @@ for i in range(nper): period_data.append((perlen[i], nstp[i], tsmult[i])) +user_time = 100.0 + def build_gwf_sim(name, gwf_ws, mf6, newton=False): sim = flopy.mf6.MFSimulation( @@ -319,6 +321,8 @@ def build_prt_sim(name, gwf, prt_ws, mf6, drape=False, dry_tracking_method=False trackcsv_filerecord=[trackcsvfile], saverecord=[("BUDGET", "ALL")], pname="oc", + ntracktimes=1 if "stay" in name else 0, + tracktimes=[(user_time,)] if "stay" in name else None, ) rel_prt_folder = os.path.relpath(gwf_ws, start=prt_ws) @@ -373,6 +377,8 @@ def check_output(idx, test, snapshot): # ignore particle 4, it terminates early with optimization=2 when built with ifort if "drop" in name: actual = actual.drop(actual[actual.irpt == 4].index) + if "stay" in name: + assert len(actual[actual.t == user_time]) == 5 assert snapshot == actual.to_records(index=False) plot_pathlines = False diff --git a/src/Solution/ParticleTracker/Method.f90 b/src/Solution/ParticleTracker/Method.f90 index 4035a0ad8a1..2369b836d6c 100644 --- a/src/Solution/ParticleTracker/Method.f90 +++ b/src/Solution/ParticleTracker/Method.f90 @@ -191,15 +191,18 @@ end subroutine save !! tracking the particle or terminate it, as well as whether to !! record any output data as per selected reporting conditions. !< - subroutine check(this, particle, cell_defn) + subroutine check(this, particle, cell_defn, tmax) ! modules use TdisModule, only: endofsimulation, totim ! dummy class(MethodType), intent(inout) :: this type(ParticleType), pointer, intent(inout) :: particle type(CellDefnType), pointer, intent(inout) :: cell_defn + real(DP), intent(in) :: tmax ! local logical(LGP) :: dry_cell, dry_particle, no_exit_face, stop_zone, weak_sink + integer(I4B) :: i + real(DP) :: t dry_cell = this%fmi%ibdgwfsat0(cell_defn%icell) == 0 dry_particle = particle%z > cell_defn%top @@ -235,6 +238,8 @@ subroutine check(this, particle, cell_defn) if (dry_cell) then if (particle%idrymeth == 0) then + ! drop to cell bottom. handled by pass + ! to bottom method, nothing to do here no_exit_face = .false. else if (particle%idrymeth == 1) then ! stop @@ -246,23 +251,38 @@ subroutine check(this, particle, cell_defn) ! stay no_exit_face = .false. particle%advancing = .false. + + ! record user tracking times + call this%tracktimes%advance() + if (this%tracktimes%any()) then + do i = this%tracktimes%selection(1), this%tracktimes%selection(2) + t = this%tracktimes%times(i) + if (t < particle%ttrack) cycle + if (t >= tmax) exit + particle%ttrack = t + call this%save(particle, reason=5) + end do + end if + + ! update tracking time to time + ! step end time and save record particle%ttrack = totim + call this%save(particle, reason=2) + ! terminate if last period/step if (endofsimulation) then particle%istatus = 5 call this%save(particle, reason=3) return end if - call this%save(particle, reason=2) end if else if (dry_particle .and. this%name /= "passtobottom") then - ! dry particle if (particle%idrymeth == 0) then ! drop to water table particle%z = cell_defn%top call this%save(particle, reason=1) else if (particle%idrymeth == 1) then - ! terminate + ! stop particle%advancing = .false. particle%istatus = 7 call this%save(particle, reason=3) @@ -271,6 +291,23 @@ subroutine check(this, particle, cell_defn) ! stay no_exit_face = .false. particle%advancing = .false. + + ! record user tracking times + call this%tracktimes%advance() + if (this%tracktimes%any()) then + do i = this%tracktimes%selection(1), this%tracktimes%selection(2) + t = this%tracktimes%times(i) + if (t < particle%ttrack) cycle + if (t >= tmax) exit + particle%ttrack = t + call this%save(particle, reason=5) + end do + end if + + ! update tracking time to time + ! step end time and save record + particle%ttrack = totim + call this%save(particle, reason=2) return end if end if diff --git a/src/Solution/ParticleTracker/MethodCellPassToBot.f90 b/src/Solution/ParticleTracker/MethodCellPassToBot.f90 index 433b3ca9650..c56dab1b104 100644 --- a/src/Solution/ParticleTracker/MethodCellPassToBot.f90 +++ b/src/Solution/ParticleTracker/MethodCellPassToBot.f90 @@ -47,7 +47,7 @@ subroutine apply_ptb(this, particle, tmax) real(DP), intent(in) :: tmax ! Check termination/reporting conditions - call this%check(particle, this%cell%defn) + call this%check(particle, this%cell%defn, tmax) if (.not. particle%advancing) return ! Pass to bottom face diff --git a/src/Solution/ParticleTracker/MethodCellPollock.f90 b/src/Solution/ParticleTracker/MethodCellPollock.f90 index 86b1d09f4ca..06790e8f7f5 100644 --- a/src/Solution/ParticleTracker/MethodCellPollock.f90 +++ b/src/Solution/ParticleTracker/MethodCellPollock.f90 @@ -130,7 +130,7 @@ subroutine apply_mcp(this, particle, tmax) select type (cell => this%cell) type is (CellRectType) ! Check termination/reporting conditions - call this%check(particle, cell%defn) + call this%check(particle, cell%defn, tmax) if (.not. particle%advancing) return ! Transform model coordinates to local cell coordinates diff --git a/src/Solution/ParticleTracker/MethodCellPollockQuad.f90 b/src/Solution/ParticleTracker/MethodCellPollockQuad.f90 index c5f61f2c221..5bd749e192d 100644 --- a/src/Solution/ParticleTracker/MethodCellPollockQuad.f90 +++ b/src/Solution/ParticleTracker/MethodCellPollockQuad.f90 @@ -209,7 +209,7 @@ subroutine apply_mcpq(this, particle, tmax) select type (cell => this%cell) type is (CellRectQuadType) ! Check termination/reporting conditions - call this%check(particle, cell%defn) + call this%check(particle, cell%defn, tmax) if (.not. particle%advancing) return ! Transform model coordinates to local cell coordinates diff --git a/src/Solution/ParticleTracker/MethodCellTernary.f90 b/src/Solution/ParticleTracker/MethodCellTernary.f90 index d6447db6638..03a401173c2 100644 --- a/src/Solution/ParticleTracker/MethodCellTernary.f90 +++ b/src/Solution/ParticleTracker/MethodCellTernary.f90 @@ -165,7 +165,7 @@ subroutine apply_mct(this, particle, tmax) select type (cell => this%cell) type is (CellPolyType) ! Check termination/reporting conditions - call this%check(particle, this%cell%defn) + call this%check(particle, this%cell%defn, tmax) if (.not. particle%advancing) return ! Number of vertices