fix(prt): sum boundary and face flows when loading cell (#1900)

Previously cell face flows were not properly summed with boundary flows, rather for cell faces with active neighbors, the face flow replaced any boundary flows. The reason this has gone unnoticed til now is probably that IFLOWFACE is typically applied to faces on the boundary of the active domain.

Also factor out a subroutine to load face flows into the cell definition — could be pulled into base MethodType but it is a small amount of duplication and implementations may diverge.

doc/ReleaseNotes/develop.tex

@@ -26,6 +26,7 @@
 		\item For ASCII input files erroneously containing a mix of line endings, MODFLOW would sometimes proceed with unexpected results. The program was corrected to stop with an error message if an input line contained both carriage returns and line feeds.
 		\item Previously the PRT model's default behavior was to track particles until termination, as with MODPATH 7 stop time option 2 (extend). Under extended tracking, the program may not halt if particles enter a cycle in the flow system. This PR changes the default to the equivalent of MP7 stop time option 1 (final), terminating at simulation end unless a new Particle Release Point (PRP) package keyword option EXTEND\_TRACKING is provided. This is meant to provide a stronger guarantee that the program halts under default settings.
 		\item A refactor of the energy storage and transfer (EST) package associated with the GWE model type results in different input requirements that breaks backward compatibility with what was required in version 6.5.0.  The PACKAGEDATA block was removed from the EST package input.  In its place, the heat capabity of water, the specified density of water, and the latent heat of vaporization are instead given default values that can be overridden by specifying alternative values in the OPTIONS block.
+		\item The PRT model's cell face flows were improperly combined with boundary flows; for cell faces with active neighbors, the face flow replaced any boundary flows (likely a rare situation because IFLOWFACE is typically applied to faces on the boundary of the active domain). The face flow calculation has been corrected.
 	%	\item xxx
 	\end{itemize}
 

src/Solution/ParticleTracker/MethodDis.f90

@@ -19,20 +19,22 @@ module MethodDisModule
   public :: create_method_dis
 
   type, extends(MethodType) :: MethodDisType
+    private
   contains
-    procedure, public :: apply => apply_dis !< apply the method
+    procedure, public :: apply => apply_dis !< apply the DIS tracking method
     procedure, public :: deallocate !< deallocate arrays and scalars
     procedure, public :: load => load_dis !< load the method
     procedure, public :: pass => pass_dis !< pass the particle to the next domain
-    procedure, private :: get_top !< get cell top elevation
-    procedure, private :: update_flowja !< load intercell mass flows
-    procedure, private :: load_particle !< load particle properties
-    procedure, private :: load_properties !< load cell properties
-    procedure, private :: load_neighbors !< load face neighbors
-    procedure, private :: load_faceflows !< load face flows
-    procedure, private :: load_boundary_flows !< load boundary flows
-    procedure, private :: load_celldefn !< load cell definition from the grid
-    procedure, private :: load_cell !< load cell geometry and flows
+    procedure :: get_top !< get cell top elevation
+    procedure :: update_flowja !< load intercell mass flows
+    procedure :: load_particle !< load particle properties
+    procedure :: load_properties !< load cell properties
+    procedure :: load_neighbors !< load cell face neighbors
+    procedure :: load_flows !< load cell face flows
+    procedure :: load_boundary_flows_to_defn !< load boundary flows to the cell definition
+    procedure :: load_face_flows_to_defn !< load face flows to the cell definition
+    procedure :: load_celldefn !< load cell definition from the grid
+    procedure :: load_cell !< load cell geometry and flows
   end type MethodDisType
 
 contains
@@ -317,7 +319,7 @@ subroutine load_celldefn(this, ic, defn)
     defn%ispv180(1:defn%npolyverts + 1) = .false.
 
     ! -- Load face flows (assumes face neighbors already loaded)
-    call this%load_faceflows(defn)
+    call this%load_flows(defn)
 
   end subroutine load_celldefn
 
@@ -413,50 +415,57 @@ subroutine load_neighbors(this, defn)
   end subroutine load_neighbors
 
   !> @brief Load flows into the cell definition.
-  !! These include face flows and net distributed flows.
+  !! These include face, boundary and net distributed flows.
   !! Assumes cell index and number of vertices are already loaded.
-  subroutine load_faceflows(this, defn)
-    ! -- dummy
+  subroutine load_flows(this, defn)
     class(MethodDisType), intent(inout) :: this
     type(CellDefnType), intent(inout) :: defn
-    ! -- local
-    integer(I4B) :: m
-    integer(I4B) :: n
-    real(DP) :: q
 
-    ! -- Load face flows, including boundary flows. As with cell verts,
-    !    the face flow array wraps around. Top and bottom flows make up
-    !    the last two elements, respectively, for size npolyverts + 3.
-    !    If there is no flow through any face, set a no-exit-face flag.
+    ! Load face flows, including boundary flows. As with cell verts,
+    ! the face flow array wraps around. Top and bottom flows make up
+    ! the last two elements, respectively, for size npolyverts + 3.
+    ! If there is no flow through any face, set a no-exit-face flag.
     defn%faceflow = DZERO
     defn%inoexitface = 1
-    call this%load_boundary_flows(defn)
-    do m = 1, defn%npolyverts + 3
-      n = defn%facenbr(m)
-      if (n > 0) then
-        q = this%fmi%gwfflowja(this%fmi%dis%con%ia(defn%icell) + n)
-        defn%faceflow(m) = q
-      end if
-      if (defn%faceflow(m) < DZERO) defn%inoexitface = 0
-    end do
+    call this%load_boundary_flows_to_defn(defn)
+    call this%load_face_flows_to_defn(defn)
 
-    ! -- Add up net distributed flow
+    ! Add up net distributed flow
     defn%distflow = this%fmi%SourceFlows(defn%icell) + &
                     this%fmi%SinkFlows(defn%icell) + &
                     this%fmi%StorageFlows(defn%icell)
 
-    ! -- Set weak sink flag
+    ! Set weak sink flag
     if (this%fmi%SinkFlows(defn%icell) .ne. DZERO) then
       defn%iweaksink = 1
     else
       defn%iweaksink = 0
     end if
 
-  end subroutine load_faceflows
+  end subroutine load_flows
+
+  subroutine load_face_flows_to_defn(this, defn)
+    ! dummy
+    class(MethodDisType), intent(inout) :: this
+    type(CellDefnType), intent(inout) :: defn
+    ! local
+    integer(I4B) :: m, n, nfaces
+    real(DP) :: q
+
+    nfaces = defn%npolyverts + 3
+    do m = 1, nfaces
+      n = defn%facenbr(m)
+      if (n > 0) then
+        q = this%fmi%gwfflowja(this%fmi%dis%con%ia(defn%icell) + n)
+        defn%faceflow(m) = defn%faceflow(m) + q
+      end if
+      if (defn%faceflow(m) < DZERO) defn%inoexitface = 0
+    end do
+  end subroutine load_face_flows_to_defn
 
   !> @brief Add boundary flows to the cell definition faceflow array.
   !! Assumes cell index and number of vertices are already loaded.
-  subroutine load_boundary_flows(this, defn)
+  subroutine load_boundary_flows_to_defn(this, defn)
     ! -- dummy
     class(MethodDisType), intent(inout) :: this
     type(CellDefnType), intent(inout) :: defn
@@ -477,6 +486,6 @@ subroutine load_boundary_flows(this, defn)
                        this%fmi%BoundaryFlows(ioffset + 9)
     defn%faceflow(7) = defn%faceflow(7) + &
                        this%fmi%BoundaryFlows(ioffset + 10)
-  end subroutine load_boundary_flows
+  end subroutine load_boundary_flows_to_defn
 
 end module MethodDisModule

src/Solution/ParticleTracker/MethodDisv.f90

@@ -21,24 +21,26 @@ module MethodDisvModule
   public :: create_method_disv
 
   type, extends(MethodType) :: MethodDisvType
+    private
     type(CellDefnType), pointer :: neighbor => null() !< ptr to a neighbor defn
   contains
-    procedure, public :: apply => apply_disv !< apply the DISV-grid method
+    procedure, public :: apply => apply_disv !< apply the DISV tracking method
     procedure, public :: deallocate !< deallocate arrays and scalars
     procedure, public :: load => load_disv !< load the cell method
     procedure, public :: load_cell_defn !< load cell definition from the grid
     procedure, public :: pass => pass_disv !< pass the particle to the next cell
-    procedure, private :: map_neighbor !< map a location on the cell face to the shared face of a neighbor
-    procedure, private :: update_flowja !< update intercell mass flows
-    procedure, private :: load_particle !< load particle properties
-    procedure, private :: load_properties !< load cell properties
-    procedure, private :: load_polygon !< load cell polygon
-    procedure, private :: load_neighbors !< load face neighbors
-    procedure, private :: load_indicators !< load 180-degree vertex indicator
-    procedure, private :: load_faceflows !< load flows
-    procedure, private :: load_boundary_flows_to_defn_rect !< load boundary flows to a rectangular cell
-    procedure, private :: load_boundary_flows_to_defn_rect_quad !< load boundary flows to a rectangular-quad cell
-    procedure, private :: load_boundary_flows_to_defn_poly !< load boundary flows to a polygonal cell
+    procedure :: map_neighbor !< map a location on the cell face to the shared face of a neighbor
+    procedure :: update_flowja !< update intercell mass flows
+    procedure :: load_particle !< load particle properties
+    procedure :: load_properties !< load cell properties
+    procedure :: load_polygon !< load cell polygon
+    procedure :: load_neighbors !< load cell face neighbors
+    procedure :: load_indicators !< load cell 180-degree vertex indicator
+    procedure :: load_flows !< load the cell's flows
+    procedure :: load_boundary_flows_to_defn_rect !< load boundary flows to a rectangular cell definition
+    procedure :: load_boundary_flows_to_defn_rect_quad !< load boundary flows to a rectangular-quad cell definition
+    procedure :: load_boundary_flows_to_defn_poly !< load boundary flows to a polygonal cell definition
+    procedure :: load_face_flows_to_defn_poly !< load face flows to a polygonal cell definition
   end type MethodDisvType
 
 contains
@@ -317,7 +319,7 @@ subroutine load_cell_defn(this, ic, defn)
     call this%load_indicators(defn)
 
     ! -- Load flows (assumes face neighbors already loaded)
-    call this%load_faceflows(defn)
+    call this%load_flows(defn)
   end subroutine load_cell_defn
 
   !> @brief Loads cell properties to cell definition from the grid.
@@ -435,53 +437,62 @@ subroutine load_neighbors(this, defn)
   end subroutine load_neighbors
 
   !> @brief Load flows into the cell definition.
-  !! These include face flows and net distributed flows.
+  !! These include face, boundary and net distributed flows.
   !! Assumes cell index and number of vertices are already loaded.
-  subroutine load_faceflows(this, defn)
-    ! -- dummy
+  subroutine load_flows(this, defn)
+    ! dummy
     class(MethodDisvType), intent(inout) :: this
     type(CellDefnType), intent(inout) :: defn
-    ! -- local
+    ! local
     integer(I4B) :: nfaces
     integer(I4B) :: nslots
-    integer(I4B) :: m
-    integer(I4B) :: n
-    real(DP) :: q
 
-    ! -- expand faceflow array if needed
+    ! expand faceflow array if needed
     nfaces = defn%npolyverts + 3
     nslots = size(defn%faceflow)
     if (nslots < nfaces) call ExpandArray(defn%faceflow, nfaces - nslots)
 
-    ! -- Load face flows, including boundary flows. As with cell verts,
-    !    the face flow array wraps around. Top and bottom flows make up
-    !    the last two elements, respectively, for size npolyverts + 3.
-    !    If there is no flow through any face, set a no-exit-face flag.
+    ! Load face flows, including boundary flows. As with cell verts,
+    ! the face flow array wraps around. Top and bottom flows make up
+    ! the last two elements, respectively, for size npolyverts + 3.
+    ! If there is no flow through any face, set a no-exit-face flag.
     defn%faceflow = DZERO
     defn%inoexitface = 1
     call this%load_boundary_flows_to_defn_poly(defn)
-    do m = 1, nfaces
-      n = defn%facenbr(m)
-      if (n > 0) then
-        q = this%fmi%gwfflowja(this%fmi%dis%con%ia(defn%icell) + n)
-        defn%faceflow(m) = q
-      end if
-      if (defn%faceflow(m) < DZERO) defn%inoexitface = 0
-    end do
+    call this%load_face_flows_to_defn_poly(defn)
 
-    ! -- Add up net distributed flow
+    ! Add up net distributed flow
     defn%distflow = this%fmi%SourceFlows(defn%icell) + &
                     this%fmi%SinkFlows(defn%icell) + &
                     this%fmi%StorageFlows(defn%icell)
 
-    ! -- Set weak sink flag
+    ! Set weak sink flag
     if (this%fmi%SinkFlows(defn%icell) .ne. DZERO) then
       defn%iweaksink = 1
     else
       defn%iweaksink = 0
     end if
 
-  end subroutine load_faceflows
+  end subroutine load_flows
+
+  subroutine load_face_flows_to_defn_poly(this, defn)
+    ! dummy
+    class(MethodDisvType), intent(inout) :: this
+    type(CellDefnType), intent(inout) :: defn
+    ! local
+    integer(I4B) :: m, n, nfaces
+    real(DP) :: q
+
+    nfaces = defn%npolyverts + 3
+    do m = 1, nfaces
+      n = defn%facenbr(m)
+      if (n > 0) then
+        q = this%fmi%gwfflowja(this%fmi%dis%con%ia(defn%icell) + n)
+        defn%faceflow(m) = defn%faceflow(m) + q
+      end if
+      if (defn%faceflow(m) < DZERO) defn%inoexitface = 0
+    end do
+  end subroutine load_face_flows_to_defn_poly
 
   !> @brief Load boundary flows from the grid into a rectangular cell.
   !! Assumes cell index and number of vertices are already loaded.