MethodCell.f90

Go to the documentation of this file.

module methodcellmodule
 
  use kindmodule, only: dp, i4b, lgp
  use errorutilmodule, only: pstop
  use constantsmodule, only: done, dzero
  use methodmodule, only: methodtype, level_feature
  use particlemodule, only: particletype, active, term_no_exits, term_boundary
  use particleeventmodule, only: particleeventtype
  use cellexiteventmodule, only: cellexiteventtype
  use celldefnmodule, only: celldefntype, saturation_dry
  use iteratormodule, only: iteratortype
  implicit none
 
  private
  public :: methodcelltype
 
  type, abstract, extends(methodtype) :: methodcelltype
  contains
    procedure, public :: assess
    procedure, public :: cellexit
    procedure, public :: forms_cycle
    procedure, public :: store_event
    procedure, public :: get_level
    procedure, public :: try_pass
    procedure :: iboundary_to_icellface
  end type methodcelltype
 
contains
  !> @brief Try passing the particle to the next subdomain.
  !! or to a boundary of this method's tracking domain.
  !!
  !! If the particle is not advancing, or if it's not at
  !! a boundary, nothing happens. Otherwise, raise a cell
  !! exit event. Then, if at a net outflow boundary face,
  !! or top face of a partially saturated cell, terminate.
  !<
  subroutine try_pass(this, particle, nextlevel, advancing)
    class(methodcelltype), intent(inout) :: this
    type(particletype), pointer, intent(inout) :: particle
    logical(LGP) :: advancing
    integer(I4B) :: nextlevel
    ! local
    integer(I4B) :: ic, iboundary, icellface
 
    if (.not. particle%advancing) then
      advancing = .false.
      particle%iboundary = 0
      return
    end if
 
    call this%pass(particle)
 
    iboundary = particle%iboundary(level_feature)
    icellface = this%iboundary_to_icellface(iboundary)
    if (icellface <= 0) return
 
    ! on a cell face, done advancing. raise an exit event
    advancing = .false.
    call this%cellexit(particle)
 
    ! assigned boundary face with net outflow? terminate
    ic = particle%itrdomain(level_feature)
    if (this%fmi%is_net_out_boundary_face(ic, icellface)) then
      call this%terminate(particle, status=term_boundary)
      return
    end if
 
  end subroutine try_pass
 
  !> @brief Check reporting/terminating conditions before tracking
  !! the particle across the cell.
  !!
  !! Check a number of conditions determining whether to continue
  !! tracking the particle or terminate it, as well as whether to
  !! record any output data as per selected reporting conditions.
  !<
  subroutine assess(this, particle, cell_defn, tmax)
    ! modules
    use tdismodule, only: endofsimulation, totimc, totim
    use particlemodule, only: term_weaksink, term_no_exits, &
                              term_stopzone, term_inactive
    use particleeventmodule, only: terminate, timestep, weaksink, usertime
    ! dummy
    class(methodcelltype), intent(inout) :: this
    type(particletype), pointer, intent(inout) :: particle
    type(celldefntype), pointer, intent(inout) :: cell_defn
    real(DP), intent(in) :: tmax
    ! local
    logical(LGP) :: dry_cell, dry_particle, no_exit_face, stop_zone, weak_sink
    integer(I4B) :: i
    real(DP) :: t, ttrackmax
 
    dry_cell = cell_defn%isatstat == saturation_dry
    dry_particle = particle%z > cell_defn%top
    no_exit_face = cell_defn%inoexitface > 0
    stop_zone = cell_defn%izone > 0 .and. particle%istopzone == cell_defn%izone
    weak_sink = cell_defn%iweaksink > 0
 
    particle%izone = cell_defn%izone
    if (stop_zone) then
      call this%terminate(particle, status=term_stopzone)
      return
    end if
 
    if (no_exit_face .and. .not. dry_cell) then
      call this%terminate(particle, status=term_no_exits)
      return
    end if
 
    if (weak_sink) then
      if (particle%istopweaksink > 0) then
        call this%terminate(particle, status=term_weaksink)
        return
      else
        call this%weaksink(particle)
      end if
    end if
 
    if (dry_cell) then
      if (particle%idrymeth == 0) then
        ! drop to cell bottom. handled by pass
        ! to bottom method, nothing to do here
        no_exit_face = .false.
      else if (particle%idrymeth == 1) then
        ! stop
        call this%terminate(particle, status=term_inactive)
        return
      else if (particle%idrymeth == 2) then
        ! stay
        particle%advancing = .false.
        no_exit_face = .false.
 
        ! we might report tracking times
        ! out of order here, but we want
        ! the particle termination event
        ! (if this is the last time step)
        ! to have the maximum tracking t,
        ! so we need to keep tabs on it.
        ttrackmax = totim
 
        ! update tracking time to time
        ! step end time and save record
        particle%ttrack = totim
        call this%timestep(particle)
 
        ! record user tracking times
        call this%tracktimes%advance()
        if (this%tracktimes%any()) then
          do i = this%tracktimes%selection(1), this%tracktimes%selection(2)
            t = this%tracktimes%times(i)
            if (t < totimc) cycle
            if (t >= tmax) exit
            particle%ttrack = t
            call this%usertime(particle)
            if (t > ttrackmax) ttrackmax = t
          end do
        end if
 
        ! terminate if last period/step
        if (endofsimulation) then
          particle%ttrack = ttrackmax
          call this%terminate(particle, status=term_no_exits)
          return
        end if
      end if
    else if (dry_particle .and. this%name /= "passtobottom") then
      if (particle%idrymeth == 0) then
        ! drop to water table
        particle%z = cell_defn%top
        call this%dropped(particle)
      else if (particle%idrymeth == 1) then
        ! stop
        call this%terminate(particle, status=term_inactive)
        return
      else if (particle%idrymeth == 2) then
        ! stay
        particle%advancing = .false.
        no_exit_face = .false.
 
        ! we might report tracking times
        ! out of order here, but we want
        ! the particle termination event
        ! (if this is the last time step)
        ! to have the maximum tracking t,
        ! so we need to keep tabs on it.
        ttrackmax = totim
 
        ! update tracking time to time
        ! step end time and save record
        particle%ttrack = totim
        call this%timestep(particle)
 
        ! record user tracking times
        call this%tracktimes%advance()
        if (this%tracktimes%any()) then
          do i = this%tracktimes%selection(1), this%tracktimes%selection(2)
            t = this%tracktimes%times(i)
            if (t < totimc) cycle
            if (t >= tmax) exit
            particle%ttrack = t
            call this%usertime(particle)
            if (t > ttrackmax) ttrackmax = t
          end do
        end if
      end if
    end if
 
    if (no_exit_face) then
      particle%advancing = .false.
      particle%istatus = term_no_exits
      call this%terminate(particle)
      return
    end if
 
  end subroutine assess
 
  !> @brief Particle exits a cell.
  subroutine cellexit(this, particle)
    class(methodcelltype), intent(inout) :: this
    type(particletype), pointer, intent(inout) :: particle
    class(particleeventtype), pointer :: event
    ! local
    integer(I4B) :: i, nhist
    class(*), pointer :: prev
 
    allocate (cellexiteventtype :: event)
    select type (event)
    type is (cellexiteventtype)
      event%exit_face = particle%iboundary(level_feature)
    end select
    call this%events%dispatch(particle, event)
    if (particle%icycwin == 0) then
      deallocate (event)
      return
    end if
    if (this%forms_cycle(particle, event)) then
      ! print event history
      print *, "Cyclic pathline detected"
      nhist = particle%history%Count()
      do i = 1, nhist
        prev => particle%history%GetItem(i)
        select type (prev)
        class is (particleeventtype)
          print *, "Back ", nhist - i + 1, ": ", prev%get_text()
        end select
      end do
      print *, "Current :", event%get_text()
      call pstop(1, 'Cyclic pathline detected, aborting')
    else
      call this%store_event(particle, event)
    end if
  end subroutine cellexit
 
  !> @brief Check if the event forms a cycle in the particle path.
  function forms_cycle(this, particle, event) result(found_cycle)
    ! dummy
    class(methodcelltype), intent(inout) :: this
    type(particletype), pointer, intent(inout) :: particle
    class(particleeventtype), pointer, intent(in) :: event
    ! local
    class(iteratortype), allocatable :: itr
    logical(LGP) :: found_cycle
 
    found_cycle = .false.
    select type (event)
    type is (cellexiteventtype)
      itr = particle%history%Iterator()
      do while (itr%has_next())
        call itr%next()
        select type (prev => itr%value())
        class is (cellexiteventtype)
          if (event%icu == prev%icu .and. &
              event%ilay == prev%ilay .and. &
              event%izone == prev%izone .and. &
              event%exit_face == prev%exit_face .and. &
              event%exit_face /= 0) then
            found_cycle = .true.
            exit
          end if
        end select
      end do
    end select
  end function forms_cycle
 
  !> @brief Save the event in the particle's history.
  !! Acts like a queue, the oldest event is removed
  !! when the event count exceeds the maximum size.
  subroutine store_event(this, particle, event)
    ! dummy
    class(methodcelltype), intent(inout) :: this
    type(particletype), pointer, intent(inout) :: particle
    class(particleeventtype), pointer, intent(in) :: event
    ! local
    class(*), pointer :: p
 
    select type (event)
    type is (cellexiteventtype)
      p => event
      call particle%history%Add(p)
      if (particle%history%Count() > particle%icycwin) &
        call particle%history%RemoveNode(1, .true.)
    end select
  end subroutine store_event
 
  !> @brief Get the cell method's level.
  function get_level(this) result(level)
    class(methodcelltype), intent(in) :: this
    integer(I4B) :: level
    level = level_feature
  end function get_level
 
  !> @brief Convert an iboundary number to an iface number
  function iboundary_to_icellface(this, iboundary) result(iface)
    class(methodcelltype), intent(inout) :: this
    integer(I4B), intent(in) :: iboundary
    integer(I4B) :: iface
    integer(I4B) :: nfaces
 
    iface = iboundary
    nfaces = this%cell%defn%npolyverts + 2
    if (iface >= nfaces) &
      ! uncompress and drop wraparound index
      iface = iface + (this%fmi%max_faces - nfaces) - 1
  end function iboundary_to_icellface
 
end module methodcellmodule