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, &
97 tracktimes=this%tracktimes)
108 integer(I4B) :: exitFace
109 integer(I4B) :: inface
114 select case (exitface)
121 if (inface .eq. 0) inface = this%nverts
125 if (isc .gt. this%nverts) isc = 1
132 if (isc .lt. 1) isc = this%nverts
138 inface = this%nverts + 2
141 inface = this%nverts + 3
144 if (inface .eq. -1)
then
146 else if (inface .eq. 0)
then
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)
168 call this%assess(particle, this%cell%defn, tmax)
169 if (.not. particle%advancing)
return
172 this%nverts = cell%defn%npolyverts
175 if (
allocated(this%xvert))
then
176 deallocate (this%xvert)
177 deallocate (this%yvert)
178 deallocate (this%vne)
179 deallocate (this%vv0x)
180 deallocate (this%vv0y)
181 deallocate (this%vv1x)
182 deallocate (this%vv1y)
183 deallocate (this%iprev)
184 deallocate (this%xvertnext)
185 deallocate (this%yvertnext)
188 allocate (this%xvert(this%nverts))
189 allocate (this%yvert(this%nverts))
190 allocate (this%vne(this%nverts))
191 allocate (this%vv0x(this%nverts))
192 allocate (this%vv0y(this%nverts))
193 allocate (this%vv1x(this%nverts))
194 allocate (this%vv1y(this%nverts))
195 allocate (this%iprev(this%nverts))
196 allocate (this%xvertnext(this%nverts))
197 allocate (this%yvertnext(this%nverts))
201 allocate (xs(this%nverts))
202 allocate (ys(this%nverts))
203 xs = cell%defn%polyvert(1, :)
204 ys = cell%defn%polyvert(2, :)
205 xo = xs(minloc(abs(xs), dim=1))
206 yo = ys(minloc(abs(ys), dim=1))
211 do i = 1, this%nverts
212 x = cell%defn%polyvert(1, i)
213 y = cell%defn%polyvert(2, i)
214 call transform(x, y,
dzero, x, y, z, xo, yo)
220 this%ztop = cell%defn%top
221 this%zbot = cell%defn%bot
222 this%dz = this%ztop - this%zbot
225 do i = 1, this%nverts
228 this%iprev = cshift(this%iprev, -1)
229 this%xvertnext = cshift(this%xvert, 1)
230 this%yvertnext = cshift(this%yvert, 1)
236 call particle%transform(xo, yo)
239 call this%track(particle, 2, tmax)
242 call particle%transform(xo, yo, invert=.true.)
243 call particle%reset_transform()
254 integer(I4B) :: ic, isc, iv0
255 real(DP) :: x0, y0, x1, y1, x2, y2, xi, yi
256 real(DP) :: x1rel, y1rel, x2rel, y2rel
257 real(DP) :: di2, d02, d12, di1, d01
258 real(DP) :: alphai, betai
260 select type (cell => this%cell)
268 do iv0 = 1, this%nverts
271 x1 = this%xvertnext(iv0)
272 y1 = this%yvertnext(iv0)
279 di2 = xi * y2rel - yi * x2rel
280 d02 = x0 * y2rel - y0 * x2rel
281 d12 = x1rel * y2rel - y1rel * x2rel
282 di1 = xi * y1rel - yi * x1rel
283 d01 = x0 * y1rel - y0 * x1rel
284 alphai = (di2 - d02) / d12
285 betai = -(di1 - d01) / d12
287 if ((alphai .ge.
dzero) .and. &
288 (alphai + betai .le.
done))
then
294 print *,
"error -- initial triangle not found in cell ", ic, &
295 " for particle at ", particle%x, particle%y, particle%z
302 subcell%isubcell = isc
306 subcell%x0 = this%xvert(iv0)
307 subcell%y0 = this%yvert(iv0)
308 subcell%x1 = this%xvertnext(iv0)
309 subcell%y1 = this%yvertnext(iv0)
310 subcell%x2 = this%xctr
311 subcell%y2 = this%yctr
312 subcell%v0x = this%vv0x(iv0)
313 subcell%v0y = this%vv0y(iv0)
314 subcell%v1x = this%vv1x(iv0)
315 subcell%v1y = this%vv1y(iv0)
316 subcell%v2x = this%vctrx
317 subcell%v2y = this%vctry
318 subcell%ztop = this%ztop
319 subcell%zbot = this%zbot
321 subcell%vzbot = this%vzbot
322 subcell%vztop = this%vztop
335 real(DP),
allocatable,
dimension(:) :: xvals
336 real(DP),
allocatable,
dimension(:) :: yvals
343 real(DP),
allocatable,
dimension(:) :: wk1
344 real(DP),
allocatable,
dimension(:) :: wk2
345 real(DP),
allocatable,
dimension(:) :: unextxnext
346 real(DP),
allocatable,
dimension(:) :: unextynext
347 real(DP),
allocatable,
dimension(:) :: le
348 real(DP),
allocatable,
dimension(:) :: unextx
349 real(DP),
allocatable,
dimension(:) :: unexty
351 real(DP),
allocatable,
dimension(:) :: areasub
353 real(DP),
allocatable,
dimension(:) :: li
354 real(DP),
allocatable,
dimension(:) :: unintx
355 real(DP),
allocatable,
dimension(:) :: uninty
356 real(DP),
allocatable,
dimension(:) :: xmid
357 real(DP),
allocatable,
dimension(:) :: ymid
358 real(DP),
allocatable,
dimension(:) :: lm
359 real(DP),
allocatable,
dimension(:) :: umx
360 real(DP),
allocatable,
dimension(:) :: umy
361 real(DP),
allocatable,
dimension(:) :: kappax
362 real(DP),
allocatable,
dimension(:) :: kappay
363 real(DP),
allocatable,
dimension(:) :: vm0x
364 real(DP),
allocatable,
dimension(:) :: vm0y
365 real(DP),
allocatable,
dimension(:) :: vm1x
366 real(DP),
allocatable,
dimension(:) :: vm1y
368 integer(I4B) :: nvert
370 select type (cell => this%cell)
374 allocate (le(this%nverts))
375 allocate (unextx(this%nverts))
376 allocate (unexty(this%nverts))
377 allocate (areasub(this%nverts))
378 allocate (li(this%nverts))
379 allocate (unintx(this%nverts))
380 allocate (uninty(this%nverts))
381 allocate (xmid(this%nverts))
382 allocate (ymid(this%nverts))
383 allocate (lm(this%nverts))
384 allocate (umx(this%nverts))
385 allocate (umy(this%nverts))
386 allocate (kappax(this%nverts))
387 allocate (kappay(this%nverts))
388 allocate (vm0x(this%nverts))
389 allocate (vm0y(this%nverts))
390 allocate (vm1x(this%nverts))
391 allocate (vm1y(this%nverts))
392 allocate (unextxnext(this%nverts))
393 allocate (unextynext(this%nverts))
394 allocate (wk1(this%nverts))
395 allocate (wk2(this%nverts))
400 wk1 = this%xvertnext - this%xvert
401 wk2 = this%yvertnext - this%yvert
402 le = dsqrt(wk1 * wk1 + wk2 * wk2)
407 areacell = area(this%xvert, this%yvert)
408 sixa = areacell * dsix
409 wk1 = -(this%xvert * this%yvertnext - this%xvertnext * this%yvert)
410 nvert =
size(this%xvert)
411 this%xctr = sum((this%xvert + this%xvertnext) * wk1) / sixa
412 this%yctr = sum((this%yvert + this%yvertnext) * wk1) / sixa
418 do i = 1, this%nverts
419 xvals(1) = this%xvert(i)
420 xvals(2) = this%xvertnext(i)
422 yvals(1) = this%yvert(i)
423 yvals(2) = this%yvertnext(i)
425 areasub(i) = area(xvals, yvals)
429 term =
done / (cell%defn%porosity * cell%defn%retfactor * this%dz)
430 do i = 1, this%nverts
431 this%vne(i) = -cell%defn%faceflow(i) * term / le(i)
435 divcell = sum(le * this%vne) / areacell
438 wk1 = this%xvert - this%xctr
439 wk2 = this%yvert - this%yctr
440 li = dsqrt(wk1 * wk1 + wk2 * wk2)
444 unextxnext = cshift(unintx, 1)
445 unextynext = cshift(uninty, 1)
448 xmid = 5.d-1 * (this%xvert + this%xctr)
449 ymid = 5.d-1 * (this%yvert + this%yctr)
452 wk1 = cshift(xmid, 1) - xmid
453 wk2 = cshift(ymid, 1) - ymid
454 lm = dsqrt(wk1 * wk1 + wk2 * wk2)
471 call this%calc_thru_hcsum(vm0i0, divcell, le, li, lm, areasub, areacell, &
472 unintx, uninty, unextx, unexty, &
473 unextxnext, unextynext, &
474 kappax, kappay, vm0x, vm0y, vm1x, vm1y, hcsum0)
476 vm0ival = vm0i0 + perturb
477 call this%calc_thru_hcsum(vm0ival, divcell, le, li, lm, areasub, areacell, &
478 unintx, uninty, unextx, unexty, &
479 unextxnext, unextynext, &
480 kappax, kappay, vm0x, vm0y, vm1x, vm1y, hcsum)
482 jac = (hcsum - hcsum0) / perturb
483 vm0ival = vm0i0 - hcsum0 / jac
484 call this%calc_thru_hcsum(vm0ival, divcell, le, li, lm, areasub, areacell, &
485 unintx, uninty, unextx, unexty, &
486 unextxnext, unextynext, &
487 kappax, kappay, vm0x, vm0y, vm1x, vm1y, hcsum)
492 this%vv0x = 2.d0 * vm0x - this%vctrx
493 this%vv0y = 2.d0 * vm0y - this%vctry
494 this%vv1x = 2.d0 * vm1x - this%vctrx
495 this%vv1y = 2.d0 * vm1y - this%vctry
498 term =
done / (cell%defn%retfactor * cell%defn%porosity * areacell)
499 this%vzbot = cell%defn%faceflow(this%nverts + 2) * term
500 this%vztop = -cell%defn%faceflow(this%nverts + 3) * term
521 deallocate (unextxnext)
522 deallocate (unextynext)
532 le, li, lm, areasub, areacell, &
533 unintx, uninty, unextx, unexty, &
534 unintxnext, unintynext, &
535 kappax, kappay, vm0x, vm0y, vm1x, vm1y, hcsum)
541 real(DP),
dimension(:) :: le
542 real(DP),
dimension(:) :: li
543 real(DP),
dimension(:) :: lm
544 real(DP),
dimension(:) :: areasub
546 real(DP),
dimension(:) :: unintx
547 real(DP),
dimension(:) :: uninty
548 real(DP),
dimension(:) :: unextx
549 real(DP),
dimension(:) :: unexty
550 real(DP),
dimension(:) :: unintxnext
551 real(DP),
dimension(:) :: unintynext
552 real(DP),
dimension(:) :: kappax
553 real(DP),
dimension(:) :: kappay
554 real(DP),
dimension(:) :: vm0x
555 real(DP),
dimension(:) :: vm0y
556 real(DP),
dimension(:) :: vm1x
557 real(DP),
dimension(:) :: vm1y
559 real(DP),
allocatable,
dimension(:) :: vm0i
560 real(DP),
allocatable,
dimension(:) :: vm0e
561 real(DP),
allocatable,
dimension(:) :: vm1i
562 real(DP),
allocatable,
dimension(:) :: vm1e
563 real(DP),
allocatable,
dimension(:) :: uprod
564 real(DP),
allocatable,
dimension(:) :: det
565 real(DP),
allocatable,
dimension(:) :: wt
566 real(DP),
allocatable,
dimension(:) :: bi0x
567 real(DP),
allocatable,
dimension(:) :: be0x
568 real(DP),
allocatable,
dimension(:) :: bi0y
569 real(DP),
allocatable,
dimension(:) :: be0y
570 real(DP),
allocatable,
dimension(:) :: bi1x
571 real(DP),
allocatable,
dimension(:) :: be1x
572 real(DP),
allocatable,
dimension(:) :: bi1y
573 real(DP),
allocatable,
dimension(:) :: be1y
574 real(DP),
allocatable,
dimension(:) :: be01x
575 real(DP),
allocatable,
dimension(:) :: be01y
587 allocate (vm0i(this%nverts))
588 allocate (vm0e(this%nverts))
589 allocate (vm1i(this%nverts))
590 allocate (vm1e(this%nverts))
591 allocate (uprod(this%nverts))
592 allocate (det(this%nverts))
593 allocate (wt(this%nverts))
594 allocate (bi0x(this%nverts))
595 allocate (be0x(this%nverts))
596 allocate (bi0y(this%nverts))
597 allocate (be0y(this%nverts))
598 allocate (bi1x(this%nverts))
599 allocate (be1x(this%nverts))
600 allocate (bi1y(this%nverts))
601 allocate (be1y(this%nverts))
602 allocate (be01x(this%nverts))
603 allocate (be01y(this%nverts))
609 do i = 2, this%nverts
611 vm0i(i) = (li(ip) * vm0i(ip) - le(ip) * this%vne(ip) &
612 + areasub(ip) * divcell) / li(i)
616 vm1i = cshift(vm0i, 1)
619 uprod = unintx * unextx + uninty * unexty
620 det =
done - uprod * uprod
621 bi0x = (unintx - unextx * uprod) / det
622 be0x = (unextx - unintx * uprod) / det
623 bi0y = (uninty - unexty * uprod) / det
624 be0y = (unexty - uninty * uprod) / det
625 uprod = unintxnext * unextx + unintynext * unexty
626 det =
done - uprod * uprod
627 bi1x = (unintxnext - unextx * uprod) / det
628 be1x = (unextx - unintxnext * uprod) / det
629 bi1y = (unintynext - unexty * uprod) / det
630 be1y = (unexty - unintynext * uprod) / det
631 be01x = 5.d-1 * (be0x + be1x)
632 be01y = 5.d-1 * (be0y + be1y)
633 wt = areasub / areacell
634 emxx =
dtwo - sum(wt * be01x * unextx)
635 emxy = -sum(wt * be01x * unexty)
636 emyx = -sum(wt * be01y * unextx)
637 emyy =
dtwo - sum(wt * be01y * unexty)
638 rx = sum(wt * (bi0x * vm0i + bi1x * vm1i + be01x * this%vne))
639 ry = sum(wt * (bi0y * vm0i + bi1y * vm1i + be01y * this%vne))
640 emdet = emxx * emyy - emxy * emyx
641 this%vctrx = (emyy * rx - emxy * ry) / emdet
642 this%vctry = (emxx * ry - emyx * rx) / emdet
645 vm0e = 5.d-1 * (this%vne + unextx * this%vctrx + unexty * this%vctry)
651 vm0x = bi0x * vm0i + be0x * vm0e
652 vm0y = bi0y * vm0i + be0y * vm0e
653 vm1x = bi1x * vm1i + be1x * vm0e
654 vm1y = bi1y * vm1i + be1y * vm0e
658 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.
@, public level_subfeature
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.