methoddisvmodule Module Reference

Data Types
type	methoddisvtype

Functions/Subroutines
subroutine, public	create_method_disv (method)
	Create a new vertex grid (DISV) tracking method. More...
subroutine	deallocate (this)
	Destroy the tracking method. More...
subroutine	load_disv (this, particle, next_level, submethod)
	Load the cell and the tracking method. More...
subroutine	load_particle (this, cell, particle)
subroutine	update_flowja (this, cell, particle)
subroutine	pass_disv (this, particle)
	Pass a particle to the next cell, if there is one. More...
subroutine	map_neighbor (this, defn, inface, z)
	Map location on cell face to shared cell face of neighbor. More...
subroutine	apply_disv (this, particle, tmax)
	Apply the DISV tracking method to a particle. More...
subroutine	load_cell_defn (this, ic, defn)
	Load cell definition from the grid. More...
subroutine	load_properties (this, ic, defn)
	Loads cell properties to cell definition from the grid. More...
subroutine	load_polygon (this, defn)
subroutine	load_neighbors (this, defn)
	Loads face neighbors to cell definition from the grid Assumes cell index and number of vertices are already loaded. More...
subroutine	load_flows (this, defn)
	Load flows into the cell definition. These include face, boundary and net distributed flows. Assumes cell index and number of vertices are already loaded. More...
subroutine	load_face_flows_to_defn_poly (this, defn)
subroutine	load_boundary_flows_to_defn_poly (this, defn)
	Load boundary flows from the grid into a polygonal cell. Assumes cell index and number of vertices are already loaded. More...
subroutine	load_indicators (this, defn)
	Load 180-degree vertex indicator array and set flags indicating how cell can be represented. Assumes cell index and number of vertices are already loaded. More...

Function/Subroutine Documentation

apply_disv()

subroutine methoddisvmodule::apply_disv ( class(methoddisvtype), intent(inout)  this, type(particletype), intent(inout), pointer  particle, real(dp), intent(in)  tmax  )

private

Definition at line 321 of file MethodDisv.f90.

    class(MethodDisvType), intent(inout) :: this
    type(ParticleType), pointer, intent(inout) :: particle
    real(DP), intent(in) :: tmax
 
    call this%track(particle, 1, tmax)

create_method_disv()

subroutine, public methoddisvmodule::create_method_disv

(

type(methoddisvtype), pointer

method

)

Definition at line 46 of file MethodDisv.f90.

    ! dummy
    type(MethodDisvType), pointer :: method
    ! local
    type(CellPolyType), pointer :: cell
 
    allocate (method)
    allocate (method%name)
    call create_cell_poly(cell)
    method%cell => cell
    method%name = "disv"
    method%delegates = .true.
    call create_defn(method%neighbor)

Here is the call graph for this function:

Here is the caller graph for this function:

deallocate()

subroutine methoddisvmodule::deallocate ( class(methoddisvtype), intent(inout)  this )

private

Definition at line 62 of file MethodDisv.f90.

    class(MethodDisvType), intent(inout) :: this
    deallocate (this%name)

load_boundary_flows_to_defn_poly()

subroutine methoddisvmodule::load_boundary_flows_to_defn_poly ( class(methoddisvtype), intent(inout)  this, type(celldefntype), intent(inout)  defn  )

private

Definition at line 519 of file MethodDisv.f90.

    ! dummy
    class(MethodDisvType), intent(inout) :: this
    type(CellDefnType), intent(inout) :: defn
 
    ! local
    integer(I4B) :: ic, iv, ioffset, npolyverts, max_faces
 
    ic = defn%icell
    npolyverts = defn%npolyverts
    max_faces = this%fmi%max_faces
    ioffset = (ic - 1) * max_faces
    do iv = 1, npolyverts
      defn%faceflow(iv) = &
        defn%faceflow(iv) + &
        this%fmi%BoundaryFlows(ioffset + iv)
    end do
    defn%faceflow(npolyverts + 1) = defn%faceflow(1)
    defn%faceflow(npolyverts + 2) = &
      defn%faceflow(npolyverts + 2) + &
      this%fmi%BoundaryFlows(ioffset + max_faces - 1)
    defn%faceflow(npolyverts + 3) = &
      defn%faceflow(npolyverts + 3) + &
      this%fmi%BoundaryFlows(ioffset + max_faces)
 

