MODFLOW 6  version 6.7.0.dev1
USGS Modular Hydrologic Model
MethodSubcellTernary.f90
Go to the documentation of this file.
2  use kindmodule, only: dp, i4b, lgp
4  use errorutilmodule, only: pstop
5  use geomutilmodule, only: clamp_bary, skew
6  use methodmodule, only: methodtype
7  use cellmodule, only: celltype
8  use subcellmodule, only: subcelltype
10  use particlemodule, only: particletype
12  use prtfmimodule, only: prtfmitype
13  use basedismodule, only: disbasetype
14  use mathutilmodule, only: is_close
15  implicit none
16 
17  private
18  public :: methodsubcellternarytype
20 
21  !> @brief Ternary triangular subcell tracking method.
23  integer(I4B), public, pointer :: zeromethod
24  contains
25  procedure, public :: apply => apply_mst
26  procedure, public :: deallocate
27  procedure, private :: track_subcell
29 
30 contains
31 
32  !> @brief Create a new ternary subcell tracking method.
33  subroutine create_method_subcell_ternary(method)
34  ! dummy
35  type(methodsubcellternarytype), pointer :: method
36  ! local
37  type(subcelltritype), pointer :: subcell
38 
39  allocate (method)
40  call create_subcell_tri(subcell)
41  method%subcell => subcell
42  method%name => method%subcell%type
43  method%delegates = .false.
44  end subroutine create_method_subcell_ternary
45 
46  !> @brief Deallocate the ternary subcell tracking method.
47  subroutine deallocate (this)
48  class(methodsubcellternarytype), intent(inout) :: this
49  deallocate (this%name)
50  end subroutine deallocate
51 
52  !> @brief Apply the ternary subcell tracking method.
53  subroutine apply_mst(this, particle, tmax)
54  class(methodsubcellternarytype), intent(inout) :: this
55  type(particletype), pointer, intent(inout) :: particle
56  real(DP), intent(in) :: tmax
57 
58  select type (subcell => this%subcell)
59  type is (subcelltritype)
60  call this%track_subcell(subcell, particle, tmax)
61  end select
62  end subroutine apply_mst
63 
64  !> @brief Track a particle across a triangular subcell.
65  subroutine track_subcell(this, subcell, particle, tmax)
67  ! dummy
68  class(methodsubcellternarytype), intent(inout) :: this
69  class(subcelltritype), intent(in) :: subcell
70  type(particletype), pointer, intent(inout) :: particle
71  real(DP), intent(in) :: tmax
72  ! local
73  integer(I4B) :: exitFace
74  real(DP) :: x0
75  real(DP) :: y0
76  real(DP) :: x1
77  real(DP) :: y1
78  real(DP) :: x2
79  real(DP) :: y2
80  real(DP) :: v0x
81  real(DP) :: v0y
82  real(DP) :: v1x
83  real(DP) :: v1y
84  real(DP) :: v2x
85  real(DP) :: v2y
86  real(DP) :: xi
87  real(DP) :: yi
88  real(DP) :: zi
89  real(DP) :: zirel
90  real(DP) :: ztop
91  real(DP) :: zbot
92  real(DP) :: dz
93  real(DP) :: rxx
94  real(DP) :: rxy
95  real(DP) :: ryx
96  real(DP) :: ryy
97  real(DP) :: sxx
98  real(DP) :: sxy
99  real(DP) :: syy
100  real(DP) :: alp0
101  real(DP) :: bet0
102  real(DP) :: alp1
103  real(DP) :: bet1
104  real(DP) :: alp2
105  real(DP) :: bet2
106  real(DP) :: alpi
107  real(DP) :: beti
108  real(DP) :: gami
109  real(DP) :: vzbot
110  real(DP) :: vztop
111  real(DP) :: vzi
112  real(DP) :: vziodz
113  real(DP) :: az
114  real(DP) :: dtexitz
115  real(DP) :: dt
116  real(DP) :: t
117  real(DP) :: t0
118  real(DP) :: dtexitxy
119  real(DP) :: texit
120  real(DP) :: x
121  real(DP) :: y
122  real(DP) :: z
123  integer(I4B) :: izstatus
124  integer(I4B) :: itopbotexit
125  integer(I4B) :: ntmax
126  integer(I4B) :: isolv
127  integer(I4B) :: itrifaceenter
128  integer(I4B) :: itrifaceexit
129  real(DP) :: tol
130  real(DP) :: dtexit
131  real(DP) :: alpexit
132  real(DP) :: betexit
133  integer(I4B) :: reason
134  integer(I4B) :: i
135 
136  ! Set solution method
137  if (particle%iexmeth == 0) then
138  isolv = 1 ! default to Brent's
139  else
140  isolv = particle%iexmeth
141  end if
142 
143  ntmax = 10000
144  tol = particle%extol
145  reason = -1
146 
147  ! Set some local variables for convenience.
148  xi = particle%x
149  yi = particle%y
150  zi = particle%z
151  x0 = subcell%x0
152  y0 = subcell%y0
153  x1 = subcell%x1
154  y1 = subcell%y1
155  x2 = subcell%x2
156  y2 = subcell%y2
157  v0x = subcell%v0x
158  v0y = subcell%v0y
159  v1x = subcell%v1x
160  v1y = subcell%v1y
161  v2x = subcell%v2x
162  v2y = subcell%v2y
163  zbot = subcell%zbot
164  ztop = subcell%ztop
165  dz = subcell%dz
166  vzbot = subcell%vzbot
167  vztop = subcell%vztop
168 
169  ! Transform coordinates to the "canonical" configuration:
170  ! barycentric in two dimensions with alpha, beta & gamma
171  ! such that at f2 alpha = 0, f0 beta = 0, f1 gamma = 0.
172  !
173  ! v2
174  ! |\
175  ! f2| \f1
176  ! |__\
177  ! v0 f0 v1
178  !
179  call canonical(x0, y0, x1, y1, x2, y2, &
180  v0x, v0y, v1x, v1y, v2x, v2y, &
181  xi, yi, &
182  rxx, rxy, ryx, ryy, &
183  sxx, sxy, syy, &
184  alp0, bet0, alp1, bet1, alp2, bet2, alpi, beti)
185 
186  ! Clamp particle coordinates to the canonical triangular
187  ! subcell and nudge it ever so slightly inside if needed.
188  call clamp_bary(alpi, beti, gami, pad=dsame * dep3)
189 
190  ! Do calculations related to the analytical z solution.
191  ! (TODO: just once for each cell? store at cell-level?)
192  ! Clamp the relative z coordinate to the unit interval.
193  zirel = (zi - zbot) / dz
194  if (zirel > done) then
195  zirel = done
196  else if (zirel < dzero) then
197  zirel = dzero
198  end if
199  call calculate_dt(vzbot, vztop, dz, zirel, vzi, &
200  az, dtexitz, izstatus, &
201  itopbotexit)
202  vziodz = vzi / dz
203 
204  ! If possible, track the particle across the subcell.
205  itrifaceenter = particle%iboundary(3) - 1
206  if (itrifaceenter == -1) itrifaceenter = 999
207  call traverse_triangle(isolv, tol, &
208  dtexitxy, alpexit, betexit, &
209  itrifaceenter, itrifaceexit, &
210  alp1, bet1, alp2, bet2, alpi, beti)
211 
212  ! If the subcell has no exit face, terminate the particle.
213  ! todo: after initial release, consider ramifications
214  if (itopbotexit == 0 .and. itrifaceexit == 0) then
215  particle%istatus = term_no_exits_sub
216  particle%advancing = .false.
217  call this%save(particle, reason=3)
218  return
219  end if
220 
221  ! -- Determine (earliest) exit face and corresponding travel time to exit
222  if (itopbotexit == 0) then
223  ! -- Exits through triangle face first
224  exitface = itrifaceexit
225  dtexit = dtexitxy
226  else if (itrifaceexit == 0 .or. dtexitz < dtexitxy) then
227  ! -- Exits through top/bottom first
228  exitface = 45
229  dtexit = dtexitz
230  else
231  ! -- Exits through triangle face first
232  exitface = itrifaceexit
233  dtexit = dtexitxy
234  end if
235  if (exitface == 45) then
236  if (itopbotexit == -1) then
237  exitface = 4
238  else
239  exitface = 5
240  end if
241  end if
242 
243  ! -- Make sure dt is positive
244  if (dtexit < dzero) then
245  particle%istatus = term_no_exits_sub
246  particle%advancing = .false.
247  call this%save(particle, reason=3)
248  return
249  end if
250 
251  texit = particle%ttrack + dtexit
252  t0 = particle%ttrack
253 
254  ! Select user tracking times to solve. If this is the last time step
255  ! in the simulation, times falling after the simulation end time are
256  ! only included if the 'extend' option is on, otherwise only times in
257  ! the time step are included.
258  call this%tracktimes%advance()
259  if (this%tracktimes%any()) then
260  do i = this%tracktimes%selection(1), this%tracktimes%selection(2)
261  t = this%tracktimes%times(i)
262  if (t < particle%ttrack) cycle
263  if (t >= texit .or. t >= tmax) exit
264  dt = t - t0
265  call calculate_xyz_position(dt, rxx, rxy, ryx, ryy, sxx, sxy, syy, &
266  izstatus, x0, y0, az, vzi, vzbot, &
267  ztop, zbot, zi, x, y, z)
268  particle%x = x
269  particle%y = y
270  particle%z = z
271  particle%ttrack = t
272  particle%istatus = active
273  call this%save(particle, reason=5)
274  end do
275  end if
276 
277  ! Compute exit time and face and update the particle's coordinates
278  ! (local, unscaled) and other properties. The particle may at this
279  ! point lie on a boundary of the subcell or may still be within it.
280  if (texit .gt. tmax) then
281  ! The computed exit time is greater than the maximum time, so set
282  ! final time for particle trajectory equal to maximum time.
283  t = tmax
284  dt = t - t0
285  exitface = 0
286  particle%istatus = active
287  particle%advancing = .false.
288  reason = 2 ! timestep end
289  else
290  ! The computed exit time is less than or equal to the maximum time,
291  ! so set final time for particle trajectory equal to exit time.
292  t = texit
293  dt = dtexit
294  reason = 1 ! (sub)cell transition
295  end if
296  call calculate_xyz_position(dt, rxx, rxy, ryx, ryy, sxx, sxy, syy, &
297  izstatus, x0, y0, az, vzi, vzbot, &
298  ztop, zbot, zi, x, y, z, exitface)
299  particle%x = x
300  particle%y = y
301  particle%z = z
302  particle%ttrack = t
303  particle%iboundary(3) = exitface
304 
305  call this%save(particle, reason=reason)
306  end subroutine track_subcell
307 
308  !> @brief Do calculations related to analytical z solution
309  !!
310  !! This subroutine consists partly or entirely of code written by
311  !! David W. Pollock of the USGS for MODPATH 7. The authors of the present
312  !! code are responsible for its appropriate application in this context
313  !! and for any modifications or errors.
314  !<
315  subroutine calculate_dt(v1, v2, dx, xL, v, dvdx, &
316  dt, status, itopbotexit)
317  real(DP) :: v1
318  real(DP) :: v2
319  real(DP) :: dx
320  real(DP) :: xL
321  real(DP) :: v
322  real(DP) :: dvdx
323  real(DP) :: dt
324  real(DP) :: v2a
325  real(DP) :: v1a
326  real(DP) :: dv
327  real(DP) :: dva
328  real(DP) :: vv
329  real(DP) :: vvv
330  real(DP) :: zro
331  real(DP) :: zrom
332  real(DP) :: x
333  real(DP) :: tol
334  real(DP) :: vr1
335  real(DP) :: vr2
336  real(DP) :: vr
337  real(DP) :: v1v2
338  integer(I4B) :: status
339  integer(I4B) :: itopbotexit
340  logical(LGP) :: noOutflow
341 
342  ! Initialize variables
343  status = -1
344  dt = 1.0d+20
345  v2a = v2
346  if (v2a .lt. dzero) v2a = -v2a
347  v1a = v1
348  if (v1a .lt. dzero) v1a = -v1a
349  dv = v2 - v1
350  dva = dv
351  if (dva .lt. dzero) dva = -dva
352 
353  ! Check for a uniform zero velocity in this direction.
354  ! If so, set status = 2 and return (dt = 1.0d+20).
355  tol = 1.0d-15
356  if ((v2a .lt. tol) .and. (v1a .lt. tol)) then
357  v = dzero
358  dvdx = dzero
359  status = 2
360  itopbotexit = 0
361  return
362  end if
363 
364  ! Check for uniform non-zero velocity in this direction.
365  ! If so, set compute dt using the constant velocity,
366  ! set status = 1 and return.
367  vv = v1a
368  if (v2a .gt. vv) vv = v2a
369  vvv = dva / vv
370  if (vvv .lt. 1.0d-4) then
371  zro = tol
372  zrom = -zro
373  v = v1
374  x = xl * dx
375  if (v1 .gt. zro) then
376  dt = (dx - x) / v1
377  itopbotexit = -2
378  end if
379  if (v1 .lt. zrom) then
380  dt = -x / v1
381  itopbotexit = -1
382  end if
383  dvdx = dzero
384  status = 1
385  return
386  end if
387 
388  ! Velocity has a linear variation.
389  ! Compute velocity corresponding to particle position
390  dvdx = dv / dx
391  v = (done - xl) * v1 + xl * v2
392 
393  ! If flow is into the cell from both sides there is no outflow.
394  ! In that case, set status = 3 and return
395  nooutflow = .true.
396  if (v1 .lt. dzero) nooutflow = .false.
397  if (v2 .gt. dzero) nooutflow = .false.
398  if (nooutflow) then
399  status = 3
400  itopbotexit = 0
401  return
402  end if
403 
404  ! If there is a divide in the cell for this flow direction, check to see if the
405  ! particle is located exactly on the divide. If it is, move it very slightly to
406  ! get it off the divide. This avoids possible numerical problems related to
407  ! stagnation points.
408  if ((v1 .le. dzero) .and. (v2 .ge. dzero)) then
409  if (abs(v) .le. dzero) then
410  v = 1.0d-20
411  if (v2 .le. dzero) v = -v
412  end if
413  end if
414 
415  ! If there is a flow divide, find out what side of the divide the particle
416  ! is on and set the value of vr appropriately to reflect that location.
417  vr1 = v1 / v
418  vr2 = v2 / v
419  vr = vr1
420  itopbotexit = -1
421  if (vr .le. dzero) then
422  vr = vr2
423  itopbotexit = -2
424  end if
425 
426  ! Check if velocity is in the same direction throughout cell (i.e. no flow divide).
427  ! Check if product v1*v2 > 0 then the velocity is in the same direction throughout
428  ! the cell (i.e. no flow divide). If so, set vr to reflect appropriate direction.
429  v1v2 = v1 * v2
430  if (v1v2 .gt. dzero) then
431  if (v .gt. dzero) then
432  vr = vr2
433  itopbotexit = -2
434  end if
435  if (v .lt. dzero) then
436  vr = vr1
437  itopbotexit = -1
438  end if
439  end if
440 
441  ! Compute travel time to exit face. Return with status = 0
442  dt = log(abs(vr)) / dvdx
443  status = 0
444  end subroutine calculate_dt
445 
446  !> @brief Calculate the particle's local unscaled xyz coordinates after dt.
447  subroutine calculate_xyz_position(dt, rxx, rxy, ryx, ryy, sxx, sxy, syy, &
448  izstatus, x0, y0, az, vzi, vzbot, &
449  ztop, zbot, zi, x, y, z, exitFace)
450  ! dummy
451  real(DP) :: dt
452  real(DP) :: rxx
453  real(DP) :: rxy
454  real(DP) :: ryx
455  real(DP) :: ryy
456  real(DP) :: sxx
457  real(DP) :: sxy
458  real(DP) :: syy
459  integer(I4B) :: izstatus
460  real(DP) :: x0
461  real(DP) :: y0
462  real(DP) :: az
463  real(DP) :: vzi
464  real(DP) :: vzbot
465  real(DP) :: ztop
466  real(DP) :: zbot
467  real(DP) :: zi
468  real(DP) :: x
469  real(DP) :: y
470  real(DP) :: z
471  integer(I4B), optional :: exitFace
472  ! local
473  integer(I4B) :: lexitface
474  real(DP) :: rot(2, 2), res(2), loc(2)
475  real(DP) :: alp
476  real(DP) :: bet
477 
478  ! process optional exit face argument
479  if (present(exitface)) then
480  lexitface = exitface
481  else
482  lexitface = 0
483  end if
484 
485  ! calculate alpha and beta
486  call step_analytical(dt, alp, bet)
487 
488  ! if exit face is known, set alpha or beta coordinate
489  ! corresponding to the exit face exactly.
490  if (lexitface .eq. 1) then
491  bet = dzero
492  else if (lexitface .eq. 2) then
493  alp = done - bet
494  else if (lexitface .eq. 3) then
495  alp = dzero
496  end if
497 
498  ! if exit face is top or bottom, set z coordinate exactly.
499  if (lexitface .eq. 4) then
500  z = zbot
501  else if (lexitface .eq. 5) then
502  z = ztop
503  else
504  ! otherwise calculate z.
505  if (izstatus .eq. 2) then
506  ! vz uniformly zero
507  z = zi
508  else if (izstatus .eq. 1) then
509  ! vz uniform, nonzero
510  z = zi + vzi * dt
511  else
512  ! vz nonuniform
513  z = zbot + (vzi * dexp(az * dt) - vzbot) / az
514  end if
515  end if
516 
517  ! transform (alp, beta) to (x, y)
518  loc = (/alp, bet/)
519  loc = skew(loc, (/sxx, sxy, syy/), invert=.true.)
520  rot = reshape((/rxx, rxy, ryx, ryy/), shape(rot))
521  res = matmul(rot, loc) ! rotate vector
522  x = res(1) + x0
523  y = res(2) + y0
524 
525  end subroutine calculate_xyz_position
526 
This module contains simulation constants.
Definition: Constants.f90:9
real(dp), parameter dsame
real constant for values that are considered the same based on machine precision
Definition: Constants.f90:122
real(dp), parameter dep3
real constant 1000
Definition: Constants.f90:88
real(dp), parameter donethird
real constant 1/3
Definition: Constants.f90:67
real(dp), parameter dhalf
real constant 1/2
Definition: Constants.f90:68
real(dp), parameter dzero
real constant zero
Definition: Constants.f90:65
real(dp), parameter dtwo
real constant 2
Definition: Constants.f90:79
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
subroutine, public clamp_bary(alpha, beta, gamma, pad)
Clamp barycentric coordinates to the interior of a triangle, with optional padding some minimum dista...
Definition: GeomUtil.f90:466
pure real(dp) function, dimension(2), public skew(v, s, invert)
Skew a 2D vector along the x-axis.
Definition: GeomUtil.f90:149
This module defines variable data types.
Definition: kind.f90:8
pure logical function, public is_close(a, b, rtol, atol, symmetric)
Check if a real value is approximately equal to another.
Definition: MathUtil.f90:46
Particle tracking strategies.
Definition: Method.f90:2
subroutine track_subcell(this, subcell, particle, tmax)
Track a particle across a triangular subcell.
subroutine apply_mst(this, particle, tmax)
Apply the ternary subcell tracking method.
subroutine, public create_method_subcell_ternary(method)
Create a new ternary subcell tracking method.
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.
subroutine calculate_dt(v1, v2, dx, xL, v, dvdx, dt, status, itopbotexit)
Do calculations related to analytical z solution.
@ term_no_exits_sub
terminated in a subcell with no exit face
Definition: Particle.f90:39
subroutine, public create_subcell_tri(subcell)
Create a new triangular subcell.
Definition: SubcellTri.f90:26
subroutine, public canonical(x0, y0, x1, y1, x2, y2, v0x, v0y, v1x, v1y, v2x, v2y, xi, yi, rxx, rxy, ryx, ryy, sxx, sxy, syy, alp0, bet0, alp1, bet1, alp2, bet2, alpi, beti)
Set coordinates to "canonical" configuration.
subroutine, public traverse_triangle(isolv, tol, texit, alpexit, betexit, itrifaceenter, itrifaceexit, alp1, bet1, alp2, bet2, alpi, beti)
Traverse triangular cell.
subroutine, public step_analytical(t, alp, bet)
Step (evaluate) analytically depending on case.
Base type for grid cells of a concrete type. Contains a cell-definition which is information shared b...
Definition: Cell.f90:13
Base type for particle tracking methods.
Definition: Method.f90:31
Ternary triangular subcell tracking method.
Particle tracked by the PRT model.
Definition: Particle.f90:78
A subcell of a cell.
Definition: Subcell.f90:9