fix(prt): record user tracking times for stationary particles (MODFLO…

…W-USGS#2177)

Fix MODFLOW-USGS#2172 and other problems affecting tracking events, in context of the DRY_TRACKING_METHOD STAY option and more generally.

1. User-selected times were not reported for stationary particles (the reported issue).
2. Particles stationary in the final time step did not record a time step ending event.
3. Particles stationary above the water table in partially saturated cells did not record a time step end event.
4. Particles stationary above the water table in partially saturated cells in the final time step did not record a termination event.
5. Particles still active at the end of the simulation did not record a termination event.

1-4 are only relevant/possible with newton on.

autotest/test_prt_dry.py

@@ -70,6 +70,17 @@
 for i in range(nper):
     period_data.append((perlen[i], nstp[i], tsmult[i]))
 
+user_time = 100.0
+
+# particles are released in layer 0
+offsets = [
+    (-1, 0, 0),
+    (-1, -1, 0),
+    (-1, 1, 0),
+    (-1, 0, -1),
+    (-1, 0, 1),
+]
+
 
 def build_gwf_sim(name, gwf_ws, mf6, newton=False):
     sim = flopy.mf6.MFSimulation(
@@ -267,15 +278,6 @@ def build_prt_sim(name, gwf, prt_ws, mf6, drape=False, dry_tracking_method=False
 
     flopy.mf6.ModflowPrtmip(prt, porosity=porosity, pname="mip")
 
-    # particles are released in layer 0
-    offsets = [
-        (-1, 0, 0),
-        (-1, -1, 0),
-        (-1, 1, 0),
-        (-1, 0, -1),
-        (-1, 0, 1),
-    ]
-
     lay = 1
     row, col = (int(nrow / 4), int(ncol / 4))
     prp_cells = [(lay + k, row + i, col + j) for k, i, j in offsets]
@@ -319,6 +321,8 @@ def build_prt_sim(name, gwf, prt_ws, mf6, drape=False, dry_tracking_method=False
         trackcsv_filerecord=[trackcsvfile],
         saverecord=[("BUDGET", "ALL")],
         pname="oc",
+        ntracktimes=1 if "stay" in name else 0,
+        tracktimes=[(user_time,)] if "stay" in name else None,
     )
 
     rel_prt_folder = os.path.relpath(gwf_ws, start=prt_ws)
@@ -368,11 +372,31 @@ def check_output(idx, test, snapshot):
     strtpts = pls[pls.ireason == 0]
 
     # compare to expected results
-    decimals = 1 if "drop" in name else 2
-    actual = pls.drop(["name", "icell"], axis=1).round(decimals).reset_index(drop=True)
-    # ignore particle 4, it terminates early with optimization=2 when built with ifort
-    if "drop" in name:
+    places = 1 if "drop" in name else 2
+    actual = pls.drop(["name", "icell"], axis=1).round(places).reset_index(drop=True)
+    nparts = len(offsets)  # number of particles
+
+    if "drape" in name:
+        assert len(actual[actual.ireason == 0]) == nparts  # release
+    elif "drop" in name:
+        # ignore particle 4, it terminates early when
+        # mf6 is built with optimization=2 with ifort
         actual = actual.drop(actual[actual.irpt == 4].index)
+        nparts -= 1
+        assert len(actual[actual.ireason == 0]) == nparts  # release
+    elif "stop" in name:
+        assert len(actual[actual.ireason == 0]) == nparts  # release
+    elif "stay" in name:
+        assert len(actual[actual.ireason == 0]) == nparts  # release
+        assert len(actual[actual.t == user_time]) == nparts  # user time
+    else:
+        # immediate termination, permanently unreleased
+        assert len(actual) == nparts
+
+    # in all cases, all particles should have a termination event
+    assert len(actual[actual.ireason == 3]) == nparts
+
+    # snapshot comparison
     assert snapshot == actual.to_records(index=False)
 
     plot_pathlines = False

doc/ReleaseNotes/develop.tex

@@ -24,6 +24,7 @@
 	\item GWT, GWE and PRT FMI Packages can now read a GWF Model's binary grid file via a new GWFGRID entry. This allows FMI Packages to perform the same grid equivalence checks as exchanges, which guarantees identical error-checking behavior whether a GWT, GWE or PRT Model is coupled to a GWF Model or running as a post-processor. The GWFGRID file entry is optional but recommended. A future version may make the GWFGRID entry mandatory.
 	\item A regression was recently introduced into the PRT model's generalized tracking method, in which a coordinate transformation was carried out prematurely. This could result in incorrect particle positions in and near quad-refined cells. This bug has been fixed.
 	\item The PRT model previously allowed particles to be released at any time. Release times falling outside the bounds of the simulation's time discretization could produce undefined behavior. Any release times occurring before the simulation begins (i.e. negative times) will now be skipped with a warning message. If EXTEND\_TRACKING is not enabled, release times occurring after the end of the simulation will now be skipped with a warning message as well. If EXTEND\_TRACKING is enabled, release times after the end of the simulation are allowed.
+	\item The PRT Model did not previously report all expected tracking events. In particular, time step end and termination events could go unreported with the DRY\_TRACKING\_METHOD options DROP and STAY (only relevant for Newton formulation models), and termination events were not always reported at the end of the simulation. Reporting has been corrected for the cases identified.
 	\item A profiling module is added for more enhanced performance diagnostics of the program. It can be activated through the PROFILE\_OPTION in the simulation name file.
 \end{itemize}
 

src/Model/ModelUtilities/TrackFile.f90

@@ -47,12 +47,13 @@ module TrackFileModule
   !!     1: active
   !!     2: terminated at boundary face
   !!     3: terminated in weak sink cell
-  !!     4: terminated in weak source cell**
+  !!     4: terminated in weak source cell
   !!     5: terminated in cell with no exit face
   !!     6: terminated in cell with specified zone number
   !!     7: terminated in inactive cell
   !!     8: permanently unreleased***
   !!     9: terminated in subcell with no exit face*****
+  !!     10: terminated upon simulation's end
   !!
   !! PRT shares the same status enumeration as MODPATH 7. However, some
   !! don't apply to PRT; for instance, MODPATH 7 distinguishes forwards

src/Model/ParticleTracking/prt-prp.f90

@@ -436,17 +436,11 @@ subroutine release(this, ip, trelease)
     type(ParticleType), pointer :: particle
 
     call this%initialize_particle(particle, ip, trelease)
-
-    ! Increment cumulative particle count
     np = this%nparticles + 1
     this%nparticles = np
-
-    ! Save the particle to the store
-    call this%particles%save_particle(particle, np)
+    call this%particles%put(particle, np)
     deallocate (particle)
-
-    ! Accumulate mass for release point
-    this%rptm(ip) = this%rptm(ip) + DONE
+    this%rptm(ip) = this%rptm(ip) + DONE ! TODO configurable mass
 
   end subroutine release
 
@@ -486,7 +480,6 @@ subroutine initialize_particle(this, particle, ip, trelease)
     end select
     particle%ilay = ilay
     particle%izone = this%rptzone(ic)
-
     particle%istatus = 0
     ! If the cell is inactive, either drape the particle
     ! to the top-most active cell beneath it if drape is
@@ -495,8 +488,9 @@ subroutine initialize_particle(this, particle, ip, trelease)
       ic_old = ic
       if (this%idrape > 0) then
         call this%dis%highest_active(ic, this%ibound)
-        if (ic == ic_old .or. this%ibound(ic) == 0) &
+        if (ic == ic_old .or. this%ibound(ic) == 0) then
           particle%istatus = 8
+        end if
       else
         particle%istatus = 8
       end if