Skip to content

refactor(prt): simplify face numbering scheme #2586

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.

Sign up for GitHub

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

Open
wants to merge 1 commit into
base: develop
Choose a base branch
from
+58 −44
Open

refactor(prt): simplify face numbering scheme #2586

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
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
78 changes: 48 additions & 30 deletions src/Model/ParticleTracking/prt-fmi.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -13,16 +13,23 @@ module PrtFmiModule
private
public :: PrtFmiType
public :: fmi_cr
public :: IFLOWFACE_TOP, IFLOWFACE_BOTTOM

character(len=LENPACKAGENAME) :: text = ' PRTFMI'

!> @brief IFLOWFACE numbers for top and bottom faces
enum, bind(c)
enumerator :: IFLOWFACE_TOP = -1
enumerator :: IFLOWFACE_BOTTOM = -2
end enum

type, extends(FlowModelInterfaceType) :: PrtFmiType
private
integer(I4B), public :: max_faces !< maximum number of faces for grid cell polygons
integer(I4B), public :: max_faces !< maximum number of 3d cell faces
real(DP), allocatable, public :: SourceFlows(:) ! cell source flows array
real(DP), allocatable, public :: SinkFlows(:) ! cell sink flows array
real(DP), allocatable, public :: StorageFlows(:) ! cell storage flows array
real(DP), allocatable, public :: BoundaryFlows(:) ! cell boundary flows array
real(DP), allocatable, public :: BoundaryFlows(:, :) ! cell boundary flows array
integer(I4B), allocatable, public :: BoundaryFaces(:) ! bitmask of assigned boundary faces

contains
Expand All @@ -34,6 +41,7 @@ module PrtFmiModule
procedure :: is_boundary_face
procedure :: is_net_out_boundary_face
procedure, private :: iflowface_to_icellface
procedure, private :: icellface_to_iflowface

end type PrtFmiType

Expand Down Expand Up @@ -144,7 +152,7 @@ subroutine prtfmi_df(this, dis, idryinactive)
allocate (this%StorageFlows(this%dis%nodes))
allocate (this%SourceFlows(this%dis%nodes))
allocate (this%SinkFlows(this%dis%nodes))
allocate (this%BoundaryFlows(this%dis%nodes * this%max_faces))
allocate (this%BoundaryFlows(this%dis%nodes, this%max_faces))
allocate (this%BoundaryFaces(this%dis%nodes))

end subroutine prtfmi_df
Expand All @@ -155,7 +163,7 @@ subroutine accumulate_flows(this)
class(PrtFmiType) :: this
! local
integer(I4B) :: j, i, ip, ib
integer(I4B) :: ioffset, iflowface, iauxiflowface, icellface
integer(I4B) :: iflowface, iauxiflowface, icellface
real(DP) :: qbnd
character(len=LENAUXNAME) :: auxname
integer(I4B) :: naux
Expand Down Expand Up @@ -196,9 +204,8 @@ subroutine accumulate_flows(this)
end if
if (icellface > 0) then
call this%mark_boundary_face(i, icellface)
ioffset = (i - 1) * this%max_faces
this%BoundaryFlows(ioffset + icellface) = &
this%BoundaryFlows(ioffset + icellface) + qbnd
this%BoundaryFlows(i, icellface) = &
this%BoundaryFlows(i, icellface) + qbnd
else if (qbnd .gt. DZERO) then
this%SourceFlows(i) = this%SourceFlows(i) + qbnd
else if (qbnd .lt. DZERO) then
Expand All @@ -210,10 +217,10 @@ subroutine accumulate_flows(this)
end subroutine accumulate_flows

!> @brief Mark a face as a boundary face.
subroutine mark_boundary_face(this, ic, iface)
subroutine mark_boundary_face(this, ic, icellface)
class(PrtFmiType) :: this
integer(I4B), intent(in) :: ic !< node number (reduced)
integer(I4B), intent(in) :: iface !< cell face number
integer(I4B), intent(in) :: icellface !< cell face number
! local
integer(I4B) :: bit_pos

