23 integer(I4B) :: nverts
24 real(dp),
allocatable,
dimension(:) :: xvert
25 real(dp),
allocatable,
dimension(:) :: yvert
26 real(dp),
allocatable,
dimension(:) :: vne
27 real(dp),
allocatable,
dimension(:) :: vv0x
28 real(dp),
allocatable,
dimension(:) :: vv0y
29 real(dp),
allocatable,
dimension(:) :: vv1x
30 real(dp),
allocatable,
dimension(:) :: vv1y
40 integer(I4B),
allocatable,
dimension(:) :: iprev
41 real(dp),
allocatable,
dimension(:) :: xvertnext
42 real(dp),
allocatable,
dimension(:) :: yvertnext
43 integer(I4B),
public,
pointer :: zeromethod
49 procedure,
public :: pass =>
pass_mct
67 method%name => method%cell%type
68 method%delegates = .true.
70 method%subcell => subcell
76 deallocate (this%name)
80 subroutine load_mct(this, particle, next_level, submethod)
84 integer(I4B),
intent(in) :: next_level
85 class(
methodtype),
pointer,
intent(inout) :: submethod
87 select type (subcell => this%subcell)
89 call this%load_subcell(particle, subcell)
95 subcell=this%subcell, &
96 trackctl=this%trackctl, &
97 tracktimes=this%tracktimes)
108 integer(I4B) :: exitFace
109 integer(I4B) :: inface
111 exitface = particle%iboundary(3)
112 isc = particle%idomain(3)
114 select case (exitface)
121 if (inface .eq. 0) inface = this%nverts
125 if (isc .gt. this%nverts) isc = 1
126 particle%idomain(3) = isc
127 particle%iboundary(3) = 3
132 if (isc .lt. 1) isc = this%nverts
133 particle%idomain(3) = isc
134 particle%iboundary(3) = 2
138 inface = this%nverts + 2
141 inface = this%nverts + 3
144 if (inface .eq. -1)
then
145 particle%iboundary(2) = 0
146 else if (inface .eq. 0)
then
147 particle%iboundary(2) = 0
149 particle%iboundary(2) = inface
159 real(DP),
intent(in) :: tmax
162 real(DP) :: x, y, z, xO, yO
163 real(DP),
allocatable :: xs(:), ys(:)
166 select type (cell => this%cell)
169 call this%check(particle, this%cell%defn, tmax)
170 if (.not. particle%advancing)
return
173 this%nverts = cell%defn%npolyverts
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)
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))
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))
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)
221 this%ztop = cell%defn%top
222 this%zbot = cell%defn%bot
223 this%dz = this%ztop - this%zbot
226 do i = 1, this%nverts
229 this%iprev = cshift(this%iprev, -1)
230 this%xvertnext = cshift(this%xvert, 1)
231 this%yvertnext = cshift(this%yvert, 1)
237 call particle%transform(xo, yo)
240 call this%track(particle, 2, tmax)
243 call particle%transform(xo, yo, invert=.true.)
244 call particle%reset_transform()
280 select type (cell => this%cell)
284 isc = particle%idomain(3)
288 do iv0 = 1, this%nverts
291 x1 = this%xvertnext(iv0)
292 y1 = this%yvertnext(iv0)
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
307 if ((alphai .ge.
dzero) .and. &
308 (alphai + betai .le.
done))
then
314 print *,
"error -- initial triangle not found in cell ", ic, &
315 " for particle at ", particle%x, particle%y, particle%z
319 particle%idomain(3) = isc
322 subcell%isubcell = 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)
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
341 subcell%vzbot = this%vzbot
342 subcell%vztop = this%vztop
355 real(DP),
allocatable,
dimension(:) :: xvals
356 real(DP),
allocatable,
dimension(:) :: yvals
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
371 real(DP),
allocatable,
dimension(:) :: areasub
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
388 integer(I4B) :: nvert
390 select type (cell => this%cell)
394 allocate (le(this%nverts))
395 allocate (unextx(this%nverts))
396 allocate (unexty(this%nverts))
397 allocate (areasub(this%nverts))
398 allocate (li(this%nverts))
399 allocate (unintx(this%nverts))
400 allocate (uninty(this%nverts))
401 allocate (xmid(this%nverts))
402 allocate (ymid(this%nverts))
403 allocate (lm(this%nverts))
404 allocate (umx(this%nverts))
405 allocate (umy(this%nverts))
406 allocate (kappax(this%nverts))
407 allocate (kappay(this%nverts))
408 allocate (vm0x(this%nverts))
409 allocate (vm0y(this%nverts))
410 allocate (vm1x(this%nverts))
411 allocate (vm1y(this%nverts))
412 allocate (unextxnext(this%nverts))
413 allocate (unextynext(this%nverts))
414 allocate (wk1(this%nverts))
415 allocate (wk2(this%nverts))
420 wk1 = this%xvertnext - this%xvert
421 wk2 = this%yvertnext - this%yvert
422 le = dsqrt(wk1 * wk1 + wk2 * wk2)
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
438 do i = 1, this%nverts
439 xvals(1) = this%xvert(i)
440 xvals(2) = this%xvertnext(i)
442 yvals(1) = this%yvert(i)
443 yvals(2) = this%yvertnext(i)
445 areasub(i) = area(xvals, yvals)
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)
455 divcell = sum(le * this%vne) / areacell
458 wk1 = this%xvert - this%xctr
459 wk2 = this%yvert - this%yctr
460 li = dsqrt(wk1 * wk1 + wk2 * wk2)
464 unextxnext = cshift(unintx, 1)
465 unextynext = cshift(uninty, 1)
468 xmid = 5.d-1 * (this%xvert + this%xctr)
469 ymid = 5.d-1 * (this%yvert + this%yctr)
472 wk1 = cshift(xmid, 1) - xmid
473 wk2 = cshift(ymid, 1) - ymid
474 lm = dsqrt(wk1 * wk1 + wk2 * wk2)
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)
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)
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)
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
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
541 deallocate (unextxnext)
542 deallocate (unextynext)
552 le, li, lm, areasub, areacell, &
553 unintx, uninty, unextx, unexty, &
554 unintxnext, unintynext, &
555 kappax, kappay, vm0x, vm0y, vm1x, vm1y, hcsum)
561 real(DP),
dimension(:) :: le
562 real(DP),
dimension(:) :: li
563 real(DP),
dimension(:) :: lm
564 real(DP),
dimension(:) :: areasub
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
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
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))
629 do i = 2, this%nverts
631 vm0i(i) = (li(ip) * vm0i(ip) - le(ip) * this%vne(ip) &
632 + areasub(ip) * divcell) / li(i)
636 vm1i = cshift(vm0i, 1)
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
665 vm0e = 5.d-1 * (this%vne + unextx * this%vctrx + unexty * this%vctry)
671 vm0x = bi0x * vm0i + be0x * vm0e
672 vm0y = bi0y * vm0i + be0y * vm0e
673 vm1x = bi1x * vm1i + be1x * vm0e
674 vm1y = bi1y * vm1i + be1y * vm0e
678 hcsum = sum(lm * (kappax * (vm0x + vm1x) + kappay * (vm0y + vm1y)))
subroutine, public create_cell_poly(cell)
Create a new polygonal cell.
This module contains simulation constants.
real(dp), parameter dhalf
real constant 1/2
real(dp), parameter dzero
real constant zero
real(dp), parameter dtwo
real constant 2
real(dp), parameter dsix
real constant 6
real(dp), parameter done
real constant 1
subroutine pstop(status, message)
Stop the program, optionally specifying an error status code.
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.
real(dp) function, public area(xv, yv, cw)
Calculate polygon area, with vertices given in CW or CCW order.
logical function, public point_in_polygon(x, y, poly)
Check if a point is within a polygon.
This module defines variable data types.
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.
Subcell-level tracking methods.
type(methodsubcellternarytype), pointer, public method_subcell_tern
subroutine, public create_subcell_tri(subcell)
Create a new triangular subcell.
Base type for particle tracking methods.
Particle tracked by the PRT model.