Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

refactor(prt): deconflict naming, use enum for tracking level #2202

Draft
wants to merge 3 commits into
base: develop
Choose a base branch
from
Draft
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 1 addition & 0 deletions doc/mf6io/prt/prt.tex
Original file line number Diff line number Diff line change
Expand Up @@ -115,6 +115,7 @@ \subsection{Particle Track Output}
\item \texttt{7}: particle terminated in an inactive cell
\item \texttt{8}: particle terminated immediately upon release into a dry cell
\item \texttt{9}: particle terminated in a subcell with no exit face
\item \texttt{10}: particle terminated at stop time or end of simulation
wpbonelli marked this conversation as resolved.
Show resolved Hide resolved
\end{description}

The IREASON field indicates the reason the particle track record was saved:
Expand Down
14 changes: 7 additions & 7 deletions src/Model/ParticleTracking/prt-prp.f90
Original file line number Diff line number Diff line change
Expand Up @@ -445,7 +445,7 @@ subroutine release(this, ip, trelease)
end subroutine release

subroutine initialize_particle(this, particle, ip, trelease)
use ParticleModule, only: TERM_UNRELEASED
use ParticleModule, only: TERM_UNRELEASED, LVL_MODEL, LVL_CELL, LVL_SUBCELL
class(PrtPrpType), intent(inout) :: this !< this instance
type(ParticleType), pointer, intent(inout) :: particle !< the particle
integer(I4B), intent(in) :: ip !< particle index
Expand Down Expand Up @@ -528,12 +528,12 @@ subroutine initialize_particle(this, particle, ip, trelease)
end if

particle%ttrack = particle%trelease
particle%idomain(1) = 0
particle%iboundary(1) = 0
particle%idomain(2) = ic
particle%iboundary(2) = 0
particle%idomain(3) = 0
particle%iboundary(3) = 0
particle%itrdomain(LVL_MODEL) = 0
particle%iboundary(LVL_MODEL) = 0
particle%itrdomain(LVL_CELL) = ic
particle%iboundary(LVL_CELL) = 0
particle%itrdomain(LVL_SUBCELL) = 0
particle%iboundary(LVL_SUBCELL) = 0
particle%ifrctrn = this%ifrctrn
particle%iexmeth = this%iexmeth
particle%iextend = this%iextend
Expand Down
3 changes: 2 additions & 1 deletion src/Model/ParticleTracking/prt.f90
Original file line number Diff line number Diff line change
Expand Up @@ -405,6 +405,7 @@ end subroutine prt_cq

!> @brief Calculate particle mass storage
subroutine prt_cq_sto(this)
use ParticleModule, only: LVL_CELL
! modules
use TdisModule, only: delt
use PrtPrpModule, only: PrtPrpType
Expand Down Expand Up @@ -436,7 +437,7 @@ subroutine prt_cq_sto(this)
istatus = packobj%particles%istatus(np)
! this may need to change if istatus flags change
if ((istatus > 0) .and. (istatus /= 8)) then
n = packobj%particles%idomain(np, 2)
n = packobj%particles%itrdomain(np, LVL_CELL)
! Each particle currently assigned unit mass
this%masssto(n) = this%masssto(n) + DONE
end if
Expand Down
4 changes: 2 additions & 2 deletions src/Solution/ParticleTracker/MethodCellPassToBot.f90
Original file line number Diff line number Diff line change
Expand Up @@ -43,7 +43,7 @@ end subroutine deallocate

!> @brief Pass particle vertically and instantaneously to the cell bottom
subroutine apply_ptb(this, particle, tmax)
use ParticleModule, only: TERM_NO_EXITS
use ParticleModule, only: TERM_NO_EXITS, LVL_CELL
! dummy
class(MethodCellPassToBotType), intent(inout) :: this
type(ParticleType), pointer, intent(inout) :: particle
Expand All @@ -57,7 +57,7 @@ subroutine apply_ptb(this, particle, tmax)

