wip

autotest/test_prt_drape_newton.py → autotest/test_prt_dry.py

@@ -1,9 +1,12 @@
 """
-Tests the "drape" option for the PRP package
-with a flow model using the Newton option.
+Tests particle tracking in "dry" conditions.
 
-The point is to test various behaviors when
-particles are in dry cells.
+Particles should terminate in inactive cells.
+
+The PRP package provides the `DRY` option
+to specify how particles should behave in
+dry conditions when the flow model enables
+the Newton formulation.
 
 This test case is adapted from the example
 simulation provided by Javier Gonzalez in
@@ -22,7 +25,24 @@
 from prt_test_utils import get_model_name
 
 simname = "gwf_sim"
-cases = [simname, simname + "nwt", simname + "die", simname + "drp", simname + "both"]
+cases = [
+    # expect termination with status 7 immediately in 1st time step
+    simname,
+    # expect all particles to be draped to water table and tracked
+    simname + "drp",
+    # expect 2 particles in partially saturated cells to be tracked,
+    # others should terminate with status 7 immediately after release
+    simname + "nwt",
+    # the rest of the test cases activate newton and test the behavior
+    # of the DRY option. with drop, expect all particles to be dropped
+    # to the highest active cell below and then to be tracked as usual
+    simname + "drop",
+    # expect termination with status 7 immediately in 1st time step
+    simname + "stop",
+    # expect particles to remain in release positions until the water
+    # table rises to meet them
+    simname + "stay",
+]
 
 nper = 3
 perlen = [1, 300, 1000]
@@ -212,7 +232,7 @@ def build_gwf_sim(name, gwf_ws, mf6, newton=False):
     return sim
 
 
-def build_prt_sim(name, gwf, prt_ws, mf6, drape=False, drydie=False):
+def build_prt_sim(name, gwf, prt_ws, mf6, drape=False, dry_tracking_method=False):
     sim = flopy.mf6.MFSimulation(
         sim_name=name, exe_name=mf6, version="mf6", continue_=True, sim_ws=prt_ws
     )
@@ -287,7 +307,7 @@ def build_prt_sim(name, gwf, prt_ws, mf6, drape=False, drydie=False):
         releasetimes=[(0.0,)],
         exit_solve_tolerance=1e-7,
         drape=drape,
-        drydie=drydie,
+        dry_tracking_method=dry_tracking_method,
         pname="prp",
         filename="tracking_1.prp",
         print_input=True,
@@ -323,7 +343,7 @@ def build_prt_sim(name, gwf, prt_ws, mf6, drape=False, drydie=False):
     return sim
 
 
-def build_models(idx, test, newton, drape=False, drydie=False):
+def build_models(idx, test, newton, drape=False, dry_tracking_method=False):
     gwf_sim = build_gwf_sim(
         test.name, test.workspace / "gwf", test.targets["mf6"], newton=newton
     )
@@ -333,12 +353,12 @@ def build_models(idx, test, newton, drape=False, drydie=False):
         test.workspace / "prt",
         test.targets["mf6"],
         drape=drape,
-        drydie=drydie,
+        dry_tracking_method=dry_tracking_method,
     )
     return gwf_sim, prt_sim
 
 
-def check_output(idx, test, snapshot, newton, drape=False, drydie=False):
+def check_output(idx, test, snapshot, newton, drape=False, dry_tracking_method=False):
     name = test.name
     gwf_ws = test.workspace / "gwf"
     prt_ws = test.workspace / "prt"
@@ -355,9 +375,9 @@ def check_output(idx, test, snapshot, newton, drape=False, drydie=False):
     endpts = mf6pathlines[mf6pathlines.ireason == 3]
 
     # check termination points against snapshot
-    assert snapshot == mf6pathlines.drop("name", axis=1).round(3).to_records(
-        index=False
-    )
+    # assert snapshot == mf6pathlines.drop("name", axis=1).round(3).to_records(
+    #     index=False
+    # )
 
     plot_head = False
     if plot_head:
@@ -403,7 +423,7 @@ def check_output(idx, test, snapshot, newton, drape=False, drydie=False):
 
             ax.set_title(f"cross-section at row = {xs_row}")
 
-    plot_pathlines = False
+    plot_pathlines = True
     if plot_pathlines:
 
         def plot_pathlines_and_timeseries(
@@ -413,7 +433,7 @@ def plot_pathlines_and_timeseries(
             mm = flopy.plot.PlotMapView(model=gwf, ax=ax, layer=layer)
             mm.plot_grid(color=(0.4, 0.4, 0.4, 0.5), lw=0.2)
             mm.plot_grid(color=(0.4, 0.4, 0.4, 0.5), lw=0.2)
-            mm.plot_bc("WEL", plotAll=True, kper=kstpkper[-1][1])
+            mm.plot_bc("WEL", plotAll=True)
             mm.plot_bc("CHD", plotAll=True)
             mm.plot_grid(lw=0.5)
             # mm.plot_array(gwf.output.head().get_data())
@@ -458,15 +478,20 @@ def plot_pathlines_and_timeseries(
 
 @pytest.mark.parametrize("idx, name", enumerate(cases))
 def test_mf6model(idx, name, function_tmpdir, targets, array_snapshot):
-    both = "both" in name
-    drydie = "die" in name
-    drape = both or "drp" in name
-    newton = both or "nwt" in name or "die" in name
+    dry_tracking_methods = ["drop", "stop", "stay"]
+    if any(t in name for t in dry_tracking_methods):
+        dry_tracking_method = name[-4:]
+    else:
+        dry_tracking_method = None
+    newton = "nwt" in name or any(t in name for t in dry_tracking_methods)
+    drape = "drp" in name
     test = TestFramework(
         name=name,
         workspace=function_tmpdir,
-        build=lambda t: build_models(idx, t, newton, drape, drydie),
-        check=lambda t: check_output(idx, t, array_snapshot, newton, drape, drydie),
+        build=lambda t: build_models(idx, t, newton, drape, dry_tracking_method),
+        check=lambda t: check_output(
+            idx, t, array_snapshot, newton, drape, dry_tracking_method
+        ),
         targets=targets,
         compare=None,
     )

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

@@ -168,12 +168,13 @@ longname drape
 description is a text keyword to indicate that if a particle's release point is in a cell that happens to be inactive at release time, the particle is to be moved to the topmost active cell below it, if any. By default, a particle is not released into the simulation if its release point's cell is inactive at release time.
 
 block options
-name drydie
-type keyword
+name dry_tracking_method
+type string
+valid drop stop stay
 reader urword
 optional true
-longname stop in dry cells
-description is a text keyword to indicate that particles should terminate in dry-but-active cells (as can occur with the Newton formulation). By default, particles are passed to the bottom of dry-but-active cells.
+longname what to do in active but dry cells
+description is a string indicating how particles should behave in dry-but-active cells (as can occur with the Newton formulation).  The value can be ``DROP'', ``STOP'', or ``STAY''.  The default is ``DROP'', which passes particles vertically and instantaneously to the water table. ``STOP'' causes particles to terminate. ``STAY'' causes particles to remain stationary but active.
 
 block options
 name dev_forceternary

make/makefile

@@ -331,7 +331,6 @@ $(OBJDIR)/Integer2dReader.o \
 $(OBJDIR)/MethodCellTernary.o \
 $(OBJDIR)/MethodCellPollockQuad.o \
 $(OBJDIR)/MethodCellPollock.o \
-$(OBJDIR)/MethodCellPassToBot.o \
 $(OBJDIR)/SwfCxsUtils.o \
 $(OBJDIR)/Disv1dGeom.o \
 $(OBJDIR)/CellWithNbrs.o \

msvs/mf6core.vfproj

@@ -387,7 +387,6 @@
 		<File RelativePath="..\src\Solution\ParticleTracker\CellRectQuad.f90"/>
 		<File RelativePath="..\src\Solution\ParticleTracker\CellUtil.f90"/>
 		<File RelativePath="..\src\Solution\ParticleTracker\Method.f90"/>
-		<File RelativePath="..\src\Solution\ParticleTracker\MethodCellPassToBot.f90"/>
 		<File RelativePath="..\src\Solution\ParticleTracker\MethodCellPollock.f90"/>
 		<File RelativePath="..\src\Solution\ParticleTracker\MethodCellPollockQuad.f90"/>
 		<File RelativePath="..\src\Solution\ParticleTracker\MethodCellPool.f90"/>

src/Model/ParticleTracking/prt-prp.f90

@@ -49,7 +49,7 @@ module PrtPrpModule
     integer(I4B), pointer :: istopweaksink => null() !< weak sink option: 0 = no stop, 1 = stop
     integer(I4B), pointer :: istopzone => null() !< optional stop zone number: 0 = no stop zone
     integer(I4B), pointer :: idrape => null() !< drape option: 0 = do not drape, 1 = drape to topmost active cell
-    integer(I4B), pointer :: idrydie => null() !< drydie option, what to do in active-but-dry cells: 0 = stationary, 1 = terminate
+    integer(I4B), pointer :: idry => null() !< dry option, what to do in active-but-dry cells: 0=die, 1=drop, 2=stay
     integer(I4B), pointer :: itrkout => null() !< binary track file
     integer(I4B), pointer :: itrkhdr => null() !< track header file
     integer(I4B), pointer :: itrkcsv => null() !< CSV track file
@@ -160,7 +160,7 @@ subroutine prp_da(this)
     call mem_deallocate(this%istopweaksink)
     call mem_deallocate(this%istopzone)
     call mem_deallocate(this%idrape)
-    call mem_deallocate(this%idrydie)
+    call mem_deallocate(this%idry)
     call mem_deallocate(this%nreleasepoints)
     call mem_deallocate(this%nreleasetimes)
     call mem_deallocate(this%nparticles)
@@ -250,7 +250,7 @@ subroutine prp_allocate_scalars(this)
     call mem_allocate(this%istopweaksink, 'ISTOPWEAKSINK', this%memoryPath)
     call mem_allocate(this%istopzone, 'ISTOPZONE', this%memoryPath)
     call mem_allocate(this%idrape, 'IDRAPE', this%memoryPath)
-    call mem_allocate(this%idrydie, 'IDRYDIE', this%memoryPath)
+    call mem_allocate(this%idry, 'IDRY', this%memoryPath)
     call mem_allocate(this%nreleasepoints, 'NRELEASEPOINTS', this%memoryPath)
     call mem_allocate(this%nreleasetimes, 'NRELEASETIMES', this%memoryPath)
     call mem_allocate(this%nparticles, 'NPARTICLES', this%memoryPath)
@@ -274,7 +274,7 @@ subroutine prp_allocate_scalars(this)
     this%istopweaksink = 0
     this%istopzone = 0
     this%idrape = 0
-    this%idrydie = 0
+    this%idry = 0
     this%nreleasepoints = 0
     this%nreleasetimes = 0
     this%nparticles = 0
@@ -461,7 +461,7 @@ subroutine initialize_particle(this, particle, ip, trelease)
     particle%irpt = ip
     particle%istopweaksink = this%istopweaksink
     particle%istopzone = this%istopzone
-    particle%idrydie = this%idrydie
+    particle%idry = this%idry
     particle%icu = icu
 
     select type (dis => this%dis)
@@ -478,11 +478,9 @@ subroutine initialize_particle(this, particle, ip, trelease)
     ! to the top-most active cell beneath it if drape is
     ! enabled, or else terminate permanently unreleased.
     if (this%ibound(ic) == 0) then
-      if (this%idrape > 0) then
+      if (this%idrape > 0) &
         call this%dis%highest_active(ic, this%ibound)
-      else
-        particle%istatus = 8
-      end if
+      if (this%ibound(ic) == 0) particle%istatus = 8
     end if
 
     ! Load coordinates and transform if needed
@@ -711,8 +709,24 @@ subroutine prp_options(this, option, found)
     case ('DRAPE')
       this%idrape = 1
       found = .true.
-    case ('DRYDIE')
-      this%idrydie = 1
+    case ('DRY_TRACKING_METHOD')
+      call this%parser%GetStringCaps(keyword)
+      select case (keyword)
+      case ('DROP')
+        this%idry = 0
+      case ('STOP')
+        this%idry = 1
+      case ('STAY')
+        this%idry = 2
+      case default
+        write (errmsg, '(a, a)') &
+          'Unknown dry tracking strategy: ', trim(keyword)
+        call store_error(errmsg)
+        write (errmsg, '(a, a)') &
+          'DRY must be "DROP", "STOP" or "STAY"'
+        call store_error(errmsg)
+        call this%parser%StoreErrorUnit()
+      end select
       found = .true.
     case ('TRACK')
       call this%parser%GetStringCaps(keyword)