MODFLOW 6  version 6.6.0.dev0
USGS Modular Hydrologic Model
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 50 of file MethodSubcellPollock.f90.

51  ! -- dummy
52  class(MethodSubcellPollockType), intent(inout) :: this
53  type(ParticleType), pointer, intent(inout) :: particle
54  real(DP), intent(in) :: tmax
55  ! -- local
56  real(DP) :: xOrigin
57  real(DP) :: yOrigin
58  real(DP) :: zOrigin
59  real(DP) :: sinrot
60  real(DP) :: cosrot
61 
62  select type (subcell => this%subcell)
63  type is (subcellrecttype)
64  ! -- Transform particle position into local subcell coordinates,
65  ! track particle across subcell, convert back to model coords
66  ! (sinrot and cosrot should be 0 and 1, respectively, i.e. no
67  ! rotation, also no z translation; only x and y translations)
68  xorigin = subcell%xOrigin
69  yorigin = subcell%yOrigin
70  zorigin = subcell%zOrigin
71  sinrot = subcell%sinrot
72  cosrot = subcell%cosrot
73  call particle%transform(xorigin, yorigin)
74  call this%track_subcell(subcell, particle, tmax)
75  call particle%transform(xorigin, yorigin, invert=.true.)
76  call particle%reset_transform()
77  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 277 of file MethodSubcellPollock.f90.

278  ! dummy
279  real(DP) :: v1
280  real(DP) :: v2
281  real(DP) :: dx
282  real(DP) :: xL
283  real(DP) :: v
284  real(DP) :: dvdx
285  real(DP) :: dt
286  ! result
287  integer(I4B) :: status
288  ! local
289  real(DP) :: v2a
290  real(DP) :: v1a
291  real(DP) :: dv
292  real(DP) :: dva
293  real(DP) :: vv
294  real(DP) :: vvv
295  real(DP) :: zro
296  real(DP) :: zrom
297  real(DP) :: x
298  real(DP) :: tol
299  real(DP) :: vr1
300  real(DP) :: vr2
301  real(DP) :: vr
302  real(DP) :: v1v2
303  logical(LGP) :: noOutflow
304 
305  ! -- Initialize variables.
306  status = -1
307  dt = 1.0d+20
308  v2a = v2
309  if (v2a .lt. dzero) v2a = -v2a
310  v1a = v1
311  if (v1a .lt. dzero) v1a = -v1a
312  dv = v2 - v1
313  dva = dv
314  if (dva .lt. dzero) dva = -dva
315 
316  ! -- Check for a uniform zero velocity in this direction.
317  ! -- If so, set status = 2 and return (dt = 1.0d+20).
318  tol = 1.0d-15
319  if ((v2a .lt. tol) .and. (v1a .lt. tol)) then
320  v = dzero
321  dvdx = dzero
322  status = 2
323  return
324  end if
325 
326  ! -- Check for uniform non-zero velocity in this direction.
327  ! -- If so, set compute dt using the constant velocity,
328  ! -- set status = 1 and return.
329  vv = v1a
330  if (v2a .gt. vv) vv = v2a
331  vvv = dva / vv
332  if (vvv .lt. 1.0d-4) then
333  zro = tol
334  zrom = -zro
335  v = v1
336  x = xl * dx
337  if (v1 .gt. zro) dt = (dx - x) / v1
338  if (v1 .lt. zrom) dt = -x / v1
339  dvdx = dzero
340  status = 1
341  return
342  end if
343 
344  ! -- Velocity has a linear variation.
345  ! -- Compute velocity corresponding to particle position.
346  dvdx = dv / dx
347  v = (done - xl) * v1 + xl * v2
348 
349  ! -- If flow is into the cell from both sides there is no outflow.
350  ! -- In that case, set status = 3 and return.
351  nooutflow = .true.
352  if (v1 .lt. dzero) nooutflow = .false.
353  if (v2 .gt. dzero) nooutflow = .false.
354  if (nooutflow) then
355  status = 3
356  return
357  end if
358 
359  ! -- If there is a divide in the cell for this flow direction, check to
360  ! -- see if the particle is located exactly on the divide. If it is, move
361  ! -- it very slightly to get it off the divide. This avoids possible
362  ! -- numerical problems related to stagnation points.
363  if ((v1 .le. dzero) .and. (v2 .ge. dzero)) then
364  if (abs(v) .le. dzero) then
365  v = 1.0d-20
366  if (v2 .le. dzero) v = -v
367  end if
368  end if
369 
370  ! -- If there is a flow divide, this check finds out what side of the
371  ! -- divide the particle is on and sets the value of vr appropriately
372  ! -- to reflect that location.
373  vr1 = v1 / v
374  vr2 = v2 / v
375  vr = vr1
376  if (vr .le. dzero) then
377  vr = vr2
378  end if
379 
380  ! -- If the product v1*v2 > 0, the velocity is in the same direction
381  ! -- throughout the cell (i.e. no flow divide). If so, set the value
382  ! -- of vr to reflect the appropriate direction.
383  v1v2 = v1 * v2
384  if (v1v2 .gt. dzero) then
385  if (v .gt. dzero) vr = vr2
386  if (v .lt. dzero) vr = vr1
387  end if
388 
389  ! -- Check if vr is (very close to) zero.
390  ! -- If so, set status = 2 and return (dt = 1.0d+20).
391  if (dabs(vr) .lt. 1.0d-10) then
392  v = dzero
393  dvdx = dzero
394  status = 2
395  return
396  end if
397 
398  ! -- Compute travel time to exit face. Return with status = 0.
399  dt = log(vr) / dvdx
400  status = 0
401 
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 30 of file MethodSubcellPollock.f90.

31  ! -- dummy
32  type(MethodSubcellPollockType), pointer :: method
33  ! -- local
34  type(SubcellRectType), pointer :: subcell
35 
36  allocate (method)
37  call create_subcell_rect(subcell)
38  method%subcell => subcell
39  method%name => method%subcell%type
40  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 44 of file MethodSubcellPollock.f90.

45  class(MethodSubcellPollockType), intent(inout) :: this
46  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 411 of file MethodSubcellPollock.f90.

412  ! dummy
413  real(DP), intent(in) :: v
414  real(DP), intent(in) :: dvdx
415  real(DP), intent(in) :: v1
416  real(DP), intent(in) :: v2
417  real(DP), intent(in) :: dt
418  real(DP), intent(in) :: x
419  real(DP), intent(in) :: dx
420  logical(LGP), intent(in), optional :: velocity_profile
421  ! result
422  real(DP) :: newx
423  logical(LGP) :: lprofile
424 
425  ! -- process optional arguments
426  if (present(velocity_profile)) then
427  lprofile = velocity_profile
428  else
429  lprofile = .false.
430  end if
431 
432  ! -- recompute coordinate
433  newx = x
434  if (lprofile) then
435  newx = newx + (v1 * dt / dx)
436  else if (v .ne. dzero) then
437  newx = newx + (v * (exp(dvdx * dt) - done) / dvdx / dx)
438  end if
439 
440  ! -- clamp to [0, 1]
441  if (newx .lt. dzero) newx = dzero
442  if (newx .gt. done) newx = done
443 
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.

