wip

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)

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

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.

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, &

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.

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