Skip to content

Commit

Permalink
wip
Browse files Browse the repository at this point in the history
  • Loading branch information
@wpbonelli
wpbonelli committed Oct 30, 2023
1 parent 327e0b6 commit 4747def
Show file tree
Hide file tree
Showing 4 changed files with 68 additions and 103 deletions.
4 changes: 2 additions & 2 deletions src/Model/ParticleTracking/prt1prp1.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down
26 changes: 26 additions & 0 deletions src/Solution/ParticleTracker/Method.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -9,6 +9,7 @@ module MethodModule

private
public :: MethodType
public :: get_iatop

type, abstract :: MethodType
! private
Expand Down Expand Up @@ -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
40 changes: 8 additions & 32 deletions src/Solution/ParticleTracker/MethodDis.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand All @@ -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)
Expand Down
101 changes: 32 additions & 69 deletions src/Solution/ParticleTracker/MethodDisv.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -376,15 +349,16 @@ 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
!! These include cell index, vertex count, (ia)top and botm,
!! porosity, retfactor, and izone.
subroutine load_cellDefn_basic(this, ic, cellDefn)
implicit none
! -- dummy
Expand All @@ -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)
Expand All @@ -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
Expand All @@ -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
Expand Down Expand Up @@ -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
Expand All @@ -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
Expand Down Expand Up @@ -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
Expand All @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand All @@ -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
Expand Down

0 comments on commit 4747def

Please sign in to comment.