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
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 the particle's exit face and travel time to exit.
222  ! The exit face is the face through which it would exit first,
223  ! considering only the velocity component in the direction of
224  ! the face. Then compute the particle's exit time.
225  if (itrifaceexit /= 0) then
226  ! Exit through lateral subcell face
227  exitface = itrifaceexit
228  dtexit = dtexitxy
229  else if (dtexitz < dtexitxy .and. dtexitz >= 0.0_dp) then
230  ! Exit through top or bottom
231  if (itopbotexit == -1) then
232  exitface = 4
233  else
234  exitface = 5
235  end if
236  dtexit = dtexitz
237  else
238  particle%istatus = term_no_exits_sub
239  particle%advancing = .false.
240  call this%save(particle, reason=3)
241  return
242  end if
243  texit = particle%ttrack + dtexit
244  t0 = particle%ttrack
245 
246  ! Select user tracking times to solve. If this is the last time step
247  ! in the simulation, times falling after the simulation end time are
248  ! only included if the 'extend' option is on, otherwise only times in
249  ! the time step are included.
250  call this%tracktimes%advance()
251  if (this%tracktimes%any()) then
252  do i = this%tracktimes%selection(1), this%tracktimes%selection(2)
253  t = this%tracktimes%times(i)
254  if (t < particle%ttrack) cycle
255  if (t >= texit .or. t >= tmax) exit
256  dt = t - t0
257  call calculate_xyz_position(dt, rxx, rxy, ryx, ryy, sxx, sxy, syy, &
258  izstatus, x0, y0, az, vzi, vzbot, &
259  ztop, zbot, zi, x, y, z)
260  particle%x = x
261  particle%y = y
262  particle%z = z
263  particle%ttrack = t
264  particle%istatus = active
265  call this%save(particle, reason=5)
266  end do
267  end if
268 
269  ! Compute exit time and face and update the particle's coordinates
270  ! (local, unscaled) and other properties. The particle may at this
271  ! point lie on a boundary of the subcell or may still be within it.
272  if (texit .gt. tmax) then
273  ! The computed exit time is greater than the maximum time, so set
274  ! final time for particle trajectory equal to maximum time.
275  t = tmax
276  dt = t - t0
277  exitface = 0
278  particle%istatus = active
279  particle%advancing = .false.
280  reason = 2 ! timestep end
281  else
282  ! The computed exit time is less than or equal to the maximum time,
283  ! so set final time for particle trajectory equal to exit time.
284  t = texit
285  dt = dtexit
286  reason = 1 ! (sub)cell transition
287  end if
288  call calculate_xyz_position(dt, rxx, rxy, ryx, ryy, sxx, sxy, syy, &
289  izstatus, x0, y0, az, vzi, vzbot, &
290  ztop, zbot, zi, x, y, z, exitface)
291  particle%x = x
292  particle%y = y
293  particle%z = z
294  particle%ttrack = t
295  particle%iboundary(3) = exitface
296 
297  call this%save(particle, reason=reason)
298  end subroutine track_subcell
299 
300  !> @brief Do calculations related to analytical z solution
301  !!
302  !! This subroutine consists partly or entirely of code written by
303  !! David W. Pollock of the USGS for MODPATH 7. The authors of the present
304  !! code are responsible for its appropriate application in this context
305  !! and for any modifications or errors.
306  !<
307  subroutine calculate_dt(v1, v2, dx, xL, v, dvdx, &
308  dt, status, itopbotexit)
309  real(DP) :: v1
310  real(DP) :: v2
311  real(DP) :: dx
312  real(DP) :: xL
313  real(DP) :: v
314  real(DP) :: dvdx
315  real(DP) :: dt
316  real(DP) :: v2a
317  real(DP) :: v1a
318  real(DP) :: dv
319  real(DP) :: dva
320  real(DP) :: vv
321  real(DP) :: vvv
322  real(DP) :: zro
323  real(DP) :: zrom
324  real(DP) :: x
325  real(DP) :: tol
326  real(DP) :: vr1
327  real(DP) :: vr2
328  real(DP) :: vr
329  real(DP) :: v1v2
330  integer(I4B) :: status
331  integer(I4B) :: itopbotexit
332  logical(LGP) :: noOutflow
333 
334  ! Initialize variables
335  status = -1
336  dt = 1.0d+20
337  v2a = v2
338  if (v2a .lt. dzero) v2a = -v2a
339  v1a = v1
340  if (v1a .lt. dzero) v1a = -v1a
341  dv = v2 - v1
342  dva = dv
343  if (dva .lt. dzero) dva = -dva
344 
345  ! Check for a uniform zero velocity in this direction.
346  ! If so, set status = 2 and return (dt = 1.0d+20).
347  tol = 1.0d-15
348  if ((v2a .lt. tol) .and. (v1a .lt. tol)) then
349  v = dzero
350  dvdx = dzero
351  status = 2
352  itopbotexit = 0
353  return
354  end if
355 
356  ! Check for uniform non-zero velocity in this direction.
357  ! If so, set compute dt using the constant velocity,
358  ! set status = 1 and return.
359  vv = v1a
360  if (v2a .gt. vv) vv = v2a
361  vvv = dva / vv
362  if (vvv .lt. 1.0d-4) then
363  zro = tol
364  zrom = -zro
365  v = v1
366  x = xl * dx
367  if (v1 .gt. zro) then
368  dt = (dx - x) / v1
369  itopbotexit = -2
370  end if
371  if (v1 .lt. zrom) then
372  dt = -x / v1
373  itopbotexit = -1
374  end if
375  dvdx = dzero
376  status = 1
377  return
378  end if
379 
380  ! Velocity has a linear variation.
381  ! Compute velocity corresponding to particle position
382  dvdx = dv / dx
383  v = (done - xl) * v1 + xl * v2
384 
385  ! If flow is into the cell from both sides there is no outflow.
386  ! In that case, set status = 3 and return
387  nooutflow = .true.
388  if (v1 .lt. dzero) nooutflow = .false.
389  if (v2 .gt. dzero) nooutflow = .false.
390  if (nooutflow) then
391  status = 3
392  itopbotexit = 0
393  return
394  end if
395 
396  ! If there is a divide in the cell for this flow direction, check to see if the
397  ! particle is located exactly on the divide. If it is, move it very slightly to
398  ! get it off the divide. This avoids possible numerical problems related to
399  ! stagnation points.
400  if ((v1 .le. dzero) .and. (v2 .ge. dzero)) then
401  if (abs(v) .le. dzero) then
402  v = 1.0d-20
403  if (v2 .le. dzero) v = -v
404  end if
405  end if
406 
407  ! If there is a flow divide, find out what side of the divide the particle
408  ! is on and set the value of vr appropriately to reflect that location.
409  vr1 = v1 / v
410  vr2 = v2 / v
411  vr = vr1
412  itopbotexit = -1
413  if (vr .le. dzero) then
414  vr = vr2
415  itopbotexit = -2
416  end if
417 
418  ! Check if velocity is in the same direction throughout cell (i.e. no flow divide).
419  ! Check if product v1*v2 > 0 then the velocity is in the same direction throughout
420  ! the cell (i.e. no flow divide). If so, set vr to reflect appropriate direction.
421  v1v2 = v1 * v2
422  if (v1v2 .gt. dzero) then
423  if (v .gt. dzero) then
424  vr = vr2
425  itopbotexit = -2
426  end if
427  if (v .lt. dzero) then
428  vr = vr1
429  itopbotexit = -1
430  end if
431  end if
432 
433  ! Compute travel time to exit face. Return with status = 0
434  dt = log(abs(vr)) / dvdx
435  status = 0
436  end subroutine calculate_dt
437 
438  !> @brief Calculate the particle's local unscaled xyz coordinates after dt.
439  subroutine calculate_xyz_position(dt, rxx, rxy, ryx, ryy, sxx, sxy, syy, &
440  izstatus, x0, y0, az, vzi, vzbot, &
441  ztop, zbot, zi, x, y, z, exitFace)
442  ! dummy
443  real(DP) :: dt
444  real(DP) :: rxx
445  real(DP) :: rxy
446  real(DP) :: ryx
447  real(DP) :: ryy
448  real(DP) :: sxx
449  real(DP) :: sxy
450  real(DP) :: syy
451  integer(I4B) :: izstatus
452  real(DP) :: x0
453  real(DP) :: y0
454  real(DP) :: az
455  real(DP) :: vzi
456  real(DP) :: vzbot
457  real(DP) :: ztop
458  real(DP) :: zbot
459  real(DP) :: zi
460  real(DP) :: x
461  real(DP) :: y
462  real(DP) :: z
463  integer(I4B), optional :: exitFace
464  ! local
465  integer(I4B) :: lexitface
466  real(DP) :: rot(2, 2), res(2), loc(2)
467  real(DP) :: alp
468  real(DP) :: bet
469 
470  ! process optional exit face argument
471  if (present(exitface)) then
472  lexitface = exitface
473  else
474  lexitface = 0
475  end if
476 
477  ! calculate alpha and beta
478  call step_analytical(dt, alp, bet)
479 
480  ! if exit face is known, set alpha or beta coordinate
481  ! corresponding to the exit face exactly.
482  if (lexitface .eq. 1) then
483  bet = dzero
484  else if (lexitface .eq. 2) then
485  alp = done - bet
486  else if (lexitface .eq. 3) then
487  alp = dzero
488  end if
489 
490  ! if exit face is top or bottom, set z coordinate exactly.
491  if (lexitface .eq. 4) then
492  z = zbot
493  else if (lexitface .eq. 5) then
494  z = ztop
495  else
496  ! otherwise calculate z.
497  if (izstatus .eq. 2) then
498  ! vz uniformly zero
499  z = zi
500  else if (izstatus .eq. 1) then
501  ! vz uniform, nonzero
502  z = zi + vzi * dt
503  else
504  ! vz nonuniform
505  z = zbot + (vzi * dexp(az * dt) - vzbot) / az
506  end if
507  end if
508 
509  ! transform (alp, beta) to (x, y)
510  loc = (/alp, bet/)
511  loc = skew(loc, (/sxx, sxy, syy/), invert=.true.)
512  rot = reshape((/rxx, rxy, ryx, ryy/), shape(rot))
513  res = matmul(rot, loc) ! rotate vector
514  x = res(1) + x0
515  y = res(2) + y0
516 
517  end subroutine calculate_xyz_position
518 
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