methodcellpollockquadmodule Module Reference

Data Types
type	methodcellpollockquadtype

Functions/Subroutines
procedure subroutine, public	create_method_cell_quad (method)
	Create a new tracking method. More...
subroutine	destroy (this)
	Destroy the tracking method. More...
subroutine	load_mcpq (this, particle, next_level, submethod)
	Load subcell into tracking method. More...
subroutine	pass_mcpq (this, particle)
	Pass particle to next subcell if there is one, or to the cell face. More...
subroutine	apply_mcpq (this, particle, tmax)
	Solve the quad-rectangular cell via Pollock's method. More...
subroutine	load_subcell (this, particle, subcell)
	Load the rectangular subcell from the rectangular cell. More...

Function/Subroutine Documentation

apply_mcpq()

subroutine methodcellpollockquadmodule::apply_mcpq ( class(methodcellpollockquadtype), intent(inout)  this, type(particletype), intent(inout), pointer  particle, real(dp), intent(in)  tmax  )

private

Definition at line 199 of file MethodCellPollockQuad.f90.

    ! -- dummy
    class(MethodCellPollockQuadType), intent(inout) :: this
    type(ParticleType), pointer, intent(inout) :: particle
    real(DP), intent(in) :: tmax
    ! -- local
    double precision :: xOrigin, yOrigin, zOrigin, sinrot, cosrot
 
    select type (cell => this%cell)
    type is (cellrectquadtype)
      ! -- Update particle state, terminate early if done advancing
      call this%update(particle, cell%defn)
      if (.not. particle%advancing) return
 
      ! -- If the particle is above the top of the cell (which is presumed to
      ! -- represent a water table above the cell bottom), pass the particle
      ! -- vertically and instantaneously to the cell top elevation and save
      ! -- the particle state to output file(s).
      if (particle%z > cell%defn%top) then
        particle%z = cell%defn%top
        call this%save(particle, reason=1) ! reason=1: cell transition
      end if
 
      ! -- Transform particle location into local cell coordinates,
      !    translated and rotated but not scaled relative to model.
      !    Then track particle, transform back to model coordinates,
      !    and reset transformation (drop accumulated roundoff error)
      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) ! kluge, hardwired to level 2
      call particle%transform(xorigin, yorigin, zorigin, &
                              sinrot, cosrot, invert=.true.)
      call particle%transform(reset=.true.)
    end select

create_method_cell_quad()

procedure subroutine, public methodcellpollockquadmodule::create_method_cell_quad

(

type(methodcellpollockquadtype), pointer

method

)

Definition at line 31 of file MethodCellPollockQuad.f90.

    ! -- dummy
    type(MethodCellPollockQuadType), pointer :: method
    ! -- local
    type(CellRectQuadType), pointer :: cell
    type(SubcellRectType), pointer :: subcell
 
    allocate (method)
    call create_cell_rect_quad(cell)
    method%cell => cell
    method%type => method%cell%type
    method%delegates = .true.
    call create_subcell_rect(subcell)
    method%subcell => subcell

Here is the call graph for this function:

Here is the caller graph for this function:

destroy()

subroutine methodcellpollockquadmodule::destroy ( class(methodcellpollockquadtype), intent(inout)  this )

private

Definition at line 48 of file MethodCellPollockQuad.f90.

    class(MethodCellPollockQuadType), intent(inout) :: this
    deallocate (this%type)

load_mcpq()

subroutine methodcellpollockquadmodule::load_mcpq ( class(methodcellpollockquadtype), intent(inout)  this, type(particletype), intent(inout), pointer  particle, integer, intent(in)  next_level, class(methodtype), intent(inout), pointer  submethod  )

private

Definition at line 54 of file MethodCellPollockQuad.f90.

    class(MethodCellPollockQuadType), 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 method_subcell_plck%init( &
      cell=this%cell, &
      subcell=this%subcell, &
      trackfilectl=this%trackfilectl, &
      tracktimes=this%tracktimes)
    submethod => method_subcell_plck

load_subcell()

subroutine methodcellpollockquadmodule::load_subcell ( class(methodcellpollockquadtype), intent(inout)  this, type(particletype), intent(inout), pointer  particle, class(subcellrecttype), intent(inout)  subcell  )

private

