refactor(prt): structure exit solution decision (#2489)

PRT instantiates a tree of tracking methods, each of which may delegate to submethods or solve particle motion itself. Leaves in the tree do the work, nodes delegate depending on various conditions.

The work consists of

- trying to obtain an exit solution from the method domain in each dimension independently
- choosing which of the candidate exit solutions, if any, is the correct (e.g., first) one
- solving particle coordinates for user-specified times up to exit time
- solving (the remaining) particle coordinates at exit time

and generally maintaining the particle's state while this happens.

This happens in the track_subcell routine in the relevant subcell methods. The four steps happen roughly in the order described, but the procedure is not structured otherwise (e.g. by subroutines, etc). So the concept of an exit solution is implicit and is so far mainly relevant for setting boundary indices (e.g. face numbers).

Up to now there has been a single criterion for selecting among candidate exit solutions: take the one that occurs first. However, #2478 is an example of a situation in which this gives unwanted results. So we may have arbitrary criteria for selecting the "best" exit solution. But this is awkward to implement right now, as the logic is interwoven across the track_subcell routine.

Separate the evaluation of candidate exit solutions and the selection of one into procedures. This makes filtering of candidate solutions cleaner.

- find_exits: generates candidates
- pick_exit: picks the best solution

Also introduce a base DomainType and module, used in the signature of find_exits. Also, miscellaneous other tidying, and clarify argument intent in some routines in TernarySolveUtils.f90.

Author: wpbonelli

make/makefile

@@ -310,11 +310,11 @@ $(OBJDIR)/Dis.o \
 $(OBJDIR)/Dis2d.o \
 $(OBJDIR)/BudgetObject.o \
 $(OBJDIR)/BoundaryPackage.o \
-$(OBJDIR)/CellDefn.o \
 $(OBJDIR)/ParticleEvent.o \
 $(OBJDIR)/sort.o \
+$(OBJDIR)/Domain.o \
+$(OBJDIR)/CellDefn.o \
 $(OBJDIR)/FlowModelInterface.o \
-$(OBJDIR)/Cell.o \
 $(OBJDIR)/WeakSinkEvent.o \
 $(OBJDIR)/UserTimeEvent.o \
 $(OBJDIR)/TimeStepEvent.o \
@@ -325,6 +325,8 @@ $(OBJDIR)/ReleaseEvent.o \
 $(OBJDIR)/prt-fmi.o \
 $(OBJDIR)/ParticleEvents.o \
 $(OBJDIR)/FeatExitEvent.o \
+$(OBJDIR)/ExitSolution.o \
+$(OBJDIR)/Cell.o \
 $(OBJDIR)/Method.o \
 $(OBJDIR)/CellExitEvent.o \
 $(OBJDIR)/SubcellExitEvent.o \

msvs/mf6core.vfproj

@@ -422,6 +422,7 @@
 		<File RelativePath="..\src\Solution\LinearMethods\ImsLinearSolver.f90"/>
 		<File RelativePath="..\src\Solution\LinearMethods\ImsReordering.f90"/></Filter>
 		<Filter Name="ParticleTracker">
+		<File RelativePath="..\src\Solution\ParticleTracker\Domain\Domain.f90"/>
 		<File RelativePath="..\src\Solution\ParticleTracker\Domain\Cell.f90"/>
 		<File RelativePath="..\src\Solution\ParticleTracker\Domain\CellDefn.f90"/>
 		<File RelativePath="..\src\Solution\ParticleTracker\Domain\CellPoly.f90"/>
@@ -431,6 +432,7 @@
 		<File RelativePath="..\src\Solution\ParticleTracker\Domain\Subcell.f90"/>
 		<File RelativePath="..\src\Solution\ParticleTracker\Domain\SubcellRect.f90"/>
 		<File RelativePath="..\src\Solution\ParticleTracker\Domain\SubcellTri.f90"/>
+		<File RelativePath="..\src\Solution\ParticleTracker\Method\ExitSolution.f90"/>
 		<File RelativePath="..\src\Solution\ParticleTracker\Method\Method.f90"/>
 		<File RelativePath="..\src\Solution\ParticleTracker\Method\MethodCell.f90"/>
 		<File RelativePath="..\src\Solution\ParticleTracker\Method\MethodCellPassToBot.f90"/>

src/Solution/ParticleTracker/Domain/Cell.f90

@@ -1,25 +1,16 @@
 module CellModule
 
   use CellDefnModule, only: CellDefnType
+  use DomainModule, only: DomainType
   use KindModule, only: I4B
   implicit none
   private
   public :: CellType
 
   !> @brief Base type for grid cells of a concrete type. Contains
   !! a cell-definition which is information shared by cell types.
-  type, abstract :: CellType
-    character(len=40), pointer :: type ! tracking domain type
+  type, abstract, extends(DomainType) :: CellType
     type(CellDefnType), pointer :: defn => null() ! cell defn
-  contains
-    procedure(destroy), deferred :: destroy !< destroy the cell
   end type CellType
 
-  abstract interface
-    subroutine destroy(this)
-      import CellType
-      class(CellType), intent(inout) :: this
-    end subroutine
-  end interface
-
 end module CellModule

src/Solution/ParticleTracker/Domain/Domain.f90

@@ -0,0 +1,21 @@
+module DomainModule
+
+  implicit none
+  private
+  public :: DomainType
+
+  !> @brief A tracking domain.
+  type, abstract :: DomainType
+    character(len=40), pointer :: type !< character string that names the tracking domain type
+  contains
+    procedure(destroy), deferred :: destroy !< destructor
+  end type DomainType
+
+  abstract interface
+    subroutine destroy(this)
+      import DomainType
+      class(DomainType), intent(inout) :: this
+    end subroutine
+  end interface
+
+end module DomainModule

src/Solution/ParticleTracker/Domain/Subcell.f90

@@ -1,27 +1,20 @@
 module SubcellModule
 
   use CellDefnModule, only: CellDefnType
+  use DomainModule, only: DomainType
   implicit none
   private
   public :: SubcellType
 
   !> @brief A subcell of a cell.
-  type, abstract :: SubcellType
-    private
-    character(len=40), pointer, public :: type !< character string that names the tracking domain type
+  type, abstract, extends(DomainType) :: SubcellType
     integer, public :: isubcell !< index of subcell in the cell
     integer, public :: icell !< index of cell in the source grid
   contains
-    procedure(destroy), deferred :: destroy !< destructor
     procedure(init), deferred :: init !< initializer
   end type SubcellType
 
   abstract interface
-    subroutine destroy(this)
-      import SubcellType
-      class(SubcellType), intent(inout) :: this
-    end subroutine
-
     subroutine init(this)
       import SubcellType
       class(SubcellType), intent(inout) :: this

src/Solution/ParticleTracker/Method/ExitSolution.f90

@@ -0,0 +1,30 @@
+module ExitSolutionModule
+
+  use KindModule, only: I4B, DP, LGP
+  use ConstantsModule, only: DZERO
+
+  implicit none
+
+  !> @brief Exit status codes
+  enum, bind(C)
+    enumerator :: OK_EXIT = 0 !< exit found using velocity interpolation
+    ! below only used for linear interp
+    enumerator :: OK_EXIT_CONSTANT = 1 !< exit found, constant velocity
+    enumerator :: NO_EXIT_STATIONARY = 2 !< no exit, zero velocity
+    enumerator :: NO_EXIT_NO_OUTFLOW = 3 !< no exit, no outflow
+  end enum
+
+  !> @brief Base type for exit solutions
+  type :: ExitSolutionType
+    integer(I4B) :: status = -1 !< domain exit status code
+    integer(I4B) :: iboundary = 0 !< boundary number
+    real(DP) :: dt = 1.0d+20 !< time to exit
+  end type ExitSolutionType
+
+  !> @brief Linear velocity interpolation exit solution
+  type, extends(ExitSolutionType) :: LinearExitSolutionType
+    real(DP) :: v = DZERO !< particle velocity
+    real(DP) :: dvdx = DZERO !< velocity gradient
+  end type LinearExitSolutionType
+
+end module ExitSolutionModule

src/Solution/ParticleTracker/Method/Method.f90

@@ -20,6 +20,9 @@ module MethodModule
   use CellDefnModule, only: CellDefnType
   use TimeSelectModule, only: TimeSelectType
   use MathUtilModule, only: is_close
+  use DomainModule, only: DomainType
+  use ExitSolutionModule, only: ExitSolutionType
+  use ListModule, only: ListType
   implicit none
 
   public :: LEVEL_MODEL, LEVEL_FEATURE, LEVEL_SUBFEATURE
@@ -68,10 +71,12 @@ module MethodModule
     procedure(apply), deferred :: apply !< apply the method to the particle
     procedure(assess), deferred :: assess !< assess conditions before tracking
     procedure(deallocate), deferred :: deallocate !< deallocate the method object
+    procedure :: get_level !< get the tracking method level
     ! Overridden in subtypes that delegate
     procedure :: pass !< pass the particle to the next subdomain
     procedure :: load !< load the subdomain tracking method
-    procedure :: get_level !< get the tracking method level
+    procedure :: find_exits !< find domain exit solutions
+    procedure :: pick_exit
     ! Implemented here
     procedure :: init
     procedure :: track
@@ -111,6 +116,7 @@ end subroutine deallocate
 
 contains
 
+  !> @brief Initialize the method with pointers to model data.
   subroutine init(this, fmi, cell, subcell, events, tracktimes, &
                   izone, flowja, porosity, retfactor)
     class(MethodType), intent(inout) :: this
@@ -178,7 +184,15 @@ subroutine try_pass(this, particle, nextlevel, advancing)
     end if
   end subroutine try_pass
 
-  !> @brief Load the subdomain tracking method (submethod).
+  !> @brief Get tracking method level.
+  function get_level(this) result(level)
+    class(MethodType), intent(in) :: this
+    integer(I4B) :: level
+    level = -1 ! suppress compiler warning
+    call pstop(1, "get_level must be overridden")
+  end function get_level
+
+  !> @brief Load subdomain tracking method (submethod).
   subroutine load(this, particle, next_level, submethod)
     class(MethodType), intent(inout) :: this
     type(ParticleType), pointer, intent(inout) :: particle
@@ -187,22 +201,35 @@ subroutine load(this, particle, next_level, submethod)
     call pstop(1, "load must be overridden")
   end subroutine load
 
-  !> @brief Get the tracking method's level.
-  function get_level(this) result(level)
-    class(MethodType), intent(in) :: this
-    integer(I4B) :: level
-    level = -1 ! suppress compiler warning
-    call pstop(1, "get_level must be overridden")
-  end function get_level
-
-  !> @brief Pass the particle to the next subdomain.
+  !> @brief Pass particle to the next subdomain or to a domain boundary.
   subroutine pass(this, particle)
     class(MethodType), intent(inout) :: this
     type(ParticleType), pointer, intent(inout) :: particle
     call pstop(1, "pass must be overridden")
   end subroutine pass
 
-  !> @brief Particle is released.
+  !> @brief Compute candidate exit solutions.
+  subroutine find_exits(this, particle, domain)
+    class(MethodType), intent(inout) :: this
+    type(ParticleType), pointer, intent(inout) :: particle
+    class(DomainType), intent(in) :: domain
+    if (.not. this%delegates) &
+      call pstop(1, "find_exits called on non-delegating method")
+    call pstop(1, "find_exits must be overridden in delegating methods")
+  end subroutine find_exits
+
+  !> @brief Choose an exit solution among candidates.
+  function pick_exit(this, particle) result(exit_soln)
+    class(MethodType), intent(inout) :: this
+    type(ParticleType), pointer, intent(inout) :: particle
+    integer(I4B) :: exit_soln
+    exit_soln = 0 ! suppress compiler warning
+    if (.not. this%delegates) &
+      call pstop(1, "pick_exit called on non-delegating method")
+    call pstop(1, "pick_exit must be overridden in delegating methods")
+  end function pick_exit
+
+  !> @brief A particle is released.
   subroutine release(this, particle)
     class(MethodType), intent(inout) :: this
     type(ParticleType), pointer, intent(inout) :: particle
@@ -212,7 +239,7 @@ subroutine release(this, particle)
     deallocate (event)
   end subroutine release
 
-  !> @brief Particle terminates.
+  !> @brief A particle terminates.
   subroutine terminate(this, particle, status)
     class(MethodType), intent(inout) :: this
     type(ParticleType), pointer, intent(inout) :: particle
@@ -225,7 +252,7 @@ subroutine terminate(this, particle, status)
     deallocate (event)
   end subroutine terminate
 
-  !> @brief Time step ends.
+  !> @brief A time step ends.
   subroutine timestep(this, particle)
     class(MethodType), intent(inout) :: this
     type(ParticleType), pointer, intent(inout) :: particle
@@ -235,7 +262,7 @@ subroutine timestep(this, particle)
     deallocate (event)
   end subroutine timestep
 
-  !> @brief Particle leaves a weak sink.
+  !> @brief A particle leaves a weak sink.
   subroutine weaksink(this, particle)
     class(MethodType), intent(inout) :: this
     type(ParticleType), pointer, intent(inout) :: particle
@@ -245,7 +272,7 @@ subroutine weaksink(this, particle)
     deallocate (event)
   end subroutine weaksink
 
-  !> @brief User-defined tracking time occurs.
+  !> @brief A user-defined tracking time occurs.
   subroutine usertime(this, particle)
     class(MethodType), intent(inout) :: this
     type(ParticleType), pointer, intent(inout) :: particle

src/Solution/ParticleTracker/Method/MethodSubcell.f90

@@ -10,6 +10,7 @@ module MethodSubcellModule
   private
   public :: MethodSubcellType
 
+  !> @brief Abstract base type for subcell tracking methods
   type, abstract, extends(MethodCellType) :: MethodSubcellType
   contains
     procedure, public :: assess
@@ -19,6 +20,7 @@ module MethodSubcellModule
 
 contains
 
+  !> @brief Assess conditions before tracking
   subroutine assess(this, particle, cell_defn, tmax)
     ! dummy
     class(MethodSubcellType), intent(inout) :: this