refactor(prt-prp): require exit_solve_tolerance option (MODFLOW-USGS#…

…1816)

* make EXIT_SOLVE_TOLERANCE required for DISV grids for all grids, motivation: one should be aware the generalized method is semi-analytical and subcells are solved numerically for grid cells we can't solve with Pollock's method
* accept any value for DIS grids, following convention in e.g. rnb in gwf-csub.dfn packagedata where a sometimes-required option can have any value when not required
* fix double comma in bib file

autotest/prt_test_utils.py

@@ -1,12 +1,13 @@
 import os
-from types import SimpleNamespace
-from typing import Tuple
 
 import flopy
 import matplotlib as mpl
 import numpy as np
 
 
+DEFAULT_EXIT_SOLVE_TOL = 1e-5
+
+
 def all_equal(series, val):
     a = series.to_numpy()
     return a[0] == val and (a[0] == a).all()

autotest/test_prt_budget.py

@@ -26,6 +26,7 @@
     get_model_name,
     get_partdata,
     has_default_boundnames,
+    DEFAULT_EXIT_SOLVE_TOL,
 )
 
 simname = "prtbud"
@@ -72,7 +73,6 @@ def build_prt_sim(name, gwf_ws, prt_ws, mf6):
     partdata = get_partdata(prt.modelgrid, HorizontalCase.releasepts_mp7)
     coords = partdata.to_coords(prt.modelgrid)
     releasepts = [(i, 0, 0, 0, c[0], c[1], c[2]) for i, c in enumerate(coords)]
-    # assert np.allclose(HorizontalCase.releasepts_prt, releasepts)
 
     # create prp package
     prp_track_file = f"{prt_name}.prp.trk"
@@ -88,6 +88,7 @@ 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,
     )
 
     # create output control package

autotest/test_prt_disv1.py

@@ -235,6 +235,7 @@ def build_prt_sim(idx, gwf_ws, prt_ws, mf6):
             boundnames=True,
             print_input=True,
             dev_forceternary=i == 1,
+            exit_solve_tolerance=1e-10,
         )
 
     # create output control package

autotest/test_prt_drape.py

@@ -19,7 +19,12 @@
 import pytest
 from flopy.utils.binaryfile import HeadFile
 from framework import TestFramework
-from prt_test_utils import all_equal, check_track_data, get_model_name
+from prt_test_utils import (
+    all_equal,
+    check_track_data,
+    get_model_name,
+    DEFAULT_EXIT_SOLVE_TOL,
+)
 
 simname = "prtdrape"
 cases = [simname, f"{simname}_drp"]
@@ -184,6 +189,7 @@ 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,
     )
 
     # create output control package

autotest/test_prt_exg.py

@@ -24,7 +24,12 @@
 from flopy.utils import PathlineFile
 from flopy.utils.binaryfile import HeadFile
 from framework import TestFramework
-from prt_test_utils import FlopyReadmeCase, check_budget_data, check_track_data
+from prt_test_utils import (
+    FlopyReadmeCase,
+    check_budget_data,
+    check_track_data,
+    DEFAULT_EXIT_SOLVE_TOL,
+)
 
 simname = "prtexg01"
 cases = [simname, f"{simname}bnms"]
@@ -73,6 +78,7 @@ def build_mf6_sim(idx, test):
         packagedata=rpts,
         perioddata={0: ["FIRST"]},
         boundnames="bnms" in name,
+        exit_solve_tolerance=DEFAULT_EXIT_SOLVE_TOL,
     )
 
     # create output control package

autotest/test_prt_fmi.py

@@ -45,6 +45,7 @@
     get_model_name,
     get_partdata,
     has_default_boundnames,
+    DEFAULT_EXIT_SOLVE_TOL,
 )
 
 simname = "prtfmi01"
@@ -121,6 +122,7 @@ 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,
     )
 
     # create output control package

autotest/test_prt_release_timing.py

@@ -42,6 +42,7 @@
     check_track_data,
     get_model_name,
     get_partdata,
+    DEFAULT_EXIT_SOLVE_TOL,
 )
 
 pytest_plugins = ["modflow_devtools.snapshots"]
@@ -169,6 +170,7 @@ def build_prt_sim(name, gwf_ws, prt_ws, mf6, fraction=None):
             else None
         ),
         print_input=True,
+        exit_solve_tolerance=DEFAULT_EXIT_SOLVE_TOL,
     )
 
     # create output control package

autotest/test_prt_stop_zones.py

@@ -39,6 +39,7 @@
     check_budget_data,
     check_track_data,
     get_model_name,
+    DEFAULT_EXIT_SOLVE_TOL,
 )
 
 simname = "prtfmi03"
@@ -129,6 +130,7 @@ 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,
     )
 
     # create output control package

autotest/test_prt_ternary_methods.py

@@ -51,7 +51,9 @@
 ]
 
 
-def build_prt_sim(idx, name, gwf_ws, prt_ws, targets, exit_solve_tolerance=1e-5):
+def build_prt_sim(
+    idx, name, gwf_ws, prt_ws, targets, exit_solve_tolerance=1e-5
+):
     prt_ws = Path(prt_ws)
     gwfname = get_model_name(name, "gwf")
     prtname = get_model_name(name, "prt")
@@ -79,9 +81,7 @@ def build_prt_sim(idx, name, gwf_ws, prt_ws, targets, exit_solve_tolerance=1e-5)
         vertices=vertices,
         cell2d=cell2d,
     )
-    flopy.mf6.ModflowPrtmip(
-        prt, pname="mip", porosity=porosity
-    )
+    flopy.mf6.ModflowPrtmip(prt, pname="mip", porosity=porosity)
     prpdata = [
         # particle index, (layer, cell index), x, y, z
         (0, (0, 88), 95, 92, 0.5),
@@ -101,7 +101,7 @@ def build_prt_sim(idx, name, gwf_ws, prt_ws, targets, exit_solve_tolerance=1e-5)
         boundnames=True,
         stop_at_weak_sink=True,  # currently required for this problem
         dev_exit_solve_method=methods[idx],
-        exit_solve_tolerance=exit_solve_tolerance
+        exit_solve_tolerance=exit_solve_tolerance,
     )
     prt_track_file = f"{prtname}.trk"
     prt_track_csv_file = f"{prtname}.trk.csv"
@@ -134,7 +134,12 @@ def build_models(idx, test, exit_solve_tolerance=1e-7):
         test.name, test.workspace, test.targets, ["left", "botm"]
     )
     prt_sim = build_prt_sim(
-        idx, test.name, test.workspace, test.workspace / "prt", test.targets, exit_solve_tolerance
+        idx,
+        test.name,
+        test.workspace,
+        test.workspace / "prt",
+        test.targets,
+        exit_solve_tolerance,
     )
     return gwf_sim, prt_sim
 
@@ -258,7 +263,9 @@ def check_output(idx, test, snapshot):
 
 
 @pytest.mark.parametrize("idx, name", enumerate(cases))
-def test_mf6model(idx, name, function_tmpdir, targets, benchmark, array_snapshot):
+def test_mf6model(
+    idx, name, function_tmpdir, targets, benchmark, array_snapshot
+):
     test = TestFramework(
         name=name,
         workspace=function_tmpdir,

autotest/test_prt_track_events.py

@@ -39,6 +39,7 @@
     check_budget_data,
     check_track_data,
     get_model_name,
+    DEFAULT_EXIT_SOLVE_TOL,
 )
 
 simname = "prtevnt"
@@ -141,6 +142,7 @@ 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,
         )
         for grp in ["a", "b"]
     ]

autotest/test_prt_triangle.py

@@ -186,6 +186,7 @@ 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,
         )
     prt_track_file = f"{prtname}.trk"
     prt_track_csv_file = f"{prtname}.trk.csv"