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 (this, particle)
	Pick the exit solution with the shortest travel time. More...
subroutine	find_exits (this, particle, domain)
	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 61 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) :: z_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)
      x_origin = subcell%xOrigin
      y_origin = subcell%yOrigin
      z_origin = subcell%zOrigin
      sinrot = subcell%sinrot
      cosrot = subcell%cosrot
      call particle%transform(x_origin, y_origin)
      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 342 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 41 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 55 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 322 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 ( class(methodsubcellpollocktype), intent(inout)  this, type(particletype), intent(inout), pointer  particle, class(domaintype), intent(in)  domain  )

private

Definition at line 281 of file MethodSubcellPollock.f90.

    class(MethodSubcellPollockType), intent(inout) :: this
    type(ParticleType), pointer, intent(inout) :: particle
    class(DomainType), intent(in) :: domain
    ! local
    real(DP) :: 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
      this%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 (this%exit_solutions(1)%v < dzero) then
        this%exit_solutions(1)%iboundary = 1
      else if (this%exit_solutions(1)%v > dzero) then
        this%exit_solutions(1)%iboundary = 2
      end if
      if (this%exit_solutions(2)%v < dzero) then
        this%exit_solutions(2)%iboundary = 3
      else if (this%exit_solutions(2)%v > dzero) then
        this%exit_solutions(2)%iboundary = 4
      end if
      if (this%exit_solutions(3)%v < dzero) then
        this%exit_solutions(3)%iboundary = 5
      else if (this%exit_solutions(3)%v > dzero) then
        this%exit_solutions(3)%iboundary = 6
      end if
    end select

Here is the call 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 476 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	(	class(methodsubcellpollocktype), intent(inout)	this,
		type(particletype), intent(inout), pointer	particle
	)

Definition at line 254 of file MethodSubcellPollock.f90.

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

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 97 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) :: 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_x, exit_y, exit_z
 
    t0 = particle%ttrack
    x0 = particle%x / subcell%dx
    y0 = particle%y / subcell%dy
    z0 = particle%z / subcell%dz
 
    ! Find exit solution in each direction
    call this%find_exits(particle, subcell)
 
    exit_x = this%exit_solutions(1)
    exit_y = this%exit_solutions(2)
    exit_z = this%exit_solutions(3)
 
    ! Subcell has no exit face, terminate the particle
    ! TODO: consider ramifications
    if (all([this%exit_solutions%status] == no_exit_no_outflow)) 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 (all([this%exit_solutions%status] == no_exit_stationary) .and. &
        particle%extend .and. endofsimulation) then
      call this%terminate(particle, status=term_timeout)
      return
    end if
 
    ! Pick exit solution, face, travel time, and time
    exit_soln = this%pick_exit(particle)
    if (exit_soln == 0) then
      exit_face = 0
      dtexit = 1.0d+30
    else
      exit_face = this%exit_solutions(exit_soln)%iboundary
      dtexit = this%exit_solutions(exit_soln)%dt
    end if
    texit = particle%ttrack + dtexit
 
    ! 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
 
    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(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.
 
      ! Set final particle location in local (unscaled) subcell coordinates,
      ! final time for particle trajectory
      particle%x = x * subcell%dx
      particle%y = y * subcell%dy
      particle%z = z * subcell%dz
      particle%ttrack = t
      particle%iboundary(level_subfeature) = exit_face
 
      ! Save particle track record
      call this%timestep(particle)
    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 ((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
 
      ! 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) = exit_face
 
      ! Save particle track record
      call this%subcellexit(particle)
    end if
 

Here is the call graph for this function: