22 integer(I4B),
public,
pointer :: zeromethod
39 allocate (method%zeromethod)
41 method%subcell => subcell
42 method%type => method%subcell%type
43 method%delegates = .false.
50 deallocate (this%type)
57 real(DP),
intent(in) :: tmax
59 select type (subcell => this%subcell)
61 call this%track_subcell(subcell, particle, tmax)
71 real(DP),
intent(in) :: tmax
73 integer(I4B) :: exitFace
101 real(DP) :: rot(2, 2), res(2), loc(2)
126 integer(I4B) :: izstatus
127 integer(I4B) :: itopbotexit
128 integer(I4B) :: ntmax
129 integer(I4B) :: nsave
130 integer(I4B) :: isolv
131 integer(I4B) :: itrifaceenter
132 integer(I4B) :: itrifaceexit
138 integer(I4B) :: ntdebug
139 integer(I4B) :: reason
141 integer(I4B) :: tslice(2)
146 isolv = this%zeromethod
170 vzbot = subcell%vzbot
171 vztop = subcell%vztop
175 v0x, v0y, v1x, v1y, v2x, v2y, &
177 rxx, rxy, ryx, ryy, &
179 alp0, bet0, alp1, bet1, alp2, bet2, alpi, beti, &
186 zirel = (zi - zbot) / dz
188 az, dtexitz, izstatus, &
194 itrifaceenter = particle%iboundary(3) - 1
195 if (itrifaceenter .eq. -1) itrifaceenter = 999
201 dtexitxy, alpexit, betexit, &
202 itrifaceenter, itrifaceexit, &
203 rxx, rxy, ryx, ryy, &
204 alp0, bet0, alp1, bet1, alp2, bet2, alpi, beti, &
208 if ((itopbotexit .eq. 0) .and. (itrifaceexit .eq. 0))
then
209 call this%save(particle, reason=9)
214 if (itopbotexit .eq. 0)
then
216 exitface = itrifaceexit
218 else if (itrifaceexit .eq. 0)
then
222 else if (dtexitz .lt. dtexitxy)
then
228 exitface = itrifaceexit
231 if (exitface .eq. 45)
then
232 if (itopbotexit .eq. -1)
then
240 texit = particle%ttrack + dtexit
247 call this%tracktimes%try_advance()
248 tslice = this%tracktimes%selection
249 if (all(tslice > 0))
then
250 do i = tslice(1), tslice(2)
251 t = this%tracktimes%times(i)
252 if (t < particle%ttrack .or. t >= texit .or. t >= tmax) cycle
256 if (lbary) loc =
skew(loc, (/sxx, sxy, syy/), invert=.true.)
257 rot = reshape((/rxx, rxy, ryx, ryy/), shape(rot))
258 res = matmul(rot, loc)
262 if (izstatus .eq. 2)
then
265 else if (izstatus .eq. 1)
then
270 z = zbot + (vzi * dexp(az * dt) - vzbot) / az
277 call this%save(particle, reason=5)
281 if (texit .gt. tmax)
then
288 particle%advancing = .false.
302 if (exitface .eq. 1)
then
304 else if (exitface .eq. 2)
then
306 else if (exitface .eq. 3)
then
310 if (lbary) loc =
skew(loc, (/sxx, sxy, syy/), invert=.true.)
311 rot = reshape((/rxx, rxy, ryx, ryy/), shape(rot))
312 res = matmul(rot, loc)
315 if (exitface .eq. 4)
then
317 else if (exitface .eq. 5)
then
320 if (izstatus .eq. 2)
then
323 else if (izstatus .eq. 1)
then
328 z = zbot + (vzi * dexp(az * dt) - vzbot) / az
338 particle%iboundary(3) = exitface
342 call this%save(particle, reason=reason)
353 dt, status, itopbotexit)
375 integer(I4B) :: status
376 integer(I4B) :: itopbotexit
377 logical(LGP) :: noOutflow
383 if (v2a .lt. 0d0) v2a = -v2a
385 if (v1a .lt. 0d0) v1a = -v1a
388 if (dva .lt. 0d0) dva = -dva
393 if ((v2a .lt. tol) .and. (v1a .lt. tol))
then
405 if (v2a .gt. vv) vv = v2a
407 if (vvv .lt. 1.0d-4)
then
412 if (v1 .gt. zro)
then
416 if (v1 .lt. zrom)
then
428 v = (1.0d0 - xl) * v1 + xl * v2
433 if (v1 .lt. 0d0) nooutflow = .false.
434 if (v2 .gt. 0d0) nooutflow = .false.
445 if ((v1 .le. 0d0) .and. (v2 .ge. 0d0))
then
446 if (abs(v) .le. 0d0)
then
448 if (v2 .le. 0d0) v = -v
458 if (vr .le. 0d0)
then
467 if (v1v2 .gt. 0d0)
then
subroutine pstop(status, message)
Stop the program, optionally specifying an error status code.
pure real(dp) function, dimension(2), public skew(v, s, invert)
Skew a 2D vector along the x-axis.
This module defines variable data types.
Particle tracking strategies.
subroutine track_subcell(this, subcell, particle, tmax)
Track a particle across a triangular subcell using the ternary method.
subroutine apply_mst(this, particle, tmax)
Apply the ternary subcell method.
subroutine, public create_method_subcell_ternary(method)
Create a new ternary subcell-method object.
subroutine calculate_dt(v1, v2, dx, xL, v, dvdx, dt, status, itopbotexit)
Do calculations related to analytical z solution.
subroutine, public create_subcell_tri(subcell)
Create a new triangular subcell.
subroutine, public traverse_triangle(isolv, tol, step, texit, alpexit, betexit, itrifaceenter, itrifaceexit, rxx, rxy, ryx, ryy, alp0, bet0, alp1, bet1, alp2, bet2, alpi, beti, vziodz, az, bary)
Traverse triangular cell.
subroutine, public canonical(x0, y0, x1, y1, x2, y2, v0x, v0y, v1x, v1y, v2x, v2y, xi, yi, rxx, rxy, ryx, ryy, sxx, sxy, syy, alp0, bet0, alp1, bet1, alp2, bet2, alpi, beti, bary)
Set coordinates to "canonical" configuration.
subroutine, public step_analytical(t, alp, bet)
Step (evaluate) analytically depending on case.
Base type for grid cells of a concrete type. Contains a cell-definition which is information shared b...
Base type for particle tracking methods.
Ternary triangular subcell tracking method.
A particle tracked by the PRT model.
Manages particle track (i.e. pathline) files.