Skip to content

Commit

Permalink
various renaming, refactoring/cleanup
Browse files Browse the repository at this point in the history
  • Loading branch information
@wpbonelli
wpbonelli committed Dec 13, 2023
1 parent 0d88824 commit 98eba93
Show file tree
Hide file tree
Showing 4 changed files with 186 additions and 192 deletions.
181 changes: 91 additions & 90 deletions src/Solution/ParticleTracker/MethodSubcellPollock.f90
View file
Original file line number Diff line number Diff line change
@@ -1,5 +1,6 @@
module MethodSubcellPollockModule
use KindModule, only: DP, I4B
use KindModule, only: DP, I4B, LGP
use ErrorUtilModule, only: pstop
use MethodModule
use SubcellRectModule
use ParticleModule
Expand All @@ -8,37 +9,35 @@ module MethodSubcellPollockModule
private
public :: MethodSubcellPollockType
public :: create_methodSubcellPollock
public :: CalculateDT
public :: calculate_dt

! -- Extend MethodType to the Pollock's subcell-method type (MethodSubcellPollockType)
!> @brief Rectangular subcell tracking method (reduces to Pollock's method for MODPATH)
type, extends(MethodType) :: MethodSubcellPollockType
private
type(SubcellRectType), pointer, public :: subcellRect => null() ! tracking domain for the method
double precision, allocatable, public :: qextl1(:), qextl2(:), qintl(:) ! external and internal subcell flows for the cell
type(SubcellRectType), pointer, public :: subcell => null() !< tracking domain
double precision, allocatable, public :: qextl1(:), qextl2(:), qintl(:) !< external and internal subcell flows
contains
procedure, public :: destroy ! destructor for the method
procedure, public :: init ! initializes the method
procedure, public :: apply => apply_mSP ! applies Pollock's subcell method (tracks particle)
procedure, private :: track_sub
procedure, public :: destroy !< destructor
procedure, public :: init !< initializer
procedure, public :: apply => apply_msp !< apply the method
procedure, private :: track_sub !< track over the subdomain
end type MethodSubcellPollockType

contains

!> @brief Create a new Pollock's subcell-method object
subroutine create_methodSubcellPollock(methodSubcellPollock)
type(MethodSubcellPollockType), pointer :: methodSubcellPollock

allocate (methodSubcellPollock)
subroutine create_methodSubcellPollock(method)
type(MethodSubcellPollockType), pointer :: method

! -- This method does not delegate tracking to a submethod; it performs
! -- actual particle-tracking calculations
methodSubcellPollock%delegatesTracking = .false.
allocate (method)
method%delegatesTracking = .false.

! -- Create tracking domain for this method and set trackingDomain pointer
call create_subcellRect(methodSubcellPollock%subcellRect)
methodSubcellPollock%trackingDomainType => &
methodSubcellPollock%subcellRect%type

call create_subcellRect(method%subcell)
method%trackingDomainType => &
method%subcell%type
end subroutine create_methodSubcellPollock

!> @brief Destructor for a Pollock's subcell-method object
Expand All @@ -48,18 +47,16 @@ subroutine destroy(this)
end subroutine destroy

!> @brief Initialize a Pollock's subcell-method object
subroutine init(this, subcellRect, trackctl)
subroutine init(this, subcell, trackctl)
class(MethodSubcellPollockType), intent(inout) :: this
type(SubcellRectType), pointer :: subcellRect
type(SubcellRectType), pointer :: subcell
type(TrackControlType), pointer :: trackctl

this%subcellRect => subcellRect
this%subcell => subcell
this%trackctl => trackctl

end subroutine init

!> @brief Apply Pollock's method to a rectangular subcell
subroutine apply_mSP(this, particle, tmax)
subroutine apply_msp(this, particle, tmax)
! -- dummy
class(MethodSubcellPollockType), intent(inout) :: this
type(ParticleType), pointer, intent(inout) :: particle
Expand All @@ -68,24 +65,24 @@ subroutine apply_mSP(this, particle, tmax)
double precision :: xOrigin, yOrigin, zOrigin, sinrot, cosrot

! -- Transform particle location into local subcell coordinates
xOrigin = this%subcellRect%xOrigin
yOrigin = this%subcellRect%yOrigin
zOrigin = this%subcellRect%zOrigin
sinrot = this%subcellRect%sinrot
cosrot = this%subcellRect%cosrot
xOrigin = this%subcell%xOrigin
yOrigin = this%subcell%yOrigin
zOrigin = this%subcell%zOrigin
sinrot = this%subcell%sinrot
cosrot = this%subcell%cosrot

! -- sinrot and cosrot should be 0. and 1., respectively, i.e., there
! -- is no rotation. Also no z translation; only x and y translations.
call particle%transf_coords(xOrigin, yOrigin)

! -- Track the particle across the subcell
! call track_sub(this%subcellRect,particle,initialTime,maximumTime,t)
call this%track_sub(this%subcellRect, particle, tmax)
call this%track_sub(this%subcell, particle, tmax)

! -- Transform particle location back to local cell coordinates
call particle%transf_coords(xOrigin, yOrigin, invert=.true.)

end subroutine apply_mSP
end subroutine apply_msp

!> @brief Track a particle across a rectangular subcell using Pollock's method
!!
Expand All @@ -94,13 +91,13 @@ end subroutine apply_mSP
!! code are responsible for its appropriate application in this context
!! and for any modifications or errors.
!<
subroutine track_sub(this, subcellRect, particle, tmax) ! kluge note: rename???
subroutine track_sub(this, subcell, particle, tmax)
! modules
use ParticleModule, only: get_particle_id
use TdisModule, only: kper, kstp
! dummy
class(MethodSubcellPollockType), intent(inout) :: this
class(SubcellRectType), intent(in) :: subcellRect
class(SubcellRectType), intent(in) :: subcell
type(ParticleType), pointer, intent(inout) :: particle
real(DP), intent(in) :: tmax
! local
Expand All @@ -116,25 +113,25 @@ subroutine track_sub(this, subcellRect, particle, tmax) ! kluge note: rename???
reason = -1

! -- Initial particle location in scaled subcell coordinates
initialX = particle%x / subcellRect%dx
initialY = particle%y / subcellRect%dy
initialZ = particle%z / subcellRect%dz
initialX = particle%x / subcell%dx
initialY = particle%y / subcell%dy
initialZ = particle%z / subcell%dz

! -- Make local copies of face velocities for convenience
vx1 = subcellRect%vx1
vx2 = subcellRect%vx2
vy1 = subcellRect%vy1
vy2 = subcellRect%vy2
vz1 = subcellRect%vz1
vz2 = subcellRect%vz2
vx1 = subcell%vx1
vx2 = subcell%vx2
vy1 = subcell%vy1
vy2 = subcell%vy2
vz1 = subcell%vz1
vz2 = subcell%vz2

! -- Compute time of travel to each possible exit face
statusVX = CalculateDT(vx1, vx2, subcellRect%dx, initialX, vx, dvxdx, &
dtexitx)
statusVY = CalculateDT(vy1, vy2, subcellRect%dy, initialY, vy, dvydy, &
dtexity)
statusVZ = CalculateDT(vz1, vz2, subcellRect%dz, initialZ, vz, dvzdz, &
dtexitz)
statusVX = calculate_dt(vx1, vx2, subcell%dx, initialX, vx, dvxdx, &
dtexitx)
statusVY = calculate_dt(vy1, vy2, subcell%dy, initialY, vy, dvydy, &
dtexity)
statusVZ = calculate_dt(vz1, vz2, subcell%dz, initialZ, vz, dvzdz, &
dtexitz)

! -- Check for no exit face
if ((statusVX .eq. 3) .and. (statusVY .eq. 3) .and. (statusVZ .eq. 3)) then
Expand All @@ -143,14 +140,10 @@ subroutine track_sub(this, subcellRect, particle, tmax) ! kluge note: rename???
! particle%istatus = 5
! return

! kluge note: contact the developer situation
print *, "======================================"
print *, "Subcell with no exit face" ! kluge
print *, "Particle ", get_particle_id(particle)
print *, "Cell ", particle%iTrackingDomain(2)
print *, "======================================"
!!pause
stop
! contact the developer situation (for now? always?)
print *, "Subcell with no exit face: particle", get_particle_id(particle), &
"cell", particle%iTrackingDomain(2)
call pstop(1)
end if

! -- Determine (earliest) exit face and corresponding travel time to exit
Expand Down Expand Up @@ -194,9 +187,9 @@ subroutine track_sub(this, subcellRect, particle, tmax) ! kluge note: rename???
! -- calculate particle location at that final time.
t = tmax
dt = t - particle%ttrack
x = NewXYZ(vx, dvxdx, vx1, vx2, dt, initialX, subcellRect%dx, statusVX)
y = NewXYZ(vy, dvydy, vy1, vy2, dt, initialY, subcellRect%dy, statusVY)
z = NewXYZ(vz, dvzdz, vz1, vz2, dt, initialZ, subcellRect%dz, statusVZ)
x = new_x(vx, dvxdx, vx1, vx2, dt, initialX, subcell%dx, statusVX == 1)
y = new_x(vy, dvydy, vy1, vy2, dt, initialY, subcell%dy, statusVY == 1)
z = new_x(vz, dvzdz, vz1, vz2, dt, initialZ, subcell%dz, statusVZ == 1)
exitFace = 0
particle%istatus = 1
particle%advancing = .false.
Expand All @@ -209,31 +202,31 @@ subroutine track_sub(this, subcellRect, particle, tmax) ! kluge note: rename???
dt = dtexit
if ((exitFace .eq. 1) .or. (exitFace .eq. 2)) then
x = 0d0
y = NewXYZ(vy, dvydy, vy1, vy2, dt, initialY, subcellRect%dy, statusVY)
z = NewXYZ(vz, dvzdz, vz1, vz2, dt, initialZ, subcellRect%dz, statusVZ)
y = new_x(vy, dvydy, vy1, vy2, dt, initialY, subcell%dy, statusVY == 1)
z = new_x(vz, dvzdz, vz1, vz2, dt, initialZ, subcell%dz, statusVZ == 1)
if (exitFace .eq. 2) x = 1.0d0
else if ((exitFace .eq. 3) .or. (exitFace .eq. 4)) then
x = NewXYZ(vx, dvxdx, vx1, vx2, dt, initialX, subcellRect%dx, statusVX)
x = new_x(vx, dvxdx, vx1, vx2, dt, initialX, subcell%dx, statusVX == 1)
y = 0d0
z = NewXYZ(vz, dvzdz, vz1, vz2, dt, initialZ, subcellRect%dz, statusVZ)
z = new_x(vz, dvzdz, vz1, vz2, dt, initialZ, subcell%dz, statusVZ == 1)
if (exitFace .eq. 4) y = 1.0d0
else if ((exitFace .eq. 5) .or. (exitFace .eq. 6)) then
x = NewXYZ(vx, dvxdx, vx1, vx2, dt, initialX, subcellRect%dx, statusVX)
y = NewXYZ(vy, dvydy, vy1, vy2, dt, initialY, subcellRect%dy, statusVY)
x = new_x(vx, dvxdx, vx1, vx2, dt, initialX, subcell%dx, statusVX == 1)
y = new_x(vy, dvydy, vy1, vy2, dt, initialY, subcell%dy, statusVY == 1)
z = 0d0
if (exitFace .eq. 6) z = 1.0d0
else
print *, "something went wrong in track_sub" ! kluge
stop ! kluge
print *, "programmer error, invalid exit face", exitFace
call pstop(1)
end if
reason = 1 ! cell transition
end if

! -- Set final particle location in local (unscaled) subcell coordinates,
! -- final time for particle trajectory, and exit face
particle%x = x * subcellRect%dx
particle%y = y * subcellRect%dy
particle%z = z * subcellRect%dz
particle%x = x * subcell%dx
particle%y = y * subcell%dy
particle%z = z * subcell%dz
particle%ttrack = t
particle%iTrackingDomainBoundary(3) = exitFace

Expand All @@ -251,7 +244,7 @@ end subroutine track_sub
!! code are responsible for its appropriate application in this context
!! and for any modifications or errors.
!<
function CalculateDT(v1, v2, dx, xL, v, dvdx, dt) result(status)
function calculate_dt(v1, v2, dx, xL, v, dvdx, dt) result(status)
! dummy
doubleprecision, intent(in) :: v1, v2, dx, xL
doubleprecision, intent(inout) :: v, dvdx, dt
Expand Down Expand Up @@ -359,34 +352,42 @@ function CalculateDT(v1, v2, dx, xL, v, dvdx, dt) result(status)
dt = log(vr) / dvdx
status = 0

end function CalculateDT
end function calculate_dt

!> @brief Update particle coordinates based on time increment
!> @brief Update a cell-local coordinate based on a time increment.
!!
!! This subroutine consists partly or entirely of code written by
!! David W. Pollock of the USGS for MODPATH 7. The authors of the present
!! code are responsible for its appropriate application in this context
!! and for any modifications or errors.
!<
function NewXYZ(v, dvdx, v1, v2, dt, x, dx, velocityProfileStatus) result(newX)
pure function new_x(v, dvdx, v1, v2, dt, x, dx, velocity_profile) result(newx)
! dummy
integer, intent(in) :: velocityProfileStatus
doubleprecision, intent(in) :: v, dvdx, v1, v2, dt, x, dx
real(DP), intent(in) :: v, dvdx, v1, v2, dt, x, dx
logical(LGP), intent(in), optional :: velocity_profile
! result
doubleprecision :: newX

newX = x
select case (velocityProfileStatus)
case (1)
newX = newX + (v1 * dt / dx)
case default
if (v .ne. 0d0) then
newX = newX + (v * (exp(dvdx * dt) - 1.0d0) / dvdx / dx)
end if
end select
if (newX .lt. 0d0) newX = 0d0
if (newX .gt. 1.0d0) newX = 1.0d0
real(DP) :: newx
logical(LGP) :: lprofile

! -- process optional arguments
if (present(velocity_profile)) then
lprofile = velocity_profile
else
lprofile = .false.
end if

! -- recompute coordinate
newx = x
if (lprofile) then
newx = newx + (v1 * dt / dx)
else if (v .ne. 0d0) then
newx = newx + (v * (exp(dvdx * dt) - 1.0d0) / dvdx / dx)
end if

! -- clamp to [0, 1]
if (newx .lt. 0d0) newx = 0d0
if (newx .gt. 1.0d0) newx = 1.0d0

end function NewXYZ
end function new_x

end module MethodSubcellPollockModule
Loading

0 comments on commit 98eba93

Please sign in to comment.