MethodCellPollock.f90

Go to the documentation of this file.

module methodcellpollockmodule
 
  use kindmodule, only: dp, i4b
  use constantsmodule, only: done, dzero
  use methodmodule, only: methodtype, level_feature, level_subfeature
  use methodcellmodule, only: methodcelltype
  use methodsubcellpollockmodule, only: methodsubcellpollocktype, &
                                        create_method_subcell_pollock
  use cellrectmodule, only: cellrecttype, create_cell_rect
  use subcellrectmodule, only: subcellrecttype, create_subcell_rect
  use particlemodule, only: particletype
  implicit none
 
  private
  public :: methodcellpollocktype
  public :: create_method_cell_pollock
 
  type, extends(methodcelltype) :: methodcellpollocktype
    type(methodsubcellpollocktype), pointer :: method_subcell_plck => null()
  contains
    procedure, public :: apply => apply_mcp
    procedure, public :: deallocate => destroy_mcp
    procedure, public :: load => load_mcp
    procedure, public :: load_subcell
    procedure, public :: pass => pass_mcp
  end type methodcellpollocktype
 
contains
 
  !> @brief Create a tracking method
  subroutine create_method_cell_pollock(method)
    ! dummy
    type(methodcellpollocktype), pointer :: method
    ! local
    type(cellrecttype), pointer :: cell
    type(subcellrecttype), pointer :: subcell
 
    allocate (method)
    call create_cell_rect(cell)
    method%cell => cell
    method%name => method%cell%type
    method%delegates = .true.
    call create_subcell_rect(subcell)
    method%subcell => subcell
 
    ! Create subcell method this method can delegate to
    call create_method_subcell_pollock(method%method_subcell_plck)
  end subroutine create_method_cell_pollock
 
  !> @brief Destroy the tracking method
  subroutine destroy_mcp(this)
    class(methodcellpollocktype), intent(inout) :: this
    deallocate (this%name)
 
    ! Deallocate owned subcell method
    call this%method_subcell_plck%deallocate()
    deallocate (this%method_subcell_plck)
  end subroutine destroy_mcp
 
  !> @brief Load subcell tracking method
  subroutine load_mcp(this, particle, next_level, submethod)
    ! modules
    use subcellmodule, only: subcelltype
    ! dummy
    class(methodcellpollocktype), intent(inout) :: this
    type(particletype), pointer, intent(inout) :: particle
    integer, intent(in) :: next_level
    class(methodtype), pointer, intent(inout) :: submethod
 
    select type (subcell => this%subcell)
    type is (subcellrecttype)
      call this%load_subcell(particle, subcell)
    end select
    call this%method_subcell_plck%init( &
      fmi=this%fmi, &
      cell=this%cell, &
      subcell=this%subcell, &
      events=this%events, &
      tracktimes=this%tracktimes)
    submethod => this%method_subcell_plck
    particle%itrdomain(next_level) = 1
  end subroutine load_mcp
 
  !> @brief Having exited the lone subcell, pass the particle to the cell face
  !! In this case the lone subcell is the cell.
  subroutine pass_mcp(this, particle)
    ! dummy
    class(methodcellpollocktype), intent(inout) :: this
    type(particletype), pointer, intent(inout) :: particle
    ! local
    integer(I4B) :: exitface
    integer(I4B) :: entryface
 
    exitface = particle%iboundary(level_subfeature)
    ! Map subcell exit face to cell face
    select case (exitface) ! note: exitFace uses Dave's iface convention
    case (0)
      entryface = -1
    case (1)
      entryface = 1
    case (2)
      entryface = 3
    case (3)
      entryface = 4
    case (4)
      entryface = 2
    case (5)
      entryface = 6 ! note: inface=5 same as inface=1 due to wraparound
    case (6)
      entryface = 7
    end select
    if (entryface .eq. -1) then
      particle%iboundary(level_feature) = 0
    else
      if ((entryface .ge. 1) .and. (entryface .le. 4)) then
        ! Account for local cell rotation
        select type (cell => this%cell)
        type is (cellrecttype)
          entryface = entryface + cell%ipvOrigin - 1
        end select
        if (entryface .gt. 4) entryface = entryface - 4
      end if
      particle%iboundary(level_feature) = entryface
    end if
  end subroutine pass_mcp
 
  !> @brief Apply Pollock's method to a rectangular cell
  subroutine apply_mcp(this, particle, tmax)
    ! dummy
    class(methodcellpollocktype), intent(inout) :: this
    type(particletype), pointer, intent(inout) :: particle
    real(DP), intent(in) :: tmax
    ! local
    real(DP) :: xOrigin
    real(DP) :: yOrigin
    real(DP) :: zOrigin
    real(DP) :: sinrot
    real(DP) :: cosrot
 
    select type (cell => this%cell)
    type is (cellrecttype)
      call this%assess(particle, this%cell%defn, tmax)
      if (.not. particle%advancing) return
 
      ! Transform model coordinates to local cell coordinates
      ! (translated/rotated but not scaled relative to model),
      ! track the particle over the cell, then transform back.
      xorigin = cell%xOrigin
      yorigin = cell%yOrigin
      zorigin = cell%zOrigin
      sinrot = cell%sinrot
      cosrot = cell%cosrot
      call particle%transform(xorigin, yorigin, zorigin, &
                              sinrot, cosrot)
      call this%track(particle, 2, tmax)
      call particle%transform(xorigin, yorigin, zorigin, &
                              sinrot, cosrot, invert=.true.)
      call particle%reset_transform()
    end select
  end subroutine apply_mcp
 
  !> @brief Loads the lone rectangular subcell from the rectangular cell
  subroutine load_subcell(this, particle, subcell) !
    ! dummy
    class(methodcellpollocktype), intent(inout) :: this
    type(particletype), pointer, intent(inout) :: particle
    type(subcellrecttype), intent(inout) :: subcell
 
    select type (cell => this%cell)
    type is (cellrecttype)
      ! Set cell/subcell numbers
      subcell%icell = cell%defn%icell
      subcell%isubcell = 1
 
      ! Subcell calculations will be done in local subcell coordinates
      subcell%dx = cell%dx
      subcell%dy = cell%dy
      subcell%dz = cell%dz
      subcell%sinrot = dzero
      subcell%cosrot = done
      subcell%xOrigin = dzero
      subcell%yOrigin = dzero
      subcell%zOrigin = dzero
 
      ! Set subcell edge velocities
      subcell%vx1 = cell%vx1 ! cell velocities already account for retfactor and porosity
      subcell%vx2 = cell%vx2
      subcell%vy1 = cell%vy1
      subcell%vy2 = cell%vy2
      subcell%vz1 = cell%vz1
      subcell%vz2 = cell%vz2
    end select
  end subroutine load_subcell
 
end module methodcellpollockmodule