From ba29c04395a48a809a2ec8f17e44663b664f201f Mon Sep 17 00:00:00 2001 From: wpbonelli Date: Fri, 21 Jun 2024 08:16:50 -0400 Subject: [PATCH] fix(prt): make extended tracking opt-in (#1888) Previously the PRT model's default behavior was to track particles until termination, as with MODPATH 7 stop time option 2 (extend). Under extended tracking, the program may not halt if particles enter a cycle in the flow system. This PR changes the default to the equivalent of MP7 stop time option 1 (final), terminating at simulation end unless a new Particle Release Point (PRP) package keyword option EXTEND is provided. This is meant to provide a stronger guarantee that the program halts under default settings. The extend option is stored as an integer rather than logical in case we need to enumerate more alternatives in future. This PR also avoids unnecessary iterations in MethodSubcellPollock/Ternary, and changes the default visibility of ParticleStore members to private, as per style conventions. --- autotest/test_prt_budget.py | 1 + autotest/test_prt_disv1.py | 1 + autotest/test_prt_drape.py | 1 + autotest/test_prt_exg.py | 1 + autotest/test_prt_fmi.py | 53 ++++++++++-------- autotest/test_prt_quad_refinement.py | 1 + autotest/test_prt_release_timing.py | 1 + autotest/test_prt_stop_zones.py | 1 + autotest/test_prt_ternary_methods.py | 1 + autotest/test_prt_track_events.py | 1 + autotest/test_prt_triangle.py | 1 + autotest/test_prt_voronoi1.py | 1 + autotest/test_prt_voronoi2.py | 1 + autotest/test_prt_weak_sinks.py | 1 + doc/ReleaseNotes/develop.tex | 1 + doc/mf6io/mf6ivar/dfn/prt-prp.dfn | 8 +++ doc/mf6io/mf6ivar/md/mf6ivar.md | 1 + doc/mf6io/mf6ivar/tex/prt-prp-desc.tex | 2 + doc/mf6io/mf6ivar/tex/prt-prp-options.dat | 1 + doc/mf6io/prt/prt.tex | 2 +- src/Model/ParticleTracking/prt-prp.f90 | 8 +++ src/Model/ParticleTracking/prt.f90 | 23 +++++--- .../ParticleTracker/MethodSubcellPollock.f90 | 7 ++- .../ParticleTracker/MethodSubcellTernary.f90 | 11 ++-- src/Solution/ParticleTracker/Particle.f90 | 56 +++++++++++-------- 25 files changed, 120 insertions(+), 66 deletions(-) diff --git a/autotest/test_prt_budget.py b/autotest/test_prt_budget.py index dc53ea29b4d..b8dc8f26809 100644 --- a/autotest/test_prt_budget.py +++ b/autotest/test_prt_budget.py @@ -89,6 +89,7 @@ def build_prt_sim(name, gwf_ws, prt_ws, mf6): stop_at_weak_sink=False, boundnames=True, exit_solve_tolerance=DEFAULT_EXIT_SOLVE_TOL, + extend_tracking=True, ) # create output control package diff --git a/autotest/test_prt_disv1.py b/autotest/test_prt_disv1.py index e35101b2b3d..8c4750a9eb9 100644 --- a/autotest/test_prt_disv1.py +++ b/autotest/test_prt_disv1.py @@ -236,6 +236,7 @@ def build_prt_sim(idx, gwf_ws, prt_ws, mf6): print_input=True, dev_forceternary=i == 1, exit_solve_tolerance=1e-10, + extend_tracking=True, ) # create output control package diff --git a/autotest/test_prt_drape.py b/autotest/test_prt_drape.py index 417f50edf79..a2aebf863c2 100644 --- a/autotest/test_prt_drape.py +++ b/autotest/test_prt_drape.py @@ -190,6 +190,7 @@ def build_prt_sim(name, gwf_ws, prt_ws, mf6): trackcsv_filerecord=[prp_track_csv_file], drape="drp" in name, exit_solve_tolerance=DEFAULT_EXIT_SOLVE_TOL, + extend_tracking=True, ) # create output control package diff --git a/autotest/test_prt_exg.py b/autotest/test_prt_exg.py index 00d6e2276d3..ecac06bda66 100644 --- a/autotest/test_prt_exg.py +++ b/autotest/test_prt_exg.py @@ -79,6 +79,7 @@ def build_mf6_sim(idx, test): perioddata={0: ["FIRST"]}, boundnames="bnms" in name, exit_solve_tolerance=DEFAULT_EXIT_SOLVE_TOL, + extend_tracking=True, ) # create output control package diff --git a/autotest/test_prt_fmi.py b/autotest/test_prt_fmi.py index 1aa15012331..736c4f108c9 100644 --- a/autotest/test_prt_fmi.py +++ b/autotest/test_prt_fmi.py @@ -48,8 +48,8 @@ DEFAULT_EXIT_SOLVE_TOL, ) -simname = "prtfmi01" -cases = [simname, f"{simname}saws", f"{simname}bprp"] +simname = "prtfmi" +cases = [simname, f"{simname}saws", f"{simname}bprp", f"{simname}noext"] def build_prt_sim(name, gwf_ws, prt_ws, mf6): @@ -123,6 +123,7 @@ def build_prt_sim(name, gwf_ws, prt_ws, mf6): stop_at_weak_sink="saws" in prt_name, boundnames=True, exit_solve_tolerance=DEFAULT_EXIT_SOLVE_TOL, + extend_tracking="noext" not in prt_name, ) # create output control package @@ -346,28 +347,32 @@ def check_output(idx, test): plt.show() plt.savefig(gwf_ws / f"test_{simname}.png") - # convert mf6 pathlines to mp7 format - mf6_pls = to_mp7_pathlines(mf6_pls) - - # sort both dataframes by particleid and time - mf6_pls.sort_values(by=["particleid", "time"], inplace=True) - mp7_pls.sort_values(by=["particleid", "time"], inplace=True) - - # drop columns for which there is no direct correspondence between mf6 and mp7 - del mf6_pls["sequencenumber"] - del mf6_pls["particleidloc"] - del mf6_pls["xloc"] - del mf6_pls["yloc"] - del mf6_pls["zloc"] - del mp7_pls["sequencenumber"] - del mp7_pls["particleidloc"] - del mp7_pls["xloc"] - del mp7_pls["yloc"] - del mp7_pls["zloc"] - - # compare mf6 / mp7 pathline data - assert mf6_pls.shape == mp7_pls.shape - assert np.allclose(mf6_pls, mp7_pls, atol=1e-3) + if "noext" in name: + # maximum tracking time should be simulation stop time + assert mf6_pls.t.max() == FlopyReadmeCase.nper * FlopyReadmeCase.perlen + else: + # convert mf6 pathlines to mp7 format + mf6_pls = to_mp7_pathlines(mf6_pls) + + # sort both dataframes by particleid and time + mf6_pls.sort_values(by=["particleid", "time"], inplace=True) + mp7_pls.sort_values(by=["particleid", "time"], inplace=True) + + # drop columns for which there is no direct correspondence between mf6 and mp7 + del mf6_pls["sequencenumber"] + del mf6_pls["particleidloc"] + del mf6_pls["xloc"] + del mf6_pls["yloc"] + del mf6_pls["zloc"] + del mp7_pls["sequencenumber"] + del mp7_pls["particleidloc"] + del mp7_pls["xloc"] + del mp7_pls["yloc"] + del mp7_pls["zloc"] + + # compare mf6 / mp7 pathline data + assert mf6_pls.shape == mp7_pls.shape + assert np.allclose(mf6_pls, mp7_pls, atol=1e-3) @pytest.mark.parametrize("idx, name", enumerate(cases)) diff --git a/autotest/test_prt_quad_refinement.py b/autotest/test_prt_quad_refinement.py index ced713c2a98..997f354811f 100644 --- a/autotest/test_prt_quad_refinement.py +++ b/autotest/test_prt_quad_refinement.py @@ -154,6 +154,7 @@ def build_mf6_sim(idx, test, **kwargs): perioddata={0: ["FIRST"]}, exit_solve_tolerance=DEFAULT_EXIT_SOLVE_TOL, dev_forceternary=tracking_method == "ternary", + extend_tracking=True, ) prt_track_file = f"{prt_name}.trk" prt_track_csv_file = f"{prt_name}.trk.csv" diff --git a/autotest/test_prt_release_timing.py b/autotest/test_prt_release_timing.py index 0977ca3bd27..f17b0a91a40 100644 --- a/autotest/test_prt_release_timing.py +++ b/autotest/test_prt_release_timing.py @@ -177,6 +177,7 @@ def build_prt_sim(name, gwf_ws, prt_ws, mf6, fraction=None): ), print_input=True, exit_solve_tolerance=DEFAULT_EXIT_SOLVE_TOL, + extend_tracking=True, ) # create output control package diff --git a/autotest/test_prt_stop_zones.py b/autotest/test_prt_stop_zones.py index 6b9fffb6393..d6cfcb0037d 100644 --- a/autotest/test_prt_stop_zones.py +++ b/autotest/test_prt_stop_zones.py @@ -131,6 +131,7 @@ def build_prt_sim(name, gwf_ws, prt_ws, mf6): perioddata={0: ["FIRST"]}, istopzone=1, exit_solve_tolerance=DEFAULT_EXIT_SOLVE_TOL, + extend_tracking=True, ) # create output control package diff --git a/autotest/test_prt_ternary_methods.py b/autotest/test_prt_ternary_methods.py index e496628fe00..9a9f706420d 100644 --- a/autotest/test_prt_ternary_methods.py +++ b/autotest/test_prt_ternary_methods.py @@ -100,6 +100,7 @@ def build_prt_sim( stop_at_weak_sink=True, # currently required for this problem dev_exit_solve_method=methods[idx], exit_solve_tolerance=exit_solve_tolerance, + extend_tracking=True, ) prt_track_file = f"{prtname}.trk" prt_track_csv_file = f"{prtname}.trk.csv" diff --git a/autotest/test_prt_track_events.py b/autotest/test_prt_track_events.py index 5912c7e679c..3f63d43cc35 100644 --- a/autotest/test_prt_track_events.py +++ b/autotest/test_prt_track_events.py @@ -143,6 +143,7 @@ def build_prt_sim(name, gwf_ws, prt_ws, mf6): packagedata=releasepts_prt[grp], perioddata={0: ["FIRST"]}, exit_solve_tolerance=DEFAULT_EXIT_SOLVE_TOL, + extend_tracking=True, ) for grp in ["a", "b"] ] diff --git a/autotest/test_prt_triangle.py b/autotest/test_prt_triangle.py index 89c3f8fcc3a..619d058cac8 100644 --- a/autotest/test_prt_triangle.py +++ b/autotest/test_prt_triangle.py @@ -185,6 +185,7 @@ def build_prt_sim(idx, name, gwf_ws, prt_ws, targets): boundnames=True, stop_at_weak_sink=True, # currently required for this problem exit_solve_tolerance=1e-5, + extend_tracking=True, ) prt_track_file = f"{prtname}.trk" prt_track_csv_file = f"{prtname}.trk.csv" diff --git a/autotest/test_prt_voronoi1.py b/autotest/test_prt_voronoi1.py index b3249abd750..800fe59a588 100644 --- a/autotest/test_prt_voronoi1.py +++ b/autotest/test_prt_voronoi1.py @@ -554,6 +554,7 @@ def build_prt_sim(idx, name, gwf_ws, prt_ws, targets, cell_ids): boundnames=True, stop_at_weak_sink=True, exit_solve_tolerance=1e-10, + extend_tracking=True, ) prt_track_file = f"{prt_name}.trk" prt_track_csv_file = f"{prt_name}.trk.csv" diff --git a/autotest/test_prt_voronoi2.py b/autotest/test_prt_voronoi2.py index 08c77babe20..59eca88d59e 100644 --- a/autotest/test_prt_voronoi2.py +++ b/autotest/test_prt_voronoi2.py @@ -173,6 +173,7 @@ def build_prt_sim(name, gwf_ws, prt_ws, targets, cell_ids): boundnames=True, stop_at_weak_sink=True, # currently required for this problem exit_solve_tolerance=1e-10, + extend_tracking=True, ) prt_track_file = f"{prt_name}.trk" prt_track_csv_file = f"{prt_name}.trk.csv" diff --git a/autotest/test_prt_weak_sinks.py b/autotest/test_prt_weak_sinks.py index 25168dedc3f..2d9b3ab0f07 100644 --- a/autotest/test_prt_weak_sinks.py +++ b/autotest/test_prt_weak_sinks.py @@ -101,6 +101,7 @@ def build_prt_sim(name, gwf_ws, prt_ws, mf6): perioddata={0: ["FIRST"]}, stop_at_weak_sink="saws" in name, exit_solve_tolerance=DEFAULT_EXIT_SOLVE_TOL, + extend_tracking=True, ) # create output control package diff --git a/doc/ReleaseNotes/develop.tex b/doc/ReleaseNotes/develop.tex index 4d310021a29..42ccf616372 100644 --- a/doc/ReleaseNotes/develop.tex +++ b/doc/ReleaseNotes/develop.tex @@ -24,6 +24,7 @@ \item The PRT model could hang upon particle termination due to no (sub)cell exit face. This occurred because a flag signaling particle advance was not set in the proper location. \item Entangled with the previous issue, the ternary method could erroneously terminate a particle and report no exit face (before now this would hang) due to precision error in the exit time/position calculation. This could happen when two conditions are both met: the particle enters the subcell very close to one of its vertices, and flow very nearly parallels one of the subcell's faces. We have encountered similar situations before, solved by nudging the particle a small distance into the interior of the subcell before applying the tracking method. This particular case is resolved by increasing the padding distance from $\sisetup{input-digits = 0123456789\epsilon} \num{\epsilon e+2}$ to $\sisetup{input-digits = 0123456789\epsilon} \num{\epsilon e+5}$, where $\epsilon$ is machine precision. \item For ASCII input files erroneously containing a mix of line endings, MODFLOW would sometimes proceed with unexpected results. The program was corrected to stop with an error message if an input line contained both carriage returns and line feeds. + \item Previously the PRT model's default behavior was to track particles until termination, as with MODPATH 7 stop time option 2 (extend). Under extended tracking, the program may not halt if particles enter a cycle in the flow system. This PR changes the default to the equivalent of MP7 stop time option 1 (final), terminating at simulation end unless a new Particle Release Point (PRP) package keyword option EXTEND\_TRACKING is provided. This is meant to provide a stronger guarantee that the program halts under default settings. % \item xxx % \item xxx \end{itemize} diff --git a/doc/mf6io/mf6ivar/dfn/prt-prp.dfn b/doc/mf6io/mf6ivar/dfn/prt-prp.dfn index 5385a42966a..7517d1c98d5 100644 --- a/doc/mf6io/mf6ivar/dfn/prt-prp.dfn +++ b/doc/mf6io/mf6ivar/dfn/prt-prp.dfn @@ -42,6 +42,14 @@ optional true longname whether to use local z coordinates description indicates that ``zrpt'' defines the local z coordinate of the release point within the cell, with value of 0 at the bottom and 1 at the top of the cell. If the cell is partially saturated at release time, the top of the cell is considered to be the water table elevation (the head in the cell) rather than the top defined by the user. +block options +name extend_tracking +type keyword +reader urword +optional true +longname whether to extend tracking beyond the end of the simulation +description indicates that particles should be tracked beyond the end of the simulation's final time step (using that time step's flows) until particles terminate or reach a specified stop time. By default, particles are terminated at the end of the simulation's final time step. + block options name track_filerecord type record track fileout trackfile diff --git a/doc/mf6io/mf6ivar/md/mf6ivar.md b/doc/mf6io/mf6ivar/md/mf6ivar.md index 8e45432fc93..f6b3a094c40 100644 --- a/doc/mf6io/mf6ivar/md/mf6ivar.md +++ b/doc/mf6io/mf6ivar/md/mf6ivar.md @@ -1619,6 +1619,7 @@ | PRT | PRP | OPTIONS | DEV_EXIT_SOLVE_METHOD | INTEGER | the method for iterative solution of particle exit location and time in the generalized Pollock's method. 1 Brent, 2 Chandrupatla. The default is Brent. | | PRT | PRP | OPTIONS | EXIT_SOLVE_TOLERANCE | DOUBLE PRECISION | the convergence tolerance for iterative solution of particle exit location and time in the generalized Pollock's method. A value of 0.00001 works well for many problems, but the value that strikes the best balance between accuracy and runtime is problem-dependent. | | PRT | PRP | OPTIONS | LOCAL_Z | KEYWORD | indicates that ``zrpt'' defines the local z coordinate of the release point within the cell, with value of 0 at the bottom and 1 at the top of the cell. If the cell is partially saturated at release time, the top of the cell is considered to be the water table elevation (the head in the cell) rather than the top defined by the user. | +| PRT | PRP | OPTIONS | EXTEND_TRACKING | KEYWORD | indicates that particles should be tracked beyond the end of the simulation's final time step (using that time step's flows) until particles terminate or reach a specified stop time. By default, particles are terminated at the end of the simulation's final time step. | | PRT | PRP | OPTIONS | TRACK | KEYWORD | keyword to specify that record corresponds to a binary track output file. Each PRP Package may have a distinct binary track output file. | | PRT | PRP | OPTIONS | FILEOUT | KEYWORD | keyword to specify that an output filename is expected next. | | PRT | PRP | OPTIONS | TRACKFILE | STRING | name of the binary output file to write tracking information. | diff --git a/doc/mf6io/mf6ivar/tex/prt-prp-desc.tex b/doc/mf6io/mf6ivar/tex/prt-prp-desc.tex index 830577fbfaf..11ae436ce4e 100644 --- a/doc/mf6io/mf6ivar/tex/prt-prp-desc.tex +++ b/doc/mf6io/mf6ivar/tex/prt-prp-desc.tex @@ -11,6 +11,8 @@ \item \texttt{LOCAL\_Z}---indicates that ``zrpt'' defines the local z coordinate of the release point within the cell, with value of 0 at the bottom and 1 at the top of the cell. If the cell is partially saturated at release time, the top of the cell is considered to be the water table elevation (the head in the cell) rather than the top defined by the user. +\item \texttt{EXTEND\_TRACKING}---indicates that particles should be tracked beyond the end of the simulation's final time step (using that time step's flows) until particles terminate or reach a specified stop time. By default, particles are terminated at the end of the simulation's final time step. + \item \texttt{TRACK}---keyword to specify that record corresponds to a binary track output file. Each PRP Package may have a distinct binary track output file. \item \texttt{FILEOUT}---keyword to specify that an output filename is expected next. diff --git a/doc/mf6io/mf6ivar/tex/prt-prp-options.dat b/doc/mf6io/mf6ivar/tex/prt-prp-options.dat index dc663bb3868..9305b16cd89 100644 --- a/doc/mf6io/mf6ivar/tex/prt-prp-options.dat +++ b/doc/mf6io/mf6ivar/tex/prt-prp-options.dat @@ -3,6 +3,7 @@ BEGIN OPTIONS [PRINT_INPUT] EXIT_SOLVE_TOLERANCE [LOCAL_Z] + [EXTEND_TRACKING] [TRACK FILEOUT ] [TRACKCSV FILEOUT ] [STOPTIME ] diff --git a/doc/mf6io/prt/prt.tex b/doc/mf6io/prt/prt.tex index 63a9a756667..196eba092c6 100644 --- a/doc/mf6io/prt/prt.tex +++ b/doc/mf6io/prt/prt.tex @@ -12,7 +12,7 @@ \subsection{Units of Length and Time} The PRT Model formulates the particle trajectory equations without using prescribed length and time units. Any consistent units of length and time can be used when specifying the input data for a simulation. This capability gives a certain amount of freedom to the user, but care must be exercised to avoid mixing units. The program cannot detect the use of inconsistent units. \subsection{Time Stepping} -In \mf time step lengths are controlled by the user and specified in the Temporal Discretization (TDIS) input file. When the flow model and particle-tracking model run in the same simulation, the time step length specified in TDIS is used for both models. If the PRT Model runs in a separate simulation, the time discretization may differ. Instructions for specifying time steps are described in the TDIS section of this user guide; additional information on GWF and PRT configurations are in the Flow Model Interface section. +In \mf time step lengths are controlled by the user and specified in the Temporal Discretization (TDIS) input file. When the flow model and particle-tracking model run in the same simulation, the time step length specified in TDIS is used for both models. If the PRT Model runs in a separate simulation, the time discretization may differ. Instructions for specifying time steps are described in the TDIS section of this user guide; additional information on GWF and PRT configurations are in the Flow Model Interface section. By default, the PRT model will terminate particles at the end of the simulation's final time step; PRT can be configured to extend tracking until particles terminate naturally (i.e. at boundary conditions, sinks, or specified stop times) via Particle Release Point (PRP) package settings. \subsection{Specifying Cell Face Flows using IFLOWFACE} diff --git a/src/Model/ParticleTracking/prt-prp.f90 b/src/Model/ParticleTracking/prt-prp.f90 index 558788f647c..2365b71d8c8 100644 --- a/src/Model/ParticleTracking/prt-prp.f90 +++ b/src/Model/ParticleTracking/prt-prp.f90 @@ -47,6 +47,7 @@ module PrtPrpModule integer(I4B), pointer :: itrkcsv => null() !< CSV track file integer(I4B), pointer :: irlstls => null() !< release time file logical(LGP), pointer :: localz => null() !< compute z coordinates local to the cell + integer(I4B), pointer :: extend => null() !< extend tracking beyond simulation's end logical(LGP), pointer :: rlsall => null() !< release in all time step logical(LGP), pointer :: rlsfirst => null() !< release in first time step logical(LGP), pointer :: rlstimelist => null() !< use global release time @@ -146,6 +147,7 @@ subroutine prp_da(this) call mem_deallocate(this%rlsfirst) call mem_deallocate(this%rlstimelist) call mem_deallocate(this%localz) + call mem_deallocate(this%extend) call mem_deallocate(this%offset) call mem_deallocate(this%stoptime) call mem_deallocate(this%stoptraveltime) @@ -246,6 +248,7 @@ subroutine prp_allocate_scalars(this) call mem_allocate(this%rlsfirst, 'RLSFIRST', this%memoryPath) call mem_allocate(this%rlstimelist, 'RELEASETIME', this%memoryPath) call mem_allocate(this%localz, 'LOCALZ', this%memoryPath) + call mem_allocate(this%extend, 'EXTEND', this%memoryPath) call mem_allocate(this%offset, 'OFFSET', this%memoryPath) call mem_allocate(this%stoptime, 'STOPTIME', this%memoryPath) call mem_allocate(this%stoptraveltime, 'STOPTRAVELTIME', this%memoryPath) @@ -268,6 +271,7 @@ subroutine prp_allocate_scalars(this) this%rlsfirst = .false. this%rlstimelist = .false. this%localz = .false. + this%extend = 0 this%offset = DZERO this%stoptime = huge(1d0) this%stoptraveltime = huge(1d0) @@ -471,6 +475,7 @@ subroutine prp_ad(this) particle%ifrctrn = this%ifrctrn particle%iexmethod = this%iexmethod particle%extol = this%extol + particle%extend = this%extend ! -- Persist particle to particle store call this%particles%save_particle(particle, np) @@ -828,6 +833,9 @@ subroutine prp_options(this, option, found) case ('LOCAL_Z') this%localz = .true. found = .true. + case ('EXTEND_TRACKING') + this%extend = 1 + found = .true. case ('DEV_FORCETERNARY') call this%parser%DevOpt() this%ifrctrn = 1 diff --git a/src/Model/ParticleTracking/prt.f90 b/src/Model/ParticleTracking/prt.f90 index e0c401b1a18..5059d9c66bf 100644 --- a/src/Model/ParticleTracking/prt.f90 +++ b/src/Model/ParticleTracking/prt.f90 @@ -933,25 +933,30 @@ subroutine prt_solve(this) if (particle%istatus == 8) & call this%method%save(particle, reason=3) ! reason=3: termination - ! -- If particle is inactive or not yet to be released, cycle + ! If particle is inactive or not yet to be released, cycle if (particle%istatus > 1) cycle - ! -- If particle released this time step, record its initial state + ! If particle released this time step, record its initial state particle%istatus = 1 if (particle%trelease >= totimc) & call this%method%save(particle, reason=0) ! reason=0: release - ! -- Maximum time is end of time step unless this is the last - ! time step in the simulation, which case it's the particle - ! stop time. - tmax = totimc + delt - if (nper == kper .and. nstp(kper) == kstp) & + ! Maximum time is the end of the time step or the particle + ! stop time, whichever comes first, unless it's the final + ! time step and the extend option is on, in which case + ! it's just the particle stop time. + if (nper == kper .and. & + nstp(kper) == kstp .and. & + particle%extend > 0) then tmax = particle%tstop + else + tmax = min(totimc + delt, particle%tstop) + end if - ! -- Get and apply the tracking method + ! Get and apply the tracking method call this%method%apply(particle, tmax) - ! -- Update particle storage + ! Update particle storage call packobj%particles%save_particle(particle, np) end do end select diff --git a/src/Solution/ParticleTracker/MethodSubcellPollock.f90 b/src/Solution/ParticleTracker/MethodSubcellPollock.f90 index 20ccd39f6c3..dcb7729597f 100644 --- a/src/Solution/ParticleTracker/MethodSubcellPollock.f90 +++ b/src/Solution/ParticleTracker/MethodSubcellPollock.f90 @@ -181,15 +181,16 @@ subroutine track_subcell(this, subcell, particle, tmax) ! -- 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. + ! last time step include all times after it ends only if the 'extend' + ! option is on, otherwise times in this 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 .or. t >= texit .or. t >= tmax) cycle + if (t < particle%ttrack) cycle + if (t >= texit .or. t >= tmax) exit dt = t - t0 x = new_x(vx, dvxdx, subcell%vx1, subcell%vx2, & dt, initialX, subcell%dx, statusVX == 1) diff --git a/src/Solution/ParticleTracker/MethodSubcellTernary.f90 b/src/Solution/ParticleTracker/MethodSubcellTernary.f90 index 09f51c91e4d..c8717519d47 100644 --- a/src/Solution/ParticleTracker/MethodSubcellTernary.f90 +++ b/src/Solution/ParticleTracker/MethodSubcellTernary.f90 @@ -226,16 +226,17 @@ subroutine track_subcell(this, subcell, particle, tmax) texit = particle%ttrack + dtexit t0 = particle%ttrack - ! Solve user-specified tracking times within the current stress period and - ! time step. If this is the last time step, solve all times after it ends, - ! as per MODPATH 7 with stop time option 'extend'. - ! todo: reconsider whether this should be default? + ! -- 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 only if the 'extend' + ! option is on, otherwise times in this 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 .or. t >= texit .or. t >= tmax) cycle + if (t < particle%ttrack) cycle + if (t >= texit .or. t >= tmax) exit dt = t - t0 call calculate_xyz_position(dt, rxx, rxy, ryx, ryy, sxx, sxy, syy, & izstatus, x0, y0, az, vzi, vzbot, & diff --git a/src/Solution/ParticleTracker/Particle.f90 b/src/Solution/ParticleTracker/Particle.f90 index fb35c29f296..21d8c0e8b25 100644 --- a/src/Solution/ParticleTracker/Particle.f90 +++ b/src/Solution/ParticleTracker/Particle.f90 @@ -62,6 +62,7 @@ module ParticleModule integer(I4B), public :: ifrctrn !< whether to force solving the particle with the ternary method integer(I4B), public :: iexmethod !< method for iterative solution of particle exit location and time in generalized Pollock's method real(DP), public :: extol !< tolerance for iterative solution of particle exit location and time in generalized Pollock's method + integer(I4B), public :: extend !< whether to extend tracking beyond the end of the simulation contains procedure, public :: get_model_coords procedure, public :: load_particle @@ -70,30 +71,32 @@ module ParticleModule !> @brief Structure of arrays to store particles. type ParticleStoreType + private ! identity - character(len=LENBOUNDNAME), dimension(:), pointer, contiguous :: name !< optional particle label - integer(I4B), dimension(:), pointer, contiguous :: imdl !< index of model particle originated in - integer(I4B), dimension(:), pointer, contiguous :: iprp !< index of release package the particle originated in - integer(I4B), dimension(:), pointer, contiguous :: irpt !< index of release point in the particle release package the particle originated in + character(len=LENBOUNDNAME), dimension(:), pointer, public, contiguous :: name !< optional particle label + integer(I4B), dimension(:), pointer, public, contiguous :: imdl !< index of model particle originated in + integer(I4B), dimension(:), pointer, public, contiguous :: iprp !< index of release package the particle originated in + integer(I4B), dimension(:), pointer, public, contiguous :: irpt !< index of release point in the particle release package the particle originated in ! stopping criteria - integer(I4B), dimension(:), pointer, contiguous :: istopweaksink !< weak sink option: 0 = do not stop, 1 = stop - integer(I4B), dimension(:), pointer, contiguous :: istopzone !< stop zone number + integer(I4B), dimension(:), pointer, public, contiguous :: istopweaksink !< weak sink option: 0 = do not stop, 1 = stop + integer(I4B), dimension(:), pointer, public, contiguous :: istopzone !< stop zone number ! state - integer(I4B), dimension(:, :), pointer, contiguous :: idomain !< array of indices for domains in the tracking domain hierarchy - integer(I4B), dimension(:, :), pointer, contiguous :: iboundary !< array of indices for tracking domain boundaries - integer(I4B), dimension(:), pointer, contiguous :: icu !< cell number (user, not reduced) - integer(I4B), dimension(:), pointer, contiguous :: ilay !< layer - integer(I4B), dimension(:), pointer, contiguous :: izone !< current zone number - integer(I4B), dimension(:), pointer, contiguous :: istatus !< particle status - real(DP), dimension(:), pointer, contiguous :: x !< model x coord of particle - real(DP), dimension(:), pointer, contiguous :: y !< model y coord of particle - real(DP), dimension(:), pointer, contiguous :: z !< model z coord of particle - real(DP), dimension(:), pointer, contiguous :: trelease !< particle release time - real(DP), dimension(:), pointer, contiguous :: tstop !< particle stop time - real(DP), dimension(:), pointer, contiguous :: ttrack !< current tracking time - integer(I4B), dimension(:), pointer, contiguous :: ifrctrn !< force ternary method - integer(I4B), dimension(:), pointer, contiguous :: iexmethod !< method for iterative solution of particle exit location and time in generalized Pollock's method - real(DP), dimension(:), pointer, contiguous :: extol !< tolerance for iterative solution of particle exit location and time in generalized Pollock's method + integer(I4B), dimension(:, :), pointer, public, contiguous :: idomain !< array of indices for domains in the tracking domain hierarchy + integer(I4B), dimension(:, :), pointer, public, contiguous :: iboundary !< array of indices for tracking domain boundaries + integer(I4B), dimension(:), pointer, public, contiguous :: icu !< cell number (user, not reduced) + integer(I4B), dimension(:), pointer, public, contiguous :: ilay !< layer + integer(I4B), dimension(:), pointer, public, contiguous :: izone !< current zone number + integer(I4B), dimension(:), pointer, public, contiguous :: istatus !< particle status + real(DP), dimension(:), pointer, public, contiguous :: x !< model x coord of particle + real(DP), dimension(:), pointer, public, contiguous :: y !< model y coord of particle + real(DP), dimension(:), pointer, public, contiguous :: z !< model z coord of particle + real(DP), dimension(:), pointer, public, contiguous :: trelease !< particle release time + real(DP), dimension(:), pointer, public, contiguous :: tstop !< particle stop time + real(DP), dimension(:), pointer, public, contiguous :: ttrack !< current tracking time + integer(I4B), dimension(:), pointer, public, contiguous :: ifrctrn !< force ternary method + integer(I4B), dimension(:), pointer, public, contiguous :: iexmethod !< method for iterative solution of particle exit location and time in generalized Pollock's method + real(DP), dimension(:), pointer, public, contiguous :: extol !< tolerance for iterative solution of particle exit location and time in generalized Pollock's method + integer(LGP), dimension(:), pointer, public, contiguous :: extend !< whether to extend tracking beyond the end of the simulation contains procedure, public :: deallocate procedure, public :: resize @@ -136,6 +139,7 @@ subroutine allocate_particle_store(this, np, mempath) call mem_allocate(this%ifrctrn, np, 'PLIFRCTRN', mempath) call mem_allocate(this%iexmethod, np, 'PLIEXMETHOD', mempath) call mem_allocate(this%extol, np, 'PLEXTOL', mempath) + call mem_allocate(this%extend, np, 'PLEXTEND', mempath) call mem_allocate(this%idomain, np, levelmax, 'PLIDOMAIN', mempath) call mem_allocate(this%iboundary, np, levelmax, 'PLIBOUNDARY', mempath) end subroutine allocate_particle_store @@ -164,6 +168,7 @@ subroutine deallocate (this, mempath) call mem_deallocate(this%ifrctrn, 'PLIFRCTRN', mempath) call mem_deallocate(this%iexmethod, 'PLIEXMETHOD', mempath) call mem_deallocate(this%extol, 'PLEXTOL', mempath) + call mem_deallocate(this%extend, 'PLEXTEND', mempath) call mem_deallocate(this%idomain, 'PLIDOMAIN', mempath) call mem_deallocate(this%iboundary, 'PLIBOUNDARY', mempath) end subroutine deallocate @@ -195,6 +200,7 @@ subroutine resize(this, np, mempath) call mem_reallocate(this%ifrctrn, np, 'PLIFRCTRN', mempath) call mem_reallocate(this%iexmethod, np, 'PLIEXMETHOD', mempath) call mem_reallocate(this%extol, np, 'PLEXTOL', mempath) + call mem_reallocate(this%extend, np, 'PLEXTEND', mempath) call mem_reallocate(this%idomain, np, levelmax, 'PLIDOMAIN', mempath) call mem_reallocate(this%iboundary, np, levelmax, 'PLIBOUNDARY', mempath) end subroutine resize @@ -238,6 +244,7 @@ subroutine load_particle(this, store, imdl, iprp, ip) this%ifrctrn = store%ifrctrn(ip) this%iexmethod = store%iexmethod(ip) this%extol = store%extol(ip) + this%extend = store%extend(ip) end subroutine load_particle !> @brief Save a particle's state to the particle store @@ -270,9 +277,10 @@ subroutine save_particle(this, particle, ip) ip, & 1:levelmax) = & particle%iboundary(1:levelmax) - this%ifrctrn = particle%ifrctrn - this%iexmethod = particle%iexmethod - this%extol = particle%extol + this%ifrctrn(ip) = particle%ifrctrn + this%iexmethod(ip) = particle%iexmethod + this%extol(ip) = particle%extol + this%extend(ip) = particle%extend end subroutine save_particle !> @brief Apply the given global-to-local transformation to the particle.