load_cell_defn()

subroutine methoddisvmodule::load_cell_defn ( class(methoddisvtype), intent(inout)  this, integer(i4b), intent(in)  ic, type(celldefntype), intent(inout), pointer  defn  )

private

Definition at line 330 of file MethodDisv.f90.

    ! dummy
    class(MethodDisvType), intent(inout) :: this
    integer(I4B), intent(in) :: ic
    type(CellDefnType), pointer, intent(inout) :: defn
 
    call this%load_properties(ic, defn)
    call this%load_polygon(defn)
    call this%load_neighbors(defn)
    call this%load_indicators(defn)
    call this%load_flows(defn)

load_disv()

subroutine methoddisvmodule::load_disv ( class(methoddisvtype), intent(inout)  this, type(particletype), intent(inout), pointer  particle, integer(i4b), intent(in)  next_level, class(methodtype), intent(inout), pointer  submethod  )

private

Definition at line 68 of file MethodDisv.f90.

    use cellmodule
    use cellrectmodule
    use cellrectquadmodule
    use cellutilmodule
    ! dummy
    class(MethodDisvType), intent(inout) :: this
    type(ParticleType), pointer, intent(inout) :: particle
    integer(I4B), intent(in) :: next_level
    class(MethodType), pointer, intent(inout) :: submethod
    ! local
    integer(I4B) :: ic
    class(CellType), pointer :: base
    type(CellRectType), pointer :: rect
    type(CellRectQuadType), pointer :: quad
 
    select type (cell => this%cell)
    type is (cellpolytype)
      ic = particle%itrdomain(next_level)
      call this%load_cell_defn(ic, cell%defn)
      if (this%fmi%ibdgwfsat0(ic) == 0) then
        ! Cell is active but dry, so select and initialize pass-to-bottom
        ! cell method and set cell method pointer
        call method_cell_ptb%init( &
          fmi=this%fmi, &
          cell=this%cell, &
          events=this%events, &
          tracktimes=this%tracktimes)
        submethod => method_cell_ptb
      else if (particle%ifrctrn > 0) then
        ! Force the ternary method
        call method_cell_tern%init( &
          fmi=this%fmi, &
          cell=this%cell, &
          events=this%events, &
          tracktimes=this%tracktimes)
        submethod => method_cell_tern
      else if (cell%defn%can_be_rect) then
        ! Cell is a rectangle, convert it to a rectangular cell type and
        ! initialize Pollock's method
        call cell_poly_to_rect(cell, rect)
        base => rect
        call method_cell_plck%init( &
          fmi=this%fmi, &
          cell=base, &
          events=this%events, &
          tracktimes=this%tracktimes)
        submethod => method_cell_plck
      else if (cell%defn%can_be_quad) then
        ! Cell is quad-refined, convert to a quad rect cell type and
        ! initialize the corresponding method
        call cell_poly_to_quad(cell, quad)
        base => quad
        call method_cell_quad%init( &
          fmi=this%fmi, &
          cell=base, &
          events=this%events, &
          tracktimes=this%tracktimes)
        submethod => method_cell_quad
      else
        ! Default to the ternary method
        call method_cell_tern%init( &
          fmi=this%fmi, &
          cell=this%cell, &
          events=this%events, &
          tracktimes=this%tracktimes)
        submethod => method_cell_tern
      end if
    end select

Here is the call graph for this function:

load_face_flows_to_defn_poly()

subroutine methoddisvmodule::load_face_flows_to_defn_poly ( class(methoddisvtype), intent(inout)  this, type(celldefntype), intent(inout)  defn  )

private

Definition at line 498 of file MethodDisv.f90.

    ! dummy
    class(MethodDisvType), intent(inout) :: this
    type(CellDefnType), intent(inout) :: defn
    ! local
    integer(I4B) :: m, n, nfaces
    real(DP) :: q
 
    nfaces = defn%npolyverts + 3
    do m = 1, nfaces
      n = defn%facenbr(m)
      if (n > 0) then
        q = this%fmi%gwfflowja(this%fmi%dis%con%ia(defn%icell) + n)
        defn%faceflow(m) = defn%faceflow(m) + q
      end if
      if (defn%faceflow(m) < dzero) defn%inoexitface = 0
    end do

load_flows()

subroutine methoddisvmodule::load_flows ( class(methoddisvtype), intent(inout)  this, type(celldefntype), intent(inout)  defn  )

private

Definition at line 464 of file MethodDisv.f90.

    ! dummy
    class(MethodDisvType), intent(inout) :: this
    type(CellDefnType), intent(inout) :: defn
    ! local
    integer(I4B) :: nfaces, nslots
 
    ! expand faceflow array if needed
    nfaces = defn%npolyverts + 3
    nslots = size(defn%faceflow)
    if (nslots < nfaces) call expandarray(defn%faceflow, nfaces - nslots)
 
    ! Load face flows, including boundary flows. As with cell verts,
    ! the face flow array wraps around. Top and bottom flows make up
    ! the last two elements, respectively, for size npolyverts + 3.
    ! If there is no flow through any face, set a no-exit-face flag.
    defn%faceflow = dzero
    defn%inoexitface = 1
    call this%load_boundary_flows_to_defn_poly(defn)
    call this%load_face_flows_to_defn_poly(defn)
 
    ! Add up net distributed flow
    defn%distflow = this%fmi%SourceFlows(defn%icell) + &
                    this%fmi%SinkFlows(defn%icell) + &
                    this%fmi%StorageFlows(defn%icell)
 
    ! Set weak sink flag
    if (this%fmi%SinkFlows(defn%icell) .ne. dzero) then
      defn%iweaksink = 1
    else
      defn%iweaksink = 0
    end if

load_indicators()

subroutine methoddisvmodule::load_indicators ( class(methoddisvtype), intent(inout)  this, type(celldefntype), intent(inout), pointer  defn  )

private

Definition at line 550 of file MethodDisv.f90.

    ! dummy
    class(MethodDisvType), intent(inout) :: this
    type(CellDefnType), pointer, intent(inout) :: defn
    ! local
    integer(I4B) :: npolyverts
    integer(I4B) :: m
    integer(I4B) :: m0
    integer(I4B) :: m1
    integer(I4B) :: m2
    integer(I4B) :: ic
    integer(I4B) :: num90
    integer(I4B) :: num180
    real(DP) :: x0
    real(DP) :: y0
    real(DP) :: x1
    real(DP) :: y1
    real(DP) :: x2
    real(DP) :: y2
    real(DP) :: epsang
    real(DP) :: s0x
    real(DP) :: s0y
    real(DP) :: &
      s0mag, s2x, s2y, s2mag, sinang, cosang, dotprod
    logical(LGP) last180
 
    ic = defn%icell
    npolyverts = defn%npolyverts
 
    if (size(defn%ispv180) < npolyverts + 3) &
      call expandarray(defn%ispv180, npolyverts + 1)
 
    defn%ispv180(1:npolyverts + 1) = .false.
    defn%can_be_rect = .false.
    defn%can_be_quad = .false.
    epsang = 1d-5
    num90 = 0
    num180 = 0
    last180 = .false.
 
    ! assumes non-self-intersecting polygon
    do m = 1, npolyverts
      m1 = m
      if (m1 .eq. 1) then
        m0 = npolyverts
        m2 = 2
      else if (m1 .eq. npolyverts) then
        m0 = npolyverts - 1
        m2 = 1
      else
        m0 = m1 - 1
        m2 = m1 + 1
      end if
      x0 = defn%polyvert(1, m0)
      y0 = defn%polyvert(2, m0)
      x1 = defn%polyvert(1, m1)
      y1 = defn%polyvert(2, m1)
      x2 = defn%polyvert(1, m2)
      y2 = defn%polyvert(2, m2)
      s0x = x0 - x1
      s0y = y0 - y1
      s0mag = dsqrt(s0x * s0x + s0y * s0y)
      s2x = x2 - x1
      s2y = y2 - y1
      s2mag = dsqrt(s2x * s2x + s2y * s2y)
      sinang = (s0x * s2y - s0y * s2x) / (s0mag * s2mag)
      cosang = dsqrt(dabs(done - (sinang * sinang)))
      if (dabs(sinang) .lt. epsang) then
        dotprod = s0x * s2x + s0y * s2y
        if (dotprod .lt. dzero) then
          num180 = num180 + 1
          last180 = .true.
          defn%ispv180(m) = .true.
        end if
      else
        if (dabs(cosang) .lt. epsang) num90 = num90 + 1
        last180 = .false.
      end if
    end do
 
    ! List of 180-degree indicators wraps around for convenience
    defn%ispv180(npolyverts + 1) = defn%ispv180(1)
 
    ! Set rect/quad flags
    if (num90 .eq. 4) then
      if (num180 .eq. 0) then
        defn%can_be_rect = .true.
      else
        defn%can_be_quad = .true.
      end if
    end if

load_neighbors()

subroutine methoddisvmodule::load_neighbors ( class(methoddisvtype), intent(inout)  this, type(celldefntype), intent(inout), pointer  defn  )

private

Definition at line 390 of file MethodDisv.f90.

    ! dummy
    class(MethodDisvType), intent(inout) :: this
    type(CellDefnType), pointer, intent(inout) :: defn
    ! local
    integer(I4B) :: ic1
    integer(I4B) :: ic2
    integer(I4B) :: icu1
    integer(I4B) :: icu2
    integer(I4B) :: j1
    integer(I4B) :: j2
    integer(I4B) :: k1
    integer(I4B) :: k2
    integer(I4B) :: iloc
    integer(I4B) :: ipos
    integer(I4B) :: istart1
    integer(I4B) :: istart2
    integer(I4B) :: istop1
    integer(I4B) :: istop2
    integer(I4B) :: isharedface
    integer(I4B) :: ncpl
    integer(I4B) :: nfaces
    integer(I4B) :: nslots
 
    ! expand facenbr array if needed
    nfaces = defn%npolyverts + 3
    nslots = size(defn%facenbr)
    if (nslots < nfaces) call expandarray(defn%facenbr, nfaces - nslots)
 
    select type (dis => this%fmi%dis)
    type is (disvtype)
      ! Load face neighbors.
      defn%facenbr = 0
      ic1 = defn%icell
      icu1 = dis%get_nodeuser(ic1)
      ncpl = dis%get_ncpl()
      call get_jk(icu1, ncpl, dis%nlay, j1, k1)
      istart1 = dis%iavert(j1)
      istop1 = dis%iavert(j1 + 1) - 1
      do iloc = 1, dis%con%ia(ic1 + 1) - dis%con%ia(ic1) - 1
        ipos = dis%con%ia(ic1) + iloc
        ! mask could become relevant if PRT uses interface model
        if (dis%con%mask(ipos) == 0) cycle
        ic2 = dis%con%ja(ipos)
        icu2 = dis%get_nodeuser(ic2)
        call get_jk(icu2, ncpl, dis%nlay, j2, k2)
        istart2 = dis%iavert(j2)
        istop2 = dis%iavert(j2 + 1) - 1
        call shared_face(dis%javert(istart1:istop1), &
                         dis%javert(istart2:istop2), &
                         isharedface)
        if (isharedface /= 0) then
          ! Edge (polygon) face neighbor
          defn%facenbr(isharedface) = int(iloc, 1)
        else
          if (k2 > k1) then
            ! Bottom face neighbor
            defn%facenbr(defn%npolyverts + 2) = int(iloc, 1)
          else if (k2 < k1) then
            ! Top face neighbor
            defn%facenbr(defn%npolyverts + 3) = int(iloc, 1)
          else
            call pstop(1, "k2 should be <> k1, since no shared edge face")
          end if
        end if
      end do
    end select
    ! List of edge (polygon) faces wraps around
    defn%facenbr(defn%npolyverts + 1) = defn%facenbr(1)

Here is the call graph for this function:

load_particle()

subroutine methoddisvmodule::load_particle	(	class(methoddisvtype), intent(inout)	this,
		type(cellpolytype), intent(inout), pointer	cell,
		type(particletype), intent(inout), pointer	particle
	)

Definition at line 139 of file MethodDisv.f90.

    ! modules
    use disvmodule, only: disvtype
    use particlemodule, only: term_boundary
    use particleeventmodule, only: terminate
    ! dummy
    class(MethodDisvType), intent(inout) :: this
    type(CellPolyType), pointer, intent(inout) :: cell
    type(ParticleType), pointer, intent(inout) :: particle
    ! local
    integer(I4B) :: inface
    integer(I4B) :: ipos
    integer(I4B) :: ic
    integer(I4B) :: icu
    integer(I4B) :: inbr
    integer(I4B) :: idiag
    integer(I4B) :: icpl
    integer(I4B) :: ilay
    real(DP) :: z
 
    select type (dis => this%fmi%dis)
    type is (disvtype)
      inface = particle%iboundary(level_feature)
      idiag = dis%con%ia(cell%defn%icell)
      inbr = cell%defn%facenbr(inface)
      ipos = idiag + inbr
      ic = dis%con%ja(ipos)
      icu = dis%get_nodeuser(ic)
      call get_jk(icu, dis%ncpl, dis%nlay, icpl, ilay)
 
      ! if returning to a cell through the bottom
      ! face after previously leaving it through
      ! that same face, we've entered a cycle
      ! as can occur e.g. in wells. terminate
      ! in the previous cell.
      if (ic == particle%icp .and. inface == 7 .and. ilay < particle%ilay) then
        particle%itrdomain(level_feature) = particle%icp
        particle%izone = particle%izp
        call this%terminate(particle, &
                            status=term_boundary)
        return
      else
        particle%icp = particle%itrdomain(level_feature)
        particle%izp = particle%izone
      end if
 
      particle%itrdomain(level_feature) = ic
      particle%icu = icu
      particle%ilay = ilay
 
      z = particle%z
      call this%map_neighbor(cell%defn, inface, z)
 
      particle%iboundary(level_feature) = inface
      particle%itrdomain(level_subfeature:) = 0
      particle%iboundary(level_subfeature:) = 0
      particle%z = z
    end select
 

Here is the call graph for this function:

load_polygon()

subroutine methoddisvmodule::load_polygon ( class(methoddisvtype), intent(inout)  this, type(celldefntype), intent(inout), pointer  defn  )

private

Definition at line 376 of file MethodDisv.f90.

    ! dummy
    class(MethodDisvType), intent(inout) :: this
    type(CellDefnType), pointer, intent(inout) :: defn
 
    call this%fmi%dis%get_polyverts( &
      defn%icell, &
      defn%polyvert, &
      closed=.true.)
    defn%npolyverts = size(defn%polyvert, dim=2) - 1

load_properties()

subroutine methoddisvmodule::load_properties ( class(methoddisvtype), intent(inout)  this, integer(i4b), intent(in)  ic, type(celldefntype), intent(inout), pointer  defn  )

private

Definition at line 344 of file MethodDisv.f90.

    ! dummy
    class(MethodDisvType), intent(inout) :: this
    integer(I4B), intent(in) :: ic
    type(CellDefnType), pointer, intent(inout) :: defn
    ! local
    real(DP) :: top
    real(DP) :: bot
    real(DP) :: sat
    integer(I4B) :: icu, icpl, ilay
 
    defn%icell = ic
    defn%iatop = get_iatop(this%fmi%dis%get_ncpl(), &
                           this%fmi%dis%get_nodeuser(ic))
    top = this%fmi%dis%top(ic)
    bot = this%fmi%dis%bot(ic)
    sat = this%fmi%gwfsat(ic)
    top = bot + sat * (top - bot)
    defn%top = top
    defn%bot = bot
    defn%sat = sat
    defn%porosity = this%porosity(ic)
    defn%retfactor = this%retfactor(ic)
    defn%izone = this%izone(ic)
    select type (dis => this%fmi%dis)
    type is (disvtype)
      icu = dis%get_nodeuser(ic)
      call get_jk(icu, dis%ncpl, dis%nlay, icpl, ilay)
      defn%ilay = ilay
    end select

Here is the call graph for this function:

map_neighbor()

subroutine methoddisvmodule::map_neighbor	(	class(methoddisvtype), intent(inout)	this,
		type(celldefntype), intent(inout), pointer	defn,
		integer(i4b), intent(inout)	inface,
		double precision, intent(inout)	z
	)

Definition at line 255 of file MethodDisv.f90.

    ! dummy
    class(MethodDisvType), intent(inout) :: this
    type(CellDefnType), pointer, intent(inout) :: defn
    integer(I4B), intent(inout) :: inface
    double precision, intent(inout) :: z
    ! local
    integer(I4B) :: icin
    integer(I4B) :: npolyvertsin
    integer(I4B) :: ic
    integer(I4B) :: npolyverts
    integer(I4B) :: inbr
    integer(I4B) :: inbrnbr
    integer(I4B) :: j
    integer(I4B) :: m
    real(DP) :: zrel
    real(DP) :: topfrom
    real(DP) :: botfrom
    real(DP) :: top
    real(DP) :: bot
    real(DP) :: sat
 
    ! Map to shared cell face of neighbor
    inbr = defn%facenbr(inface)
    if (inbr .eq. 0) then
      ! Exterior face; no neighbor to map to
      ! todo AMP: reconsider when multiple models allowed
      inface = -1
    else
      ! Load definition for neighbor cell (neighbor with shared face)
      icin = defn%icell
      j = this%fmi%dis%con%ia(icin)
      ic = this%fmi%dis%con%ja(j + inbr)
      call this%load_cell_defn(ic, this%neighbor)
 
      npolyvertsin = defn%npolyverts
      npolyverts = this%neighbor%npolyverts
      if (inface .eq. npolyvertsin + 2) then
        ! Exits through bot, enters through top
        inface = npolyverts + 3
      else if (inface .eq. npolyvertsin + 3) then
        ! Exits through top, enters through bot
        inface = npolyverts + 2
      else
        ! Exits and enters through shared polygon face
        j = this%fmi%dis%con%ia(ic)
        do m = 1, npolyverts + 3
          inbrnbr = this%neighbor%facenbr(m)
          if (this%fmi%dis%con%ja(j + inbrnbr) .eq. icin) then
            inface = m
            exit
          end if
        end do
        ! Map z between cells
        topfrom = defn%top
        botfrom = defn%bot
        zrel = (z - botfrom) / (topfrom - botfrom)
        top = this%fmi%dis%top(ic)
        bot = this%fmi%dis%bot(ic)
        sat = this%fmi%gwfsat(ic)
        z = bot + zrel * sat * (top - bot)
      end if
    end if

pass_disv()

subroutine methoddisvmodule::pass_disv ( class(methoddisvtype), intent(inout)  this, type(particletype), intent(inout), pointer  particle  )

private

Definition at line 224 of file MethodDisv.f90.

    use particlemodule, only: term_boundary
    use particleeventmodule, only: terminate
    ! dummy
    class(MethodDisvType), intent(inout) :: this
    type(ParticleType), pointer, intent(inout) :: particle
    ! local
    type(CellPolyType), pointer :: cell
 
    select type (c => this%cell)
    type is (cellpolytype)
      cell => c
      ! If the entry face has no neighbors it's a
      ! boundary face, so terminate the particle.
      ! todo AMP: reconsider when multiple models supported
      if (cell%defn%facenbr(particle%iboundary(level_feature)) .eq. 0) then
        call this%terminate(particle, &
                            status=term_boundary)
      else
        ! Otherwise, load cell properties into the
        ! particle. It may be marked to terminate.
        call this%load_particle(cell, particle)
        if (.not. particle%advancing) return
 
        ! Update intercell mass flows
        call this%update_flowja(cell, particle)
      end if
    end select

update_flowja()

subroutine methoddisvmodule::update_flowja	(	class(methoddisvtype), intent(inout)	this,
		type(cellpolytype), intent(inout), pointer	cell,
		type(particletype), intent(inout), pointer	particle
	)

Definition at line 200 of file MethodDisv.f90.

    ! dummy
    class(MethodDisvType), intent(inout) :: this
    type(CellPolyType), pointer, intent(inout) :: cell
    type(ParticleType), pointer, intent(inout) :: particle
    ! local
    integer(I4B) :: inbr
    integer(I4B) :: idiag
    integer(I4B) :: ipos
 
    idiag = this%fmi%dis%con%ia(cell%defn%icell)
    inbr = cell%defn%facenbr(particle%iboundary(level_feature))
    ipos = idiag + inbr
 
    ! leaving old cell
    this%flowja(ipos) = this%flowja(ipos) - done
 
    ! entering new cell
    this%flowja(this%fmi%dis%con%isym(ipos)) &
      = this%flowja(this%fmi%dis%con%isym(ipos)) + done