! Pass to bottom face
particle%z = this%cell%defn%bot
particle%iboundary(2) = this%cell%defn%npolyverts + 2
particle%iboundary(LVL_CELL) = this%cell%defn%npolyverts + 2

! Terminate if in bottom layer
select type (dis => this%fmi%dis)
Expand Down
9 changes: 5 additions & 4 deletions src/Solution/ParticleTracker/MethodCellPollock.f90
Original file line number Diff line number Diff line change
Expand Up @@ -68,20 +68,21 @@ subroutine load_mcp(this, particle, next_level, submethod)
trackctl=this%trackctl, &
tracktimes=this%tracktimes)
submethod => method_subcell_plck
particle%idomain(next_level) = 1
particle%itrdomain(next_level) = 1
end subroutine load_mcp

!> @brief Having exited the lone subcell, pass the particle to the cell face
!! In this case the lone subcell is the cell.
subroutine pass_mcp(this, particle)
use ParticleModule, only: LVL_CELL, LVL_SUBCELL
! dummy
class(MethodCellPollockType), intent(inout) :: this
type(ParticleType), pointer, intent(inout) :: particle
! local
integer(I4B) :: exitface
integer(I4B) :: entryface

exitface = particle%iboundary(3)
exitface = particle%iboundary(LVL_SUBCELL)
! Map subcell exit face to cell face
select case (exitface) ! note: exitFace uses Dave's iface convention
case (0)
Expand All @@ -100,7 +101,7 @@ subroutine pass_mcp(this, particle)
entryface = 7
end select
if (entryface .eq. -1) then
particle%iboundary(2) = 0
particle%iboundary(LVL_CELL) = 0
else
if ((entryface .ge. 1) .and. (entryface .le. 4)) then
! Account for local cell rotation
Expand All @@ -110,7 +111,7 @@ subroutine pass_mcp(this, particle)
end select
if (entryface .gt. 4) entryface = entryface - 4
end if
particle%iboundary(2) = entryface
particle%iboundary(LVL_CELL) = entryface
end if
end subroutine pass_mcp

Expand Down
48 changes: 25 additions & 23 deletions src/Solution/ParticleTracker/MethodCellPollockQuad.f90
Original file line number Diff line number Diff line change
Expand Up @@ -71,6 +71,7 @@ end subroutine load_mcpq

!> @brief Pass particle to next subcell if there is one, or to the cell face
subroutine pass_mcpq(this, particle)
use ParticleModule, only: LVL_CELL, LVL_SUBCELL
! dummy
class(MethodCellPollockQuadType), intent(inout) :: this
type(ParticleType), pointer, intent(inout) :: particle
Expand All @@ -79,8 +80,8 @@ subroutine pass_mcpq(this, particle)

select type (cell => this%cell)
type is (CellRectQuadType)
exitFace = particle%iboundary(3)
isc = particle%idomain(3)
exitFace = particle%iboundary(LVL_SUBCELL)
isc = particle%itrdomain(LVL_SUBCELL)
npolyverts = cell%defn%npolyverts

! exitFace uses MODPATH 7 iface convention here
Expand All @@ -92,13 +93,13 @@ subroutine pass_mcpq(this, particle)
select case (isc)
case (1)
! W face, subcell 1 --> E face, subcell 4 (cell interior)
particle%idomain(3) = 4
particle%iboundary(3) = 2
particle%itrdomain(LVL_SUBCELL) = 4
particle%iboundary(LVL_SUBCELL) = 2
inface = 0 ! want Domain(2) unchanged; Boundary(2) = 0
case (2)
! W face, subcell 2 --> E face, subcell 3 (cell interior)
particle%idomain(3) = 3
particle%iboundary(3) = 2
particle%itrdomain(LVL_SUBCELL) = 3
particle%iboundary(LVL_SUBCELL) = 2
inface = 0 ! want Domain(2) unchanged; Boundary(2) = 0
case (3)
! W face, subcell 3 (cell face)
Expand All @@ -121,21 +122,21 @@ subroutine pass_mcpq(this, particle)
infaceoff = -1
case (3)
! E face, subcell 3 --> W face, subcell 2 (cell interior)
particle%idomain(3) = 2
particle%iboundary(3) = 1
particle%itrdomain(LVL_SUBCELL) = 2
particle%iboundary(LVL_SUBCELL) = 1
inface = 0 ! want Domain(2) unchanged; Boundary(2) = 0
case (4)
! E face, subcell 4 --> W face subcell 1 (cell interior)
particle%idomain(3) = 1
particle%iboundary(3) = 1
particle%itrdomain(LVL_SUBCELL) = 1
particle%iboundary(LVL_SUBCELL) = 1
inface = 0 ! want Domain(2) unchanged; Boundary(2) = 0
end select
case (3)
select case (isc)
case (1)
! S face, subcell 1 --> N face, subcell 2 (cell interior)
particle%idomain(3) = 2
particle%iboundary(3) = 4
particle%itrdomain(LVL_SUBCELL) = 2
particle%iboundary(LVL_SUBCELL) = 4
inface = 0 ! want Domain(2) unchanged; Boundary(2) = 0
case (2)
! S face, subcell 2 (cell face)
Expand All @@ -147,8 +148,8 @@ subroutine pass_mcpq(this, particle)
infaceoff = -1
case (4)
! S face, subcell 4 --> N face, subcell 3 (cell interior)
particle%idomain(3) = 3
particle%iboundary(3) = 4
particle%itrdomain(LVL_SUBCELL) = 3
particle%iboundary(LVL_SUBCELL) = 4
inface = 0 ! want Domain(2) unchanged; Boundary(2) = 0
end select
case (4)
Expand All @@ -159,13 +160,13 @@ subroutine pass_mcpq(this, particle)
infaceoff = -1
case (2)
! N face, subcell 2 --> S face, subcell 1 (cell interior)
particle%idomain(3) = 1
particle%iboundary(3) = 3
particle%itrdomain(LVL_SUBCELL) = 1
particle%iboundary(LVL_SUBCELL) = 3
inface = 0 ! want Domain(2) unchanged; Boundary(2) = 0
case (3)
! N face, subcell 3 --> S face, subcell 4 (cell interior)
particle%idomain(3) = 4
particle%iboundary(3) = 3
particle%itrdomain(LVL_SUBCELL) = 4
particle%iboundary(LVL_SUBCELL) = 3
inface = 0 ! want Domain(2) unchanged; Boundary(2) = 0
case (4)
! N face, subcell 4 (cell face)
Expand All @@ -181,9 +182,9 @@ subroutine pass_mcpq(this, particle)
end select

if (inface .eq. -1) then
particle%iboundary(2) = 0
particle%iboundary(LVL_CELL) = 0
else if (inface .eq. 0) then
particle%iboundary(2) = 0
particle%iboundary(LVL_CELL) = 0
else
if ((inface .ge. 1) .and. (inface .le. 4)) then
! Account for local cell rotation
Expand All @@ -192,7 +193,7 @@ subroutine pass_mcpq(this, particle)
inface = cell%irectvert(inface) + infaceoff
if (inface .lt. 1) inface = inface + npolyverts
end if
particle%iboundary(2) = inface
particle%iboundary(LVL_CELL) = inface
end if
end select
end subroutine pass_mcpq
Expand Down Expand Up @@ -234,6 +235,7 @@ end subroutine apply_mcpq

!> @brief Load the rectangular subcell from the rectangular cell
subroutine load_subcell(this, particle, subcell)
use ParticleModule, only: LVL_SUBCELL
! dummy
class(MethodCellPollockQuadType), intent(inout) :: this
type(ParticleType), pointer, intent(inout) :: particle
Expand All @@ -251,7 +253,7 @@ subroutine load_subcell(this, particle, subcell)
factor = factor / cell%defn%porosity
npolyverts = cell%defn%npolyverts

isc = particle%idomain(3)
isc = particle%itrdomain(LVL_SUBCELL)
! Subcells 1, 2, 3, and 4 are Pollock's subcells A, B, C, and D,
! respectively

Expand All @@ -277,7 +279,7 @@ subroutine load_subcell(this, particle, subcell)
end if

subcell%isubcell = isc
particle%idomain(3) = isc
particle%itrdomain(LVL_SUBCELL) = isc
end if
dx = 5d-1 * dx
dy = 5d-1 * dy
Expand Down
24 changes: 13 additions & 11 deletions src/Solution/ParticleTracker/MethodCellTernary.f90
Original file line number Diff line number Diff line change
Expand Up @@ -99,6 +99,7 @@ end subroutine load_mct

!> @brief Pass particle to next subcell if there is one, or to the cell face
subroutine pass_mct(this, particle)
use ParticleModule, only: LVL_CELL, LVL_SUBCELL
! dummy
class(MethodCellTernaryType), intent(inout) :: this
type(ParticleType), pointer, intent(inout) :: particle
Expand All @@ -107,8 +108,8 @@ subroutine pass_mct(this, particle)
integer(I4B) :: exitFace
integer(I4B) :: inface

exitFace = particle%iboundary(3)
isc = particle%idomain(3)
exitFace = particle%iboundary(LVL_SUBCELL)
isc = particle%itrdomain(LVL_SUBCELL)

select case (exitFace)
case (0)
Expand All @@ -122,15 +123,15 @@ subroutine pass_mct(this, particle)
! Subcell face --> next subcell in "cycle" (cell interior)
isc = isc + 1
if (isc .gt. this%nverts) isc = 1
particle%idomain(3) = isc
particle%iboundary(3) = 3
particle%itrdomain(LVL_SUBCELL) = isc
particle%iboundary(LVL_SUBCELL) = 3
inface = 0
case (3)
! Subcell face --> preceding subcell in "cycle" (cell interior)
isc = isc - 1
if (isc .lt. 1) isc = this%nverts
particle%idomain(3) = isc
particle%iboundary(3) = 2
particle%itrdomain(LVL_SUBCELL) = isc
particle%iboundary(LVL_SUBCELL) = 2
inface = 0
case (4)
! Subcell bottom (cell bottom)
Expand All @@ -141,11 +142,11 @@ subroutine pass_mct(this, particle)
end select

if (inface .eq. -1) then
particle%iboundary(2) = 0
particle%iboundary(LVL_CELL) = 0
else if (inface .eq. 0) then
particle%iboundary(2) = 0
particle%iboundary(LVL_CELL) = 0
else
particle%iboundary(2) = inface
particle%iboundary(LVL_CELL) = inface
end if
end subroutine pass_mct

Expand Down Expand Up @@ -248,6 +249,7 @@ end subroutine apply_mct

!> @brief Loads a triangular subcell from the polygonal cell
subroutine load_subcell(this, particle, subcell)
use ParticleModule, only: LVL_SUBCELL
! dummy
class(MethodCellTernaryType), intent(inout) :: this
type(ParticleType), pointer, intent(inout) :: particle
Expand Down Expand Up @@ -280,7 +282,7 @@ subroutine load_subcell(this, particle, subcell)
type is (CellPolyType)
ic = cell%defn%icell
subcell%icell = ic
isc = particle%idomain(3)
isc = particle%itrdomain(LVL_SUBCELL)
if (isc .le. 0) then
xi = particle%x
yi = particle%y
Expand Down Expand Up @@ -315,7 +317,7 @@ subroutine load_subcell(this, particle, subcell)

call pstop(1)
else
particle%idomain(3) = isc
particle%itrdomain(LVL_SUBCELL) = isc
end if
end if
subcell%isubcell = isc
Expand Down
Loading
Loading