Method.f90

Go to the documentation of this file.

!> @brief Particle tracking strategies
module methodmodule
 
  use kindmodule, only: dp, i4b, lgp
  use constantsmodule, only: dzero
  use errorutilmodule, only: pstop
  use subcellmodule, only: subcelltype
  use particlemodule, only: particletype
  use particleeventsmodule, only: particleeventdispatchertype
  use particleeventmodule, only: particleeventtype
  use releaseeventmodule, only: releaseeventtype
  use timestepeventmodule, only: timestepeventtype
  use terminationeventmodule, only: terminationeventtype
  use weaksinkeventmodule, only: weaksinkeventtype
  use usertimeeventmodule, only: usertimeeventtype
  use featexiteventmodule, only: featexiteventtype
  use droppedeventmodule, only: droppedeventtype
  use basedismodule, only: disbasetype
  use prtfmimodule, only: prtfmitype
  use cellmodule, only: celltype
  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
 
  !> @brief Tracking method level enumeration.
  !!
  !> Tracking levels: 1: model, 2: grid feature, 3: grid subfeature.
  !! A tracking level identifies the domain through which a tracking
  !! method is responsible for moving a particle. Methods operate on
  !! a particular level and delegate to submethods for levels higher
  !! than (i.e. below the scope of) their own.
  !<
  enum, bind(C)
    enumerator :: level_model = 1
    enumerator :: level_feature = 2
    enumerator :: level_subfeature = 3
  end enum
 
  private
  public :: methodtype
 
  !> @brief Base type for particle tracking methods.
  !!
  !! The PRT tracking algorithm invokes a "tracking method" for each
  !! domain. A domain can be a model, cell in a model, or subcell in
  !! a cell. Tracking proceeds recursively, delegating to a possibly
  !! arbitrary number of subdomains (currently, only the three above
  !! are recognized). A tracking method is responsible for advancing
  !! a particle through a domain, delegating to subdomains as needed
  !! depending on cell geometry (implementing the strategy pattern).
  !<
  type, abstract :: methodtype
    character(len=40), pointer, public :: name !< method name
    logical(LGP), public :: delegates !< whether the method delegates
    type(prtfmitype), pointer, public :: fmi => null() !< ptr to fmi
    class(celltype), pointer, public :: cell => null() !< ptr to the current cell
    class(subcelltype), pointer, public :: subcell => null() !< ptr to the current subcell
    type(particleeventdispatchertype), pointer, public :: events => null() !< ptr to event dispatcher
    type(timeselecttype), pointer, public :: tracktimes => null() !< ptr to user-defined tracking times
    integer(I4B), dimension(:), pointer, contiguous, public :: izone => null() !< pointer to zone numbers
    real(dp), dimension(:), pointer, contiguous, public :: flowja => null() !< pointer to intercell flows
    real(dp), dimension(:), pointer, contiguous, public :: porosity => null() !< pointer to aquifer porosity
    real(dp), dimension(:), pointer, contiguous, public :: retfactor => null() !< pointer to retardation factor
  contains
    ! Implemented in all subtypes
    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 :: find_exits !< find domain exit solutions
    procedure :: pick_exit
    ! Implemented here
    procedure :: init
    procedure :: track
    procedure :: try_pass
    ! Event firing methods
    procedure :: release
    procedure :: terminate
    procedure :: timestep
    procedure :: weaksink
    procedure :: usertime
    procedure :: dropped
  end type methodtype
 
  abstract interface
    subroutine apply(this, particle, tmax)
      import dp
      import methodtype
      import particletype
      class(methodtype), intent(inout) :: this
      type(particletype), pointer, intent(inout) :: particle
      real(DP), intent(in) :: tmax
    end subroutine apply
    subroutine assess(this, particle, cell_defn, tmax)
      import dp
      import methodtype
      import particletype
      import celldefntype
      class(methodtype), intent(inout) :: this
      type(particletype), pointer, intent(inout) :: particle
      type(celldefntype), pointer, intent(inout) :: cell_defn
      real(DP), intent(in) :: tmax
    end subroutine assess
    subroutine deallocate (this)
      import methodtype
      class(methodtype), intent(inout) :: this
    end subroutine deallocate
  end interface
 
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
    type(prtfmitype), intent(in), pointer, optional :: fmi
    class(celltype), intent(in), pointer, optional :: cell
    class(subcelltype), intent(in), pointer, optional :: subcell
    type(particleeventdispatchertype), intent(in), pointer, optional :: events
    type(timeselecttype), intent(in), pointer, optional :: tracktimes
    integer(I4B), intent(in), pointer, optional :: izone(:)
    real(DP), intent(in), pointer, optional :: flowja(:)
    real(DP), intent(in), pointer, optional :: porosity(:)
    real(DP), intent(in), pointer, optional :: retfactor(:)
 
    if (present(fmi)) this%fmi => fmi
    if (present(cell)) this%cell => cell
    if (present(subcell)) this%subcell => subcell
    if (present(events)) this%events => events
    if (present(tracktimes)) this%tracktimes => tracktimes
    if (present(izone)) this%izone => izone
    if (present(flowja)) this%flowja => flowja
    if (present(porosity)) this%porosity => porosity
    if (present(retfactor)) this%retfactor => retfactor
  end subroutine init
 
  !> @brief Track the particle over subdomains of the given
  ! level until the particle terminates or tmax is reached.
  recursive subroutine track(this, particle, level, tmax)
    ! dummy
    class(methodtype), intent(inout) :: this
    type(particletype), pointer, intent(inout) :: particle
    integer(I4B) :: level
    real(dp), intent(in) :: tmax
    ! local
    logical(LGP) :: advancing
    integer(I4B) :: nextlevel
    class(methodtype), pointer :: submethod
 
    advancing = .true.
    nextlevel = level + 1
    do while (advancing)
      call this%load(particle, nextlevel, submethod)
      call submethod%apply(particle, tmax)
      call this%try_pass(particle, nextlevel, advancing)
    end do
  end subroutine track
 
  !> @brief Try passing the particle to the next subdomain.
  subroutine try_pass(this, particle, nextlevel, advancing)
    class(methodtype), intent(inout) :: this
    type(particletype), pointer, intent(inout) :: particle
    integer(I4B) :: nextlevel
    logical(LGP) :: advancing
 
    if (particle%advancing) then
      ! if still advancing, pass to the next subdomain.
      ! if that puts us on a boundary, then we're done.
      call this%pass(particle)
      if (particle%iboundary(nextlevel - 1) .ne. 0) &
        advancing = .false.
    else
      ! otherwise we're already done so
      ! reset the domain boundary value.
      advancing = .false.
      particle%iboundary = 0
    end if
  end subroutine try_pass
 
  !> @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
    integer, intent(in) :: next_level
    class(methodtype), pointer, intent(inout) :: submethod
    call pstop(1, "load must be overridden")
  end subroutine load
 
  !> @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 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
    class(particleeventtype), pointer :: event
    allocate (releaseeventtype :: event)
    call this%events%dispatch(particle, event)
    deallocate (event)
  end subroutine release
 
  !> @brief A particle terminates.
  subroutine terminate(this, particle, status)
    class(methodtype), intent(inout) :: this
    type(particletype), pointer, intent(inout) :: particle
    integer(I4B), intent(in), optional :: status
    class(particleeventtype), pointer :: event
    particle%advancing = .false.
    if (present(status)) particle%istatus = status
    allocate (terminationeventtype :: event)
    call this%events%dispatch(particle, event)
    deallocate (event)
  end subroutine terminate
 
  !> @brief A time step ends.
  subroutine timestep(this, particle)
    class(methodtype), intent(inout) :: this
    type(particletype), pointer, intent(inout) :: particle
    class(particleeventtype), pointer :: event
    allocate (timestepeventtype :: event)
    call this%events%dispatch(particle, event)
    deallocate (event)
  end subroutine timestep
 
  !> @brief A particle leaves a weak sink.
  subroutine weaksink(this, particle)
    class(methodtype), intent(inout) :: this
    type(particletype), pointer, intent(inout) :: particle
    class(particleeventtype), pointer :: event
    allocate (weaksinkeventtype :: event)
    call this%events%dispatch(particle, event)
    deallocate (event)
  end subroutine weaksink
 
  !> @brief A user-defined tracking time occurs.
  subroutine usertime(this, particle)
    class(methodtype), intent(inout) :: this
    type(particletype), pointer, intent(inout) :: particle
    class(particleeventtype), pointer :: event
    allocate (usertimeeventtype :: event)
    call this%events%dispatch(particle, event)
    deallocate (event)
  end subroutine usertime
 
  !> @brief A particle drops to the water table.
  subroutine dropped(this, particle)
    class(methodtype), intent(inout) :: this
    type(particletype), pointer, intent(inout) :: particle
    class(particleeventtype), pointer :: event
    allocate (droppedeventtype :: event)
    call this%events%dispatch(particle, event)
    deallocate (event)
  end subroutine dropped
 
end module methodmodule