feat(prt): make release coordinate checks optional (#2424)

wpbonelli · web-flow · commit c94c7a353de9 · 2025-07-30T11:55:34.000-04:00
While it's on the user to provide both release coordinates and cell IDs, for some time we have checked at release time that release coordinates fall within the specified cell. Alternative ways to specify release points are in the works in #1901 e.g. local coordinates, which would remove the redundancy, but so long as release points are in global coordinates and we don't do cell identification ourselves, the coordinate checks are a performance bottleneck when particle count is high. Add a COORDINATE_CHECK_METHOD option to allow opting out. Checks on (by default) with eager, trust me I know what I'm doing mode with none. Later, consider adding a lazy mode where validation happens at tracking time, this might save some time if we can reuse terms of other calculations for the checks, but contradicts "fail early" philosophy so should probably not become default. Need to test and get some performance numbers. If in future we did start computing cell IDs (e.g. with something fast like the celltree algorithm) this option could just be ignored when cell IDs are not provided by the user. So one could decide whether to do it in preprocessing or let PRT handle it.
diff --git a/autotest/test_prt_release_coord_check.py b/autotest/test_prt_release_coord_check.py
@@ -0,0 +1,217 @@
+"""
+Test coordinate checking functionality for PRT release points.
+
+Tests the COORDINATE_CHECK_METHOD option which can be set to:
+- 'eager' (default): check at release time
+- 'none': skip release coordinate validation
+
+Two test cases:
+1. checks on ('eager'): Should catch mismatching points/cellids with informative error
+2. no checks ('none'): Should allow mismatches and let tracking proceed. This produces
+incorrect pathlines where particles jump to the specified cells after reporting initial
+positions at the release coordinates.
+"""
+
+import flopy
+import matplotlib.cm as cm
+import matplotlib.pyplot as plt
+import pandas as pd
+import pytest
+from flopy.utils import HeadFile
+from framework import TestFramework
+from prt_test_utils import FlopyReadmeCase, get_model_name
+
+simname = "prtchk"
+cases = [
+    f"{simname}_eager",
+    f"{simname}_none",
+]
+
+
+def build_prt_sim(name, gwf_ws, prt_ws, mf6):
+    # create simulation
+    sim = flopy.mf6.MFSimulation(
+        sim_name=name,
+        exe_name=mf6,
+        version="mf6",
+        sim_ws=prt_ws,
+    )
+
+    # create tdis package
+    flopy.mf6.modflow.mftdis.ModflowTdis(
+        sim,
+        pname="tdis",
+        time_units="DAYS",
+        nper=FlopyReadmeCase.nper,
+        perioddata=[
+            (
+                FlopyReadmeCase.perlen,
+                FlopyReadmeCase.nstp,
+                FlopyReadmeCase.tsmult,
+            )
+        ],
+    )
+
+    # create prt model
+    prt_name = get_model_name(name, "prt")
+    prt = flopy.mf6.ModflowPrt(sim, modelname=prt_name)
+
+    # create prt discretization
+    flopy.mf6.modflow.mfgwfdis.ModflowGwfdis(
+        prt,
+        pname="dis",
+        nlay=FlopyReadmeCase.nlay,
+        nrow=FlopyReadmeCase.nrow,
+        ncol=FlopyReadmeCase.ncol,
+        top=FlopyReadmeCase.top,
+        botm=FlopyReadmeCase.botm,
+    )
+
+    # create mip package
+    flopy.mf6.ModflowPrtmip(prt, pname="mip", porosity=FlopyReadmeCase.porosity)
+
+    releasepts = [
+        # particle index, k, i, j, x, y, z (wrong coordinates for cell 0,0,0)
+        [0, 0, 0, 0, 5.5, 5.5, 0.5],
+    ]
+
+    # create prp package
+    prp_track_file = f"{prt_name}.prp.trk"
+    prp_track_csv_file = f"{prt_name}.prp.trk.csv"
+
+    if "eager" in name:
+        coordinate_check_method = "eager"
+    else:
+        coordinate_check_method = "none"
+
+    flopy.mf6.ModflowPrtprp(
+        prt,
+        pname="prp1",
+        filename=f"{prt_name}_1.prp",
+        nreleasepts=len(releasepts),
+        packagedata=releasepts,
+        perioddata={0: [("FIRST",)]},
+        track_filerecord=[prp_track_file],
+        trackcsv_filerecord=[prp_track_csv_file],
+        coordinate_check_method=coordinate_check_method,
+        print_input=True,
+        extend_tracking=True,
+    )
+
+    # create output control package
+    prt_track_file = f"{prt_name}.trk"
+    prt_track_csv_file = f"{prt_name}.trk.csv"
+    flopy.mf6.ModflowPrtoc(
+        prt,
+        pname="oc",
+        track_filerecord=[prt_track_file],
+        trackcsv_filerecord=[prt_track_csv_file],
+    )
+
+    # create the flow model interface
+    gwf_name = get_model_name(name, "gwf")
+    gwf_budget_file = gwf_ws / f"{gwf_name}.bud"
+    gwf_head_file = gwf_ws / f"{gwf_name}.hds"
+    flopy.mf6.ModflowPrtfmi(
+        prt,
+        packagedata=[
+            ("GWFHEAD", gwf_head_file),
+            ("GWFBUDGET", gwf_budget_file),
+        ],
+    )
+
+    # add explicit model solution
+    ems = flopy.mf6.ModflowEms(
+        sim,
+        pname="ems",
+        filename=f"{prt_name}.ems",
+    )
+    sim.register_solution_package(ems, [prt.name])
+
+    return sim
+
+
+def build_models(test):
+    gwf_sim = FlopyReadmeCase.get_gwf_sim(
+        test.name, test.workspace, test.targets["mf6"]
+    )
+    prt_sim = build_prt_sim(
+        test.name,
+        test.workspace,
+        test.workspace / "prt",
+        test.targets["mf6"],
+    )
+    return gwf_sim, prt_sim
+
+
+def check_output(test):
+    name = test.name
+
+    if "eager" in name:
+        buff = test.buffs[1]
+        assert any("Error: release point" in l for l in buff)
+
+
+def plot_output(test):
+    name = test.name
+    gwf_ws = test.workspace
+    prt_ws = test.workspace / "prt"
+    gwf_name = get_model_name(name, "gwf")
+    prt_name = get_model_name(name, "prt")
+    gwf_sim = test.sims[0]
+    gwf = gwf_sim.get_model(gwf_name)
+    mg = gwf.modelgrid
+    drape = "drp" in name
+
+    # check mf6 output files exist
+    gwf_head_file = f"{gwf_name}.hds"
+    prt_track_csv_file = f"{prt_name}.trk.csv"
+
+    # extract head, budget, and specific discharge results from GWF model
+    hds = HeadFile(gwf_ws / gwf_head_file).get_data()
+    bud = gwf.output.budget()
+    spdis = bud.get_data(text="DATA-SPDIS")[0]
+    qx, qy, qz = flopy.utils.postprocessing.get_specific_discharge(spdis, gwf)
+
+    # load mf6 pathline results
+    mf6_pls = pd.read_csv(prt_ws / prt_track_csv_file, na_filter=False)
+
+    # set up plot
+    fig, ax = plt.subplots(nrows=1, ncols=1, figsize=(10, 10))
+    ax.set_aspect("equal")
+
+    # plot mf6 pathlines in map view
+    pmv = flopy.plot.PlotMapView(modelgrid=mg, ax=ax)
+    pmv.plot_grid()
+    pmv.plot_array(hds[0], alpha=0.1)
+    pmv.plot_vector(qx, qy, normalize=True, color="white")
+    mf6_plines = mf6_pls.groupby(["iprp", "irpt", "trelease"])
+    for ipl, ((iprp, irpt, trelease), pl) in enumerate(mf6_plines):
+        pl.plot(
+            title=f"MF6 pathlines{' (drape)' if drape else ''}",
+            kind="line",
+            x="x",
+            y="y",
+            ax=ax,
+            legend=False,
+            color=cm.plasma(ipl / len(mf6_plines)),
+        )
+
+    # view/save plot
+    plt.show()
+    plt.savefig(prt_ws / f"{name}.png")
+
+
+@pytest.mark.parametrize("name", cases)
+def test_mf6model(name, function_tmpdir, targets, plot):
+    test = TestFramework(
+        name=name,
+        workspace=function_tmpdir,
+        build=build_models,
+        check=check_output,
+        plot=plot_output if plot else None,
+        targets=targets,
+        compare=None,
+        xfail=[False, "eager" in name],
+    )
+    test.run()
diff --git a/doc/ReleaseNotes/develop.toml b/doc/ReleaseNotes/develop.toml
@@ -32,7 +32,6 @@ section = "fixes"
 subsection = "stress"
 description = "Fixed viscosity ratio assignment in the GWT MAW package when the GWF VSC package is active."
 
-
 [[items]]
 section = "fixes"
 subsection = "internal"
@@ -46,4 +45,9 @@ description = "With the PRT model's Particle Release Package's EXTEND\\_TRACKING
 [[items]]
 section = "fixes"
 subsection = "parallel"
-description = "Fixed a memory access exception that can occur when writing convergence information to a CSV file (through the CSV\\_OUTER\\_OUTPUT option in IMS) for a parallel solution."
+description = "Fixed a memory access exception that can occur when writing convergence information to a CSV file (through the CSV\\_OUTER\\_OUTPUT option in IMS) for a parallel solution."
+
+[[items]]
+section = "features"
+subsection = "solution"
+description = "Added an optional new COORDINATE\\_CHECK\\_METHOD option to the PRT model's Particle Release Package, accepting values 'EAGER' (default) or 'NONE'. 'NONE' disables release-time verification that release point coordinates and release point cell IDs match, which can significantly reduce runtime when releasing many particles."
diff --git a/doc/mf6io/mf6ivar/dfn/prt-prp.dfn b/doc/mf6io/mf6ivar/dfn/prt-prp.dfn
@@ -271,6 +271,16 @@ optional true
 longname release time frequency
 description real number indicating the time frequency at which to release particles. This option can be used to schedule releases at a regular interval for the duration of the simulation, starting at the simulation start time. The release schedule is the union of this option, the RELEASETIMES block, and PERIOD block RELEASESETTING selections. If none of these are provided, a single release time is configured at the beginning of the first time step of the simulation's first stress period.
 
+block options
+name coordinate_check_method
+type string
+valid none eager
+reader urword
+optional true
+longname coordinate checking method
+description approach for verifying that release point coordinates are in the cell with the specified ID. By default, release point coordinates are checked at release time.
+default eager
+
 # --------------------- prt prp dimensions ---------------------
 
 block dimensions
diff --git a/src/Model/ParticleTracking/prt-prp.f90 b/src/Model/ParticleTracking/prt-prp.f90
@@ -59,6 +59,7 @@ module PrtPrpModule
     integer(I4B), pointer :: iextend => null() !< extend tracking beyond simulation's end
     integer(I4B), pointer :: ifrctrn => null() !< force ternary solution for quad grids
     integer(I4B), pointer :: iexmeth => null() !< method for iterative solution of particle exit location and time in generalized Pollock's method
+    integer(I4B), pointer :: ichkmeth => null() !< method for checking particle release coordinates are in the specified cells, 0 = none, 1 = eager
     real(DP), pointer :: extol => null() !< tolerance for iterative solution of particle exit location and time in generalized Pollock's method
     real(DP), pointer :: rttol => null() !< tolerance for coincident particle release times
     real(DP), pointer :: rtfreq => null() !< frequency for regularly spaced release times
@@ -170,6 +171,7 @@ subroutine prp_da(this)
     call mem_deallocate(this%irlstls)
     call mem_deallocate(this%ifrctrn)
     call mem_deallocate(this%iexmeth)
+    call mem_deallocate(this%ichkmeth)
     call mem_deallocate(this%extol)
     call mem_deallocate(this%rttol)
     call mem_deallocate(this%rtfreq)
@@ -260,6 +262,7 @@ subroutine prp_allocate_scalars(this)
     call mem_allocate(this%irlstls, 'IRLSTLS', this%memoryPath)
     call mem_allocate(this%ifrctrn, 'IFRCTRN', this%memoryPath)
     call mem_allocate(this%iexmeth, 'IEXMETH', this%memoryPath)
+    call mem_allocate(this%ichkmeth, 'ICHKMETH', this%memoryPath)
     call mem_allocate(this%extol, 'EXTOL', this%memoryPath)
     call mem_allocate(this%rttol, 'RTTOL', this%memoryPath)
     call mem_allocate(this%rtfreq, 'RTFREQ', this%memoryPath)
@@ -283,6 +286,7 @@ subroutine prp_allocate_scalars(this)
     this%irlstls = 0
     this%ifrctrn = 0
     this%iexmeth = 0
+    this%ichkmeth = 1
     this%extol = DEM5
     this%rttol = DSAME * DEP9
     this%rtfreq = DZERO
@@ -514,7 +518,8 @@ subroutine initialize_particle(this, particle, ip, trelease)
       z = this%rptz(ip)
     end if
 
-    call this%validate_release_point(ic, x, y, z)
+    if (this%ichkmeth > 0) &
+      call this%validate_release_point(ic, x, y, z)
 
     particle%x = x
     particle%y = y
@@ -818,6 +823,23 @@ subroutine prp_options(this, option, found)
         call store_error('DEV_EXIT_SOLVE_METHOD MUST BE &
           &1 (BRENT) OR 2 (CHANDRUPATLA)')
       found = .true.
+    case ('COORDINATE_CHECK_METHOD')
+      call this%parser%GetStringCaps(keyword)
+      select case (keyword)
+      case ('NONE')
+        this%ichkmeth = 0
+      case ('EAGER')
+        this%ichkmeth = 1
+      case default
+        write (errmsg, '(a, a)') &
+          'Unsupported coordinate check method: ', trim(keyword)
+        call store_error(errmsg)
+        write (errmsg, '(a, a)') &
+          'COORDINATE_CHECK_METHOD must be "NONE" or "EAGER"'
+        call store_error(errmsg)
+        call this%parser%StoreErrorUnit()
+      end select
+      found = .true.
     case default
       found = .false.
     end select
