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	pick_exit (exit_solutions)
	Pick the exit solution with the shortest travel time. More...
subroutine	find_exits (particle, domain, exit_solutions, x0, y0, z0)
	Compute candidate exit solutions. More...
type(linearexitsolutiontype) function	find_exit (v1, v2, dx, xL)
	Find an exit solution for one dimension. 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 58 of file MethodSubcellPollock.f90.

    ! dummy
    class(MethodSubcellPollockType), intent(inout) :: this
    type(ParticleType), pointer, intent(inout) :: particle
    real(DP), intent(in) :: tmax
    ! local
    real(DP) :: x_origin
    real(DP) :: y_origin
    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,
      ! since subcell%zOrigin is always 0 -- z is translated once,
      ! at the cell level, not again here)
      x_origin = subcell%xOrigin
      y_origin = subcell%yOrigin
      sinrot = subcell%sinrot
      cosrot = subcell%cosrot
      call particle%transform(x_origin, y_origin)
 
      ! Clamp the particle into the subcell in case roundoff in the
      ! model-to-local coordinate transform (e.g. under grid rotation)
      ! left it just outside, analogous to the nudge applied to the
      ! ternary method's barycentric coordinates.
      particle%x = min(max(particle%x, dzero), subcell%dx)
      particle%y = min(max(particle%y, dzero), subcell%dy)
      particle%z = min(max(particle%z, dzero), subcell%dz)
 
      call this%track_subcell(subcell, particle, tmax)
      call particle%transform(x_origin, y_origin, 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 of code written by and/or adapted from 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 335 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 = no_exit_stationary
      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 = ok_exit_constant
      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 = no_exit_no_outflow
      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 = no_exit_stationary
      return
    end if
 
    ! Compute travel time to exit face. Return with status = 0.
    dt = log(vr) / dvdx
    status = ok_exit
 

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 38 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 52 of file MethodSubcellPollock.f90.

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

find_exit()

type(linearexitsolutiontype) function methodsubcellpollockmodule::find_exit ( real(dp), intent(in)  v1, real(dp), intent(in)  v2, real(dp), intent(in)  dx, real(dp), intent(in)  xL  )

private

Definition at line 315 of file MethodSubcellPollock.f90.

    ! dummy
    real(DP), intent(in) :: v1
    real(DP), intent(in) :: v2
    real(DP), intent(in) :: dx
    real(DP), intent(in) :: xL
    type(LinearExitSolutionType) :: solution
 
    solution = linearexitsolutiontype()
    solution%status = calculate_dt(v1, v2, dx, xl, &
                                   solution%v, solution%dvdx, solution%dt)

Here is the call graph for this function:

Here is the caller graph for this function:

find_exits()

subroutine methodsubcellpollockmodule::find_exits ( type(particletype), intent(inout), pointer  particle, class(domaintype), intent(in)  domain, type(linearexitsolutiontype), dimension(3), intent(out)  exit_solutions, real(dp), intent(out)  x0, real(dp), intent(out)  y0, real(dp), intent(out)  z0  )

private

Also returns the local subcell coordinates used (x0, y0, z0), so track_subcell doesn't need to recompute them from the particle's position.

Definition at line 275 of file MethodSubcellPollock.f90.

    type(ParticleType), pointer, intent(inout) :: particle
    class(DomainType), intent(in) :: domain
    type(LinearExitSolutionType), intent(out) :: exit_solutions(3)
    real(DP), intent(out) :: x0, y0, z0
 
    select type (domain)
    type is (subcellrecttype)
      ! Initial particle location in scaled subcell coordinates
      x0 = particle%x / domain%dx
      y0 = particle%y / domain%dy
      z0 = particle%z / domain%dz
 
      ! Calculate exit solutions for each coordinate direction
      exit_solutions = [ &
                       find_exit(domain%vx1, domain%vx2, domain%dx, x0), &
                       find_exit(domain%vy1, domain%vy2, domain%dy, y0), &
                       find_exit(domain%vz1, domain%vz2, domain%dz, z0) &
                       ]
 
      ! Set exit faces
      if (exit_solutions(1)%v < dzero) then
        exit_solutions(1)%iboundary = 1
      else if (exit_solutions(1)%v > dzero) then
        exit_solutions(1)%iboundary = 2
      end if
      if (exit_solutions(2)%v < dzero) then
        exit_solutions(2)%iboundary = 3
      else if (exit_solutions(2)%v > dzero) then
        exit_solutions(2)%iboundary = 4
      end if
      if (exit_solutions(3)%v < dzero) then
        exit_solutions(3)%iboundary = 5
      else if (exit_solutions(3)%v > dzero) then
        exit_solutions(3)%iboundary = 6
      end if
    end select

Here is the call graph for this function:

Here is the caller graph for this function:

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 469 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:

pick_exit()

integer(i4b) function methodsubcellpollockmodule::pick_exit

(

type(linearexitsolutiontype), dimension(3), intent(in)

exit_solutions

)

Definition at line 244 of file MethodSubcellPollock.f90.

    type(LinearExitSolutionType), intent(in) :: exit_solutions(3)
    integer(I4B) :: exit_soln
    ! local
    real(DP) :: dtmin
 
    exit_soln = 0
    dtmin = 1.0d+30
 
    if (exit_solutions(1)%status < 2) then
      exit_soln = 1 ! x
      dtmin = exit_solutions(1)%dt
    end if
    if (exit_solutions(2)%status < 2 .and. &
        exit_solutions(2)%dt < dtmin) then
      exit_soln = 2 ! y
      dtmin = exit_solutions(2)%dt
    end if
    if (exit_solutions(3)%status < 2 .and. &
        exit_solutions(3)%dt < dtmin) then
      exit_soln = 3 ! z
    end if
 

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 103 of file MethodSubcellPollock.f90.

    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) :: dt, dtexit, texit
    real(DP) :: t, x, y, z
    real(DP) :: t0, x0, y0, z0
    integer(I4B) :: i, exit_face, exit_soln
    type(LinearExitSolutionType) :: exit_solutions(3)
    type(LinearExitSolutionType) :: exit_x, exit_y, exit_z
 
    t0 = particle%ttrack
 
    ! Find exit solution in each direction; also returns the
    ! local subcell coordinates used, so we don't recompute them
    call find_exits(particle, subcell, exit_solutions, x0, y0, z0)
    exit_x = exit_solutions(1)
    exit_y = exit_solutions(2)
    exit_z = exit_solutions(3)
 
    ! Set solution, face, & travel time
    exit_soln = pick_exit(exit_solutions)
    if (exit_soln == 0) then
      exit_face = 0
      dtexit = 1.0d+30
    else
      exit_face = exit_solutions(exit_soln)%iboundary
      dtexit = exit_solutions(exit_soln)%dt
    end if
    texit = particle%ttrack + dtexit
 
    ! Terminate if no valid exit solution was found.
    ! MP7 is more nuanced in determining what to do
    ! here. It considers whether this stress period
    ! is steady state or transient, and whether the
    ! flow is stationary, in determining whether to
    ! terminate or allow the particle to survive to
    ! the next time step. It does not compute paths
    ! within the subcell even for particles it lets
    ! remain active under this circumstance, though.
    ! While we may consider that someday, we simply
    ! terminate and sidestep the complexity for now.
    if (all([exit_solutions%status] >= no_exit_stationary)) then
      call this%terminate(particle, status=term_no_exits_sub)
      return
    end if
 
    ! 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(exit_x%v, exit_x%dvdx, subcell%vx1, subcell%vx2, &
                  dt, x0, subcell%dx, exit_x%status == 1)
        y = new_x(exit_y%v, exit_y%dvdx, subcell%vy1, subcell%vy2, &
                  dt, y0, subcell%dy, exit_y%status == 1)
        z = new_x(exit_z%v, exit_z%dvdx, subcell%vz1, subcell%vz2, &
                  dt, z0, subcell%dz, exit_z%status == 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
 
    ! Computed exit time greater than the maximum time? Set the
    ! tracking time to tmax and calculate the particle location.
    if (texit .gt. tmax) then
      t = tmax
      dt = t - t0
      x = new_x(exit_x%v, exit_x%dvdx, subcell%vx1, subcell%vx2, &
                dt, x0, subcell%dx, exit_x%status == 1)
      y = new_x(exit_y%v, exit_y%dvdx, subcell%vy1, subcell%vy2, &
                dt, y0, subcell%dy, exit_y%status == 1)
      z = new_x(exit_z%v, exit_z%dvdx, subcell%vz1, subcell%vz2, &
                dt, z0, subcell%dz, exit_z%status == 1)
      exit_face = 0
      particle%istatus = active
      particle%advancing = .false.
      particle%x = x * subcell%dx
      particle%y = y * subcell%dy
      particle%z = z * subcell%dz
      particle%ttrack = t
      particle%iboundary(level_subfeature) = exit_face
      call this%timestep(particle)
      return
    end if
 
    ! If we get to here, the particle is exiting the subcell.
    ! Its exit time is less than or equal to the maximum time.
    ! Set tracking time to the exit time and set exit location.
    t = texit
    dt = dtexit
    if ((exit_face .eq. 1) .or. (exit_face .eq. 2)) then
      x = dzero
      y = new_x(exit_y%v, exit_y%dvdx, subcell%vy1, subcell%vy2, &
                dt, y0, subcell%dy, exit_y%status == 1)
      z = new_x(exit_z%v, exit_z%dvdx, subcell%vz1, subcell%vz2, &
                dt, z0, subcell%dz, exit_z%status == 1)
      if (exit_face .eq. 2) x = done
    else if ((exit_face .eq. 3) .or. (exit_face .eq. 4)) then
      x = new_x(exit_x%v, exit_x%dvdx, subcell%vx1, subcell%vx2, dt, &
                x0, subcell%dx, exit_x%status == 1)
      y = dzero
      z = new_x(exit_z%v, exit_z%dvdx, subcell%vz1, subcell%vz2, dt, &
                z0, subcell%dz, exit_z%status == 1)
      if (exit_face .eq. 4) y = done
    else if ((exit_face .eq. 5) .or. (exit_face .eq. 6)) then
      x = new_x(exit_x%v, exit_x%dvdx, subcell%vx1, subcell%vx2, &
                dt, x0, subcell%dx, exit_x%status == 1)
      y = new_x(exit_y%v, exit_y%dvdx, subcell%vy1, subcell%vy2, &
                dt, y0, subcell%dy, exit_y%status == 1)
      z = dzero
      if (exit_face .eq. 6) z = done
    else
      print *, "programmer error, invalid exit face", exit_face
      call pstop(1)
    end if
    particle%x = x * subcell%dx
    particle%y = y * subcell%dy
    particle%z = z * subcell%dz
    particle%ttrack = t
    particle%iboundary(level_subfeature) = exit_face
    call this%subcellexit(particle)
 

Here is the call graph for this function: