41 method%type => method%cell%type
42 method%delegates = .true.
44 method%subcell => subcell
50 deallocate (this%type)
54 subroutine load_mcpq(this, particle, next_level, submethod)
57 integer,
intent(in) :: next_level
58 class(
methodtype),
pointer,
intent(inout) :: submethod
60 select type (subcell => this%subcell)
62 call this%load_subcell(particle, subcell)
66 subcell=this%subcell, &
67 trackfilectl=this%trackfilectl, &
68 tracktimes=this%tracktimes)
78 integer :: isc, exitFace, npolyverts, inface, infaceoff
80 select type (cell => this%cell)
82 exitface = particle%iboundary(3)
83 isc = particle%idomain(3)
84 npolyverts = cell%defn%npolyverts
86 select case (exitface)
94 particle%idomain(3) = 4
95 particle%iboundary(3) = 2
99 particle%idomain(3) = 3
100 particle%iboundary(3) = 2
123 particle%idomain(3) = 2
124 particle%iboundary(3) = 1
128 particle%idomain(3) = 1
129 particle%iboundary(3) = 1
136 particle%idomain(3) = 2
137 particle%iboundary(3) = 4
149 particle%idomain(3) = 3
150 particle%iboundary(3) = 4
161 particle%idomain(3) = 1
162 particle%iboundary(3) = 3
166 particle%idomain(3) = 4
167 particle%iboundary(3) = 3
176 inface = npolyverts + 2
179 inface = npolyverts + 3
181 if (inface .eq. -1)
then
182 particle%iboundary(2) = 0
183 else if (inface .eq. 0)
then
184 particle%iboundary(2) = 0
186 if ((inface .ge. 1) .and. (inface .le. 4))
then
188 inface = inface + cell%irvOrigin - 1
189 if (inface .gt. 4) inface = inface - 4
190 inface = cell%irectvert(inface) + infaceoff
191 if (inface .lt. 1) inface = inface + npolyverts
193 particle%iboundary(2) = inface
203 real(DP),
intent(in) :: tmax
205 double precision :: xOrigin, yOrigin, zOrigin, sinrot, cosrot
207 select type (cell => this%cell)
210 call this%update(particle, cell%defn)
211 if (.not. particle%advancing)
return
217 if (particle%z > cell%defn%top)
then
218 particle%z = cell%defn%top
219 call this%save(particle, reason=1)
226 xorigin = cell%xOrigin
227 yorigin = cell%yOrigin
228 zorigin = cell%zOrigin
231 call particle%transform(xorigin, yorigin, zorigin, &
233 call this%track(particle, 2, tmax)
234 call particle%transform(xorigin, yorigin, zorigin, &
235 sinrot, cosrot, invert=.true.)
236 call particle%transform(reset=.true.)
247 double precision :: dx, dy, dz, areax, areay, areaz
248 double precision :: dxprel, dyprel
249 integer :: isc, npolyverts, m1, m2
250 double precision :: qextl1, qextl2, qintl1, qintl2
251 double precision :: factor, term
253 select type (cell => this%cell)
255 factor =
done / cell%defn%retfactor
256 factor = factor / cell%defn%porosity
257 npolyverts = cell%defn%npolyverts
259 isc = particle%idomain(3)
267 dxprel = particle%x / dx
268 dyprel = particle%y / dy
269 if (dxprel .lt. 5d-1)
then
270 if (dyprel .lt. 5d-1)
then
272 else if (dyprel .gt. 5d-1)
then
276 call pstop(1,
"particle initially on shared subcell edge")
278 else if (dxprel .gt. 5d-1)
then
279 if (dyprel .lt. 5d-1)
then
281 else if (dyprel .gt. 5d-1)
then
285 call pstop(1,
"particle initially on shared subcell edge")
289 call pstop(1,
"particle initially on shared subcell edge")
291 subcell%isubcell = isc
294 particle%idomain(3) = isc
299 dz = cell%defn%top - &
304 qintl1 = cell%qintl(isc)
306 qintl2 = cell%qintl(isc + 1)
307 qextl1 = cell%qextl1(isc)
308 qextl2 = cell%qextl2(isc)
315 subcell%zOrigin = 0d0
320 term = factor / areax
321 subcell%vx1 = qintl1 * term
322 subcell%vx2 = -qextl2 * term
323 term = factor / areay
324 subcell%vy1 = -qintl2 * term
325 subcell%vy2 = -qextl1 * term
328 subcell%yOrigin = 0d0
329 term = factor / areax
330 subcell%vx1 = -qintl2 * term
331 subcell%vx2 = -qextl1 * term
332 term = factor / areay
333 subcell%vy1 = qextl2 * term
334 subcell%vy2 = -qintl1 * term
336 subcell%xOrigin = 0d0
337 subcell%yOrigin = 0d0
338 term = factor / areax
339 subcell%vx1 = qextl2 * term
340 subcell%vx2 = -qintl1 * term
341 term = factor / areay
342 subcell%vy1 = qextl1 * term
343 subcell%vy2 = qintl2 * term
345 subcell%xOrigin = 0d0
347 term = factor / areax
348 subcell%vx1 = qextl1 * term
349 subcell%vx2 = qintl2 * term
350 term = factor / areay
351 subcell%vy1 = qintl1 * term
352 subcell%vy2 = -qextl2 * term
356 term = factor / areaz
357 subcell%vz1 = 2.5d-1 * cell%defn%faceflow(m1) * term
358 subcell%vz2 = -2.5d-1 * cell%defn%faceflow(m2) * term
subroutine, public create_cell_rect_quad(cell)
Create a new rectangular-quad cell.
This module contains simulation constants.
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.
procedure subroutine, public create_method_cell_quad(method)
Create a new tracking method.
subroutine load_mcpq(this, particle, next_level, submethod)
Load subcell into tracking method.
subroutine pass_mcpq(this, particle)
Pass particle to next subcell if there is one, or to the cell face.
subroutine load_subcell(this, particle, subcell)
Load the rectangular subcell from the rectangular cell.
subroutine apply_mcpq(this, particle, tmax)
Solve the quad-rectangular cell via Pollock's method.
Particle tracking strategies.
Subcell-level tracking methods.
type(methodsubcellpollocktype), pointer, public method_subcell_plck
subroutine, public create_subcell_rect(subcell)
Create a new rectangular subcell.
Base grid cell definition.
Base type for particle tracking methods.
A particle tracked by the PRT model.
Manages particle track (i.e. pathline) files.