diff --git a/src/Model/ParticleTracking/prt1prp1.f90 b/src/Model/ParticleTracking/prt1prp1.f90 index 35b88afdc82..fc830e3578e 100644 --- a/src/Model/ParticleTracking/prt1prp1.f90 +++ b/src/Model/ParticleTracking/prt1prp1.f90 @@ -470,14 +470,14 @@ subroutine prp_ad(this) if (size(this%boundname) /= 0) & bndName = this%boundname(nps) - ! -- Add particle to particle list (todo: factor out a routine?) + ! -- Add particle to particle list (todo: factor out a routine? ! The routine should branch based on whether this is an exchange ! PRP or a normal PRP. If exchange PRP, particle identity info ! should have been injected into the PRP object by the exchange ! (e.g. this%imdl, this%iprp, this%irpt, this%trelease arrays). ! If normal PRP, imdl and iprp should be set from the pointers ! set on this PRP from the PRT model, and irpt and trelease set - ! as below. + ! as below.) this%partlist%x(np) = this%x(nps) this%partlist%y(np) = this%y(nps) this%partlist%z(np) = this%z(nps) diff --git a/src/Solution/ParticleTracker/Method.f90 b/src/Solution/ParticleTracker/Method.f90 index c566556cace..7a449f15eb6 100644 --- a/src/Solution/ParticleTracker/Method.f90 +++ b/src/Solution/ParticleTracker/Method.f90 @@ -9,6 +9,7 @@ module MethodModule private public :: MethodType + public :: get_iatop type, abstract :: MethodType ! private @@ -179,4 +180,29 @@ subroutine update(this, particle, celldefn) end subroutine update + !> @brief Get the index used to locate top elevation of a cell in the grid. + !! This is -1 if the cell is in the top layer and 1 otherwise. (kluge?) + function get_iatop(dis, ic) result(iatop) + use BaseDisModule, only: DisBaseType + implicit none + ! -- dummy + class(DisBaseType), intent(in) :: dis + integer, intent(in) :: ic + ! -- result + integer :: iatop + ! -- local + integer :: ncpl, icu + ! + ncpl = dis%get_ncpl() + icu = dis%get_nodeuser(ic) + ! + if (icu .le. ncpl) then + iatop = -1 + else + iatop = 1 + end if + ! + return + end function get_iatop + end module MethodModule diff --git a/src/Solution/ParticleTracker/MethodDis.f90 b/src/Solution/ParticleTracker/MethodDis.f90 index e975767a9bd..5a5c4c047b0 100644 --- a/src/Solution/ParticleTracker/MethodDis.f90 +++ b/src/Solution/ParticleTracker/MethodDis.f90 @@ -2,9 +2,8 @@ module MethodDisModule use KindModule, only: DP, I4B use ConstantsModule, only: DONE - use MethodModule + use MethodModule, only: MethodType, get_iatop use MethodCellPoolModule - ! use CellModule use CellDefnModule use CellRectModule use ParticleModule @@ -32,15 +31,16 @@ 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_iatop ! returns index used to locate top elevation of cell in the grid procedure, private :: get_top ! returns top elevation based on index iatop + + ! todo: maybe separate LoadCellDefn module? procedure, public :: load_cellDefn ! loads cell definition from the 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 + procedure, private :: addBoundaryFlows_cellRect ! adds BoundaryFlows from the grid to the faceflow array of a rectangular cell end type MethodDisType contains @@ -226,8 +226,6 @@ subroutine pass_mGD(this, particle) particle%iTrackingDomainBoundary(2) = inface if (inface < 5) then ! -- Map z between cells - ! iatop = this%get_iatop(ic) - ! top = this%get_top(iatop) topfrom = this%cellRect%cellDefn%top botfrom = this%cellRect%cellDefn%bot zrel = (z - botfrom) / (topfrom - botfrom) @@ -264,29 +262,6 @@ subroutine apply_mGD(this, particle, tmax) ! end subroutine apply_mGD - !> @brief Get the index used to locate top elevation of a cell in the grid - function get_iatop(this, ic) result(iatop) - use GwfDisModule - implicit none - ! -- dummy - class(MethodDisType), intent(inout) :: this - integer, intent(in) :: ic - ! -- result - integer :: iatop - ! -- local - integer :: ncpl, icu - ! - ncpl = this%fmi%dis%get_ncpl() - icu = this%fmi%dis%get_nodeuser(ic) - if (icu .le. ncpl) then - iatop = -1 ! kluge note: just returns -1 if in top layer and 1 otherwise - else - iatop = 1 - end if - ! - return - end function get_iatop - !> @brief Returns a top elevation based on index iatop function get_top(this, iatop) result(top) implicit none @@ -333,7 +308,8 @@ subroutine load_cellDefn(this, ic, cellDefn) end subroutine load_cellDefn !> @brief Loads basic cell definition components from the grid. - !! These include cell index, vertex count, and (ia)top and botm. + !! These include cell index, vertex count, (ia)top and botm, + !! porosity, retfactor, and izone. subroutine load_cellDefn_basic(this, ic, cellDefn) implicit none ! -- dummy @@ -347,11 +323,11 @@ subroutine load_cellDefn_basic(this, ic, cellDefn) ! -- cell index cellDefn%icell = ic ! - ! -- number of vertices + ! -- number of polygon vertices cellDefn%npolyverts = 4 ! rectangular cell always has 4 vertices ! ! -- iatop, top, and bot - iatop = this%get_iatop(ic) + iatop = get_iatop(this%fmi%dis, ic) cellDefn%iatop = iatop top = this%fmi%dis%top(ic) bot = this%fmi%dis%bot(ic) diff --git a/src/Solution/ParticleTracker/MethodDisv.f90 b/src/Solution/ParticleTracker/MethodDisv.f90 index a78ed5d08c5..f7fab69afb7 100644 --- a/src/Solution/ParticleTracker/MethodDisv.f90 +++ b/src/Solution/ParticleTracker/MethodDisv.f90 @@ -2,7 +2,7 @@ module MethodDisvModule use KindModule, only: DP, I4B use ConstantsModule, only: DONE - use MethodModule + use MethodModule, only: MethodType, get_iatop use MethodCellPoolModule use CellDefnModule use CellPolyModule @@ -33,20 +33,18 @@ module MethodDisvModule procedure, public :: loadsub => loadsub_mGDv ! loads the cell method procedure, public :: mapToNbrCell ! maps a location on the cell face to the shared face of a neighbor 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 + + ! todo: maybe separate LoadCellDefn module? 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 - ! adds BoundaryFlows from grid to faceflow array of a rectangular cell - procedure, public :: addBoundaryFlows_cellRect - ! adds BoundaryFlows from the grid to the faceflow array of a rectangular-quad cell - procedure, public :: addBoundaryFlows_cellRectQuad - ! adds BoundaryFlows from the grid to the faceflow array of a polygonal cell - procedure, public :: addBoundaryFlows_cellPoly + 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, private :: addBoundaryFlows_cellRect ! adds BoundaryFlows from the grid to the faceflow array of a rectangular cell + procedure, private :: addBoundaryFlows_cellRectQuad ! adds BoundaryFlows from the grid to the faceflow array of a rectangular-quad cell + procedure, private :: addBoundaryFlows_cellPoly ! adds BoundaryFlows from the grid to the faceflow array of a polygonal cell end type MethodDisvType contains @@ -256,8 +254,6 @@ subroutine mapToNbrCell(this, cellDefnin, inface, z) end if end do ! -- Map z between cells - ! iatop = this%get_iatop(ic) - ! top = this%get_top(iatop) topfrom = cellDefnin%top botfrom = cellDefnin%bot zrel = (z - botfrom) / (topfrom - botfrom) @@ -314,29 +310,6 @@ function get_npolyverts(this, ic) result(npolyverts) return end function get_npolyverts - !> @brief Get index used to locate top elevation of a cell in the grid - function get_iatop(this, ic) result(iatop) - use GwfDisvModule ! kluge??? - implicit none - ! -- dummy - class(MethodDisvType), intent(inout) :: this - integer, intent(in) :: ic - ! -- result - integer :: iatop - ! -- local - integer :: ncpl, icu - ! - ncpl = this%fmi%dis%get_ncpl() - icu = this%fmi%dis%get_nodeuser(ic) - if (icu .le. ncpl) then - iatop = -1 ! kluge note: just returns -1 if in top layer and 1 otherwise - else - iatop = 1 - end if - ! - return - end function get_iatop - !> @brief Get top elevation based on index iatop !! kluge note: not needed??? function get_top(this, iatop) result(top) @@ -376,8 +349,7 @@ 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 @@ -385,6 +357,8 @@ subroutine load_cellDefn(this, ic, cellDefn) end subroutine load_cellDefn !> @brief Loads basic cell definition components from the grid + !! These include cell index, vertex count, (ia)top and botm, + !! porosity, retfactor, and izone. subroutine load_cellDefn_basic(this, ic, cellDefn) implicit none ! -- dummy @@ -396,13 +370,14 @@ subroutine load_cellDefn_basic(this, ic, cellDefn) integer :: iatop real(DP) :: top, bot, sat ! + ! -- cell index cellDefn%icell = ic ! - ! -- Load npolyverts + ! -- number of polygon vertices cellDefn%npolyverts = this%get_npolyverts(ic) ! - ! -- Load iatop, top, and bot - iatop = this%get_iatop(ic) + ! -- iatop, top, and bot + iatop = get_iatop(this%fmi%dis, ic) cellDefn%iatop = iatop top = this%fmi%dis%top(ic) bot = this%fmi%dis%bot(ic) @@ -422,6 +397,7 @@ subroutine load_cellDefn_basic(this, ic, cellDefn) end subroutine load_cellDefn_basic !> @brief Loads polygon vertices to cell definition from the grid + !! Assumes cell index and number of vertices are already loaded. subroutine load_cellDefn_polyverts(this, cellDefn) use GwfDisvModule implicit none @@ -433,8 +409,6 @@ subroutine load_cellDefn_polyverts(this, cellDefn) integer :: ncpl ! kluge??? ! ic = cellDefn%icell - ! - if (cellDefn%npolyverts .eq. 0) cellDefn%npolyverts = this%get_npolyverts(ic) npolyverts = cellDefn%npolyverts ! ! -- Load polygon vertices. Note that the polyverts array @@ -466,6 +440,7 @@ subroutine load_cellDefn_polyverts(this, cellDefn) end subroutine load_cellDefn_polyverts !> @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 GwfDisvModule @@ -480,8 +455,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 @@ -568,12 +541,9 @@ subroutine shared_edgeface(ivlist1, ivlist2, iedgeface) end do outerloop end subroutine shared_edgeface - !> @brief Load flows into the cell from the grid - !! - !! Load polygon face, top and bottom, and net distributed - !! flows to cell definition from the grid. - !! - !< + !> @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 @@ -583,8 +553,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 @@ -624,7 +592,8 @@ subroutine load_cellDefn_flows(this, cellDefn) ! end subroutine load_cellDefn_flows - !> @brief Load boundary flows from the grid into a rectangular cell + !> @brief Load boundary flows from the grid into a rectangular cell. + !! Assumes cell index and number of vertices are already loaded. subroutine addBoundaryFlows_cellRect(this, cellDefn) implicit none ! -- dummy @@ -660,7 +629,8 @@ subroutine addBoundaryFlows_cellRect(this, cellDefn) ! end subroutine addBoundaryFlows_cellRect - !> @brief Load boundary flows from the grid into rectangular quadcell + !> @brief Load boundary flows from the grid into rectangular quadcell. + !! Assumes cell index and number of vertices are already loaded. subroutine addBoundaryFlows_cellRectQuad(this, cellDefn) implicit none ! -- dummy @@ -727,7 +697,8 @@ subroutine addBoundaryFlows_cellRectQuad(this, cellDefn) ! end subroutine addBoundaryFlows_cellRectQuad - !> @brief Load boundary flows from the grid into a polygonal cell + !> @brief Load boundary flows from the grid into a polygonal cell. + !! Assumes cell index and number of vertices are already loaded. subroutine addBoundaryFlows_cellPoly(this, cellDefn) implicit none ! -- dummy @@ -760,20 +731,15 @@ subroutine addBoundaryFlows_cellPoly(this, cellDefn) ! end subroutine addBoundaryFlows_cellPoly - !> @brief Load 180-degree vertex indicator array into cell. - !! - !! Loads 180-degree vertex indicator array to cell - !! definition and sets flags that indicate how cell - !! can be represented - !! - !< + !> @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) ! kluge note: rename??? implicit none ! -- dummy class(MethodDisvType), intent(inout) :: this type(CellDefnType), pointer, intent(inout) :: cellDefn ! -- local - ! integer :: ig,ic,npolyverts,iapnbr,iapv180 integer :: npolyverts, m, m0, m1, m2 integer :: ic integer :: num90, num180, numacute @@ -783,9 +749,6 @@ subroutine load_cellDefn_ispv180(this, cellDefn) ! kluge note: rename??? logical last180 ! ic = cellDefn%icell - ! - ! -- 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