fix(PRT): fix conditionals controlling particle advance (#1681)

* conditions to determine whether particle should continue advancing were not independent
* set cell when initializing subcell methods, not currently needed but might be in future
* minor cleanup in test_prt_budget.py

autotest/test_prt_budget.py

@@ -87,7 +87,7 @@ def build_prt_sim(name, gwf_ws, prt_ws, mf6):
         perioddata={0: ["FIRST"]},
         track_filerecord=[prp_track_file],
         trackcsv_filerecord=[prp_track_csv_file],
-        stop_at_weak_sink="saws" in prt_name,
+        stop_at_weak_sink=False,
         boundnames=True,
     )
 

src/Solution/ParticleTracker/Method.f90

@@ -173,21 +173,27 @@ subroutine update(this, particle, cell_defn)
         particle%istatus = 6
         call this%trackfilectl%save(particle, kper=kper, &
                                     kstp=kstp, reason=3) ! reason=3: termination
+        return
       end if
-    else if (cell_defn%inoexitface .ne. 0) then
+    end if
+    if (cell_defn%inoexitface .ne. 0) then
       particle%advancing = .false.
       particle%istatus = 5
       call this%trackfilectl%save(particle, kper=kper, &
                                   kstp=kstp, reason=3) ! reason=3: termination
-    else if (cell_defn%iweaksink .ne. 0) then
+      return
+    end if
+    if (cell_defn%iweaksink .ne. 0) then
       if (particle%istopweaksink .ne. 0) then
         particle%advancing = .false.
         particle%istatus = 3
         call this%trackfilectl%save(particle, kper=kper, &
                                     kstp=kstp, reason=3) ! reason=3: termination
+        return
       else
         call this%trackfilectl%save(particle, kper=kper, &
                                     kstp=kstp, reason=4) ! reason=4: exited weak sink
+        return
       end if
     end if
   end subroutine update

src/Solution/ParticleTracker/MethodCellPollock.f90

@@ -64,6 +64,7 @@ subroutine load_mcp(this, particle, next_level, submethod)
       call this%load_subcell(particle, subcell)
     end select
     call method_subcell_plck%init( &
+      cell=this%cell, &
       subcell=this%subcell, &
       trackfilectl=this%trackfilectl, &
       tracktimes=this%tracktimes)

src/Solution/ParticleTracker/MethodCellPollockQuad.f90

@@ -62,6 +62,7 @@ subroutine load_mcpq(this, particle, next_level, submethod)
       call this%load_subcell(particle, subcell)
     end select
     call method_subcell_plck%init( &
+      cell=this%cell, &
       subcell=this%subcell, &
       trackfilectl=this%trackfilectl, &
       tracktimes=this%tracktimes)

src/Solution/ParticleTracker/MethodCellTernary.f90

@@ -80,6 +80,7 @@ subroutine load_mct(this, particle, next_level, submethod)
       call this%load_subcell(particle, subcell)
     end select
     call method_subcell_tern%init( &
+      cell=this%cell, &
       subcell=this%subcell, &
       trackfilectl=this%trackfilectl, &
       tracktimes=this%tracktimes)

src/Solution/ParticleTracker/MethodDis.f90

@@ -383,9 +383,7 @@ subroutine load_flows_to_defn(this, defn)
     ! -- allocate faceflow array
     call ExpandArray(defn%faceflow, npolyverts + 3)
 
-    ! -- Load face flows. Note that the faceflow array
-    ! -- does not get reallocated if it is already allocated
-    ! -- to a size greater than or equal to npolyverts+3.
+    ! -- Load face flows.
     defn%faceflow = 0d0 ! kluge note: eventually use DZERO for 0d0 throughout
     ! -- As with polygon nbrs, polygon face flows wrap around for
     ! -- convenience at position npolyverts+1, and bot and top flows
@@ -394,8 +392,8 @@ subroutine load_flows_to_defn(this, defn)
       n = defn%facenbr(m)
       if (n > 0) &
         defn%faceflow(m) = this%fmi%gwfflowja(this%fmi%dis%con%ia(ic) + n)
-      ! if (cellDefn%faceflow(m) < 0d0) defn%inoexitface = 0
     end do
+
     ! -- Add BoundaryFlows to face flows
     call this%load_boundary_flows_to_defn(defn)
     ! -- Set inoexitface flag
@@ -414,6 +412,7 @@ subroutine load_flows_to_defn(this, defn)
     else
       defn%iweaksink = 0
     end if
+
   end subroutine load_flows_to_defn
 
   !> @brief Add boundary flows to the cell definition faceflow array.
@@ -423,9 +422,11 @@ subroutine load_boundary_flows_to_defn(this, defn)
     class(MethodDisType), intent(inout) :: this
     type(CellDefnType), intent(inout) :: defn
     ! -- local
+    integer(I4B) :: ic
     integer(I4B) :: ioffset
 
-    ioffset = (defn%icell - 1) * 10
+    ic = defn%icell
+    ioffset = (ic - 1) * 10
     defn%faceflow(1) = defn%faceflow(1) + &
                        this%fmi%BoundaryFlows(ioffset + 1) ! kluge note: should these be additive (seems so)???
     defn%faceflow(2) = defn%faceflow(2) + &