From 7fc574959a6bacb5a54357c3431015c03f6176e7 Mon Sep 17 00:00:00 2001 From: wpbonelli Date: Tue, 11 Jun 2024 17:01:16 -0400 Subject: [PATCH] fix(prt): calculate global from local z before checking release point * previously checked that z coordinate was within grid before local->global conversion --- src/Model/ParticleTracking/prt-prp.f90 | 27 ++++++++++--------- src/Model/ParticleTracking/prt.f90 | 6 ++--- src/Solution/ParticleTracker/Particle.f90 | 33 +++++++++++------------ 3 files changed, 34 insertions(+), 32 deletions(-) diff --git a/src/Model/ParticleTracking/prt-prp.f90 b/src/Model/ParticleTracking/prt-prp.f90 index 45c6417f7d3..558788f647c 100644 --- a/src/Model/ParticleTracking/prt-prp.f90 +++ b/src/Model/ParticleTracking/prt-prp.f90 @@ -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)') & @@ -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 @@ -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 diff --git a/src/Model/ParticleTracking/prt.f90 b/src/Model/ParticleTracking/prt.f90 index 5263d88c7b6..3dbdab59937 100644 --- a/src/Model/ParticleTracking/prt.f90 +++ b/src/Model/ParticleTracking/prt.f90 @@ -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) & @@ -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 diff --git a/src/Solution/ParticleTracker/Particle.f90 b/src/Solution/ParticleTracker/Particle.f90 index 0d13dac96be..fb35c29f296 100644 --- a/src/Solution/ParticleTracker/Particle.f90 +++ b/src/Solution/ParticleTracker/Particle.f90 @@ -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 @@ -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 @@ -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 @@ -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 @@ -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 @@ -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 @@ -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, &