MODFLOW 6  version 6.7.0.dev1
USGS Modular Hydrologic Model
MethodCellTernary.f90
Go to the documentation of this file.
2 
3  use kindmodule, only: dp, i4b, lgp
4  use errorutilmodule, only: pstop
5  use methodmodule, only: methodtype
11  use particlemodule
12  use geomutilmodule, only: area, transform
13  use constantsmodule, only: dzero, done, dtwo, dsix
15  implicit none
16 
17  private
18  public :: methodcellternarytype
20 
22  private
23  integer(I4B) :: nverts !< number of vertices
24  real(dp), allocatable, dimension(:) :: xvert
25  real(dp), allocatable, dimension(:) :: yvert !< cell vertex coordinates
26  real(dp), allocatable, dimension(:) :: vne !< cell edge normal velocities
27  real(dp), allocatable, dimension(:) :: vv0x
28  real(dp), allocatable, dimension(:) :: vv0y
29  real(dp), allocatable, dimension(:) :: vv1x
30  real(dp), allocatable, dimension(:) :: vv1y !< cell vertex velocities
31  real(dp) :: xctr
32  real(dp) :: yctr !< cell center coordinates
33  real(dp) :: vctrx
34  real(dp) :: vctry !< cell center velocities
35  real(dp) :: ztop
36  real(dp) :: zbot !< cell top and bottom elevations
37  real(dp) :: dz !< cell thickness
38  real(dp) :: vztop
39  real(dp) :: vzbot !< cell top and bottom velocities
40  integer(I4B), allocatable, dimension(:) :: iprev !< array of shifted indices
41  real(dp), allocatable, dimension(:) :: xvertnext
42  real(dp), allocatable, dimension(:) :: yvertnext !< arrays of shifted cell vertex coordinates
43  integer(I4B), public, pointer :: zeromethod
44  contains
45  procedure, public :: apply => apply_mct
46  procedure, public :: deallocate => destroy_mct
47  procedure, public :: load => load_mct
48  procedure, public :: load_subcell
49  procedure, public :: pass => pass_mct
50  procedure :: vertvelo
51  procedure :: calc_thru_hcsum
52  end type methodcellternarytype
53 
54 contains
55 
56  !> @brief Create a tracking method
57  subroutine create_method_cell_ternary(method)
58  ! dummy
59  type(methodcellternarytype), pointer :: method
60  ! local
61  type(cellpolytype), pointer :: cell
62  type(subcelltritype), pointer :: subcell
63 
64  allocate (method)
65  call create_cell_poly(cell)
66  method%cell => cell
67  method%name => method%cell%type
68  method%delegates = .true.
69  call create_subcell_tri(subcell)
70  method%subcell => subcell
71  end subroutine create_method_cell_ternary
72 
73  !> @brief Destroy the tracking method
74  subroutine destroy_mct(this)
75  class(methodcellternarytype), intent(inout) :: this
76  deallocate (this%name)
77  end subroutine destroy_mct
78 
79  !> @brief Load subcell into tracking method
80  subroutine load_mct(this, particle, next_level, submethod)
81  ! dummy
82  class(methodcellternarytype), intent(inout) :: this
83  type(particletype), pointer, intent(inout) :: particle
84  integer(I4B), intent(in) :: next_level
85  class(methodtype), pointer, intent(inout) :: submethod
86 
87  select type (subcell => this%subcell)
88  type is (subcelltritype)
89  call this%load_subcell(particle, subcell)
90  end select
91 
92  call method_subcell_tern%init( &
93  fmi=this%fmi, &
94  cell=this%cell, &
95  subcell=this%subcell, &
96  trackctl=this%trackctl, &
97  tracktimes=this%tracktimes)
98  submethod => method_subcell_tern
99  end subroutine load_mct
100 
101  !> @brief Pass particle to next subcell if there is one, or to the cell face
102  subroutine pass_mct(this, particle)
103  ! dummy
104  class(methodcellternarytype), intent(inout) :: this
105  type(particletype), pointer, intent(inout) :: particle
106  ! local
107  integer(I4B) :: isc
108  integer(I4B) :: exitFace
109  integer(I4B) :: inface
110 
111  exitface = particle%iboundary(3)
112  isc = particle%idomain(3)
113 
114  select case (exitface)
115  case (0)
116  ! Subcell interior (cell interior)
117  inface = -1
118  case (1)
119  ! Subcell face 1 (cell face)
120  inface = isc
121  if (inface .eq. 0) inface = this%nverts
122  case (2)
123  ! Subcell face --> next subcell in "cycle" (cell interior)
124  isc = isc + 1
125  if (isc .gt. this%nverts) isc = 1
126  particle%idomain(3) = isc
127  particle%iboundary(3) = 3
128  inface = 0
129  case (3)
130  ! Subcell face --> preceding subcell in "cycle" (cell interior)
131  isc = isc - 1
132  if (isc .lt. 1) isc = this%nverts
133  particle%idomain(3) = isc
134  particle%iboundary(3) = 2
135  inface = 0
136  case (4)
137  ! Subcell bottom (cell bottom)
138  inface = this%nverts + 2
139  case (5)
140  ! Subcell top (cell top)
141  inface = this%nverts + 3
142  end select
143 
144  if (inface .eq. -1) then
145  particle%iboundary(2) = 0
146  else if (inface .eq. 0) then
147  particle%iboundary(2) = 0
148  else
149  particle%iboundary(2) = inface
150  end if
151  end subroutine pass_mct
152 
153  !> @brief Apply the ternary method to a polygonal cell
154  subroutine apply_mct(this, particle, tmax)
155  use constantsmodule, only: dzero, done, dhalf
156  ! dummy
157  class(methodcellternarytype), intent(inout) :: this
158  type(particletype), pointer, intent(inout) :: particle
159  real(DP), intent(in) :: tmax
160  ! local
161  integer(I4B) :: i
162  real(DP) :: x, y, z, xO, yO
163  real(DP), allocatable :: xs(:), ys(:)
164 
165  ! (Re)allocate type-bound arrays
166  select type (cell => this%cell)
167  type is (cellpolytype)
168  ! Check termination/reporting conditions
169  call this%check(particle, this%cell%defn, tmax)
170  if (.not. particle%advancing) return
171 
172  ! Number of vertices
173  this%nverts = cell%defn%npolyverts
174 
175  ! (Re)allocate type-bound arrays
176  if (allocated(this%xvert)) then
177  deallocate (this%xvert)
178  deallocate (this%yvert)
179  deallocate (this%vne)
180  deallocate (this%vv0x)
181  deallocate (this%vv0y)
182  deallocate (this%vv1x)
183  deallocate (this%vv1y)
184  deallocate (this%iprev)
185  deallocate (this%xvertnext)
186  deallocate (this%yvertnext)
187  end if
188 
189  allocate (this%xvert(this%nverts))
190  allocate (this%yvert(this%nverts))
191  allocate (this%vne(this%nverts))
192  allocate (this%vv0x(this%nverts))
193  allocate (this%vv0y(this%nverts))
194  allocate (this%vv1x(this%nverts))
195  allocate (this%vv1y(this%nverts))
196  allocate (this%iprev(this%nverts))
197  allocate (this%xvertnext(this%nverts))
198  allocate (this%yvertnext(this%nverts))
199 
200  ! New origin is the corner of the cell's
201  ! bounding box closest to the old origin
202  allocate (xs(this%nverts))
203  allocate (ys(this%nverts))
204  xs = cell%defn%polyvert(1, :)
205  ys = cell%defn%polyvert(2, :)
206  xo = xs(minloc(abs(xs), dim=1))
207  yo = ys(minloc(abs(ys), dim=1))
208  deallocate (xs)
209  deallocate (ys)
210 
211  ! Cell vertices
212  do i = 1, this%nverts
213  x = cell%defn%polyvert(1, i)
214  y = cell%defn%polyvert(2, i)
215  call transform(x, y, dzero, x, y, z, xo, yo)
216  this%xvert(i) = x
217  this%yvert(i) = y
218  end do
219 
220  ! Top, bottom, and thickness
221  this%ztop = cell%defn%top
222  this%zbot = cell%defn%bot
223  this%dz = this%ztop - this%zbot
224 
225  ! Shifted arrays
226  do i = 1, this%nverts
227  this%iprev(i) = i
228  end do
229  this%iprev = cshift(this%iprev, -1)
230  this%xvertnext = cshift(this%xvert, 1)
231  this%yvertnext = cshift(this%yvert, 1)
232 
233  ! Calculate vertex velocities.
234  call this%vertvelo()
235 
236  ! Transform particle coordinates
237  call particle%transform(xo, yo)
238 
239  ! Track the particle across the cell.
240  call this%track(particle, 2, tmax)
241 
242  ! Transform particle coordinates back
243  call particle%transform(xo, yo, invert=.true.)
244  call particle%reset_transform()
245 
246  end select
247 
248  end subroutine apply_mct
249 
250  !> @brief Loads a triangular subcell from the polygonal cell
251  subroutine load_subcell(this, particle, subcell)
252  ! dummy
253  class(methodcellternarytype), intent(inout) :: this
254  type(particletype), pointer, intent(inout) :: particle
255  class(subcelltritype), intent(inout) :: subcell
256  ! local
257  integer(I4B) :: ic
258  integer(I4B) :: isc
259  integer(I4B) :: iv0
260  real(DP) :: x0
261  real(DP) :: y0
262  real(DP) :: x1
263  real(DP) :: y1
264  real(DP) :: x2
265  real(DP) :: y2
266  real(DP) :: x1rel
267  real(DP) :: y1rel
268  real(DP) :: x2rel
269  real(DP) :: y2rel
270  real(DP) :: xi
271  real(DP) :: yi
272  real(DP) :: di2
273  real(DP) :: d02
274  real(DP) :: d12
275  real(DP) :: di1
276  real(DP) :: d01
277  real(DP) :: alphai
278  real(DP) :: betai
279 
280  select type (cell => this%cell)
281  type is (cellpolytype)
282  ic = cell%defn%icell
283  subcell%icell = ic
284  isc = particle%idomain(3)
285  if (isc .le. 0) then
286  xi = particle%x
287  yi = particle%y
288  do iv0 = 1, this%nverts
289  x0 = this%xvert(iv0)
290  y0 = this%yvert(iv0)
291  x1 = this%xvertnext(iv0)
292  y1 = this%yvertnext(iv0)
293  x2 = this%xctr
294  y2 = this%yctr
295  x1rel = x1 - x0
296  y1rel = y1 - y0
297  x2rel = x2 - x0
298  y2rel = y2 - y0
299  di2 = xi * y2rel - yi * x2rel
300  d02 = x0 * y2rel - y0 * x2rel
301  d12 = x1rel * y2rel - y1rel * x2rel
302  di1 = xi * y1rel - yi * x1rel
303  d01 = x0 * y1rel - y0 * x1rel
304  alphai = (di2 - d02) / d12
305  betai = -(di1 - d01) / d12
306  ! assumes particle is in the cell, so no check needed for beta
307  if ((alphai .ge. dzero) .and. &
308  (alphai + betai .le. done)) then
309  isc = iv0
310  exit
311  end if
312  end do
313  if (isc .le. 0) then
314  print *, "error -- initial triangle not found in cell ", ic, &
315  " for particle at ", particle%x, particle%y, particle%z
316 
317  call pstop(1)
318  else
319  particle%idomain(3) = isc
320  end if
321  end if
322  subcell%isubcell = isc
323 
324  ! Set coordinates and velocities at vertices of triangular subcell
325  iv0 = isc
326  subcell%x0 = this%xvert(iv0)
327  subcell%y0 = this%yvert(iv0)
328  subcell%x1 = this%xvertnext(iv0)
329  subcell%y1 = this%yvertnext(iv0)
330  subcell%x2 = this%xctr
331  subcell%y2 = this%yctr
332  subcell%v0x = this%vv0x(iv0)
333  subcell%v0y = this%vv0y(iv0)
334  subcell%v1x = this%vv1x(iv0) ! the indices here actually refer to subcells, not vertices
335  subcell%v1y = this%vv1y(iv0)
336  subcell%v2x = this%vctrx
337  subcell%v2y = this%vctry
338  subcell%ztop = this%ztop
339  subcell%zbot = this%zbot
340  subcell%dz = this%dz
341  subcell%vzbot = this%vzbot
342  subcell%vztop = this%vztop
343  end select
344  end subroutine load_subcell
345 
346  !> @brief Calculate vertex velocities
347  subroutine vertvelo(this)
348  use constantsmodule, only: dzero, done, dhalf
349  ! dummy
350  class(methodcellternarytype), intent(inout) :: this
351  ! local
352  real(DP) :: term
353  integer(I4B) :: i
354  real(DP) :: perturb
355  real(DP), allocatable, dimension(:) :: xvals
356  real(DP), allocatable, dimension(:) :: yvals
357  real(DP) :: sixa
358  real(DP) :: vm0i0
359  real(DP) :: vm0ival
360  real(DP) :: hcsum0
361  real(DP) :: hcsum
362  real(DP) :: jac
363  real(DP), allocatable, dimension(:) :: wk1
364  real(DP), allocatable, dimension(:) :: wk2
365  real(DP), allocatable, dimension(:) :: unextxnext
366  real(DP), allocatable, dimension(:) :: unextynext
367  real(DP), allocatable, dimension(:) :: le
368  real(DP), allocatable, dimension(:) :: unextx
369  real(DP), allocatable, dimension(:) :: unexty
370  real(DP) :: areacell
371  real(DP), allocatable, dimension(:) :: areasub
372  real(DP) :: divcell
373  real(DP), allocatable, dimension(:) :: li
374  real(DP), allocatable, dimension(:) :: unintx
375  real(DP), allocatable, dimension(:) :: uninty
376  real(DP), allocatable, dimension(:) :: xmid
377  real(DP), allocatable, dimension(:) :: ymid
378  real(DP), allocatable, dimension(:) :: lm
379  real(DP), allocatable, dimension(:) :: umx
380  real(DP), allocatable, dimension(:) :: umy
381  real(DP), allocatable, dimension(:) :: kappax
382  real(DP), allocatable, dimension(:) :: kappay
383  real(DP), allocatable, dimension(:) :: vm0x
384  real(DP), allocatable, dimension(:) :: vm0y
385  real(DP), allocatable, dimension(:) :: vm1x
386  real(DP), allocatable, dimension(:) :: vm1y
387  ! real(DP), allocatable, dimension(:, :) :: poly
388  integer(I4B) :: nvert
389 
390  select type (cell => this%cell)
391  type is (cellpolytype)
392 
393  ! Allocate local arrays
394  allocate (le(this%nverts)) ! lengths of exterior (cell) edges
395  allocate (unextx(this%nverts)) ! x components of unit normals to exterior edges
396  allocate (unexty(this%nverts)) ! y components of unit normals to exterior edges
397  allocate (areasub(this%nverts)) ! subcell areas
398  allocate (li(this%nverts)) ! lengths of interior edges ("spokes")
399  allocate (unintx(this%nverts)) ! x components of unit normals to interior edges
400  allocate (uninty(this%nverts)) ! y components of unit normals to interior edges
401  allocate (xmid(this%nverts)) ! x coordinates of midpoints
402  allocate (ymid(this%nverts)) ! y coordinates of midpoints
403  allocate (lm(this%nverts)) ! lengths of midpoint connectors
404  allocate (umx(this%nverts)) ! x components of midpoint-connector (cw) unit vectors
405  allocate (umy(this%nverts)) ! y components of midpoint-connector (cw) unit vectors
406  allocate (kappax(this%nverts)) ! x components of kappa vectors
407  allocate (kappay(this%nverts)) ! y components of kappa vectors
408  allocate (vm0x(this%nverts)) ! x component of vm0
409  allocate (vm0y(this%nverts)) ! y component of vm0
410  allocate (vm1x(this%nverts)) ! x component of vm1
411  allocate (vm1y(this%nverts)) ! y component of vm1
412  allocate (unextxnext(this%nverts)) ! vector of "next" interior unit-normal x coordinates defined for convenience
413  allocate (unextynext(this%nverts)) ! vector of "next" interior unit-normal y coordinates defined for convenience
414  allocate (wk1(this%nverts))
415  allocate (wk2(this%nverts))
416  allocate (xvals(3))
417  allocate (yvals(3))
418 
419  ! Exterior edge unit normals (outward) and lengths
420  wk1 = this%xvertnext - this%xvert
421  wk2 = this%yvertnext - this%yvert
422  le = dsqrt(wk1 * wk1 + wk2 * wk2)
423  unextx = -wk2 / le
424  unexty = wk1 / le
425 
426  ! Cell centroid (in general, this is NOT the average of the vertex coordinates)
427  areacell = area(this%xvert, this%yvert)
428  sixa = areacell * dsix
429  wk1 = -(this%xvert * this%yvertnext - this%xvertnext * this%yvert)
430  nvert = size(this%xvert)
431  this%xctr = sum((this%xvert + this%xvertnext) * wk1) / sixa
432  this%yctr = sum((this%yvert + this%yvertnext) * wk1) / sixa
433 
434  ! TODO: can we use some of the terms of the centroid calculation
435  ! to do a cheap point in polygon check?
436 
437  ! Subcell areas
438  do i = 1, this%nverts
439  xvals(1) = this%xvert(i)
440  xvals(2) = this%xvertnext(i)
441  xvals(3) = this%xctr
442  yvals(1) = this%yvert(i)
443  yvals(2) = this%yvertnext(i)
444  yvals(3) = this%yctr
445  areasub(i) = area(xvals, yvals)
446  end do
447 
448  ! Cell-edge normal velocities (outward)
449  term = done / (cell%defn%porosity * cell%defn%retfactor * this%dz)
450  do i = 1, this%nverts
451  this%vne(i) = -cell%defn%faceflow(i) * term / le(i)
452  end do
453 
454  ! Cell divergence (2D)
455  divcell = sum(le * this%vne) / areacell
456 
457  ! Interior edge (cw) unit normals and lengths
458  wk1 = this%xvert - this%xctr
459  wk2 = this%yvert - this%yctr
460  li = dsqrt(wk1 * wk1 + wk2 * wk2)
461  unintx = wk2 / li
462  uninty = -wk1 / li
463  ! Shifted arrays for convenience
464  unextxnext = cshift(unintx, 1)
465  unextynext = cshift(uninty, 1)
466 
467  ! Midpoints of interior edges
468  xmid = 5.d-1 * (this%xvert + this%xctr)
469  ymid = 5.d-1 * (this%yvert + this%yctr)
470 
471  ! Unit midpoint-connector (cw) vectors and lengths
472  wk1 = cshift(xmid, 1) - xmid
473  wk2 = cshift(ymid, 1) - ymid
474  lm = dsqrt(wk1 * wk1 + wk2 * wk2)
475  umx = wk1 / lm
476  umy = wk2 / lm
477 
478  ! Kappa vectors (K tensor times unit midpoint-connector vectors) do not
479  ! account for anisotropy, which is consistent with the way internal face
480  ! flow calculations are done in MP7. The isotropic value of K does not
481  ! matter in this case because it cancels out of the calculations, so
482  ! K = 1 is assumed for simplicity.
483  kappax = umx
484  kappay = umy
485 
486  ! Use linearity to find vm0i[0] such that curl of the head gradient
487  ! is zero
488  perturb = 1.d-2
489  ! Calculations at base value
490  vm0i0 = 0.d0
491  call this%calc_thru_hcsum(vm0i0, divcell, le, li, lm, areasub, areacell, &
492  unintx, uninty, unextx, unexty, &
493  unextxnext, unextynext, &
494  kappax, kappay, vm0x, vm0y, vm1x, vm1y, hcsum0)
495  ! Calculations at perturbed value
496  vm0ival = vm0i0 + perturb
497  call this%calc_thru_hcsum(vm0ival, divcell, le, li, lm, areasub, areacell, &
498  unintx, uninty, unextx, unexty, &
499  unextxnext, unextynext, &
500  kappax, kappay, vm0x, vm0y, vm1x, vm1y, hcsum)
501  ! Calculations at root value
502  jac = (hcsum - hcsum0) / perturb
503  vm0ival = vm0i0 - hcsum0 / jac
504  call this%calc_thru_hcsum(vm0ival, divcell, le, li, lm, areasub, areacell, &
505  unintx, uninty, unextx, unexty, &
506  unextxnext, unextynext, &
507  kappax, kappay, vm0x, vm0y, vm1x, vm1y, hcsum)
508 
509  ! Project linearly to get corner (vertex) velocities. Note that velocity
510  ! vv1 is at the next vertex cw from vv0, so vv0(i) and vv1(i) are the
511  ! two vertex velocities used by triangular subcell i.
512  this%vv0x = 2.d0 * vm0x - this%vctrx
513  this%vv0y = 2.d0 * vm0y - this%vctry
514  this%vv1x = 2.d0 * vm1x - this%vctrx
515  this%vv1y = 2.d0 * vm1y - this%vctry
516 
517  ! Set top and bottom velocities
518  term = done / (cell%defn%retfactor * cell%defn%porosity * areacell)
519  this%vzbot = cell%defn%faceflow(this%nverts + 2) * term
520  this%vztop = -cell%defn%faceflow(this%nverts + 3) * term
521 
522  ! Deallocate local arrays
523  deallocate (le)
524  deallocate (unextx)
525  deallocate (unexty)
526  deallocate (areasub)
527  deallocate (li)
528  deallocate (unintx)
529  deallocate (uninty)
530  deallocate (xmid)
531  deallocate (ymid)
532  deallocate (lm)
533  deallocate (umx)
534  deallocate (umy)
535  deallocate (kappax)
536  deallocate (kappay)
537  deallocate (vm0x)
538  deallocate (vm0y)
539  deallocate (vm1x)
540  deallocate (vm1y)
541  deallocate (unextxnext)
542  deallocate (unextynext)
543  deallocate (wk1)
544  deallocate (wk2)
545  deallocate (xvals)
546  deallocate (yvals)
547 
548  end select
549  end subroutine vertvelo
550 
551  subroutine calc_thru_hcsum(this, vm0ival, divcell, &
552  le, li, lm, areasub, areacell, &
553  unintx, uninty, unextx, unexty, &
554  unintxnext, unintynext, &
555  kappax, kappay, vm0x, vm0y, vm1x, vm1y, hcsum)
556  ! dummy
557  class(methodcellternarytype), intent(inout) :: this
558  real(DP) :: vm0ival
559  real(DP) :: divcell
560  real(DP) :: hcsum
561  real(DP), dimension(:) :: le
562  real(DP), dimension(:) :: li
563  real(DP), dimension(:) :: lm
564  real(DP), dimension(:) :: areasub
565  real(DP) :: areacell
566  real(DP), dimension(:) :: unintx
567  real(DP), dimension(:) :: uninty
568  real(DP), dimension(:) :: unextx
569  real(DP), dimension(:) :: unexty
570  real(DP), dimension(:) :: unintxnext
571  real(DP), dimension(:) :: unintynext
572  real(DP), dimension(:) :: kappax
573  real(DP), dimension(:) :: kappay
574  real(DP), dimension(:) :: vm0x
575  real(DP), dimension(:) :: vm0y
576  real(DP), dimension(:) :: vm1x
577  real(DP), dimension(:) :: vm1y
578  ! local
579  real(DP), allocatable, dimension(:) :: vm0i
580  real(DP), allocatable, dimension(:) :: vm0e
581  real(DP), allocatable, dimension(:) :: vm1i
582  real(DP), allocatable, dimension(:) :: vm1e
583  real(DP), allocatable, dimension(:) :: uprod
584  real(DP), allocatable, dimension(:) :: det
585  real(DP), allocatable, dimension(:) :: wt
586  real(DP), allocatable, dimension(:) :: bi0x
587  real(DP), allocatable, dimension(:) :: be0x
588  real(DP), allocatable, dimension(:) :: bi0y
589  real(DP), allocatable, dimension(:) :: be0y
590  real(DP), allocatable, dimension(:) :: bi1x
591  real(DP), allocatable, dimension(:) :: be1x
592  real(DP), allocatable, dimension(:) :: bi1y
593  real(DP), allocatable, dimension(:) :: be1y
594  real(DP), allocatable, dimension(:) :: be01x
595  real(DP), allocatable, dimension(:) :: be01y
596  real(DP) :: emxx
597  real(DP) :: emxy
598  real(DP) :: emyx
599  real(DP) :: emyy
600  real(DP) :: rx
601  real(DP) :: ry
602  real(DP) :: emdet
603  integer(I4B) :: i
604  integer(I4B) :: ip
605 
606  ! Allocate local arrays
607  allocate (vm0i(this%nverts))
608  allocate (vm0e(this%nverts))
609  allocate (vm1i(this%nverts))
610  allocate (vm1e(this%nverts))
611  allocate (uprod(this%nverts))
612  allocate (det(this%nverts))
613  allocate (wt(this%nverts))
614  allocate (bi0x(this%nverts))
615  allocate (be0x(this%nverts))
616  allocate (bi0y(this%nverts))
617  allocate (be0y(this%nverts))
618  allocate (bi1x(this%nverts))
619  allocate (be1x(this%nverts))
620  allocate (bi1y(this%nverts))
621  allocate (be1y(this%nverts))
622  allocate (be01x(this%nverts))
623  allocate (be01y(this%nverts))
624 
625  ! Set vm0i(1)
626  vm0i(1) = vm0ival
627 
628  ! Get remaining vm0i's sequentially using divergence conditions
629  do i = 2, this%nverts
630  ip = this%iprev(i)
631  vm0i(i) = (li(ip) * vm0i(ip) - le(ip) * this%vne(ip) &
632  + areasub(ip) * divcell) / li(i)
633  end do
634 
635  ! Get vm1i's from vm0i's using continuity conditions
636  vm1i = cshift(vm0i, 1)
637 
638  ! Get centroid velocity by setting up and solving 2x2 linear system
639  uprod = unintx * unextx + uninty * unexty
640  det = done - uprod * uprod
641  bi0x = (unintx - unextx * uprod) / det
642  be0x = (unextx - unintx * uprod) / det
643  bi0y = (uninty - unexty * uprod) / det
644  be0y = (unexty - uninty * uprod) / det
645  uprod = unintxnext * unextx + unintynext * unexty
646  det = done - uprod * uprod
647  bi1x = (unintxnext - unextx * uprod) / det
648  be1x = (unextx - unintxnext * uprod) / det
649  bi1y = (unintynext - unexty * uprod) / det
650  be1y = (unexty - unintynext * uprod) / det
651  be01x = 5.d-1 * (be0x + be1x)
652  be01y = 5.d-1 * (be0y + be1y)
653  wt = areasub / areacell
654  emxx = dtwo - sum(wt * be01x * unextx)
655  emxy = -sum(wt * be01x * unexty)
656  emyx = -sum(wt * be01y * unextx)
657  emyy = dtwo - sum(wt * be01y * unexty)
658  rx = sum(wt * (bi0x * vm0i + bi1x * vm1i + be01x * this%vne))
659  ry = sum(wt * (bi0y * vm0i + bi1y * vm1i + be01y * this%vne))
660  emdet = emxx * emyy - emxy * emyx
661  this%vctrx = (emyy * rx - emxy * ry) / emdet
662  this%vctry = (emxx * ry - emyx * rx) / emdet
663 
664  ! Get vm0e's using "known" conditions
665  vm0e = 5.d-1 * (this%vne + unextx * this%vctrx + unexty * this%vctry)
666 
667  ! Get vm1e's from uniformity along exterior edges
668  vm1e = vm0e
669 
670  ! Transform vm0 and vm1 to (x, y) coordinates
671  vm0x = bi0x * vm0i + be0x * vm0e
672  vm0y = bi0y * vm0i + be0y * vm0e
673  vm1x = bi1x * vm1i + be1x * vm0e
674  vm1y = bi1y * vm1i + be1y * vm0e
675 
676  ! Calculate head-cycle summation (which is proportional to
677  ! the curl of the head gradient)
678  hcsum = sum(lm * (kappax * (vm0x + vm1x) + kappay * (vm0y + vm1y)))
679 
680  ! Deallocate local arrays
681  deallocate (vm0i)
682  deallocate (vm0e)
683  deallocate (vm1i)
684  deallocate (vm1e)
685  deallocate (uprod)
686  deallocate (det)
687  deallocate (wt)
688  deallocate (bi0x)
689  deallocate (be0x)
690  deallocate (bi0y)
691  deallocate (be0y)
692  deallocate (bi1x)
693  deallocate (be1x)
694  deallocate (bi1y)
695  deallocate (be1y)
696  deallocate (be01x)
697  deallocate (be01y)
698 
699  end subroutine calc_thru_hcsum
700 
701 end module methodcellternarymodule
702 
subroutine, public create_cell_poly(cell)
Create a new polygonal cell.
Definition: CellPoly.f90:20
This module contains simulation constants.
Definition: Constants.f90:9
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 dsix
real constant 6
Definition: Constants.f90:82
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 transform(xin, yin, zin, xout, yout, zout, xorigin, yorigin, zorigin, sinrot, cosrot, invert)
Apply a 3D translation and optional 2D rotation to coordinates.
Definition: GeomUtil.f90:183
real(dp) function, public area(xv, yv, cw)
Calculate polygon area, with vertices given in CW or CCW order.
Definition: GeomUtil.f90:370
logical function, public point_in_polygon(x, y, poly)
Check if a point is within a polygon.
Definition: GeomUtil.f90:27
This module defines variable data types.
Definition: kind.f90:8
subroutine apply_mct(this, particle, tmax)
Apply the ternary method to a polygonal cell.
subroutine load_mct(this, particle, next_level, submethod)
Load subcell into tracking method.
subroutine destroy_mct(this)
Destroy the tracking method.
subroutine, public create_method_cell_ternary(method)
Create a tracking method.
subroutine load_subcell(this, particle, subcell)
Loads a triangular subcell from the polygonal cell.
subroutine calc_thru_hcsum(this, vm0ival, divcell, le, li, lm, areasub, areacell, unintx, uninty, unextx, unexty, unintxnext, unintynext, kappax, kappay, vm0x, vm0y, vm1x, vm1y, hcsum)
subroutine vertvelo(this)
Calculate vertex velocities.
subroutine pass_mct(this, particle)
Pass particle to next subcell if there is one, or to the cell face.
Particle tracking strategies.
Definition: Method.f90:2
Subcell-level tracking methods.
type(methodsubcellternarytype), pointer, public method_subcell_tern
subroutine, public create_subcell_tri(subcell)
Create a new triangular subcell.
Definition: SubcellTri.f90:26
Base type for particle tracking methods.
Definition: Method.f90:31
Particle tracked by the PRT model.
Definition: Particle.f90:78