Definition at line 241 of file MethodCellPollockQuad.f90.

    ! -- dummy
    class(MethodCellPollockQuadType), intent(inout) :: this
    type(ParticleType), pointer, intent(inout) :: particle
    class(SubcellRectType), intent(inout) :: subcell
    ! -- local
    double precision :: dx, dy, dz, areax, areay, areaz
    double precision :: dxprel, dyprel
    integer :: isc, npolyverts, m1, m2
    double precision :: qextl1, qextl2, qintl1, qintl2
    double precision :: factor, term
 
    select type (cell => this%cell)
    type is (cellrectquadtype)
      factor = done / cell%defn%retfactor
      factor = factor / cell%defn%porosity
      npolyverts = cell%defn%npolyverts
 
      isc = particle%idomain(3)
      ! -- Subcells 1, 2, 3, and 4 are Pollock's subcells A, B, C, and D,
      ! -- respectively
 
      dx = cell%dx
      dy = cell%dy
      ! -- If not already known, determine subcell number
      if (isc .le. 0) then
        dxprel = particle%x / dx
        dyprel = particle%y / dy
        if (dxprel .lt. 5d-1) then
          if (dyprel .lt. 5d-1) then
            isc = 3
          else if (dyprel .gt. 5d-1) then
            isc = 4
          else
            ! kluge note: need to resolve this ambiguity based on flow direction
            call pstop(1, "particle initially on shared subcell edge")
          end if
        else if (dxprel .gt. 5d-1) then
          if (dyprel .lt. 5d-1) then
            isc = 2
          else if (dyprel .gt. 5d-1) then
            isc = 1
          else
            ! kluge note: need to resolve this ambiguity based on flow direction
            call pstop(1, "particle initially on shared subcell edge")
          end if
        else
          ! kluge note: need to resolve this ambiguity based on flow direction
          call pstop(1, "particle initially on shared subcell edge")
        end if
        subcell%isubcell = isc
        ! kluge note: as a matter of form, do we want to allow
        ! this subroutine to modify the particle???
        particle%idomain(3) = isc
        ! kluge note: initial insubface is not currently being determined
      end if
      dx = 5d-1 * dx
      dy = 5d-1 * dy
      dz = cell%defn%top - &
           cell%defn%bot ! kluge note: need to account for partial saturation
      areax = dy * dz
      areay = dx * dz
      areaz = dx * dy
      qintl1 = cell%qintl(isc)
      ! qintl list wraps around, so isc+1=5 is ok
      qintl2 = cell%qintl(isc + 1)
      qextl1 = cell%qextl1(isc)
      qextl2 = cell%qextl2(isc)
      !
      subcell%dx = dx
      subcell%dy = dy
      subcell%dz = dz
      subcell%sinrot = 0d0
      subcell%cosrot = 1d0
      subcell%zOrigin = 0d0
      select case (isc)
      case (1)
        subcell%xOrigin = dx
        subcell%yOrigin = dy
        term = factor / areax
        subcell%vx1 = qintl1 * term
        subcell%vx2 = -qextl2 * term
        term = factor / areay
        subcell%vy1 = -qintl2 * term
        subcell%vy2 = -qextl1 * term
      case (2)
        subcell%xOrigin = dx
        subcell%yOrigin = 0d0
        term = factor / areax
        subcell%vx1 = -qintl2 * term
        subcell%vx2 = -qextl1 * term
        term = factor / areay
        subcell%vy1 = qextl2 * term
        subcell%vy2 = -qintl1 * term
      case (3)
        subcell%xOrigin = 0d0
        subcell%yOrigin = 0d0
        term = factor / areax
        subcell%vx1 = qextl2 * term
        subcell%vx2 = -qintl1 * term
        term = factor / areay
        subcell%vy1 = qextl1 * term
        subcell%vy2 = qintl2 * term
      case (4)
        subcell%xOrigin = 0d0
        subcell%yOrigin = dy
        term = factor / areax
        subcell%vx1 = qextl1 * term
        subcell%vx2 = qintl2 * term
        term = factor / areay
        subcell%vy1 = qintl1 * term
        subcell%vy2 = -qextl2 * term
      end select
      m1 = npolyverts + 2
      m2 = m1 + 1
      term = factor / areaz
      subcell%vz1 = 2.5d-1 * cell%defn%faceflow(m1) * term
      subcell%vz2 = -2.5d-1 * cell%defn%faceflow(m2) * term
    end select

Here is the call graph for this function:

pass_mcpq()

subroutine methodcellpollockquadmodule::pass_mcpq ( class(methodcellpollockquadtype), intent(inout)  this, type(particletype), intent(inout), pointer  particle  )

private

Definition at line 73 of file MethodCellPollockQuad.f90.

    ! -- dummy
    class(MethodCellPollockQuadType), intent(inout) :: this
    type(ParticleType), pointer, intent(inout) :: particle
    ! -- local
    integer :: isc, exitFace, npolyverts, inface, infaceoff
 
    select type (cell => this%cell)
    type is (cellrectquadtype)
      exitface = particle%iboundary(3)
      isc = particle%idomain(3)
      npolyverts = cell%defn%npolyverts
 
      select case (exitface) ! kluge note: exitFace uses Dave's iface convention
      case (0)
        ! -- Subcell interior (cell interior)
        inface = -1
      case (1)
        select case (isc)
        case (1)
          ! -- W face, subcell 1 --> E face, subcell 4  (cell interior)
          particle%idomain(3) = 4
          particle%iboundary(3) = 2
          inface = 0 ! kluge note: want Domain(2) unchanged; Boundary(2) = 0
        case (2)
          ! -- W face, subcell 2 --> E face, subcell 3 (cell interior)
          particle%idomain(3) = 3
          particle%iboundary(3) = 2
          inface = 0 ! kluge note: want Domain(2) unchanged; Boundary(2) = 0
        case (3)
          ! -- W face, subcell 3 (cell face)
          inface = 1 ! kluge note: want Domain(2) = -Domain(2); Boundary(2) = inface
          infaceoff = 0
        case (4)
          ! -- W face, subcell 4 (cell face)
          inface = 2 ! kluge note: want Domain(2) = -Domain(2); Boundary(2) = inface
          infaceoff = -1
        end select
      case (2)
        select case (isc)
        case (1)
          ! -- E face, subcell 1 (cell face)
          inface = 3 ! kluge note: want Domain(2) = -Domain(2); Boundary(2) = inface
          infaceoff = 0
        case (2)
          ! -- E face, subcell 2 (cell face)
          inface = 4 ! kluge note: want Domain(2) = -Domain(2); Boundary(2) = inface
          infaceoff = -1
        case (3)
          ! -- E face, subcell 3 --> W face, subcell 2 (cell interior)
          particle%idomain(3) = 2
          particle%iboundary(3) = 1
          inface = 0 ! kluge note: want Domain(2) unchanged; Boundary(2) = 0
        case (4)
          ! -- E face, subcell 4 --> W face subcell 1 (cell interior)
          particle%idomain(3) = 1
          particle%iboundary(3) = 1
          inface = 0 ! kluge note: want Domain(2) unchanged; Boundary(2) = 0
        end select
      case (3)
        select case (isc)
        case (1)
          ! -- S face, subcell 1 --> N face, subcell 2 (cell interior)
          particle%idomain(3) = 2
          particle%iboundary(3) = 4
          inface = 0 ! kluge note: want Domain(2) unchanged; Boundary(2) = 0
        case (2)
          ! -- S face, subcell 2 (cell face)
          inface = 4 ! kluge note: want Domain(2) = -Domain(2); Boundary(2) = inface
          infaceoff = 0
        case (3)
          ! -- S face, subcell 3 (cell face)
          inface = 1 ! kluge note: want Domain(2) = -Domain(2); Boundary(2) = inface
          infaceoff = -1
        case (4)
          ! -- S face, subcell 4 --> N face, subcell 3 (cell interior)
          particle%idomain(3) = 3
          particle%iboundary(3) = 4
          inface = 0 ! kluge note: want Domain(2) unchanged; Boundary(2) = 0
        end select
      case (4)
        select case (isc)
        case (1)
          ! -- N face, subcell 1 (cell face)
          inface = 3 ! kluge note: want Domain(2) = -Domain(2); Boundary(2) = inface
          infaceoff = -1
        case (2)
          ! -- N face, subcell 2 --> S face, subcell 1 (cell interior)
          particle%idomain(3) = 1
          particle%iboundary(3) = 3
          inface = 0 ! kluge note: want Domain(2) unchanged; Boundary(2) = 0
        case (3)
          ! -- N face, subcell 3 --> S face, subcell 4 (cell interior)
          particle%idomain(3) = 4
          particle%iboundary(3) = 3
          inface = 0 ! kluge note: want Domain(2) unchanged; Boundary(2) = 0
        case (4)
          ! -- N face, subcell 4 (cell face)
          inface = 2 ! kluge note: want Domain(2) = -Domain(2); Boundary(2) = inface
          infaceoff = 0
        end select
      case (5)
        ! -- Subcell bottom (cell bottom)
        inface = npolyverts + 2 ! kluge note: want Domain(2) = -Domain(2); Boundary(2) = inface
      case (6)
        ! -- Subcell top (cell top)
        inface = npolyverts + 3 ! kluge note: want Domain(2) = -Domain(2); Boundary(2) = inface
      end select
      if (inface .eq. -1) then
        particle%iboundary(2) = 0
      else if (inface .eq. 0) then
        particle%iboundary(2) = 0
      else
        if ((inface .ge. 1) .and. (inface .le. 4)) then
          ! -- Account for local cell rotation
          inface = inface + cell%irvOrigin - 1
          if (inface .gt. 4) inface = inface - 4
          inface = cell%irectvert(inface) + infaceoff
          if (inface .lt. 1) inface = inface + npolyverts
        end if
        particle%iboundary(2) = inface
      end if
    end select