diff --git a/src/Model/ParticleTracking/prt-prp.f90 b/src/Model/ParticleTracking/prt-prp.f90 index 7c747915763..08d3bd5384d 100644 --- a/src/Model/ParticleTracking/prt-prp.f90 +++ b/src/Model/ParticleTracking/prt-prp.f90 @@ -430,6 +430,8 @@ subroutine release(this, ip, trelease) ! Accumulate mass for release point this%rptm(ip) = this%rptm(ip) + DONE + print *, "released particle: ", ip + end subroutine release subroutine initialize_particle(this, particle, ip, trelease) diff --git a/src/Model/ParticleTracking/prt.f90 b/src/Model/ParticleTracking/prt.f90 index 0730d8f6b02..d4e7420a83d 100644 --- a/src/Model/ParticleTracking/prt.f90 +++ b/src/Model/ParticleTracking/prt.f90 @@ -911,17 +911,19 @@ subroutine prt_solve(this) ! -- Update PRP index iprp = iprp + 1 - ! -- Initialize PRP-specific track files, if enabled - if (packobj%itrkout > 0) then - call this%trackctl%init_track_file( & - packobj%itrkout, & - iprp=iprp) - end if - if (packobj%itrkcsv > 0) then - call this%trackctl%init_track_file( & - packobj%itrkcsv, & - csv=.true., & - iprp=iprp) + ! -- Initialize PRP-specific track files + if (kper == 1 .and. kstp == 1) then + if (packobj%itrkout > 0) then + call this%trackctl%init_track_file( & + packobj%itrkout, & + iprp=iprp) + end if + if (packobj%itrkcsv > 0) then + call this%trackctl%init_track_file( & + packobj%itrkcsv, & + csv=.true., & + iprp=iprp) + end if end if ! -- Loop over particles in package @@ -957,8 +959,9 @@ subroutine prt_solve(this) ! Get and apply the tracking method call this%method%apply(particle, tmax) - ! Reset previous cell, exit face, and zone numbers + ! Reset previous cell, subcell, and zone numbers particle%icp = 0 + particle%iscp = 0 particle%izp = 0 ! Update particle storage diff --git a/src/Solution/ParticleTracker/Method.f90 b/src/Solution/ParticleTracker/Method.f90 index 061a627e1ad..1325622746c 100644 --- a/src/Solution/ParticleTracker/Method.f90 +++ b/src/Solution/ParticleTracker/Method.f90 @@ -184,6 +184,9 @@ subroutine save(this, particle, reason) end if call this%trackctl%save(particle, kper=per, kstp=stp, reason=reason) + print *, "saving" + ! if (kper == 3) stop + end subroutine save !> @brief Check reporting/terminating conditions before tracking. diff --git a/src/Solution/ParticleTracker/MethodCellTernary.f90 b/src/Solution/ParticleTracker/MethodCellTernary.f90 index a04a65d231f..eef80b27cdb 100644 --- a/src/Solution/ParticleTracker/MethodCellTernary.f90 +++ b/src/Solution/ParticleTracker/MethodCellTernary.f90 @@ -85,7 +85,20 @@ subroutine load_mct(this, particle, next_level, submethod) select type (subcell => this%subcell) type is (SubcellTriType) call this%load_subcell(particle, subcell) + + ! don't allow a particle to return to the + ! subcell it was just in + if (subcell%isubcell == particle%iscp) then + particle%advancing = .false. + particle%idomain(3) = particle%iscp + particle%istatus = 2 + call this%save(particle, reason=3) + return + else + particle%iscp = particle%idomain(3) + end if end select + call method_subcell_tern%init( & fmi=this%fmi, & cell=this%cell, & diff --git a/src/Solution/ParticleTracker/MethodSubcellTernary.f90 b/src/Solution/ParticleTracker/MethodSubcellTernary.f90 index 97637b4c59b..69c4381833e 100644 --- a/src/Solution/ParticleTracker/MethodSubcellTernary.f90 +++ b/src/Solution/ParticleTracker/MethodSubcellTernary.f90 @@ -265,6 +265,7 @@ subroutine track_subcell(this, subcell, particle, tmax) particle%istatus = 1 particle%advancing = .false. reason = 2 ! timestep end + print *, "end of time step" else ! -- The computed exit time is less than or equal to the maximum time, ! -- so set final time for particle trajectory equal to exit time. diff --git a/src/Solution/ParticleTracker/Particle.f90 b/src/Solution/ParticleTracker/Particle.f90 index 347ae4b1a2b..014c21897e7 100644 --- a/src/Solution/ParticleTracker/Particle.f90 +++ b/src/Solution/ParticleTracker/Particle.f90 @@ -45,6 +45,7 @@ module ParticleModule integer(I4B), allocatable, public :: idomain(:) !< tracking domain hierarchy ! TODO: rename to itdomain? idomain integer(I4B), allocatable, public :: iboundary(:) !< tracking domain boundaries integer(I4B), public :: icp !< previous cell number (reduced) + integer(I4B), public :: iscp !< previous subcell number integer(I4B), public :: icu !< user cell number integer(I4B), public :: ilay !< grid layer integer(I4B), public :: izone !< current zone number @@ -90,6 +91,7 @@ module ParticleModule integer(I4B), dimension(:, :), pointer, public, contiguous :: idomain !< array of indices for domains in the tracking domain hierarchy integer(I4B), dimension(:, :), pointer, public, contiguous :: iboundary !< array of indices for tracking domain boundaries integer(I4B), dimension(:), pointer, public, contiguous :: icp !< previous cell number (reduced) + integer(I4B), dimension(:), pointer, public, contiguous :: iscp !< previous subcell number integer(I4B), dimension(:), pointer, public, contiguous :: icu !< cell number (user) integer(I4B), dimension(:), pointer, public, contiguous :: ilay !< layer integer(I4B), dimension(:), pointer, public, contiguous :: izone !< current zone number @@ -134,6 +136,7 @@ subroutine allocate_particle_store(this, np, mempath) call mem_allocate(this%iprp, np, 'PLIPRP', mempath) call mem_allocate(this%name, LENBOUNDNAME, np, 'PLNAME', mempath) call mem_allocate(this%icp, np, 'PLICP', mempath) + call mem_allocate(this%iscp, np, 'PLISCP', mempath) call mem_allocate(this%icu, np, 'PLICU', mempath) call mem_allocate(this%ilay, np, 'PLILAY', mempath) call mem_allocate(this%izone, np, 'PLIZONE', mempath) @@ -166,6 +169,7 @@ subroutine deallocate (this, mempath) call mem_deallocate(this%irpt, 'PLIRPT', mempath) call mem_deallocate(this%name, 'PLNAME', mempath) call mem_deallocate(this%icp, 'PLICP', mempath) + call mem_deallocate(this%iscp, 'PLISCP', mempath) call mem_deallocate(this%icu, 'PLICU', mempath) call mem_deallocate(this%ilay, 'PLILAY', mempath) call mem_deallocate(this%izone, 'PLIZONE', mempath) @@ -201,6 +205,7 @@ subroutine resize(this, np, mempath) call mem_reallocate(this%irpt, np, 'PLIRPT', mempath) call mem_reallocate(this%name, LENBOUNDNAME, np, 'PLNAME', mempath) call mem_reallocate(this%icp, np, 'PLICP', mempath) + call mem_reallocate(this%iscp, np, 'PLISCP', mempath) call mem_reallocate(this%icu, np, 'PLICU', mempath) call mem_reallocate(this%ilay, np, 'PLILAY', mempath) call mem_reallocate(this%izone, np, 'PLIZONE', mempath) @@ -245,6 +250,7 @@ subroutine load_particle(this, store, imdl, iprp, ip) this%istopzone = store%istopzone(ip) this%idrymeth = store%idrymeth(ip) this%icp = store%icp(ip) + this%iscp = store%iscp(ip) this%icu = store%icu(ip) this%ilay = store%ilay(ip) this%izone = store%izone(ip) @@ -282,6 +288,7 @@ subroutine save_particle(this, particle, ip) this%istopzone(ip) = particle%istopzone this%idrymeth(ip) = particle%idrymeth this%icp(ip) = particle%icp + this%iscp(ip) = particle%iscp this%icu(ip) = particle%icu this%ilay(ip) = particle%ilay this%izone(ip) = particle%izone