fix termination reporting

autotest/test_prt_dry.py

@@ -72,6 +72,15 @@
 
 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(
@@ -269,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]
@@ -372,13 +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)
-    if "stay" in name:
-        assert len(actual[actual.t == user_time]) == 5
+        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

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