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

fix(prt): calculate global from local z before checking release point #1865

Merged
merged 1 commit into from
Jun 11, 2024
Merged
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
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
Loading