methodsubcellternarymodule Module Reference

Data Types
type	barycentricexitsolutiontype
	Barycentric velocity interpolation exit solution. Inherit from LinearExitSolutionType to get around array polymorphism limitations in Fortran; the exit_solutions array below needs to be of one type for convenient use. More...
type	methodsubcellternarytype
	Ternary triangular subcell tracking method. More...

Functions/Subroutines
subroutine, public	create_method_subcell_ternary (method)
	Create a new ternary subcell tracking method. More...
subroutine	deallocate (this)
	Deallocate the ternary subcell tracking method. More...
subroutine	apply_mst (this, particle, tmax)
	Apply the ternary subcell tracking method. More...
subroutine	track_subcell (this, subcell, particle, tmax)
	Track a particle across a triangular subcell. More...
integer(i4b) function	pick_exit (exit_solutions)
	Determine earliest exit face. More...
subroutine	find_exits (particle, domain, exit_solutions, z0)
	Calculate exit solutions for each coordinate direction. More...
type(barycentricexitsolutiontype) function	find_lateral_exit (isolv, tol, itrifaceenter, alp1, bet1, alp2, bet2, alpi, beti)
type(barycentricexitsolutiontype) function	find_vertical_exit (v1, v2, dx, xL)
subroutine	calculate_dt (v1, v2, dx, xL, v, dvdx, dt, status, itopbotexit)
	Do calculations related to analytical z solution. More...
subroutine	calculate_xyz_position (dt, rxx, rxy, ryx, ryy, sxx, sxy, syy, izstatus, x0, y0, az, vzi, vzbot, ztop, zbot, zi, x, y, z, exitface)
	Calculate the particle's local unscaled xyz coordinates after dt. More...

Function/Subroutine Documentation

apply_mst()

subroutine methodsubcellternarymodule::apply_mst ( class(methodsubcellternarytype), intent(inout)  this, type(particletype), intent(inout), pointer  particle, real(dp), intent(in)  tmax  )

private

Definition at line 77 of file MethodSubcellTernary.f90.

    class(MethodSubcellTernaryType), intent(inout) :: this
    type(ParticleType), pointer, intent(inout) :: particle
    real(DP), intent(in) :: tmax
 
    select type (subcell => this%subcell)
    type is (subcelltritype)
      call this%track_subcell(subcell, particle, tmax)
    end select

calculate_dt()

subroutine methodsubcellternarymodule::calculate_dt ( real(dp)  v1, real(dp)  v2, real(dp)  dx, real(dp)  xL, real(dp)  v, real(dp)  dvdx, real(dp)  dt, integer(i4b)  status, integer(i4b)  itopbotexit  )

private

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 379 of file MethodSubcellTernary.f90.

    real(DP) :: v1
    real(DP) :: v2
    real(DP) :: dx
    real(DP) :: xL
    real(DP) :: v
    real(DP) :: dvdx
    real(DP) :: dt
    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
    integer(I4B) :: status
    integer(I4B) :: itopbotexit
    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
      itopbotexit = 0
      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) then
        dt = (dx - x) / v1
        itopbotexit = -2
      end if
      if (v1 .lt. zrom) then
        dt = -x / v1
        itopbotexit = -1
      end if
      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
      itopbotexit = 0
      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, find out what side of the divide the particle
    ! is on and set the value of vr appropriately to reflect that location.
    vr1 = v1 / v
    vr2 = v2 / v
    vr = vr1
    itopbotexit = -1
    if (vr .le. dzero) then
      vr = vr2
      itopbotexit = -2
    end if
 
    ! Check if velocity is in the same direction throughout cell (i.e. no flow divide).
    ! Check if product v1*v2 > 0 then the velocity is in the same direction throughout
    ! the cell (i.e. no flow divide). If so, set vr to reflect appropriate direction.
    v1v2 = v1 * v2
    if (v1v2 .gt. dzero) then
      if (v .gt. dzero) then
        vr = vr2
        itopbotexit = -2
      end if
      if (v .lt. dzero) then
        vr = vr1
        itopbotexit = -1
      end if
    end if
 
    ! Compute travel time to exit face. Return with status = 0
    dt = log(abs(vr)) / dvdx
    status = ok_exit

Here is the caller graph for this function:

calculate_xyz_position()

subroutine methodsubcellternarymodule::calculate_xyz_position ( real(dp)  dt, real(dp)  rxx, real(dp)  rxy, real(dp)  ryx, real(dp)  ryy, real(dp)  sxx, real(dp)  sxy, real(dp)  syy, integer(i4b)  izstatus, real(dp)  x0, real(dp)  y0, real(dp)  az, real(dp)  vzi, real(dp)  vzbot, real(dp)  ztop, real(dp)  zbot, real(dp)  zi, real(dp)  x, real(dp)  y, real(dp)  z, integer(i4b), optional  exitface  )

private

Definition at line 511 of file MethodSubcellTernary.f90.

    ! dummy
    real(DP) :: dt
    real(DP) :: rxx
    real(DP) :: rxy
    real(DP) :: ryx
    real(DP) :: ryy
    real(DP) :: sxx
    real(DP) :: sxy
    real(DP) :: syy
    integer(I4B) :: izstatus
    real(DP) :: x0
    real(DP) :: y0
    real(DP) :: az
    real(DP) :: vzi
    real(DP) :: vzbot
    real(DP) :: ztop
    real(DP) :: zbot
    real(DP) :: zi
    real(DP) :: x
    real(DP) :: y
    real(DP) :: z
    integer(I4B), optional :: exitface
    ! local
    integer(I4B) :: lexitface
    real(DP) :: rot(2, 2), res(2), loc(2)
    real(DP) :: alp
    real(DP) :: bet
 
    ! process optional exit face argument
    if (present(exitface)) then
      lexitface = exitface
    else
      lexitface = 0
    end if
 
    ! calculate alpha and beta
    call step_analytical(dt, alp, bet)
 
    ! if exit face is known, set alpha or beta coordinate
    ! corresponding to the exit face exactly.
    if (lexitface .eq. 1) then
      bet = dzero
    else if (lexitface .eq. 2) then
      alp = done - bet
    else if (lexitface .eq. 3) then
      alp = dzero
    end if
 
    ! if exit face is top or bottom, set z coordinate exactly.
    if (lexitface .eq. 4) then
      z = zbot
    else if (lexitface .eq. 5) then
      z = ztop
    else
      ! otherwise calculate z.
      if (izstatus .eq. 2) then
        ! vz uniformly zero
        z = zi
      else if (izstatus .eq. 1) then
        ! vz uniform, nonzero
        z = zi + vzi * dt
      else
        ! vz nonuniform
        z = zbot + (vzi * dexp(az * dt) - vzbot) / az
      end if
    end if
 
    ! transform (alp, beta) to (x, y)
    loc = (/alp, bet/)
    loc = skew(loc, (/sxx, sxy, syy/), invert=.true.)
    rot = reshape((/rxx, rxy, ryx, ryy/), shape(rot))
    res = matmul(rot, loc) ! rotate vector
    x = res(1) + x0
    y = res(2) + y0
 

Here is the call graph for this function:

Here is the caller graph for this function:

create_method_subcell_ternary()

subroutine, public methodsubcellternarymodule::create_method_subcell_ternary

(

type(methodsubcellternarytype), pointer

method

)

Definition at line 57 of file MethodSubcellTernary.f90.

    ! dummy
    type(MethodSubcellTernaryType), pointer :: method
    ! local
    type(SubcellTriType), pointer :: subcell
 
    allocate (method)
    call create_subcell_tri(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 methodsubcellternarymodule::deallocate ( class(methodsubcellternarytype), intent(inout)  this )

private

Definition at line 71 of file MethodSubcellTernary.f90.

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

find_exits()

subroutine methodsubcellternarymodule::find_exits ( type(particletype), intent(inout), pointer  particle, class(domaintype), intent(in)  domain, type(barycentricexitsolutiontype), dimension(2), intent(out)  exit_solutions, real(dp), intent(out)  z0  )

private

Definition at line 227 of file MethodSubcellTernary.f90.

    type(ParticleType), pointer, intent(inout) :: particle
    class(DomainType), intent(in) :: domain
    type(BarycentricExitSolutionType), intent(out) :: exit_solutions(2)
    real(DP), intent(out) :: z0
    ! local
    integer(I4B) :: ntmax
    real(DP) :: tol
    real(DP) :: zirel
    real(DP) :: rxx
    real(DP) :: rxy
    real(DP) :: ryx
    real(DP) :: ryy
    real(DP) :: sxx
    real(DP) :: sxy
    real(DP) :: syy
    real(DP) :: alp0
    real(DP) :: bet0
    real(DP) :: alp1
    real(DP) :: bet1
    real(DP) :: alp2
    real(DP) :: bet2
    real(DP) :: alpi
    real(DP) :: beti
    real(DP) :: gami
    integer(I4B) :: isolv
    integer(I4B) :: itrifaceenter
 
    ntmax = 10000
    tol = particle%extol
 
    ! Set lateral solution method
    if (particle%iexmeth == 0) then
      isolv = 1 ! default to Brent's
    else
      isolv = particle%iexmeth
    end if
 
    select type (subcell => domain)
    type is (subcelltritype)
      ! Transform coordinates to the "canonical" configuration:
      ! barycentric in two dimensions with alpha, beta & gamma
      ! such that at f2 alpha = 0, f0 beta = 0, f1 gamma = 0.
      !
      !     v2
      !     |\
      !   f2| \f1
      !     |__\
      !   v0 f0 v1
      !
      call canonical(subcell%x0, subcell%y0, &
                     subcell%x1, subcell%y1, &
                     subcell%x2, subcell%y2, &
                     subcell%v0x, subcell%v0y, &
                     subcell%v1x, subcell%v1y, &
                     subcell%v2x, subcell%v2y, &
                     particle%x, particle%y, &
                     rxx, rxy, ryx, ryy, &
                     sxx, sxy, syy, &
                     alp0, bet0, alp1, bet1, alp2, bet2, alpi, beti)
 
      ! Clamp particle coordinates to the canonical triangular
      ! subcell and nudge it ever so slightly inside if needed.
      call clamp_bary(alpi, beti, gami, pad=dsame * dep3)
 
      ! Do calculations related to the analytical z solution.
      ! (TODO: just once for each cell? store at cell-level?)
      ! Clamp the relative z coordinate to the unit interval, and nudge
      ! the particle's z ever so slightly inside the subcell if needed,
      ! analogous to the lateral clamp_bary nudge above.
      zirel = (particle%z - subcell%zbot) / subcell%dz
      if (zirel > done) then
        zirel = done
      else if (zirel < dzero) then
        zirel = dzero
      end if
      z0 = subcell%zbot + zirel * subcell%dz
      exit_solutions(1) = find_vertical_exit(subcell%vzbot, subcell%vztop, &
                                             subcell%dz, zirel)
 
      ! Calculate a semi-analytical lateral exit solution
      itrifaceenter = particle%iboundary(level_subfeature) - 1
      if (itrifaceenter == -1) itrifaceenter = 999
      exit_solutions(2) = find_lateral_exit(isolv, tol, &
                                            itrifaceenter, &
                                            alp1, bet1, alp2, &
                                            bet2, alpi, beti)
      exit_solutions(2)%rxx = rxx
      exit_solutions(2)%rxy = rxy
      exit_solutions(2)%ryx = ryx
      exit_solutions(2)%ryy = ryy
      exit_solutions(2)%sxx = sxx
      exit_solutions(2)%sxy = sxy
      exit_solutions(2)%syy = syy
 
      ! Set vertical solution exit face
      if (exit_solutions(1)%itopbotexit /= 0) then
        if (exit_solutions(1)%itopbotexit == -1) then
          exit_solutions(1)%iboundary = 4
        else
          exit_solutions(1)%iboundary = 5
        end if
      end if
 
      ! Set lateral solution exit face
      if (exit_solutions(2)%itrifaceexit /= 0) &
        exit_solutions(2)%iboundary = exit_solutions(2)%itrifaceexit
    end select

Here is the call graph for this function:

Here is the caller graph for this function:

find_lateral_exit()

type(barycentricexitsolutiontype) function methodsubcellternarymodule::find_lateral_exit ( integer(i4b)  isolv, real(dp)  tol, integer(i4b)  itrifaceenter, real(dp)  alp1, real(dp)  bet1, real(dp)  alp2, real(dp)  bet2, real(dp)  alpi, real(dp)  beti  )

private

Definition at line 337 of file MethodSubcellTernary.f90.

    integer(I4B) :: isolv
    real(DP) :: tol
    integer(I4B) :: itrifaceenter
    real(DP) :: alp1
    real(DP) :: bet1
    real(DP) :: alp2
    real(DP) :: bet2
    real(DP) :: alpi
    real(DP) :: beti
    type(BarycentricExitSolutionType) :: solution
 
    solution = barycentricexitsolutiontype()
    call traverse_triangle(isolv, tol, &
                           solution%dt, solution%alpexit, solution%betexit, &
                           itrifaceenter, solution%itrifaceexit, &
                           alp1, bet1, alp2, bet2, alpi, beti)
    if (solution%itrifaceexit > 0) solution%status = ok_exit

Here is the call graph for this function:

Here is the caller graph for this function:

find_vertical_exit()

type(barycentricexitsolutiontype) function methodsubcellternarymodule::find_vertical_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 360 of file MethodSubcellTernary.f90.

    real(DP), intent(in) :: v1
    real(DP), intent(in) :: v2
    real(DP), intent(in) :: dx
    real(DP), intent(in) :: xL
    type(BarycentricExitSolutionType) :: solution
 
    solution = barycentricexitsolutiontype()
    call calculate_dt(v1, v2, dx, xl, solution%v, solution%dvdx, &
                      solution%dt, solution%status, solution%itopbotexit)

Here is the call graph for this function:

Here is the caller graph for this function:

pick_exit()

integer(i4b) function methodsubcellternarymodule::pick_exit

(

type(barycentricexitsolutiontype), dimension(2), intent(in)

exit_solutions

)

Definition at line 211 of file MethodSubcellTernary.f90.

    type(BarycentricExitSolutionType), intent(in) :: exit_solutions(2)
    integer(I4B) :: exit_soln
 
    if (exit_solutions(1)%itopbotexit == 0) then
      exit_soln = 2 ! lateral
    else if (exit_solutions(2)%itrifaceexit == 0 .or. &
             exit_solutions(1)%dt < exit_solutions(2)%dt) then
      exit_soln = 1 ! top/bottom
    else
      exit_soln = 2 ! lateral
    end if
 

Here is the caller graph for this function:

track_subcell()

subroutine methodsubcellternarymodule::track_subcell ( class(methodsubcellternarytype), intent(inout)  this, class(subcelltritype), intent(in)  subcell, type(particletype), intent(inout), pointer  particle, real(dp), intent(in)  tmax  )

private

Definition at line 89 of file MethodSubcellTernary.f90.

    use particlemodule, only: active, term_no_exits_sub
    use particleeventmodule, only: featexit, terminate, timestep, usertime
    ! dummy
    class(MethodSubcellTernaryType), intent(inout) :: this
    class(SubcellTriType), intent(in) :: subcell
    type(ParticleType), pointer, intent(inout) :: particle
    real(DP), intent(in) :: tmax
    ! local
    real(DP) :: dt, dtexit, texit
    real(DP) :: t0, t, x, y, z0, z
    integer(I4B) :: exit_face, exit_soln, event_code, i, isolv
    type(BarycentricExitSolutionType) :: exit_solutions(2)
    type(BarycentricExitSolutionType) :: exit_z, exit_lateral
 
    event_code = -1
    if (particle%iexmeth == 0) then
      isolv = 1 ! default to Brent's solution method
    else
      isolv = particle%iexmeth
    end if
    t0 = particle%ttrack
 
    ! Find exit solutions in lateral and vertical directions. z0 comes
    ! back clamped to [zbot, ztop] (mirroring the nudge clamp_bary applies
    ! to the lateral coordinates), so it can safely anchor the position
    ! update even if the raw particle%z lies outside the subcell.
    call find_exits(particle, subcell, exit_solutions, z0)
 
    exit_z = exit_solutions(1)
    exit_lateral = exit_solutions(2)
 
    ! If the subcell has no exit face, terminate the particle.
    ! TODO: consider ramifications
    if (exit_z%itopbotexit == 0 .and. &
        exit_lateral%itrifaceexit == 0) then
      call this%terminate(particle, status=term_no_exits_sub)
      return
    end if
 
    ! Determine exit solution, face, travel time, and time
    exit_soln = pick_exit(exit_solutions)
    exit_face = exit_solutions(exit_soln)%iboundary
    dtexit = exit_solutions(exit_soln)%dt
    if (dtexit < dzero) then
      call this%terminate(particle, status=term_no_exits_sub)
      return
    end if
    texit = t0 + dtexit
 
    ! Select user tracking times to solve. If this is the last time step
    ! in the simulation, times falling after the simulation end time are
    ! only included if the 'extend' option is on, otherwise only times in
    ! the time step are included.
    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 < t0) cycle
        if (t > texit .or. t > tmax) exit
        dt = t - t0
        call calculate_xyz_position(dt, &
                                    exit_lateral%rxx, exit_lateral%rxy, &
                                    exit_lateral%ryx, exit_lateral%ryy, &
                                    exit_lateral%sxx, exit_lateral%sxy, &
                                    exit_lateral%syy, exit_z%status, &
                                    subcell%x0, subcell%y0, &
                                    exit_z%dvdx, exit_z%v, &
                                    subcell%vzbot, subcell%ztop, subcell%zbot, &
                                    z0, x, y, z)
        particle%x = x
        particle%y = y
        particle%z = z
        particle%ttrack = t
        particle%istatus = active
        call this%usertime(particle)
      end do
    end if
 
    ! Compute exit time and face and update the particle's coordinates
    ! (local, unscaled) and other properties. The particle may at this
    ! point lie on a boundary of the subcell or may still be within it.
    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.
      t = tmax
      dt = t - t0
      exit_face = 0
      particle%istatus = active
      particle%advancing = .false.
      event_code = timestep
    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.
      t = texit
      dt = dtexit
      event_code = featexit
    end if
    call calculate_xyz_position(dt, &
                                exit_lateral%rxx, exit_lateral%rxy, &
                                exit_lateral%ryx, exit_lateral%ryy, &
                                exit_lateral%sxx, exit_lateral%sxy, &
                                exit_lateral%syy, exit_z%status, &
                                subcell%x0, subcell%y0, &
                                exit_z%dvdx, exit_z%v, &
                                subcell%vzbot, subcell%ztop, subcell%zbot, &
                                z0, x, y, z, exit_face)
    particle%x = x
    particle%y = y
    particle%z = z
    particle%ttrack = t
    particle%iboundary(level_subfeature) = exit_face
 
    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: