methodsubcellpollockmodule Module Reference

Data Types
type	methodsubcellpollocktype
	Rectangular subcell tracking method. More...

Functions/Subroutines
subroutine, public	create_method_subcell_pollock (method)
	Create a new Pollock's subcell method. More...
subroutine	deallocate (this)
	Deallocate the Pollock's subcell method. More...
subroutine	apply_msp (this, particle, tmax)
	Apply Pollock's method to a rectangular subcell. More...
subroutine	track_subcell (this, subcell, particle, tmax)
	Track a particle across a rectangular subcell using Pollock's method. More...
integer(i4b) function, public	calculate_dt (v1, v2, dx, xL, v, dvdx, dt)
	Calculate particle travel time to exit and exit status. More...
pure real(dp) function	new_x (v, dvdx, v1, v2, dt, x, dx, velocity_profile)
	Update a cell-local coordinate based on a time increment. More...

Function/Subroutine Documentation

apply_msp()

subroutine methodsubcellpollockmodule::apply_msp ( class(methodsubcellpollocktype), intent(inout)  this, type(particletype), intent(inout), pointer  particle, real(dp), intent(in)  tmax  )

private

Definition at line 51 of file MethodSubcellPollock.f90.

    ! dummy
    class(MethodSubcellPollockType), 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 (subcell => this%subcell)
    type is (subcellrecttype)
      ! Transform particle position into local subcell coordinates,
      !    track particle across subcell, convert back to model coords
      !    (sinrot and cosrot should be 0 and 1, respectively, i.e. no
      !    rotation, also no z translation; only x and y translations)
      xorigin = subcell%xOrigin
      yorigin = subcell%yOrigin
      zorigin = subcell%zOrigin
      sinrot = subcell%sinrot
      cosrot = subcell%cosrot
      call particle%transform(xorigin, yorigin)
      call this%track_subcell(subcell, particle, tmax)
      call particle%transform(xorigin, yorigin, invert=.true.)
    end select

calculate_dt()

integer(i4b) function, public methodsubcellpollockmodule::calculate_dt	(	real(dp)	v1,
		real(dp)	v2,
		real(dp)	dx,
		real(dp)	xL,
		real(dp)	v,
		real(dp)	dvdx,
		real(dp)	dt
	)

This subroutine consists partly or entirely of code written by David W. Pollock of the USGS for MODPATH 7. The authors of the present code are responsible for its appropriate application in this context and for any modifications or errors.

Definition at line 294 of file MethodSubcellPollock.f90.

    ! dummy
    real(DP) :: v1
    real(DP) :: v2
    real(DP) :: dx
    real(DP) :: xL
    real(DP) :: v
    real(DP) :: dvdx
    real(DP) :: dt
    ! result
    integer(I4B) :: status
    ! local
    real(DP) :: v2a
    real(DP) :: v1a
    real(DP) :: dv
    real(DP) :: dva
    real(DP) :: vv
    real(DP) :: vvv
    real(DP) :: zro
    real(DP) :: zrom
    real(DP) :: x
    real(DP) :: tol
    real(DP) :: vr1
    real(DP) :: vr2
    real(DP) :: vr
    real(DP) :: v1v2
    logical(LGP) :: noOutflow
 
    ! Initialize variables.
    status = -1
    dt = 1.0d+20
    v2a = v2
    if (v2a .lt. dzero) v2a = -v2a
    v1a = v1
    if (v1a .lt. dzero) v1a = -v1a
    dv = v2 - v1
    dva = dv
    if (dva .lt. dzero) dva = -dva
 
    ! Check for a uniform zero velocity in this direction.
    ! If so, set status = 2 and return (dt = 1.0d+20).
    tol = 1.0d-15
    if ((v2a .lt. tol) .and. (v1a .lt. tol)) then
      v = dzero
      dvdx = dzero
      status = 2
      return
    end if
 
    ! Check for uniform non-zero velocity in this direction.
    ! If so, set compute dt using the constant velocity,
    ! set status = 1 and return.
    vv = v1a
    if (v2a .gt. vv) vv = v2a
    vvv = dva / vv
    if (vvv .lt. 1.0d-4) then
      zro = tol
      zrom = -zro
      v = v1
      x = xl * dx
      if (v1 .gt. zro) dt = (dx - x) / v1
      if (v1 .lt. zrom) dt = -x / v1
      dvdx = dzero
      status = 1
      return
    end if
 
    ! Velocity has a linear variation.
    ! Compute velocity corresponding to particle position.
    dvdx = dv / dx
    v = (done - xl) * v1 + xl * v2
 
    ! If flow is into the cell from both sides there is no outflow.
    ! In that case, set status = 3 and return.
    nooutflow = .true.
    if (v1 .lt. dzero) nooutflow = .false.
    if (v2 .gt. dzero) nooutflow = .false.
    if (nooutflow) then
      status = 3
      return
    end if
 
    ! If there is a divide in the cell for this flow direction, check to
    ! see if the particle is located exactly on the divide. If it is, move
    ! it very slightly to get it off the divide. This avoids possible
    ! numerical problems related to stagnation points.
    if ((v1 .le. dzero) .and. (v2 .ge. dzero)) then
      if (abs(v) .le. dzero) then
        v = 1.0d-20
        if (v2 .le. dzero) v = -v
      end if
    end if
 
    ! If there is a flow divide, this check finds out what side of the
    ! divide the particle is on and sets the value of vr appropriately
    ! to reflect that location.
    vr1 = v1 / v
    vr2 = v2 / v
    vr = vr1
    if (vr .le. dzero) then
      vr = vr2
    end if
 
    ! If the product v1*v2 > 0, the velocity is in the same direction
    ! throughout the cell (i.e. no flow divide). If so, set the value
    ! of vr to reflect the appropriate direction.
    v1v2 = v1 * v2
    if (v1v2 .gt. dzero) then
      if (v .gt. dzero) vr = vr2
      if (v .lt. dzero) vr = vr1
    end if
 
    ! Check if vr is (very close to) zero.
    ! If so, set status = 2 and return (dt = 1.0d+20).
    if (dabs(vr) .lt. 1.0d-10) then
      v = dzero
      dvdx = dzero
      status = 2
      return
    end if
 
    ! Compute travel time to exit face. Return with status = 0.
    dt = log(vr) / dvdx
    status = 0
 

Here is the caller graph for this function:

create_method_subcell_pollock()

subroutine, public methodsubcellpollockmodule::create_method_subcell_pollock

(

type(methodsubcellpollocktype), pointer

method

)

Definition at line 31 of file MethodSubcellPollock.f90.

    ! dummy
    type(MethodSubcellPollockType), pointer :: method
    ! local
    type(SubcellRectType), pointer :: subcell
 
    allocate (method)
    call create_subcell_rect(subcell)
    method%subcell => subcell
    method%name => method%subcell%type
    method%delegates = .false.

Here is the call graph for this function:

Here is the caller graph for this function:

deallocate()

subroutine methodsubcellpollockmodule::deallocate ( class(methodsubcellpollocktype), intent(inout)  this )

private

Definition at line 45 of file MethodSubcellPollock.f90.

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

new_x()

pure real(dp) function methodsubcellpollockmodule::new_x ( real(dp), intent(in)  v, real(dp), intent(in)  dvdx, real(dp), intent(in)  v1, real(dp), intent(in)  v2, real(dp), intent(in)  dt, real(dp), intent(in)  x, real(dp), intent(in)  dx, logical(lgp), intent(in), optional  velocity_profile  )

private

This subroutine consists partly or entirely of code written by David W. Pollock of the USGS for MODPATH 7. The authors of the present code are responsible for its appropriate application in this context and for any modifications or errors.

Definition at line 428 of file MethodSubcellPollock.f90.

    ! dummy
    real(DP), intent(in) :: v
    real(DP), intent(in) :: dvdx
    real(DP), intent(in) :: v1
    real(DP), intent(in) :: v2
    real(DP), intent(in) :: dt
    real(DP), intent(in) :: x
    real(DP), intent(in) :: dx
    logical(LGP), intent(in), optional :: velocity_profile
    ! result
    real(DP) :: newx
    logical(LGP) :: lprofile
 
    ! process optional arguments
    if (present(velocity_profile)) then
      lprofile = velocity_profile
    else
      lprofile = .false.
    end if
 
    ! recompute coordinate
    newx = x
    if (lprofile) then
      newx = newx + (v1 * dt / dx)
    else if (v .ne. dzero) then
      newx = newx + (v * (exp(dvdx * dt) - done) / dvdx / dx)
    end if
 
    ! clamp to [0, 1]
    if (newx .lt. dzero) newx = dzero
    if (newx .gt. done) newx = done
 

Here is the caller graph for this function:

track_subcell()

subroutine methodsubcellpollockmodule::track_subcell ( class(methodsubcellpollocktype), intent(inout)  this, class(subcellrecttype), intent(in)  subcell, type(particletype), intent(inout), pointer  particle, real(dp), intent(in)  tmax  )

private

This subroutine consists partly of code written by David W. Pollock of the USGS for MODPATH 7. PRT's authors take responsibility for its application in this context and for any modifications or errors.

Definition at line 87 of file MethodSubcellPollock.f90.

    use tdismodule, only: endofsimulation
    use particlemodule, only: active, term_no_exits_sub, term_timeout
    use particleeventmodule, only: timestep, featexit
    ! dummy
    class(MethodSubcellPollockType), intent(inout) :: this
    class(SubcellRectType), intent(in) :: subcell
    type(ParticleType), pointer, intent(inout) :: particle
    real(DP), intent(in) :: tmax
    ! local
    real(DP) :: vx
    real(DP) :: dvxdx
    real(DP) :: vy
    real(DP) :: dvydy
    real(DP) :: vz
    real(DP) :: dvzdz
    real(DP) :: dtexitx
    real(DP) :: dtexity
    real(DP) :: dtexitz
    real(DP) :: dtexit
    real(DP) :: texit
    real(DP) :: dt
    real(DP) :: t
    real(DP) :: t0
    real(DP) :: x
    real(DP) :: y
    real(DP) :: z
    integer(I4B) :: statusVX
    integer(I4B) :: statusVY
    integer(I4B) :: statusVZ
    integer(I4B) :: i
    real(DP) :: initialX
    real(DP) :: initialY
    real(DP) :: initialZ
    integer(I4B) :: exitFace
    integer(I4B) :: event_code
 
    event_code = -1
 
    ! Initial particle location in scaled subcell coordinates
    initialx = particle%x / subcell%dx
    initialy = particle%y / subcell%dy
    initialz = particle%z / subcell%dz
 
    ! Compute time of travel to each possible exit face
    statusvx = calculate_dt(subcell%vx1, subcell%vx2, subcell%dx, &
                            initialx, vx, dvxdx, dtexitx)
    statusvy = calculate_dt(subcell%vy1, subcell%vy2, subcell%dy, &
                            initialy, vy, dvydy, dtexity)
    statusvz = calculate_dt(subcell%vz1, subcell%vz2, subcell%dz, &
                            initialz, vz, dvzdz, dtexitz)
 
    ! Subcell has no exit face, terminate the particle
    ! todo: after initial release, consider ramifications
    if ((statusvx .eq. 3) .and. (statusvy .eq. 3) .and. (statusvz .eq. 3)) then
      call this%terminate(particle, status=term_no_exits_sub)
      return
    end if
 
    ! If particle stationary and extended tracking is on, terminate here if it's
    ! the last timestep. TODO: temporary solution, consider where to catch this?
    ! Should we really have to special case this here? We do diverge from MP7 in
    ! guaranteeing that every particle terminates at the end of the simulation..
    ! ideally that would be handled at a higher scope but with extended tracking
    ! tmax is not the end of the simulation, it's just a wildly high upper bound.
    if ((statusvx .eq. 2) .and. (statusvy .eq. 2) .and. (statusvz .eq. 2) .and. &
        particle%iextend > 0 .and. endofsimulation) then
      call this%terminate(particle, status=term_timeout)
      return
    end if
 
    ! Determine (earliest) exit face and corresponding travel time to exit
    exitface = 0
    dtexit = 1.0d+30
    if ((statusvx .lt. 2) .or. (statusvy .lt. 2) .or. (statusvz .lt. 2)) then
      ! Consider x-oriented faces
      dtexit = dtexitx
      if (vx .lt. dzero) then
        exitface = 1
      else if (vx .gt. 0) then
        exitface = 2
      end if
      ! Consider y-oriented faces
      if (dtexity .lt. dtexit) then
        dtexit = dtexity
        if (vy .lt. dzero) then
          exitface = 3
        else if (vy .gt. dzero) then
          exitface = 4
        end if
      end if
      ! Consider z-oriented faces
      if (dtexitz .lt. dtexit) then
        dtexit = dtexitz
        if (vz .lt. dzero) then
          exitface = 5
        else if (vz .gt. dzero) then
          exitface = 6
        end if
      end if
    else
    end if
 
    ! Compute exit time
    texit = particle%ttrack + dtexit
    t0 = particle%ttrack
 
    ! Select user tracking times to solve. If this is the first time step
    ! of the simulation, include all times before it begins; if it is the
    ! last time step include all times after it ends only if the 'extend'
    ! option is on, otherwise times in this period and time step only.
    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 < particle%ttrack) cycle
        if (t >= texit .or. t >= tmax) exit
        dt = t - t0
        x = new_x(vx, dvxdx, subcell%vx1, subcell%vx2, &
                  dt, initialx, subcell%dx, statusvx == 1)
        y = new_x(vy, dvydy, subcell%vy1, subcell%vy2, &
                  dt, initialy, subcell%dy, statusvy == 1)
        z = new_x(vz, dvzdz, subcell%vz1, subcell%vz2, &
                  dt, initialz, subcell%dz, statusvz == 1)
        particle%x = x * subcell%dx
        particle%y = y * subcell%dy
        particle%z = z * subcell%dz
        particle%ttrack = t
        particle%istatus = active
        call this%usertime(particle)
      end do
    end if
 
    if (texit .gt. tmax) then
      ! The computed exit time is greater than the maximum time, so set
      ! final time for particle trajectory equal to maximum time and
      ! calculate particle location at that final time.
      t = tmax
      dt = t - t0
      x = new_x(vx, dvxdx, subcell%vx1, subcell%vx2, &
                dt, initialx, subcell%dx, statusvx == 1)
      y = new_x(vy, dvydy, subcell%vy1, subcell%vy2, &
                dt, initialy, subcell%dy, statusvy == 1)
      z = new_x(vz, dvzdz, subcell%vz1, subcell%vz2, &
                dt, initialz, subcell%dz, statusvz == 1)
      exitface = 0
      particle%istatus = active
      particle%advancing = .false.
      event_code = timestep ! timestep end
    else
      ! The computed exit time is less than or equal to the maximum time,
      ! so set final time for particle trajectory equal to exit time and
      ! calculate exit location.
      t = texit
      dt = dtexit
      if ((exitface .eq. 1) .or. (exitface .eq. 2)) then
        x = dzero
        y = new_x(vy, dvydy, subcell%vy1, subcell%vy2, &
                  dt, initialy, subcell%dy, statusvy == 1)
        z = new_x(vz, dvzdz, subcell%vz1, subcell%vz2, &
                  dt, initialz, subcell%dz, statusvz == 1)
        if (exitface .eq. 2) x = done
      else if ((exitface .eq. 3) .or. (exitface .eq. 4)) then
        x = new_x(vx, dvxdx, subcell%vx1, subcell%vx2, dt, &
                  initialx, subcell%dx, statusvx == 1)
        y = dzero
        z = new_x(vz, dvzdz, subcell%vz1, subcell%vz2, dt, &
                  initialz, subcell%dz, statusvz == 1)
        if (exitface .eq. 4) y = done
      else if ((exitface .eq. 5) .or. (exitface .eq. 6)) then
        x = new_x(vx, dvxdx, subcell%vx1, subcell%vx2, &
                  dt, initialx, subcell%dx, statusvx == 1)
        y = new_x(vy, dvydy, subcell%vy1, subcell%vy2, &
                  dt, initialy, subcell%dy, statusvy == 1)
        z = dzero
        if (exitface .eq. 6) z = done
      else
        print *, "programmer error, invalid exit face", exitface
        call pstop(1)
      end if
      event_code = featexit
    end if
 
    ! Set final particle location in local (unscaled) subcell coordinates,
    ! final time for particle trajectory, and exit face
    particle%x = x * subcell%dx
    particle%y = y * subcell%dy
    particle%z = z * subcell%dz
    particle%ttrack = t
    particle%iboundary(level_subfeature) = exitface
 
    ! Save particle track record
    if (event_code == timestep) then
      call this%timestep(particle)
    else if (event_code == featexit) then
      call this%subcellexit(particle)
    end if
 

Here is the call graph for this function: