20 real(dp),
allocatable :: x_vert(:)
21 real(dp),
allocatable :: y_vert(:)
22 real(dp),
allocatable :: vx_vert_polygon(:)
23 real(dp),
allocatable :: vy_vert_polygon(:)
33 integer(I4B),
public,
pointer :: zeromethod
53 allocate (method%zeromethod)
56 method%type => method%cell%type
57 method%delegates = .true.
59 method%subcell => subcell
66 deallocate (this%type)
70 subroutine load_mct(this, particle, next_level, submethod)
74 integer(I4B),
intent(in) :: next_level
75 class(
methodtype),
pointer,
intent(inout) :: submethod
77 select type (subcell => this%subcell)
79 call this%load_subcell(particle, subcell)
83 subcell=this%subcell, &
84 trackfilectl=this%trackfilectl, &
85 tracktimes=this%tracktimes)
97 integer(I4B) :: exitFace
98 integer(I4B) :: inface
99 integer(I4B) :: npolyverts
101 exitface = particle%iboundary(3)
102 isc = particle%idomain(3)
103 select type (cell => this%cell)
105 npolyverts = cell%defn%npolyverts
108 select case (exitface)
115 if (inface .eq. 0) inface = npolyverts
119 if (isc .gt. npolyverts) isc = 1
120 particle%idomain(3) = isc
121 particle%iboundary(3) = 3
126 if (isc .lt. 1) isc = npolyverts
127 particle%idomain(3) = isc
128 particle%iboundary(3) = 2
132 inface = npolyverts + 2
135 inface = npolyverts + 3
137 if (inface .eq. -1)
then
138 particle%iboundary(2) = 0
139 else if (inface .eq. 0)
then
140 particle%iboundary(2) = 0
142 particle%iboundary(2) = inface
152 real(DP),
intent(in) :: tmax
154 integer(I4B) :: npolyverts
158 real(DP) :: retfactor
181 select type (cell => this%cell)
185 call this%update(particle, cell%defn)
188 if (.not. particle%advancing)
return
193 if (particle%z > cell%defn%top)
then
194 particle%z = cell%defn%top
195 call this%save(particle, reason=1)
198 npolyverts = cell%defn%npolyverts
199 if (
allocated(this%x_vert))
then
200 deallocate (this%x_vert)
201 deallocate (this%y_vert)
202 deallocate (this%vx_vert_polygon)
203 deallocate (this%vy_vert_polygon)
205 allocate (this%x_vert(npolyverts))
206 allocate (this%y_vert(npolyverts))
207 allocate (this%vx_vert_polygon(npolyverts))
208 allocate (this%vy_vert_polygon(npolyverts))
215 this%ztop = cell%defn%top
216 this%zbot = cell%defn%bot
217 this%dz = this%ztop - this%zbot
218 do iv = 1, npolyverts
220 if (ivp1 .gt. npolyverts) ivp1 = 1
222 if (ivm1 .lt. 1) ivm1 = npolyverts
223 x0 = cell%defn%polyvert(1, iv)
224 y0 = cell%defn%polyvert(2, iv)
225 x2 = cell%defn%polyvert(1, ivp1)
226 y2 = cell%defn%polyvert(2, ivp1)
227 x1 = cell%defn%polyvert(1, ivm1)
228 y1 = cell%defn%polyvert(2, ivm1)
229 term =
done / (cell%defn%porosity * this%dz)
230 flow0 = cell%defn%faceflow(iv) * term
231 flow1 = cell%defn%faceflow(ivm1) * term
242 det = d01y * d02x - d02y * d01x
243 retfactor = cell%defn%retfactor
247 term =
done / (retfactor * det)
249 v0x = -term * (d02x * flow1 + d01x * flow0)
251 v0y = -term * (d02y * flow1 + d01y * flow0)
252 this%vx_vert_polygon(iv) = v0x
253 this%vy_vert_polygon(iv) = v0y
260 area = area + x0 * y1 - x1 * y0
263 term =
done / (retfactor * cell%defn%porosity * area)
264 this%vzbot = cell%defn%faceflow(npolyverts + 2) * term
265 this%vztop = -cell%defn%faceflow(npolyverts + 3) * term
266 this%xctr = xsum / dble(npolyverts)
267 this%yctr = ysum / dble(npolyverts)
268 this%vxctr = vxsum / dble(npolyverts)
269 this%vyctr = vysum / dble(npolyverts)
272 call this%track(particle, 2, tmax)
287 integer(I4B) :: npolyverts
314 select type (cell => this%cell)
318 isc = particle%idomain(3)
319 npolyverts = cell%defn%npolyverts
323 do iv = 1, npolyverts
326 if (iv1 .gt. npolyverts) iv1 = 1
329 x0 = this%x_vert(ipv0)
330 y0 = this%y_vert(ipv0)
331 x1 = this%x_vert(ipv1)
332 y1 = this%y_vert(ipv1)
339 di2 = xi * y2rel - yi * x2rel
340 d02 = x0 * y2rel - y0 * x2rel
341 d12 = x1rel * y2rel - y1rel * x2rel
342 di1 = xi * y1rel - yi * x1rel
343 d01 = x0 * y1rel - y0 * x1rel
344 alphai = (di2 - d02) / d12
345 betai = -(di1 - d01) / d12
349 if ((alphai .ge. 0d0) .and. &
350 (betai .ge. betatol) .and. &
351 (alphai + betai .le. 1d0))
then
357 print *,
"error -- initial triangle not found for particle ", &
364 particle%idomain(3) = isc
367 subcell%isubcell = isc
372 if (iv1 .gt. npolyverts) iv1 = 1
375 subcell%x0 = this%x_vert(ipv0)
376 subcell%y0 = this%y_vert(ipv0)
377 subcell%x1 = this%x_vert(ipv1)
378 subcell%y1 = this%y_vert(ipv1)
379 subcell%x2 = this%xctr
380 subcell%y2 = this%yctr
381 subcell%v0x = this%vx_vert_polygon(iv0)
382 subcell%v0y = this%vy_vert_polygon(iv0)
383 subcell%v1x = this%vx_vert_polygon(iv1)
384 subcell%v1y = this%vy_vert_polygon(iv1)
385 subcell%v2x = this%vxctr
386 subcell%v2y = this%vyctr
387 subcell%ztop = this%ztop
388 subcell%zbot = this%zbot
390 subcell%vzbot = this%vzbot
391 subcell%vztop = this%vztop
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 done
real constant 1
subroutine pstop(status, message)
Stop the program, optionally specifying an error status code.
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 load_subcell(this, particle, subcell)
Loads a triangular subcell from the polygonal cell.
procedure subroutine, public create_method_cell_ternary(method)
Create a tracking method.
subroutine pass_mct(this, particle)
Pass particle to next subcell if there is one, or to the cell face.
Particle tracking strategies.
subroutine load(this, particle, next_level, submethod)
Load subdomain tracking method (submethod)
subroutine pass(this, particle)
Pass a particle to the next subdomain, internal use only.
Subcell-level tracking methods.
type(methodsubcellternarytype), pointer, public method_subcell_tern
pure character(len=lenmempath) function, public get_particle_id(particle)
Return the particle's composite ID.
subroutine, public create_subcell_tri(subcell)
Create a new triangular subcell.
Base type for particle tracking methods.
A particle tracked by the PRT model.
Manages particle track (i.e. pathline) files.