From f598f4c03a612ed79d936438da9104592e64a463 Mon Sep 17 00:00:00 2001 From: wpbonelli Date: Wed, 26 Jun 2024 08:45:13 -0400 Subject: [PATCH] fix(prt): sum boundary and face flows when loading cell (#1900) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit 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 | 1 + src/Solution/ParticleTracker/MethodDis.f90 | 81 +++++++++++--------- src/Solution/ParticleTracker/MethodDisv.f90 | 83 ++++++++++++--------- 3 files changed, 93 insertions(+), 72 deletions(-) diff --git a/doc/ReleaseNotes/develop.tex b/doc/ReleaseNotes/develop.tex index b17638a8615..2a33fc88bc3 100644 --- a/doc/ReleaseNotes/develop.tex +++ b/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} diff --git a/src/Solution/ParticleTracker/MethodDis.f90 b/src/Solution/ParticleTracker/MethodDis.f90 index e8224edc4e1..997e58842f7 100644 --- a/src/Solution/ParticleTracker/MethodDis.f90 +++ b/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 diff --git a/src/Solution/ParticleTracker/MethodDisv.f90 b/src/Solution/ParticleTracker/MethodDisv.f90 index c3170f38619..d8acddc4379 100644 --- a/src/Solution/ParticleTracker/MethodDisv.f90 +++ b/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.