refactor(prt): make exit_solve_tolerance optional, update release not…

…es (#2107)

Mention the EXIT_SOLVE_TOLERANCE option now in the PRP package in the release notes, I missed this before. This was previously a required option, but need not be — set a default of 1e-5, and simplify some tests in light of this.

Also draft a PRT migration guide, building on the development notes added in #2066, which could be distributed with the release. I'm not sure whether this is something to version, and if so, where it should go, but figured it can't hurt to have it.

autotest/test_prt_budget.py

@@ -19,7 +19,6 @@
 from flopy.utils.binaryfile import HeadFile
 from framework import TestFramework
 from prt_test_utils import (
-    DEFAULT_EXIT_SOLVE_TOL,
     HorizontalCase,
     all_equal,
     check_budget_data,
@@ -88,7 +87,6 @@ def build_prt_sim(name, gwf_ws, prt_ws, mf6):
         trackcsv_filerecord=[prp_track_csv_file],
         stop_at_weak_sink=False,
         boundnames=True,
-        exit_solve_tolerance=DEFAULT_EXIT_SOLVE_TOL,
         extend_tracking=True,
     )
 

autotest/test_prt_drape.py

@@ -20,7 +20,6 @@
 from flopy.utils.binaryfile import HeadFile
 from framework import TestFramework
 from prt_test_utils import (
-    DEFAULT_EXIT_SOLVE_TOL,
     all_equal,
     check_track_data,
     get_model_name,
@@ -182,7 +181,6 @@ def build_prt_sim(name, gwf_ws, prt_ws, mf6):
         track_filerecord=[prp_track_file],
         trackcsv_filerecord=[prp_track_csv_file],
         drape="drp" in name,
-        exit_solve_tolerance=DEFAULT_EXIT_SOLVE_TOL,
         extend_tracking=True,
     )
 

autotest/test_prt_dry.py

@@ -301,7 +301,6 @@ def build_prt_sim(name, gwf, prt_ws, mf6, drape=False, dry_tracking_method=False
         packagedata=prp_data,
         nreleasetimes=1,
         releasetimes=[(0.0,)],
-        exit_solve_tolerance=1e-7,
         drape=drape,
         dry_tracking_method=dry_tracking_method,
         pname="prp",

autotest/test_prt_exg.py

@@ -24,7 +24,6 @@
 from flopy.utils.binaryfile import HeadFile
 from framework import TestFramework
 from prt_test_utils import (
-    DEFAULT_EXIT_SOLVE_TOL,
     FlopyReadmeCase,
     check_budget_data,
     check_track_data,
@@ -73,7 +72,6 @@ def build_mf6_sim(idx, test):
         packagedata=rpts,
         perioddata={0: ["FIRST"]},
         boundnames="bnms" in name,
-        exit_solve_tolerance=DEFAULT_EXIT_SOLVE_TOL,
         extend_tracking=True,
     )
 

autotest/test_prt_fmi.py

@@ -38,7 +38,6 @@
 from flopy.utils.binaryfile import HeadFile
 from framework import TestFramework
 from prt_test_utils import (
-    DEFAULT_EXIT_SOLVE_TOL,
     FlopyReadmeCase,
     all_equal,
     check_budget_data,
@@ -116,7 +115,6 @@ def build_prt_sim(name, gwf_ws, prt_ws, mf6):
         trackcsv_filerecord=[prp_track_csv_file],
         stop_at_weak_sink="saws" in prt_name,
         boundnames=True,
-        exit_solve_tolerance=DEFAULT_EXIT_SOLVE_TOL,
         extend_tracking="noext" not in prt_name,
     )
 

autotest/test_prt_quad_refinement.py

@@ -18,7 +18,6 @@
 from flopy.utils.gridgen import Gridgen
 from framework import TestFramework
 from prt_test_utils import (
-    DEFAULT_EXIT_SOLVE_TOL,
     FlopyReadmeCase,
     check_budget_data,
     check_track_data,
@@ -144,7 +143,6 @@ def build_mf6_sim(idx, test, **kwargs):
         nreleasepts=len(rpts),
         packagedata=rpts,
         perioddata={0: ["FIRST"]},
-        exit_solve_tolerance=DEFAULT_EXIT_SOLVE_TOL,
         dev_forceternary=tracking_method == "ternary",
         extend_tracking=True,
     )

autotest/test_prt_release_timing.py

@@ -24,7 +24,6 @@
 from framework import TestFramework
 from modflow_devtools.markers import requires_pkg
 from prt_test_utils import (
-    DEFAULT_EXIT_SOLVE_TOL,
     FlopyReadmeCase,
     all_equal,
     check_budget_data,
@@ -176,7 +175,6 @@ def build_prt_sim(name, gwf_ws, prt_ws, mf6):
         else None,
         release_time_frequency=0.2 if "freq" in name else None,
         print_input=True,
-        exit_solve_tolerance=DEFAULT_EXIT_SOLVE_TOL,
         extend_tracking=True,
         release_time_tolerance=0.2 if "tol" in name else None,
     )

autotest/test_prt_stop_zones.py

@@ -35,7 +35,6 @@
 from framework import TestFramework
 from matplotlib.collections import LineCollection
 from prt_test_utils import (
-    DEFAULT_EXIT_SOLVE_TOL,
     FlopyReadmeCase,
     check_budget_data,
     check_track_data,
@@ -127,7 +126,6 @@ def build_prt_sim(name, gwf_ws, prt_ws, mf6):
         packagedata=FlopyReadmeCase.releasepts_prt,
         perioddata={0: ["FIRST"]},
         istopzone=1,
-        exit_solve_tolerance=DEFAULT_EXIT_SOLVE_TOL,
         extend_tracking=True,
     )
 

autotest/test_prt_track_events.py

@@ -34,7 +34,6 @@
 from flopy.utils.binaryfile import HeadFile
 from framework import TestFramework
 from prt_test_utils import (
-    DEFAULT_EXIT_SOLVE_TOL,
     FlopyReadmeCase,
     check_budget_data,
     check_track_data,
@@ -140,7 +139,6 @@ def build_prt_sim(name, gwf_ws, prt_ws, mf6):
             nreleasepts=len(releasepts_prt[grp]),
             packagedata=releasepts_prt[grp],
             perioddata={0: ["FIRST"]},
-            exit_solve_tolerance=DEFAULT_EXIT_SOLVE_TOL,
             extend_tracking=True,
         )
         for grp in ["a", "b"]

autotest/test_prt_triangle.py

@@ -177,7 +177,6 @@ def build_prt_sim(idx, name, gwf_ws, prt_ws, targets):
             perioddata={0: ["FIRST"]},
             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"

autotest/test_prt_weak_sinks.py

@@ -34,7 +34,6 @@
 from flopy.utils.binaryfile import HeadFile
 from framework import TestFramework
 from prt_test_utils import (
-    DEFAULT_EXIT_SOLVE_TOL,
     FlopyReadmeCase,
     check_budget_data,
     check_track_data,
@@ -97,7 +96,6 @@ def build_prt_sim(name, gwf_ws, prt_ws, mf6):
         packagedata=FlopyReadmeCase.releasepts_prt,
         perioddata={0: ["FIRST"]},
         stop_at_weak_sink="saws" in name,
-        exit_solve_tolerance=DEFAULT_EXIT_SOLVE_TOL,
         extend_tracking=True,
     )
 

doc/ReleaseNotes/develop.tex

@@ -54,6 +54,7 @@
 		\item With a flow model using the Newton formulation, the PRT model could crash upon a particle's entry into a dry cell. This has been fixed.
 		\item With a flow model using the Newton formulation, the PRT model could enter an endless loop upon a particle's entry to a dry cell if that cell contains a boundary package (e.g. a pumping well). Where the particle should be captured and terminate, it would instead be passed back and forth between the cell bottom and the top of the cell below. To avoid this, particles are forbidden from backtracking (reentering the previous cell) within the same time step.
 		\item The PRT model now allows more control over vertical particle motion in dry conditions. In addition to the existing DRAPE option, which controls release-time behavior, the PRP package now provides a DRY\_TRACKING\_METHOD option which configures how dry particles (particles in dry cells, or above the water table in partially saturated cells) behave at tracking time. This option is relevant only when the Newton formulation is used, in which case dry cells remain active; otherwise, dry cells are inactive and particles will terminate. See the MF6IO document for a detailed explanation of DRY\_TRACKING\_METHOD.
+		\item The PRT model's Particle Release Point (PRP) package now provides an option EXIT\_SOLVE\_TOLERANCE which configures the tolerance to use when solving for a particle's exit location from an unstructured grid cell. This value is only used for the generalized (ternary) tracking method on vertex grids. A value of 0.00001 is set by default. This value works well for many problems, but the value that strikes the best balance between accuracy and runtime is problem-dependent.
 	\end{itemize}
 
 	%\underline{INTERNAL FLOW PACKAGES}

doc/mf6io/mf6ivar/dfn/prt-prp.dfn

@@ -30,9 +30,10 @@ block options
 name exit_solve_tolerance
 type double precision
 reader urword
-optional false
+optional true
 longname exit solve tolerance
 description 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.
+default_value 1e-5
 
 block options
 name local_z

src/Model/ParticleTracking/prt-prp.f90

@@ -1,6 +1,6 @@
 module PrtPrpModule
   use KindModule, only: DP, I4B, LGP
-  use ConstantsModule, only: DZERO, DEM1, DONE, LENFTYPE, LINELENGTH, &
+  use ConstantsModule, only: DZERO, DEM1, DEM5, DONE, LENFTYPE, LINELENGTH, &
                              LENBOUNDNAME, LENPAKLOC, TABLEFT, TABCENTER, &
                              MNORMAL, DSAME, DEP3, DEP9
   use BndModule, only: BndType
@@ -64,7 +64,6 @@ module PrtPrpModule
     real(DP), pointer :: offset => null() !< release time offset
     real(DP), pointer :: stoptime => null() !< stop time for all release points
     real(DP), pointer :: stoptraveltime => null() !< stop travel time for all points
-    logical(LGP), pointer :: foundtol => null() !< whether tolerance option was found
     integer(I4B), pointer, contiguous :: rptnode(:) => null() !< release point reduced nns
     integer(I4B), pointer, contiguous :: rptzone(:) => null() !< release point zone numbers
     real(DP), pointer, contiguous :: rptx(:) => null() !< release point x coordinates
@@ -173,7 +172,6 @@ subroutine prp_da(this)
     call mem_deallocate(this%extol)
     call mem_deallocate(this%rttol)
     call mem_deallocate(this%rtfreq)
-    call mem_deallocate(this%foundtol)
 
     ! Deallocate arrays
     call mem_deallocate(this%rptx)
@@ -263,7 +261,6 @@ subroutine prp_allocate_scalars(this)
     call mem_allocate(this%extol, 'EXTOL', this%memoryPath)
     call mem_allocate(this%rttol, 'RTTOL', this%memoryPath)
     call mem_allocate(this%rtfreq, 'RTFREQ', this%memoryPath)
-    call mem_allocate(this%foundtol, 'FOUNDTOL', this%memoryPath)
 
     ! Set values
     this%ilocalz = 0
@@ -284,10 +281,9 @@ subroutine prp_allocate_scalars(this)
     this%irlstls = 0
     this%ifrctrn = 0
     this%iexmeth = 0
-    this%extol = DZERO
+    this%extol = DEM5
     this%rttol = DSAME * DEP9
     this%rtfreq = DZERO
-    this%foundtol = .false.
 
   end subroutine prp_allocate_scalars
 
@@ -782,7 +778,6 @@ subroutine prp_options(this, option, found)
       if (this%extol <= DZERO) &
         call store_error('EXIT_SOLVE_TOLERANCE MUST BE POSITIVE')
       found = .true.
-      this%foundtol = .true.
     case ('RELEASE_TIME_TOLERANCE')
       this%rttol = this%parser%GetDouble()
       if (this%rttol <= DZERO) &
@@ -832,9 +827,6 @@ subroutine prp_read_dimensions(this)
     integer(I4B) :: ierr
     logical :: isfound, endOfBlock
 
-    if (.not. this%foundtol) &
-      call store_error('EXIT_SOLVE_TOLERANCE MISSING, VALUE REQUIRED')
-
     ! get dimension block
     call this%parser%GetBlock('DIMENSIONS', isfound, ierr, &
                               supportOpenClose=.true.)