Skip to content

Commit

Permalink
fix(prt): calculate global from local z before checking release point (
Browse files Browse the repository at this point in the history
…#1865)

* previously checked z coord in grid before local->global conversion
* unrelated: better routine names for particle and particle store
  • Loading branch information
wpbonelli authored Jun 11, 2024
1 parent f0a13d2 commit 1e294fb
Show file tree
Hide file tree
Showing 3 changed files with 34 additions and 32 deletions.
27 changes: 15 additions & 12 deletions src/Model/ParticleTracking/prt-prp.f90
Original file line number Diff line number Diff line change
Expand Up @@ -379,11 +379,20 @@ subroutine prp_ad(this)
np = this%nparticles + 1
this%nparticles = np

! -- Check release point is within the specified cell
! and not above/below grid top/bottom respectively
! -- Particle release point. Calculate global from local z if needed.
x = this%rptx(nps)
y = this%rpty(nps)
z = this%rptz(nps)
if (this%localz) then
top = this%fmi%dis%top(ic)
bot = this%fmi%dis%bot(ic)
hds = this%fmi%gwfhead(ic)
z = bot + this%rptz(nps) * (hds - bot)
else
z = this%rptz(nps)
end if

! -- Check release point is within the specified cell
! and not above/below grid top/bottom respectively
call this%dis%get_polyverts(ic, polyverts)
if (.not. point_in_polygon(x, y, polyverts)) then
write (errmsg, '(a,g0,a,g0,a,i0)') &
Expand Down Expand Up @@ -443,14 +452,7 @@ subroutine prp_ad(this)
end if
particle%x = x
particle%y = y
if (this%localz) then
top = this%fmi%dis%top(ic)
bot = this%fmi%dis%bot(ic)
hds = this%fmi%gwfhead(ic)
particle%z = bot + this%rptz(nps) * (hds - bot)
else
particle%z = this%rptz(nps)
end if
particle%z = z
particle%trelease = trelease
! Set stopping time to earlier of times specified by STOPTIME and STOPTRAVELTIME
if (this%stoptraveltime == huge(1d0)) then
Expand All @@ -470,7 +472,8 @@ subroutine prp_ad(this)
particle%iexmethod = this%iexmethod
particle%extol = this%extol

call this%particles%load_from_particle(particle, np)
! -- Persist particle to particle store
call this%particles%save_particle(particle, np)

! -- Accumulate mass release from this point
this%rptmass(nps) = this%rptmass(nps) + DONE
Expand Down
6 changes: 3 additions & 3 deletions src/Model/ParticleTracking/prt.f90
Original file line number Diff line number Diff line change
Expand Up @@ -926,8 +926,8 @@ subroutine prt_solve(this)
! -- Loop over particles in package
do np = 1, packobj%nparticles
! -- Load particle from storage
call particle%load_from_store(packobj%particles, &
this%id, iprp, np)
call particle%load_particle(packobj%particles, &
this%id, iprp, np)

! -- If particle is permanently unreleased, record its initial/terminal state
if (particle%istatus == 8) &
Expand All @@ -952,7 +952,7 @@ subroutine prt_solve(this)
call this%method%apply(particle, tmax)

! -- Update particle storage
call packobj%particles%load_from_particle(particle, np)
call packobj%particles%save_particle(particle, np)
end do
end select
end do
Expand Down
33 changes: 16 additions & 17 deletions src/Solution/ParticleTracker/Particle.f90
Original file line number Diff line number Diff line change
Expand Up @@ -64,7 +64,7 @@ module ParticleModule
real(DP), public :: extol !< tolerance for iterative solution of particle exit location and time in generalized Pollock's method
contains
procedure, public :: get_model_coords
procedure, public :: load_from_store
procedure, public :: load_particle
procedure, public :: transform => transform_coords
end type ParticleType

Expand Down Expand Up @@ -95,9 +95,9 @@ module ParticleModule
integer(I4B), dimension(:), pointer, contiguous :: iexmethod !< method for iterative solution of particle exit location and time in generalized Pollock's method
real(DP), dimension(:), pointer, contiguous :: extol !< tolerance for iterative solution of particle exit location and time in generalized Pollock's method
contains
procedure, public :: deallocate => deallocate_particle_store
procedure, public :: resize => resize_store
procedure, public :: load_from_particle
procedure, public :: deallocate
procedure, public :: resize
procedure, public :: save_particle
end type ParticleStoreType

contains
Expand Down Expand Up @@ -141,7 +141,7 @@ subroutine allocate_particle_store(this, np, mempath)
end subroutine allocate_particle_store

!> @brief Deallocate particle arrays
subroutine deallocate_particle_store(this, mempath)
subroutine deallocate (this, mempath)
class(ParticleStoreType), intent(inout) :: this !< store
character(*), intent(in) :: mempath !< path to memory

Expand All @@ -166,10 +166,10 @@ subroutine deallocate_particle_store(this, mempath)
call mem_deallocate(this%extol, 'PLEXTOL', mempath)
call mem_deallocate(this%idomain, 'PLIDOMAIN', mempath)
call mem_deallocate(this%iboundary, 'PLIBOUNDARY', mempath)
end subroutine deallocate_particle_store
end subroutine deallocate

!> @brief Reallocate particle arrays
subroutine resize_store(this, np, mempath)
subroutine resize(this, np, mempath)
! -- dummy
class(ParticleStoreType), intent(inout) :: this !< particle store
integer(I4B), intent(in) :: np !< number of particles
Expand Down Expand Up @@ -197,15 +197,14 @@ subroutine resize_store(this, np, mempath)
call mem_reallocate(this%extol, np, 'PLEXTOL', mempath)
call mem_reallocate(this%idomain, np, levelmax, 'PLIDOMAIN', mempath)
call mem_reallocate(this%iboundary, np, levelmax, 'PLIBOUNDARY', mempath)
end subroutine resize_store
end subroutine resize

!> @brief Initialize particle from particle list.
!> @brief Load a particle from the particle store.
!!
!! This routine is used to initialize a particle from the list
!! so it can be tracked by prt_solve. The particle's advancing
!! flag is set and local coordinate transformations are reset.
!! This routine is used to initialize a particle for tracking.
!! The advancing flag and coordinate transformation are reset.
!<
subroutine load_from_store(this, store, imdl, iprp, ip)
subroutine load_particle(this, store, imdl, iprp, ip)
class(ParticleType), intent(inout) :: this !< particle
type(ParticleStoreType), intent(in) :: store !< particle storage
integer(I4B), intent(in) :: imdl !< index of model particle originated in
Expand Down Expand Up @@ -239,10 +238,10 @@ subroutine load_from_store(this, store, imdl, iprp, ip)
this%ifrctrn = store%ifrctrn(ip)
this%iexmethod = store%iexmethod(ip)
this%extol = store%extol(ip)
end subroutine load_from_store
end subroutine load_particle

!> @brief Update particle store from particle
subroutine load_from_particle(this, particle, ip)
!> @brief Save a particle's state to the particle store
subroutine save_particle(this, particle, ip)
class(ParticleStoreType), intent(inout) :: this !< particle storage
type(ParticleType), intent(in) :: particle !< particle
integer(I4B), intent(in) :: ip !< particle index
Expand Down Expand Up @@ -274,7 +273,7 @@ subroutine load_from_particle(this, particle, ip)
this%ifrctrn = particle%ifrctrn
this%iexmethod = particle%iexmethod
this%extol = particle%extol
end subroutine load_from_particle
end subroutine save_particle

!> @brief Apply the given global-to-local transformation to the particle.
subroutine transform_coords(this, xorigin, yorigin, zorigin, &
Expand Down

0 comments on commit 1e294fb

Please sign in to comment.