88  ! dummy
89  class(MethodSubcellPollockType), intent(inout) :: this
90  class(SubcellRectType), intent(in) :: subcell
91  type(ParticleType), pointer, intent(inout) :: particle
92  real(DP), intent(in) :: tmax
93  ! local
94  real(DP) :: vx
95  real(DP) :: dvxdx
96  real(DP) :: vy
97  real(DP) :: dvydy
98  real(DP) :: vz
99  real(DP) :: dvzdz
100  real(DP) :: dtexitx
101  real(DP) :: dtexity
102  real(DP) :: dtexitz
103  real(DP) :: dtexit
104  real(DP) :: texit
105  real(DP) :: dt
106  real(DP) :: t
107  real(DP) :: t0
108  real(DP) :: x
109  real(DP) :: y
110  real(DP) :: z
111  integer(I4B) :: statusVX
112  integer(I4B) :: statusVY
113  integer(I4B) :: statusVZ
114  integer(I4B) :: i
115  real(DP) :: initialX
116  real(DP) :: initialY
117  real(DP) :: initialZ
118  integer(I4B) :: exitFace
119  integer(I4B) :: reason
120 
121  reason = -1
122 
123  ! Initial particle location in scaled subcell coordinates
124  initialx = particle%x / subcell%dx
125  initialy = particle%y / subcell%dy
126  initialz = particle%z / subcell%dz
127 
128  ! Compute time of travel to each possible exit face
129  statusvx = calculate_dt(subcell%vx1, subcell%vx2, subcell%dx, &
130  initialx, vx, dvxdx, dtexitx)
131  statusvy = calculate_dt(subcell%vy1, subcell%vy2, subcell%dy, &
132  initialy, vy, dvydy, dtexity)
133  statusvz = calculate_dt(subcell%vz1, subcell%vz2, subcell%dz, &
134  initialz, vz, dvzdz, dtexitz)
135 
136  ! Subcell has no exit face, terminate the particle
137  ! todo: after initial release, consider ramifications
138  if ((statusvx .eq. 3) .and. (statusvy .eq. 3) .and. (statusvz .eq. 3)) then
139  particle%istatus = 9
140  particle%advancing = .false.
141  call this%save(particle, reason=3)
142  return
143  end if
144 
145  ! Determine (earliest) exit face and corresponding travel time to exit
146  exitface = 0
147  dtexit = 1.0d+30
148  if ((statusvx .lt. 2) .or. (statusvy .lt. 2) .or. (statusvz .lt. 2)) then
149  ! Consider x-oriented faces
150  dtexit = dtexitx
151  if (vx .lt. dzero) then
152  exitface = 1
153  else if (vx .gt. 0) then
154  exitface = 2
155  end if
156  ! Consider y-oriented faces
157  if (dtexity .lt. dtexit) then
158  dtexit = dtexity
159  if (vy .lt. dzero) then
160  exitface = 3
161  else if (vy .gt. dzero) then
162  exitface = 4
163  end if
164  end if
165  ! Consider z-oriented faces
166  if (dtexitz .lt. dtexit) then
167  dtexit = dtexitz
168  if (vz .lt. dzero) then
169  exitface = 5
170  else if (vz .gt. dzero) then
171  exitface = 6
172  end if
173  end if
174  else
175  end if
176 
177  ! Compute exit time
178  texit = particle%ttrack + dtexit
179  t0 = particle%ttrack
180 
181  ! Select user tracking times to solve. If this is the first time step
182  ! of the simulation, include all times before it begins; if it is the
183  ! last time step include all times after it ends only if the 'extend'
184  ! option is on, otherwise times in this period and time step only.
185  call this%tracktimes%advance()
186  if (this%tracktimes%any()) then
187  do i = this%tracktimes%selection(1), this%tracktimes%selection(2)
188  t = this%tracktimes%times(i)
189  if (t < particle%ttrack) cycle
190  if (t >= texit .or. t >= tmax) exit
191  dt = t - t0
192  x = new_x(vx, dvxdx, subcell%vx1, subcell%vx2, &
193  dt, initialx, subcell%dx, statusvx == 1)
194  y = new_x(vy, dvydy, subcell%vy1, subcell%vy2, &
195  dt, initialy, subcell%dy, statusvy == 1)
196  z = new_x(vz, dvzdz, subcell%vz1, subcell%vz2, &
197  dt, initialz, subcell%dz, statusvz == 1)
198  particle%x = x * subcell%dx
199  particle%y = y * subcell%dy
200  particle%z = z * subcell%dz
201  particle%ttrack = t
202  particle%istatus = 1
203  call this%save(particle, reason=5)
204  end do
205  end if
206 
207  if (texit .gt. tmax) then
208  ! The computed exit time is greater than the maximum time, so set
209  ! final time for particle trajectory equal to maximum time and
210  ! calculate particle location at that final time.
211  t = tmax
212  dt = t - t0
213  x = new_x(vx, dvxdx, subcell%vx1, subcell%vx2, &
214  dt, initialx, subcell%dx, statusvx == 1)
215  y = new_x(vy, dvydy, subcell%vy1, subcell%vy2, &
216  dt, initialy, subcell%dy, statusvy == 1)
217  z = new_x(vz, dvzdz, subcell%vz1, subcell%vz2, &
218  dt, initialz, subcell%dz, statusvz == 1)
219  exitface = 0
220  particle%istatus = 1
221  particle%advancing = .false.
222  reason = 2 ! timestep end
223  else
224  ! The computed exit time is less than or equal to the maximum time,
225  ! so set final time for particle trajectory equal to exit time and
226  ! calculate exit location.
227  t = texit
228  dt = dtexit
229  if ((exitface .eq. 1) .or. (exitface .eq. 2)) then
230  x = dzero
231  y = new_x(vy, dvydy, subcell%vy1, subcell%vy2, &
232  dt, initialy, subcell%dy, statusvy == 1)
233  z = new_x(vz, dvzdz, subcell%vz1, subcell%vz2, &
234  dt, initialz, subcell%dz, statusvz == 1)
235  if (exitface .eq. 2) x = done
236  else if ((exitface .eq. 3) .or. (exitface .eq. 4)) then
237  x = new_x(vx, dvxdx, subcell%vx1, subcell%vx2, dt, &
238  initialx, subcell%dx, statusvx == 1)
239  y = dzero
240  z = new_x(vz, dvzdz, subcell%vz1, subcell%vz2, dt, &
241  initialz, subcell%dz, statusvz == 1)
242  if (exitface .eq. 4) y = done
243  else if ((exitface .eq. 5) .or. (exitface .eq. 6)) then
244  x = new_x(vx, dvxdx, subcell%vx1, subcell%vx2, &
245  dt, initialx, subcell%dx, statusvx == 1)
246  y = new_x(vy, dvydy, subcell%vy1, subcell%vy2, &
247  dt, initialy, subcell%dy, statusvy == 1)
248  z = dzero
249  if (exitface .eq. 6) z = done
250  else
251  print *, "programmer error, invalid exit face", exitface
252  call pstop(1)
253  end if
254  reason = 1 ! cell transition
255  end if
256 
257  ! Set final particle location in local (unscaled) subcell coordinates,
258  ! final time for particle trajectory, and exit face
259  particle%x = x * subcell%dx
260  particle%y = y * subcell%dy
261  particle%z = z * subcell%dz
262  particle%ttrack = t
263  particle%iboundary(3) = exitface
264 
265  ! Save particle track record
266  if (reason >= 0) call this%save(particle, reason=reason)
267 
Here is the call graph for this function: