MODFLOW 6  version 6.7.0.dev3
USGS Modular Hydrologic Model
MethodSubcellPollock.f90
Go to the documentation of this file.
2  use kindmodule, only: dp, i4b, lgp
3  use errorutilmodule, only: pstop
8  use prtfmimodule, only: prtfmitype
9  use basedismodule, only: disbasetype
10  use cellmodule, only: celltype
11  use constantsmodule, only: dzero, done
12  use domainmodule, only: domaintype
13  use subcellmodule, only: subcelltype
14  use listmodule, only: listtype
19  implicit none
20  private
21  public :: methodsubcellpollocktype
23  public :: calculate_dt
24 
25  !> @brief Rectangular subcell tracking method
27  private
28  real(dp), allocatable, public :: qextl1(:), qextl2(:), qintl(:) !< external and internal subcell flows
29  type(linearexitsolutiontype), public :: exit_solutions(3) !< candidate exit solutions
30  contains
31  procedure, public :: find_exits
32  procedure, public :: pick_exit
33  procedure, public :: apply => apply_msp
34  procedure, public :: deallocate
35  procedure, private :: track_subcell
37 
38 contains
39 
40  !> @brief Create a new Pollock's subcell method
41  subroutine create_method_subcell_pollock(method)
42  ! dummy
43  type(methodsubcellpollocktype), pointer :: method
44  ! local
45  type(subcellrecttype), pointer :: subcell
46 
47  allocate (method)
48  call create_subcell_rect(subcell)
49  method%subcell => subcell
50  method%name => method%subcell%type
51  method%delegates = .false.
52  end subroutine create_method_subcell_pollock
53 
54  !> @brief Deallocate the Pollock's subcell method
55  subroutine deallocate (this)
56  class(methodsubcellpollocktype), intent(inout) :: this
57  deallocate (this%name)
58  end subroutine deallocate
59 
60  !> @brief Apply Pollock's method to a rectangular subcell
61  subroutine apply_msp(this, particle, tmax)
62  ! dummy
63  class(methodsubcellpollocktype), intent(inout) :: this
64  type(particletype), pointer, intent(inout) :: particle
65  real(DP), intent(in) :: tmax
66  ! local
67  real(DP) :: x_origin
68  real(DP) :: y_origin
69  real(DP) :: z_origin
70  real(DP) :: sinrot
71  real(DP) :: cosrot
72 
73  select type (subcell => this%subcell)
74  type is (subcellrecttype)
75  ! Transform particle position into local subcell coordinates,
76  ! track particle across subcell, convert back to model coords
77  ! (sinrot and cosrot should be 0 and 1, respectively, i.e. no
78  ! rotation, also no z translation; only x and y translations)
79  x_origin = subcell%xOrigin
80  y_origin = subcell%yOrigin
81  z_origin = subcell%zOrigin
82  sinrot = subcell%sinrot
83  cosrot = subcell%cosrot
84  call particle%transform(x_origin, y_origin)
85  call this%track_subcell(subcell, particle, tmax)
86  call particle%transform(x_origin, y_origin, invert=.true.)
87  end select
88  end subroutine apply_msp
89 
90  !> @brief Track a particle across a rectangular subcell using Pollock's method
91  !!
92  !! This subroutine consists partly of code written by
93  !! David W. Pollock of the USGS for MODPATH 7. PRT's
94  !! authors take responsibility for its application in
95  !! this context and for any modifications or errors.
96  !<
97  subroutine track_subcell(this, subcell, particle, tmax)
98  use tdismodule, only: endofsimulation
101  ! dummy
102  class(methodsubcellpollocktype), intent(inout) :: this
103  class(subcellrecttype), intent(in) :: subcell
104  type(particletype), pointer, intent(inout) :: particle
105  real(DP), intent(in) :: tmax
106  ! local
107  real(DP) :: dt, dtexit, texit
108  real(DP) :: t, x, y, z
109  real(DP) :: t0, x0, y0, z0
110  integer(I4B) :: i, exit_face, exit_soln
111  type(linearexitsolutiontype) :: exit_x, exit_y, exit_z
112 
113  t0 = particle%ttrack
114  x0 = particle%x / subcell%dx
115  y0 = particle%y / subcell%dy
116  z0 = particle%z / subcell%dz
117 
118  ! Find exit solution in each direction
119  call this%find_exits(particle, subcell)
120 
121  exit_x = this%exit_solutions(1)
122  exit_y = this%exit_solutions(2)
123  exit_z = this%exit_solutions(3)
124 
125  ! Subcell has no exit face, terminate the particle
126  ! TODO: consider ramifications
127  if (all([this%exit_solutions%status] == no_exit_no_outflow)) then
128  call this%terminate(particle, status=term_no_exits_sub)
129  return
130  end if
131 
132  ! If particle stationary and extended tracking is on, terminate here if it's
133  ! the last timestep. TODO: temporary solution, consider where to catch this?
134  ! Should we really have to special case this here? We do diverge from MP7 in
135  ! guaranteeing that every particle terminates at the end of the simulation..
136  ! ideally that would be handled at a higher scope but with extended tracking
137  ! tmax is not the end of the simulation, it's just a wildly high upper bound.
138  if (all([this%exit_solutions%status] == no_exit_stationary) .and. &
139  particle%extend .and. endofsimulation) then
140  call this%terminate(particle, status=term_timeout)
141  return
142  end if
143 
144  ! Pick exit solution, face, travel time, and time
145  exit_soln = this%pick_exit(particle)
146  if (exit_soln == 0) then
147  exit_face = 0
148  dtexit = 1.0d+30
149  else
150  exit_face = this%exit_solutions(exit_soln)%iboundary
151  dtexit = this%exit_solutions(exit_soln)%dt
152  end if
153  texit = particle%ttrack + dtexit
154 
155  ! Select user tracking times to solve. If this is the first time step
156  ! of the simulation, include all times before it begins; if it is the
157  ! last time step include all times after it ends only if the 'extend'
158  ! option is on, otherwise times in this period and time step only.
159  call this%tracktimes%advance()
160  if (this%tracktimes%any()) then
161  do i = this%tracktimes%selection(1), this%tracktimes%selection(2)
162  t = this%tracktimes%times(i)
163  if (t < particle%ttrack) cycle
164  if (t >= texit .or. t >= tmax) exit
165  dt = t - t0
166  x = new_x(exit_x%v, exit_x%dvdx, subcell%vx1, subcell%vx2, &
167  dt, x0, subcell%dx, exit_x%status == 1)
168  y = new_x(exit_y%v, exit_y%dvdx, subcell%vy1, subcell%vy2, &
169  dt, y0, subcell%dy, exit_y%status == 1)
170  z = new_x(exit_z%v, exit_z%dvdx, subcell%vz1, subcell%vz2, &
171  dt, z0, subcell%dz, exit_z%status == 1)
172  particle%x = x * subcell%dx
173  particle%y = y * subcell%dy
174  particle%z = z * subcell%dz
175  particle%ttrack = t
176  particle%istatus = active
177  call this%usertime(particle)
178  end do
179  end if
180 
181  if (texit .gt. tmax) then
182  ! The computed exit time is greater than the maximum time, so set
183  ! final time for particle trajectory equal to maximum time and
184  ! calculate particle location at that final time.
185  t = tmax
186  dt = t - t0
187  x = new_x(exit_x%v, exit_x%dvdx, subcell%vx1, subcell%vx2, &
188  dt, x0, subcell%dx, exit_x%status == 1)
189  y = new_x(exit_y%v, exit_y%dvdx, subcell%vy1, subcell%vy2, &
190  dt, y0, subcell%dy, exit_y%status == 1)
191  z = new_x(exit_z%v, exit_z%dvdx, subcell%vz1, subcell%vz2, &
192  dt, z0, subcell%dz, exit_z%status == 1)
193  exit_face = 0
194  particle%istatus = active
195  particle%advancing = .false.
196 
197  ! Set final particle location in local (unscaled) subcell coordinates,
198  ! final time for particle trajectory
199  particle%x = x * subcell%dx
200  particle%y = y * subcell%dy
201  particle%z = z * subcell%dz
202  particle%ttrack = t
203  particle%iboundary(level_subfeature) = exit_face
204 
205  ! Save particle track record
206  call this%timestep(particle)
207  else
208  ! The computed exit time is less than or equal to the maximum time,
209  ! so set final time for particle trajectory equal to exit time and
210  ! calculate exit location.
211  t = texit
212  dt = dtexit
213  if ((exit_face .eq. 1) .or. (exit_face .eq. 2)) then
214  x = dzero
215  y = new_x(exit_y%v, exit_y%dvdx, subcell%vy1, subcell%vy2, &
216  dt, y0, subcell%dy, exit_y%status == 1)
217  z = new_x(exit_z%v, exit_z%dvdx, subcell%vz1, subcell%vz2, &
218  dt, z0, subcell%dz, exit_z%status == 1)
219  if (exit_face .eq. 2) x = done
220  else if ((exit_face .eq. 3) .or. (exit_face .eq. 4)) then
221  x = new_x(exit_x%v, exit_x%dvdx, subcell%vx1, subcell%vx2, dt, &
222  x0, subcell%dx, exit_x%status == 1)
223  y = dzero
224  z = new_x(exit_z%v, exit_z%dvdx, subcell%vz1, subcell%vz2, dt, &
225  z0, subcell%dz, exit_z%status == 1)
226  if (exit_face .eq. 4) y = done
227  else if ((exit_face .eq. 5) .or. (exit_face .eq. 6)) then
228  x = new_x(exit_x%v, exit_x%dvdx, subcell%vx1, subcell%vx2, &
229  dt, x0, subcell%dx, exit_x%status == 1)
230  y = new_x(exit_y%v, exit_y%dvdx, subcell%vy1, subcell%vy2, &
231  dt, y0, subcell%dy, exit_y%status == 1)
232  z = dzero
233  if (exit_face .eq. 6) z = done
234  else
235  print *, "programmer error, invalid exit face", exit_face
236  call pstop(1)
237  end if
238 
239  ! Set final particle location in local (unscaled) subcell coordinates,
240  ! final time for particle trajectory, and exit face
241  particle%x = x * subcell%dx
242  particle%y = y * subcell%dy
243  particle%z = z * subcell%dz
244  particle%ttrack = t
245  particle%iboundary(level_subfeature) = exit_face
246 
247  ! Save particle track record
248  call this%subcellexit(particle)
249  end if
250 
251  end subroutine track_subcell
252 
253  !> @brief Pick the exit solution with the shortest travel time
254  function pick_exit(this, particle) result(exit_soln)
255  class(methodsubcellpollocktype), intent(inout) :: this
256  type(particletype), pointer, intent(inout) :: particle
257  integer(I4B) :: exit_soln
258  ! local
259  real(dp) :: dtmin
260 
261  exit_soln = 0
262  dtmin = 1.0d+30
263 
264  if (this%exit_solutions(1)%status < 2) then
265  exit_soln = 1 ! x
266  dtmin = this%exit_solutions(1)%dt
267  end if
268  if (this%exit_solutions(2)%status < 2 .and. &
269  this%exit_solutions(2)%dt < dtmin) then
270  exit_soln = 2 ! y
271  dtmin = this%exit_solutions(2)%dt
272  end if
273  if (this%exit_solutions(3)%status < 2 .and. &
274  this%exit_solutions(3)%dt < dtmin) then
275  exit_soln = 3 ! z
276  end if
277 
278  end function pick_exit
279 
280  !> @brief Compute candidate exit solutions
281  subroutine find_exits(this, particle, domain)
282  class(methodsubcellpollocktype), intent(inout) :: this
283  type(particletype), pointer, intent(inout) :: particle
284  class(domaintype), intent(in) :: domain
285  ! local
286  real(DP) :: x0, y0, z0
287 
288  select type (domain)
289  type is (subcellrecttype)
290  ! Initial particle location in scaled subcell coordinates
291  x0 = particle%x / domain%dx
292  y0 = particle%y / domain%dy
293  z0 = particle%z / domain%dz
294 
295  ! Calculate exit solutions for each coordinate direction
296  this%exit_solutions = [ &
297  find_exit(domain%vx1, domain%vx2, domain%dx, x0), &
298  find_exit(domain%vy1, domain%vy2, domain%dy, y0), &
299  find_exit(domain%vz1, domain%vz2, domain%dz, z0) &
300  ]
301 
302  ! Set exit faces
303  if (this%exit_solutions(1)%v < dzero) then
304  this%exit_solutions(1)%iboundary = 1
305  else if (this%exit_solutions(1)%v > dzero) then
306  this%exit_solutions(1)%iboundary = 2
307  end if
308  if (this%exit_solutions(2)%v < dzero) then
309  this%exit_solutions(2)%iboundary = 3
310  else if (this%exit_solutions(2)%v > dzero) then
311  this%exit_solutions(2)%iboundary = 4
312  end if
313  if (this%exit_solutions(3)%v < dzero) then
314  this%exit_solutions(3)%iboundary = 5
315  else if (this%exit_solutions(3)%v > dzero) then
316  this%exit_solutions(3)%iboundary = 6
317  end if
318  end select
319  end subroutine find_exits
320 
321  !> @brief Find an exit solution for one dimension
322  function find_exit(v1, v2, dx, xL) result(solution)
323  ! dummy
324  real(dp), intent(in) :: v1
325  real(dp), intent(in) :: v2
326  real(dp), intent(in) :: dx
327  real(dp), intent(in) :: xl
328  type(linearexitsolutiontype) :: solution
329 
330  solution = linearexitsolutiontype()
331  solution%status = calculate_dt(v1, v2, dx, xl, &
332  solution%v, solution%dvdx, solution%dt)
333  end function find_exit
334 
335  !> @brief Calculate particle travel time to exit and exit status.
336  !!
337  !! This subroutine consists partly of code written by and/or adapted from
338  !! David W. Pollock of the USGS for MODPATH 7. The authors of the present
339  !! code are responsible for its appropriate application in this context
340  !! and for any modifications or errors.
341  !<
342  function calculate_dt(v1, v2, dx, xL, v, dvdx, dt) result(status)
343  ! dummy
344  real(dp) :: v1
345  real(dp) :: v2
346  real(dp) :: dx
347  real(dp) :: xl
348  real(dp) :: v
349  real(dp) :: dvdx
350  real(dp) :: dt
351  ! result
352  integer(I4B) :: status
353  ! local
354  real(dp) :: v2a
355  real(dp) :: v1a
356  real(dp) :: dv
357  real(dp) :: dva
358  real(dp) :: vv
359  real(dp) :: vvv
360  real(dp) :: zro
361  real(dp) :: zrom
362  real(dp) :: x
363  real(dp) :: tol
364  real(dp) :: vr1
365  real(dp) :: vr2
366  real(dp) :: vr
367  real(dp) :: v1v2
368  logical(LGP) :: nooutflow
369 
370  ! Initialize variables.
371  status = -1
372  dt = 1.0d+20
373  v2a = v2
374  if (v2a .lt. dzero) v2a = -v2a
375  v1a = v1
376  if (v1a .lt. dzero) v1a = -v1a
377  dv = v2 - v1
378  dva = dv
379  if (dva .lt. dzero) dva = -dva
380 
381  ! Check for a uniform zero velocity in this direction.
382  ! If so, set status = 2 and return (dt = 1.0d+20).
383  tol = 1.0d-15
384  if ((v2a .lt. tol) .and. (v1a .lt. tol)) then
385  v = dzero
386  dvdx = dzero
387  status = no_exit_stationary
388  return
389  end if
390 
391  ! Check for uniform non-zero velocity in this direction.
392  ! If so, set compute dt using the constant velocity,
393  ! set status = 1 and return.
394  vv = v1a
395  if (v2a .gt. vv) vv = v2a
396  vvv = dva / vv
397  if (vvv .lt. 1.0d-4) then
398  zro = tol
399  zrom = -zro
400  v = v1
401  x = xl * dx
402  if (v1 .gt. zro) dt = (dx - x) / v1
403  if (v1 .lt. zrom) dt = -x / v1
404  dvdx = dzero
405  status = ok_exit_constant
406  return
407  end if
408 
409  ! Velocity has a linear variation.
410  ! Compute velocity corresponding to particle position.
411  dvdx = dv / dx
412  v = (done - xl) * v1 + xl * v2
413 
414  ! If flow is into the cell from both sides there is no outflow.
415  ! In that case, set status = 3 and return.
416  nooutflow = .true.
417  if (v1 .lt. dzero) nooutflow = .false.
418  if (v2 .gt. dzero) nooutflow = .false.
419  if (nooutflow) then
420  status = no_exit_no_outflow
421  return
422  end if
423 
424  ! If there is a divide in the cell for this flow direction, check to
425  ! see if the particle is located exactly on the divide. If it is, move
426  ! it very slightly to get it off the divide. This avoids possible
427  ! numerical problems related to stagnation points.
428  if ((v1 .le. dzero) .and. (v2 .ge. dzero)) then
429  if (abs(v) .le. dzero) then
430  v = 1.0d-20
431  if (v2 .le. dzero) v = -v
432  end if
433  end if
434 
435  ! If there is a flow divide, this check finds out what side of the
436  ! divide the particle is on and sets the value of vr appropriately
437  ! to reflect that location.
438  vr1 = v1 / v
439  vr2 = v2 / v
440  vr = vr1
441  if (vr .le. dzero) then
442  vr = vr2
443  end if
444 
445  ! If the product v1*v2 > 0, the velocity is in the same direction
446  ! throughout the cell (i.e. no flow divide). If so, set the value
447  ! of vr to reflect the appropriate direction.
448  v1v2 = v1 * v2
449  if (v1v2 .gt. dzero) then
450  if (v .gt. dzero) vr = vr2
451  if (v .lt. dzero) vr = vr1
452  end if
453 
454  ! Check if vr is (very close to) zero.
455  ! If so, set status = 2 and return (dt = 1.0d+20).
456  if (dabs(vr) .lt. 1.0d-10) then
457  v = dzero
458  dvdx = dzero
459  status = no_exit_stationary
460  return
461  end if
462 
463  ! Compute travel time to exit face. Return with status = 0.
464  dt = log(vr) / dvdx
465  status = ok_exit
466 
467  end function calculate_dt
468 
469  !> @brief Update a cell-local coordinate based on a time increment.
470  !!
471  !! This subroutine consists partly or entirely of code written by
472  !! David W. Pollock of the USGS for MODPATH 7. The authors of the present
473  !! code are responsible for its appropriate application in this context
474  !! and for any modifications or errors.
475  !<
476  pure function new_x(v, dvdx, v1, v2, dt, x, dx, velocity_profile) result(newx)
477  ! dummy
478  real(dp), intent(in) :: v
479  real(dp), intent(in) :: dvdx
480  real(dp), intent(in) :: v1
481  real(dp), intent(in) :: v2
482  real(dp), intent(in) :: dt
483  real(dp), intent(in) :: x
484  real(dp), intent(in) :: dx
485  logical(LGP), intent(in), optional :: velocity_profile
486  ! result
487  real(dp) :: newx
488  logical(LGP) :: lprofile
489 
490  ! process optional arguments
491  if (present(velocity_profile)) then
492  lprofile = velocity_profile
493  else
494  lprofile = .false.
495  end if
496 
497  ! recompute coordinate
498  newx = x
499  if (lprofile) then
500  newx = newx + (v1 * dt / dx)
501  else if (v .ne. dzero) then
502  newx = newx + (v * (exp(dvdx * dt) - done) / dvdx / dx)
503  end if
504 
505  ! clamp to [0, 1]
506  if (newx .lt. dzero) newx = dzero
507  if (newx .gt. done) newx = done
508 
509  end function new_x
510 
This module contains simulation constants.
Definition: Constants.f90:9
real(dp), parameter dzero
real constant zero
Definition: Constants.f90:65
real(dp), parameter done
real constant 1
Definition: Constants.f90:76
subroutine pstop(status, message)
Stop the program, optionally specifying an error status code.
Definition: ErrorUtil.f90:24
@ ok_exit_constant
exit found, constant velocity
@ ok_exit
exit found using velocity interpolation
@ no_exit_stationary
no exit, zero velocity
@ no_exit_no_outflow
no exit, no outflow
This module defines variable data types.
Definition: kind.f90:8
Particle tracking strategies.
Definition: Method.f90:2
@, public level_subfeature
Definition: Method.f90:42
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.
integer(i4b) function pick_exit(this, particle)
Pick the exit solution with the shortest travel time.
integer(i4b) function, public calculate_dt(v1, v2, dx, xL, v, dvdx, dt)
Calculate particle travel time to exit and exit status.
type(linearexitsolutiontype) function find_exit(v1, v2, dx, xL)
Find an exit solution for one dimension.
subroutine track_subcell(this, subcell, particle, tmax)
Track a particle across a rectangular subcell using Pollock's method.
subroutine, public create_method_subcell_pollock(method)
Create a new Pollock's subcell method.
subroutine find_exits(this, particle, domain)
Compute candidate exit solutions.
subroutine apply_msp(this, particle, tmax)
Apply Pollock's method to a rectangular subcell.
@, public featexit
particle exited a grid feature
@, public timestep
time step ended
@ term_timeout
terminated at stop time or end of simulation
Definition: Particle.f90:37
@ term_no_exits_sub
terminated in a subcell with no exit face
Definition: Particle.f90:36
subroutine, public create_subcell_rect(subcell)
Create a new rectangular subcell.
Definition: SubcellRect.f90:28
logical(lgp), pointer, public endofsimulation
flag indicating end of simulation
Definition: tdis.f90:28
Base type for grid cells of a concrete type. Contains a cell-definition which is information shared b...
Definition: Cell.f90:12
A tracking domain.
Definition: Domain.f90:8
Base type for exit solutions.
Linear velocity interpolation exit solution.
A generic heterogeneous doubly-linked list.
Definition: List.f90:14
Abstract base type for subcell tracking methods.
Particle tracked by the PRT model.
Definition: Particle.f90:56
A subcell of a cell.
Definition: Subcell.f90:10