fix(prt): fix drape behavior, simplify localz calculation

Author: wpbonelli

autotest/__snapshots__/test_prt_dry/test_mf6model[1-prtdry_drape].npy

@@ -0,0 +1,238 @@
+▗ ▗   ▖ ▖  Claude Code v2.1.139
+           Sonnet 4.6 · Claude Pro
+  ▘▘ ▝▝    ~/dev/modflow6/autotest
+
+✻ Conversation compacted (ctrl+o for history)
+
+
+❯ /compact                                                                      
+  ⎿  Compacted (ctrl+o to see full summary)
+  ⎿  Read test_prt_subcell_cycle.py (314 lines)
+  ⎿  Read temp/synthetic_valley/run_prt_debug.py (198 lines)
+  ⎿  Read ../src/Solution/ParticleTracker/Method/MethodSubcell.f90 (79 lines)
+  ⎿  Read ../src/Solution/ParticleTracker/Method/MethodCell.f90 (349 lines)
+
+❯ Ok, good. The plot works. But the test still passes. How could this be? What  
+  did we do differently than in the Synthetic Valley example? Would the choice  
+  of WEL instead of SFR matter? I doubt it, but you tell me. I  figure we may   
+  need to closely compare the original example with our new test case to see    
+  what we're missing.                                                           
+
+⏺ Good question. WEL vs. SFR almost certainly doesn't matter — both are
+  distributed sinks that PRT treats the same way through the FMI. Let me look at
+   the ternary method and cell defn code to understand what's really different.
+
+  Searched for 2 patterns, read 2 files (ctrl+o to expand)
+
+⏺ Now let me find where faceflow is populated for PRT cells, and what
+  inoexitface depends on:
+
+  Searched for 1 pattern (ctrl+o to expand)
+
+⏺ Bash(find /Users/wbonelli/dev/modflow6/src/Solution/ParticleTracker -name
+      "*.f90" | xargs grep -l "faceflow\|inoexitface\|distflow" | head -10)
+  ⎿  /Users/wbonelli/dev/modflow6/src/Solution/ParticleTracker/Domain/CellD
+     efn.f90
+     /Users/wbonelli/dev/modflow6/src/Solution/ParticleTracker/Domain/CellR
+     … +15 lines (ctrl+o to expand)
+
+  Read 1 file (ctrl+o to expand)
+
+⏺ Now I need to see how inoexitface and faceflow are set for DISV cells:
+
+  Read 1 file (ctrl+o to expand)
+
+⏺ Now let me find load_cell_no_exit_face which sets inoexitface:
+
+  Searched for 1 pattern, read 1 file (ctrl+o to expand)
+
+⏺ Let me look at PrtFmiModule to understand how BoundaryFlows vs SinkFlows are
+  assigned:
+
+  Searched for 1 pattern (ctrl+o to expand)
+
+⏺ Bash(find /Users/wbonelli/dev/modflow6/src -name "*.f90" | xargs grep -l
+      "BoundaryFlows\|SinkFlows\|SourceFlows" 2>/dev/null | head -10)
+  ⎿  /Users/wbonelli/dev/modflow6/src/Solution/ParticleTracker/Method/Metho
+     dDisv.f90
+     /Users/wbonelli/dev/modflow6/src/Solution/ParticleTracker/Method/Metho
+     … +2 lines (ctrl+o to expand)
+
+  Read 1 file (ctrl+o to expand)
+
+⏺ Now I have the full picture. Let me check if there's a tracking CSV from a
+  recent run, and also look at the particle starting position logic:
+
+  Searched for 1 pattern, read 1 file (ctrl+o to expand)
+
+⏺ The particle terminates immediately (barely moves). Let me check status codes
+  and look for the actual issue:
+
+  Searched for 3 patterns (ctrl+o to expand)
+
+⏺ Now I have the full picture. Let me read the actual face flows from the CBC
+  file to confirm what the GWF is producing:
+
+  Listed 1 directory (ctrl+o to expand)
+
+⏺ Bash(cd "/private/var/folders/r9/yy7jb_8x0_s7zdkp6flb42980000gp/T/pytest-of-wb
+      onelli/pytest-79/test_mf6model_0-prtscyc_0/" && python3 -c "
+      import flopy…)
+  ⎿  Error: Exit code 1
+     Traceback (most recent call last):
+       File "<string>", line 2, in <module>
+     ModuleNotFoundError: No module named 'flopy'  ⎿  Interrupted · What should Claude do instead?
+
+❯ Use the pixi `rtd` env                                                        
+
+⏺ Bash(cd "/private/var/folders/r9/yy7jb_8x0_s7zdkp6flb42980000gp/T/pytest-of-wb
+      onelli/pytest-79/test_mf6model_0-prtscyc_0/" && pixi run -e rtd python3 -c
+       "
+      import flo…)
+  ⎿  Error: Exit code 1
+     Error:   × could not find pixi.toml or pyproject.toml with tool.pixi at 
+     directory /
+       │ private/var/folders/r9/yy7jb_8x0_s7zdkp6flb42980000gp/T/pytest-of-
+       │ wbonelli/pytest-79/test_mf6model_0-prtscyc_0
+  ⎿  Interrupted · What should Claude do instead?
+
+❯ you'll have to specify --manifest-path if you're running it from the temp dir 
+
+⏺ Bash(cd "/private/var/folders/r9/yy7jb_8x0_s7zdkp6flb42980000gp/T/pytest-of-wb
+      onelli/pytest-79/test_mf6model_0-prtscyc_0/" && pixi run --manifest-path
+      /Users/wbonell…)
+  ⎿  Cell 0 face flows (positive=inflow to cell 0):
+       neighbor cell 1: q = 0.02370000 m3/day
+       neighbor cell 2: q = -0.00100000 m3/day
+     … +6 lines (ctrl+o to expand)
+  ⎿  Shell cwd was reset to /Users/wbonelli/dev/modflow6
+
+⏺ Good — the face flows exactly match FACE_Q. So inoexitface = 0, the ternary
+  method IS invoked. The problem is now clear from the tracking output.
+
+  The particle terminates with istatus=9 = TERM_NO_EXITS_SUB — that's a
+  subcell-level no-exit, not a cell-level one. The particle is released at (0, 
+  0), which is the geometric centroid — vertex 2 of every ternary subcell
+  simultaneously. The code assigns it to subcell 1 (the first one satisfying the
+   barycentric condition). But it starts there at exactly gamma=1 (the centroid
+  vertex), which is a degenerate position: the ternary subcell solver can't
+  compute a valid exit direction from a corner vertex, so it immediately fires
+  TERM_NO_EXITS_SUB. The particle never gets near subcell 2 — the one with the
+  tiny outflow face — at all.
+
+  WEL vs SFR is not the issue. Confirmed in code: WEL (no IFLOWFACE) goes to
+  SinkFlows → distflow, which is stored on CellDefnType but is never used in
+  vertvelo(). divcell is computed purely from faceflow (lateral GWF-GWF flows),
+  so the convergent field from the five inflow faces is fully in play. That part
+   of the setup is correct.
+
+  What differs from Synthetic Valley: particles there enter cell 798 through an
+  inflow face from a neighboring cell — a non-degenerate interior point. They
+  then travel through the cell before hitting the oscillation zone near face 1.
+  Our release at the centroid vertex short-circuits all of that.
+
+  Fix: release the particle inside subcell 2 (VERTS[1]→VERTS[2]→centroid), the
+  one whose exterior is the tiny outflow face. The centroid of that triangle is
+  a safe, non-degenerate starting point:
+
+  Read 1 file (ctrl+o to expand)
+
+⏺ Update(test_prt_subcell_cycle.py)
+  ⎿  Added 7 lines, removed 4 lines
+      183      )
+      184      flopy.mf6.ModflowPrtmip(prt, pname="mip", porosity=POROSIT
+           Y)
+      185  
+      186 -    # Release at the centroid of the central cell. The converg
+          -ent lateral                                                   
+      187 -    # velocity field drives the particle toward the short-face
+          - subcells                                                     
+      188 -    # where the oscillation occurs.                           
+      186 +    # Release inside subcell 1 (VERTS[1]->VERTS[2]->centroid),
+          + the subcell                                                  
+      187 +    # whose exterior is the tiny outflow face. Starting at the
+          + cell centroid                                                
+      188 +    # (shared vertex of all subcells) is degenerate for the te
+          +rnary solver                                                  
+      189 +    # and causes immediate TERM_NO_EXITS_SUB. Use the triangle
+          + centroid instead.                                            
+      190      z = 0.5 * (TOP + BOTM)
+      190 -    prpdata = [(0, (0, 0), 0.0, 0.0, z)]                      
+      191 +    sc1_cx = (VERTS[1][0] + VERTS[2][0] + 0.0) / 3.0          
+      192 +    sc1_cy = (VERTS[1][1] + VERTS[2][1] + 0.0) / 3.0          
+      193 +    prpdata = [(0, (0, 0), sc1_cx, sc1_cy, z)]                
+      194      prt_csv = f"{prtname}.trk.csv"
+      195      flopy.mf6.ModflowPrtprp(
+      196          prt,
+
+⏺ Now let's run the test and watch for the cycle detection message:
+
+⏺ Bash(cd /Users/wbonelli/dev/modflow6/autotest && pixi run --manifest-path
+      /Users/wbonelli/dev/modflow6/pixi.toml -e rtd pytest
+      test_prt_subcell_cycle.py -v -s 2>&1 |…)
+  ⎿  by the USGS or the U.S. Government as to the functionality of the
+     software and related material nor shall the fact of release
+     constitute any such warranty. The software is provided on the
+     … +23 lines (ctrl+o to expand)
+  ⎿  (timeout 2m)
+
+⏺ Still passing cleanly. Let me check the new tracking CSV to see what subcells
+  the particle visits:
+
+⏺ Bash(LATEST=$(find /private/var/folders/r9/-name "prtscyc_prt.trk.csv" -newer
+      /private/var/folders/r9/yy7jb_8x0_s7zdkp6flb42980000gp/T/pytest-of-wbonell
+      i/pytest-79/t…)
+  ⎿  Error: Exit code 1
+
+     cat: : No such file or directory
+
+  Searched for 1 pattern, read 1 file (ctrl+o to expand)
+
+⏺ Progress — the particle now travels for 323 days before terminating with
+  TERM_NO_EXITS_SUB inside subcell 2, near face 1. But with FACE_Q[1] = -0.001,
+  the outflow there is so tiny (70× smaller than the SV-scaled value) that the
+  convergent field overwhelms it, creating a velocity sink inside the subcell
+  instead of a spoke-exit oscillation.
+
+  Let me restore the original SV-extracted value and keep the new release point:
+
+⏺ Update(test_prt_subcell_cycle.py)
+  ⎿  Added 1 line, removed 1 line
+      46  # synthetic_valley model (layer 0), scaled by 0.1 for this mini
+          -model.
+      47  # Face indices correspond to edges VERTS[i]→VERTS[(i+1)%6].
+      48  # Face 1 (v1→v2, short upper-right edge) is the only exit (outf
+          low).
+      49 -FACE_Q = np.array([0.0237, -0.001, 0.1280, 1.3583, 0.0768, 0.96
+         -73])                                                           
+      49 +FACE_Q = np.array([0.0237, -0.0703, 0.1280, 1.3583, 0.0768, 0.9
+         +673])                                                          
+      50  
+      51  
+      52  def _outward_normal(vi, vj):
+
+⏺ Bash(pixi run --manifest-path /Users/wbonelli/dev/modflow6/pixi.toml -e rtd
+      pytest test_prt_subcell_cycle.py -v -s 2>&1 | grep -E
+      "PASSED|FAILED|ERROR|Cyclic|subcell…)
+  ⎿  test_prt_subcell_cycle.py::test_mf6model[0-prtscyc] writing simulation
+     ...
+      Normal termination of simulation.
+     … +3 lines (ctrl+o to expand)
+
+⏺ Bash(find /private/var/folders/r9/yy7jb_8x0_s7zdkp6flb42980000gp/T/pytest-of-w
+      bonelli/ -name "prtscyc_prt.trk.csv" | sort -t- -k2 -n | tail -1 | xargs
+      cat)
+  ⎿  kper,kstp,imdl,iprp,irpt,ilay,icell,izone,istatus,ireason,trelease,t,x
+     ,y,z,name
+     1,1,1,1,1,1,1,0,1,0,0.0000000000000000,0.0000000000000000,13.835446670
+     … +7 lines (ctrl+o to expand)
+
+⏺ Good progress. The particle now travels for 172 days and exits cleanly through
+   face 1. The oscillation sweet spot requires FACE_Q[1] to be small enough that
+   the spoke velocity wins over the face-1 exit velocity, but not so small it
+  creates a velocity sink. Let me try a few intermediate values to locate that
+  regime:
+
+⏺ Update(test_prt_subcell_cycle.py)
+  ⎿  (No output)
+  ⎿  Interrupted · What should Claude do instead?

autotest/synthetic-synthetic-valley-cycle-repro.txt

@@ -276,7 +276,8 @@ def check_output(idx, test):
         )
 
     if drape:
-        assert mf6_pls.shape[0] == 36
+        assert mf6_pls.shape[0] == 27
+        assert mf6_pls[mf6_pls.ireason == 6].empty
     else:
         # expect no movement without drape
         assert mf6_pls.shape[0] == 9