From 307e66133425929a5ae22682d7862f9505082df3 Mon Sep 17 00:00:00 2001 From: Wes Bonelli Date: Sun, 29 Oct 2023 16:22:07 -0400 Subject: [PATCH] update docstrings, simplify and make private MethodDis cell load routines --- src/Model/ParticleTracking/prt1prp1.f90 | 3 +- src/Solution/ParticleTracker/MethodDis.f90 | 96 ++++++++-------------- 2 files changed, 35 insertions(+), 64 deletions(-) diff --git a/src/Model/ParticleTracking/prt1prp1.f90 b/src/Model/ParticleTracking/prt1prp1.f90 index baf49cd06fb..35b88afdc82 100644 --- a/src/Model/ParticleTracking/prt1prp1.f90 +++ b/src/Model/ParticleTracking/prt1prp1.f90 @@ -19,6 +19,7 @@ module PrtPrpModule use ArrayHandlersModule, only: expandarray use GlobalDataModule use TrackModule, only: TrackControlType + use GeomUtilModule, only: point_in_polygon implicit none @@ -463,7 +464,7 @@ subroutine prp_ad(this) call get_jk(icu, dis%ncpl, dis%nlay, icpl, ilay) end select - ! -- todo: check that location is within the specified cell + ! -- todo: get cell polygon and check that location is within it ! -- Set particle boundname if (size(this%boundname) /= 0) & diff --git a/src/Solution/ParticleTracker/MethodDis.f90 b/src/Solution/ParticleTracker/MethodDis.f90 index f7c8e4b534f..e975767a9bd 100644 --- a/src/Solution/ParticleTracker/MethodDis.f90 +++ b/src/Solution/ParticleTracker/MethodDis.f90 @@ -32,15 +32,14 @@ module MethodDisModule procedure, public :: apply => apply_mGD ! applies the DIS-grid method procedure, public :: pass => pass_mGD ! passes the particle to the next cell procedure, public :: loadsub => loadsub_mGD ! loads the cell method - procedure, private :: get_npolyverts ! returns the number of polygon vertices for a cell in the grid procedure, private :: get_iatop ! returns index used to locate top elevation of cell in the grid procedure, private :: get_top ! returns top elevation based on index iatop procedure, public :: load_cellDefn ! loads cell definition from the grid - procedure, public :: load_cellDefn_facenbr ! loads face neighbors to a cell object - procedure, public :: load_cellDefn_polyverts ! loads polygon vertices to a cell object - procedure, public :: load_cellDefn_flows ! loads flows to a cell object - procedure, public :: load_cellDefn_ispv180 ! loads 180-degree vertex indicator to a cell object - procedure, public :: load_cellDefn_basic ! loads basic components to a cell object from its grid + procedure, private :: load_cellDefn_basic ! loads basic components to a cell object from its grid + procedure, private :: load_cellDefn_polyverts ! loads polygon vertices to a cell object + procedure, private :: load_cellDefn_facenbr ! loads face neighbors to a cell object + procedure, private :: load_cellDefn_ispv180 ! loads 180-degree vertex indicator to a cell object + procedure, private :: load_cellDefn_flows ! loads flows to a cell object procedure, public :: addBoundaryFlows_cellRect ! adds BoundaryFlows from the grid to the faceflow array of a rectangular cell end type MethodDisType @@ -117,14 +116,15 @@ subroutine loadsub_mGD(this, particle, levelNext, submethod) ic = particle%iTrackingDomain(levelNext) ! kluge note: is cell number always known coming in? ! -- load cellDefn call this%load_cellDefn(ic, this%cellRect%cellDefn) + + ! -- If cell is active but dry, select and initialize + ! -- pass-to-bottom method and set cell method pointer if (this%fmi%ibdgwfsat0(ic) == 0) then ! kluge note: use cellDefn%sat == DZERO here instead? - ! -- Cell is active but dry, so select and initialize pass-to-bottom - ! -- cell method and set cell method pointer call methodCellPassToBot%init(particle, this%cellRect%cellDefn, & this%trackctl) submethod => methodCellPassToBot else - ! -- load rectangular cell + ! -- load rectangular cell (todo: refactor into separate routine) select type (dis => this%fmi%dis) ! kluge??? type is (GwfDisType) icu = dis%get_nodeuser(ic) @@ -156,8 +156,8 @@ subroutine loadsub_mGD(this, particle, levelNext, submethod) term = factor / areaz this%cellRect%vz1 = this%cellRect%cellDefn%faceflow(6) * term this%cellRect%vz2 = -this%cellRect%cellDefn%faceflow(7) * term - ! -- Select and initialize Pollock's cell method and set cell method - ! -- pointer + + ! -- Select and initialize Pollock's method and set method pointer call methodCellPollock%init(particle, this%cellRect, this%trackctl) submethod => methodCellPollock end if @@ -264,21 +264,6 @@ subroutine apply_mGD(this, particle, tmax) ! end subroutine apply_mGD - !> @brief Return the number of polygon vertices for a cell in the grid - function get_npolyverts(this, ic) result(npolyverts) - use GwfDisModule - implicit none - ! -- dummy - class(MethodDisType), intent(inout) :: this - integer, intent(in) :: ic - ! -- result - integer :: npolyverts - ! - npolyverts = 4 - ! - return - end function get_npolyverts - !> @brief Get the index used to locate top elevation of a cell in the grid function get_iatop(this, ic) result(iatop) use GwfDisModule @@ -340,15 +325,15 @@ subroutine load_cellDefn(this, ic, cellDefn) ! -- Load 180-degree indicator call this%load_cellDefn_ispv180(cellDefn) ! - ! -- Load flows - ! -- (assumes face neighbors already loaded) + ! -- Load flows (assumes face neighbors already loaded) call this%load_cellDefn_flows(cellDefn) ! return ! end subroutine load_cellDefn - !> @brief Loads basic cell definition components from the grid + !> @brief Loads basic cell definition components from the grid. + !! These include cell index, vertex count, and (ia)top and botm. subroutine load_cellDefn_basic(this, ic, cellDefn) implicit none ! -- dummy @@ -359,15 +344,15 @@ subroutine load_cellDefn_basic(this, ic, cellDefn) integer :: iatop real(DP) :: top, bot, sat ! + ! -- cell index cellDefn%icell = ic ! - ! -- Load npolyverts - cellDefn%npolyverts = this%get_npolyverts(ic) + ! -- number of vertices + cellDefn%npolyverts = 4 ! rectangular cell always has 4 vertices ! - ! -- Load iatop, top, and bot + ! -- iatop, top, and bot iatop = this%get_iatop(ic) cellDefn%iatop = iatop - ! cellDefn%top = this%get_top(iatop) top = this%fmi%dis%top(ic) bot = this%fmi%dis%bot(ic) sat = this%fmi%gwfsat(ic) @@ -376,7 +361,7 @@ subroutine load_cellDefn_basic(this, ic, cellDefn) cellDefn%bot = bot cellDefn%sat = sat ! - ! -- Load porosity, retfactor, and izone + ! -- porosity, retfactor, and izone cellDefn%porosity = this%porosity(ic) cellDefn%retfactor = this%retfactor(ic) cellDefn%izone = this%izone(ic) @@ -386,10 +371,7 @@ subroutine load_cellDefn_basic(this, ic, cellDefn) end subroutine load_cellDefn_basic !> @brief Load polygon vertices into cell definition from the grid - !! - !! kluge note: are polyverts even needed for MethodDis??? - !! - !< + !! Assumes cell index and number of vertices are already loaded. subroutine load_cellDefn_polyverts(this, cellDefn) use InputOutputModule ! kluge use GwfDisModule ! kluge??? @@ -404,8 +386,6 @@ subroutine load_cellDefn_polyverts(this, cellDefn) double precision :: cellx, celly, dxhalf, dyhalf ! ic = cellDefn%icell - ! - if (cellDefn%npolyverts .eq. 0) cellDefn%npolyverts = this%get_npolyverts(ic) ! kluge note: just set it to 4 ??? npolyverts = cellDefn%npolyverts ! ! -- Load polygon vertices. Note that the polyverts array @@ -440,13 +420,14 @@ subroutine load_cellDefn_polyverts(this, cellDefn) cellDefn%polyvert(5)%ivert = cellDefn%polyvert(1)%ivert cellDefn%polyvert(5)%x = cellDefn%polyvert(1)%x cellDefn%polyvert(5)%y = cellDefn%polyvert(1)%y - end select ! ... kluge??? + end select ! return ! end subroutine load_cellDefn_polyverts - !> @brief Loads face neighbors to cell definition from the grid + !> @brief Loads face neighbors to cell definition from the grid. + !! Assumes cell index and number of vertices are already loaded. subroutine load_cellDefn_facenbr(this, cellDefn) use InputOutputModule use GwfDisModule @@ -462,8 +443,6 @@ subroutine load_cellDefn_facenbr(this, cellDefn) integer :: ncpl ! ic = cellDefn%icell - ! - if (cellDefn%npolyverts .eq. 0) cellDefn%npolyverts = this%get_npolyverts(ic) npolyverts = cellDefn%npolyverts ! ! -- Load face neighbors. Note that the facenbr array @@ -516,12 +495,9 @@ subroutine load_cellDefn_facenbr(this, cellDefn) ! end subroutine load_cellDefn_facenbr - !> @brief Load flows into the cell definition - !! - !! Loads polygon face, top and bottom, and net distributed - !! flows from the grid into the cell definition. - !! - !> + !> @brief Load flows into the cell definition. + !! These include face flows and net distributed flows. + !! Assumes cell index and number of vertices are already loaded. subroutine load_cellDefn_flows(this, cellDefn) implicit none ! -- dummy @@ -531,8 +507,6 @@ subroutine load_cellDefn_flows(this, cellDefn) integer :: ic, npolyverts, m, n ! ic = cellDefn%icell - ! - if (cellDefn%npolyverts .eq. 0) cellDefn%npolyverts = this%get_npolyverts(ic) npolyverts = cellDefn%npolyverts ! ! -- Load face flows. Note that the faceflow array @@ -574,7 +548,8 @@ subroutine load_cellDefn_flows(this, cellDefn) ! end subroutine load_cellDefn_flows - !> @brief Add BoundaryFlows to the cell definition faceflow array + !> @brief Add boundary flows to the cell definition faceflow array. + !! Assumes cell index and number of vertices are already loaded. subroutine addBoundaryFlows_cellRect(this, cellDefn) implicit none ! -- dummy @@ -606,13 +581,9 @@ subroutine addBoundaryFlows_cellRect(this, cellDefn) ! end subroutine addBoundaryFlows_cellRect - !> @brief Load vertex 180-degree indicator array into cell definition - !! - !! Loads 180-degree vertex indicator array to cell - !! definition and sets flags that indicate how cell - !! can be represented. Todo: rename? - !! - !< + !> @brief Load 180-degree vertex indicator array and set flags + !! indicating how cell can be represented (kluge: latter needed?). + !! Assumes cell index and number of vertices are already loaded. subroutine load_cellDefn_ispv180(this, cellDefn) implicit none ! -- dummy @@ -621,16 +592,15 @@ subroutine load_cellDefn_ispv180(this, cellDefn) ! -- local integer :: npolyverts ! - ! -- Set up polyverts if not done previously - if (cellDefn%npolyverts .eq. 0) call this%load_cellDefn_polyverts(cellDefn) npolyverts = cellDefn%npolyverts ! ! -- Load 180-degree indicator. Note that the ispv180 array ! -- does not get reallocated if it is already allocated ! -- to a size greater than or equal to npolyverts+1. - ! -- Also, set flags that indicate how cell can be represented. call allocate_as_needed(cellDefn%ispv180, npolyverts + 1) cellDefn%ispv180(1:npolyverts + 1) = .false. + ! + ! -- Set flags that indicate how cell can be represented. cellDefn%canBeCellRect = .true. cellDefn%canBeCellRectQuad = .false. !