Expand All @@ -222,25 +229,25 @@ subroutine mark_boundary_face(this, ic, iface)
print *, 'Expected a value in range [1, ', this%dis%nodes, ']'
call pstop(1)
end if
if (iface <= 0) then
print *, 'Invalid face number: ', iface
if (icellface <= 0) then
print *, 'Invalid face number: ', icellface
print *, 'Expected a value in range [1, ', this%max_faces, ']'
call pstop(1)
end if
bit_pos = iface - 1 ! bit position 0-based
bit_pos = icellface - 1 ! bit position 0-based
if (bit_pos < 0 .or. bit_pos > 31) then
print *, 'Invalid bitmask position: ', iface
print *, 'Invalid bitmask position: ', bit_pos
print *, 'Expected a value in range [0, 31]'
call pstop(1)
end if
this%BoundaryFaces(ic) = ibset(this%BoundaryFaces(ic), bit_pos)
end subroutine mark_boundary_face

!> @brief Check if a face is assigned to a boundary package.
function is_boundary_face(this, ic, iface) result(is_boundary)
function is_boundary_face(this, ic, icellface) result(is_boundary)
class(PrtFmiType) :: this
integer(I4B), intent(in) :: ic !< node number (reduced)
integer(I4B), intent(in) :: iface !< cell face number
integer(I4B), intent(in) :: icellface !< cell face number
logical(LGP) :: is_boundary
! local
integer(I4B) :: bit_pos
Expand All @@ -251,44 +258,55 @@ function is_boundary_face(this, ic, iface) result(is_boundary)
print *, 'Expected a value in range [1, ', this%dis%nodes, ']'
call pstop(1)
end if
if (iface <= 0) then
print *, 'Invalid face number: ', iface
if (icellface <= 0) then
print *, 'Invalid face number: ', icellface
print *, 'Expected a value in range [1, ', this%max_faces, ']'
call pstop(1)
end if
bit_pos = iface - 1 ! bit position 0-based
bit_pos = icellface - 1 ! bit position 0-based
if (bit_pos < 0 .or. bit_pos > 31) then
print *, 'Invalid bitmask position: ', iface
print *, 'Invalid bitmask position: ', bit_pos
print *, 'Expected a value in range [0, 31]'
call pstop(1)
end if
is_boundary = btest(this%BoundaryFaces(ic), bit_pos)
end function is_boundary_face

!> @brief Check if a face is an assigned boundary with net outflow.
function is_net_out_boundary_face(this, ic, iface) result(is_net_out_boundary)
function is_net_out_boundary_face(this, ic, icellface) &
result(is_net_out_boundary)
class(PrtFmiType) :: this
integer(I4B), intent(in) :: ic !< node number (reduced)
integer(I4B), intent(in) :: iface !< cell face number
integer(I4B), intent(in) :: icellface !< cell face number
logical(LGP) :: is_net_out_boundary
! local
integer(I4B) :: ioffset

is_net_out_boundary = .false.
if (.not. this%is_boundary_face(ic, iface)) return
ioffset = (ic - 1) * this%max_faces
if (this%BoundaryFlows(ioffset + iface) < DZERO) is_net_out_boundary = .true.
if (.not. this%is_boundary_face(ic, icellface)) return
if (this%BoundaryFlows(ic, icellface) < DZERO) &
is_net_out_boundary = .true.
end function is_net_out_boundary_face

!> @brief Convert an iflowface number to a cell face number.
!! Maps bottom (-2) -> max_faces - 1, top (-1) -> max_faces.
function iflowface_to_icellface(this, iflowface) result(iface)
function iflowface_to_icellface(this, iflowface) result(icellface)
class(PrtFmiType), intent(inout) :: this
integer(I4B), intent(in) :: iflowface
integer(I4B) :: iface
integer(I4B) :: icellface

iface = iflowface
if (iface < 0) iface = iface + this%max_faces + 1
icellface = iflowface
if (icellface < 0) icellface = icellface + this%max_faces - IFLOWFACE_TOP
end function iflowface_to_icellface

!> @brief Convert a cell face number to an iflowface number.
!! Maps bottom (max_faces - 1) -> -2, top (max_faces) -> -1.
function icellface_to_iflowface(this, icellface) result(iflowface)
class(PrtFmiType), intent(inout) :: this
integer(I4B), intent(in) :: icellface
integer(I4B) :: iflowface

iflowface = icellface
if (iflowface > this%max_faces + IFLOWFACE_BOTTOM) &
iflowface = iflowface - this%max_faces + IFLOWFACE_TOP
end function icellface_to_iflowface

end module PrtFmiModule
14 changes: 6 additions & 8 deletions src/Solution/ParticleTracker/Method/MethodDis.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -490,22 +490,20 @@ subroutine load_cell_boundary_flows(this, defn)
class(MethodDisType), intent(inout) :: this
type(CellDefnType), pointer, intent(inout) :: defn
! local
integer(I4B) :: ioffset

ioffset = (defn%icell - 1) * this%fmi%max_faces
defn%faceflow(1) = defn%faceflow(1) + &
this%fmi%BoundaryFlows(ioffset + 1)
this%fmi%BoundaryFlows(defn%icell, 1)
defn%faceflow(2) = defn%faceflow(2) + &
this%fmi%BoundaryFlows(ioffset + 2)
this%fmi%BoundaryFlows(defn%icell, 2)
defn%faceflow(3) = defn%faceflow(3) + &
this%fmi%BoundaryFlows(ioffset + 3)
this%fmi%BoundaryFlows(defn%icell, 3)
defn%faceflow(4) = defn%faceflow(4) + &
this%fmi%BoundaryFlows(ioffset + 4)
this%fmi%BoundaryFlows(defn%icell, 4)
defn%faceflow(5) = defn%faceflow(1)
defn%faceflow(6) = defn%faceflow(6) + &
this%fmi%BoundaryFlows(ioffset + this%fmi%max_faces - 1)
this%fmi%BoundaryFlows(defn%icell, this%fmi%max_faces - 1)
defn%faceflow(7) = defn%faceflow(7) + &
this%fmi%BoundaryFlows(ioffset + this%fmi%max_faces)
this%fmi%BoundaryFlows(defn%icell, this%fmi%max_faces)
end subroutine load_cell_boundary_flows

end module MethodDisModule
10 changes: 4 additions & 6 deletions src/Solution/ParticleTracker/Method/MethodDisv.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -506,24 +506,22 @@ subroutine load_cell_boundary_flows(this, defn)
type(CellDefnType), pointer, intent(inout) :: defn

! local
integer(I4B) :: ic, iv, ioffset, npolyverts, max_faces
integer(I4B) :: ic, iv, npolyverts

ic = defn%icell
npolyverts = defn%npolyverts
max_faces = this%fmi%max_faces
ioffset = (ic - 1) * max_faces
do iv = 1, npolyverts
defn%faceflow(iv) = &
defn%faceflow(iv) + &
this%fmi%BoundaryFlows(ioffset + iv)
this%fmi%BoundaryFlows(ic, iv)
end do
defn%faceflow(npolyverts + 1) = defn%faceflow(1)
defn%faceflow(npolyverts + 2) = &
defn%faceflow(npolyverts + 2) + &
this%fmi%BoundaryFlows(ioffset + max_faces - 1)
this%fmi%BoundaryFlows(ic, this%fmi%max_faces - 1)
defn%faceflow(npolyverts + 3) = &
defn%faceflow(npolyverts + 3) + &
this%fmi%BoundaryFlows(ioffset + max_faces)
this%fmi%BoundaryFlows(ic, this%fmi%max_faces)

end subroutine load_cell_boundary_flows

Expand Down
Loading