wip

Author: wpbonelli

src/Model/ParticleTracking/prt.f90

@@ -962,6 +962,7 @@ subroutine prt_solve(this)
           ! Reset previous cell, subcell, and zone numbers
           particle%icp = 0
           particle%iscp = 0
+          particle%iscefp = 0
           particle%izp = 0
 
           ! Update particle storage

src/Solution/ParticleTracker/Method.f90

@@ -267,6 +267,7 @@ subroutine check(this, particle, cell_defn)
       ! dry particle
       if (particle%idrymeth == 0) then
         ! drop to water table
+        print *, "dropping"
         particle%z = cell_defn%top
         call this%save(particle, reason=1)
       else if (particle%idrymeth == 1) then

src/Solution/ParticleTracker/MethodCellTernary.f90

@@ -1,6 +1,6 @@
 module MethodCellTernaryModule
 
-  use KindModule, only: DP, I4B
+  use KindModule, only: DP, I4B, LGP
   use ErrorUtilModule, only: pstop
   use MethodModule
   use MethodSubcellPoolModule
@@ -85,18 +85,6 @@ 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 /= 0 .and. subcell%isubcell == particle%iscp) then
-      !   print *, '---', subcell%isubcell
-      !   particle%advancing = .false.
-      !   particle%idomain(3) = particle%iscp
-      !   particle%istatus = 2
-      !   call this%save(particle, reason=3)
-      ! else
-      !   particle%iscp = particle%idomain(3)
-      ! end if
     end select
 
     call method_subcell_tern%init( &
@@ -150,6 +138,19 @@ subroutine pass_mct(this, particle)
       ! Subcell top (cell top)
       inface = this%nverts + 3
     end select
+
+    print *, "--------", isc, particle%iscp, exitFace, particle%iscefp
+    ! if (inface /= exitFace .and. isc == particle%iscp .and. exitFace == particle%iscefp) then
+    !   print *, 'subcell cycle, terminating'
+    !   particle%advancing = .false.
+    !   particle%idomain(3) = particle%iscp
+    !   particle%istatus = 2
+    !   call this%save(particle, reason=3)
+    ! else
+    !   particle%iscp = particle%idomain(3)
+    !   particle%iscefp = exitFace
+    ! end if
+
     if (inface .eq. -1) then
       particle%iboundary(2) = 0
     else if (inface .eq. 0) then

src/Solution/ParticleTracker/MethodSubcellTernary.f90

@@ -215,11 +215,13 @@ subroutine track_subcell(this, subcell, particle, tmax)
     ! considering only the velocity component in the direction of
     ! the face. Then compute the particle's exit time.
     if (itrifaceexit /= 0) then
+      print *, "exit thru side"
       ! Exit through lateral subcell face
       exitFace = itrifaceexit
       dtexit = dtexitxy
     else if (dtexitz < dtexitxy) then
       ! Exit through top or bottom
+      print *, "exit thru top or bottom"
       if (itopbotexit == -1) then
         exitFace = 4
       else
@@ -230,6 +232,16 @@ subroutine track_subcell(this, subcell, particle, tmax)
     texit = particle%ttrack + dtexit
     t0 = particle%ttrack
 
+    ! ! don't allow a particle to return to the subcell was just in
+    ! if (particle%idomain(2) == particle%icp .and. particle%idomain(3) == particle%iscp) then
+    !   print *, "isc: ", particle%idomain(3)
+    !   particle%idomain(3) = particle%iscp
+    !   particle%istatus = 5
+    !   particle%advancing = .false.
+    !   call this%save(particle, reason=3)
+    !   return
+    ! end if
+
     ! Select user tracking times to solve. If this is the first time step
     ! of the simulation, include all times before it begins; if it is the
     ! last time step include all times after it ends only if the 'extend'
@@ -253,6 +265,12 @@ subroutine track_subcell(this, subcell, particle, tmax)
       end do
     end if
 
+    ! if (izstatus > 1) then
+    !   particle%istatus = izstatus
+    !   particle%advancing = .false.
+    !   reason = 3
+    ! end if
+
     ! Compute exit time and face and update the particle's coordinates
     ! (local, unscaled) and other properties. The particle may at this
     ! point lie on a boundary of the subcell or may still be within it.
@@ -265,7 +283,6 @@ 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.
@@ -281,6 +298,13 @@ subroutine track_subcell(this, subcell, particle, tmax)
     particle%z = z
     particle%ttrack = t
     particle%iboundary(3) = exitFace
+
+    ! save previous subcell id and exit face
+    particle%iscp = particle%idomain(3)
+    particle%iscefp = particle%iboundary(3)
+    print *, "exit face: ", particle%iboundary(3)
+    print *, "subcell: ", particle%idomain(3)
+
     call this%save(particle, reason=reason)
   end subroutine track_subcell
 

src/Solution/ParticleTracker/Particle.f90

@@ -46,6 +46,7 @@ module ParticleModule
     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 :: iscefp !< previous subcell exit face
     integer(I4B), public :: icu !< user cell number
     integer(I4B), public :: ilay !< grid layer
     integer(I4B), public :: izone !< current zone number
@@ -90,8 +91,6 @@ module ParticleModule
     ! state
     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
@@ -135,8 +134,6 @@ subroutine allocate_particle_store(this, np, mempath)
     call mem_allocate(this%irpt, np, 'PLIRPT', 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)
@@ -168,8 +165,6 @@ subroutine deallocate (this, mempath)
     call mem_deallocate(this%iprp, 'PLIPRP', 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)
@@ -204,8 +199,6 @@ subroutine resize(this, np, mempath)
     call mem_reallocate(this%iprp, np, 'PLIPRP', 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)
@@ -249,8 +242,9 @@ subroutine load_particle(this, store, imdl, iprp, ip)
     this%istopweaksink = store%istopweaksink(ip)
     this%istopzone = store%istopzone(ip)
     this%idrymeth = store%idrymeth(ip)
-    this%icp = store%icp(ip)
-    this%iscp = store%iscp(ip)
+    this%icp = 0
+    this%iscp = 0
+    this%iscefp = 0
     this%icu = store%icu(ip)
     this%ilay = store%ilay(ip)
     this%izone = store%izone(ip)
@@ -287,8 +281,6 @@ subroutine save_particle(this, particle, ip)
     this%istopweaksink(ip) = particle%istopweaksink
     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