23 integer(I4B),
public,
pointer :: zeromethod
41 method%subcell => subcell
42 method%name => method%subcell%type
43 method%delegates = .false.
49 deallocate (this%name)
56 real(DP),
intent(in) :: tmax
58 select type (subcell => this%subcell)
60 call this%track_subcell(subcell, particle, tmax)
71 real(DP),
intent(in) :: tmax
73 integer(I4B) :: exitFace
123 integer(I4B) :: izstatus
124 integer(I4B) :: itopbotexit
125 integer(I4B) :: ntmax
126 integer(I4B) :: isolv
127 integer(I4B) :: itrifaceenter
128 integer(I4B) :: itrifaceexit
133 integer(I4B) :: reason
137 if (particle%iexmeth == 0)
then
140 isolv = particle%iexmeth
166 vzbot = subcell%vzbot
167 vztop = subcell%vztop
180 v0x, v0y, v1x, v1y, v2x, v2y, &
182 rxx, rxy, ryx, ryy, &
184 alp0, bet0, alp1, bet1, alp2, bet2, alpi, beti)
193 zirel = (zi - zbot) / dz
194 if (zirel >
done)
then
196 else if (zirel <
dzero)
then
200 az, dtexitz, izstatus, &
205 itrifaceenter = particle%iboundary(3) - 1
206 if (itrifaceenter == -1) itrifaceenter = 999
208 dtexitxy, alpexit, betexit, &
209 itrifaceenter, itrifaceexit, &
210 alp1, bet1, alp2, bet2, alpi, beti)
214 if (itopbotexit == 0 .and. itrifaceexit == 0)
then
216 particle%advancing = .false.
217 call this%save(particle, reason=3)
222 if (itopbotexit == 0)
then
224 exitface = itrifaceexit
226 else if (itrifaceexit == 0 .or. dtexitz < dtexitxy)
then
232 exitface = itrifaceexit
235 if (exitface == 45)
then
236 if (itopbotexit == -1)
then
244 if (dtexit <
dzero)
then
246 particle%advancing = .false.
247 call this%save(particle, reason=3)
251 texit = particle%ttrack + dtexit
258 call this%tracktimes%advance()
259 if (this%tracktimes%any())
then
260 do i = this%tracktimes%selection(1), this%tracktimes%selection(2)
261 t = this%tracktimes%times(i)
262 if (t < particle%ttrack) cycle
263 if (t >= texit .or. t >= tmax)
exit
266 izstatus, x0, y0, az, vzi, vzbot, &
267 ztop, zbot, zi, x, y, z)
273 call this%save(particle, reason=5)
280 if (texit .gt. tmax)
then
287 particle%advancing = .false.
297 izstatus, x0, y0, az, vzi, vzbot, &
298 ztop, zbot, zi, x, y, z, exitface)
303 particle%iboundary(3) = exitface
305 call this%save(particle, reason=reason)
316 dt, status, itopbotexit)
338 integer(I4B) :: status
339 integer(I4B) :: itopbotexit
340 logical(LGP) :: noOutflow
346 if (v2a .lt.
dzero) v2a = -v2a
348 if (v1a .lt.
dzero) v1a = -v1a
351 if (dva .lt.
dzero) dva = -dva
356 if ((v2a .lt. tol) .and. (v1a .lt. tol))
then
368 if (v2a .gt. vv) vv = v2a
370 if (vvv .lt. 1.0d-4)
then
375 if (v1 .gt. zro)
then
379 if (v1 .lt. zrom)
then
391 v = (
done - xl) * v1 + xl * v2
396 if (v1 .lt.
dzero) nooutflow = .false.
397 if (v2 .gt.
dzero) nooutflow = .false.
408 if ((v1 .le.
dzero) .and. (v2 .ge.
dzero))
then
409 if (abs(v) .le.
dzero)
then
411 if (v2 .le.
dzero) v = -v
421 if (vr .le.
dzero)
then
430 if (v1v2 .gt.
dzero)
then
431 if (v .gt.
dzero)
then
435 if (v .lt.
dzero)
then
442 dt = log(abs(vr)) / dvdx
448 izstatus, x0, y0, az, vzi, vzbot, &
449 ztop, zbot, zi, x, y, z, exitFace)
459 integer(I4B) :: izstatus
471 integer(I4B),
optional :: exitFace
473 integer(I4B) :: lexitface
474 real(DP) :: rot(2, 2), res(2), loc(2)
479 if (
present(exitface))
then
490 if (lexitface .eq. 1)
then
492 else if (lexitface .eq. 2)
then
494 else if (lexitface .eq. 3)
then
499 if (lexitface .eq. 4)
then
501 else if (lexitface .eq. 5)
then
505 if (izstatus .eq. 2)
then
508 else if (izstatus .eq. 1)
then
513 z = zbot + (vzi * dexp(az * dt) - vzbot) / az
519 loc =
skew(loc, (/sxx, sxy, syy/), invert=.true.)
520 rot = reshape((/rxx, rxy, ryx, ryy/), shape(rot))
521 res = matmul(rot, loc)
This module contains simulation constants.
real(dp), parameter dsame
real constant for values that are considered the same based on machine precision
real(dp), parameter dep3
real constant 1000
real(dp), parameter donethird
real constant 1/3
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 done
real constant 1
subroutine pstop(status, message)
Stop the program, optionally specifying an error status code.
subroutine, public clamp_bary(alpha, beta, gamma, pad)
Clamp barycentric coordinates to the interior of a triangle, with optional padding some minimum dista...
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.
pure logical function, public is_close(a, b, rtol, atol, symmetric)
Check if a real value is approximately equal to another.
Particle tracking strategies.
subroutine track_subcell(this, subcell, particle, tmax)
Track a particle across a triangular subcell.
subroutine apply_mst(this, particle, tmax)
Apply the ternary subcell tracking method.
subroutine, public create_method_subcell_ternary(method)
Create a new ternary subcell tracking method.
subroutine calculate_xyz_position(dt, rxx, rxy, ryx, ryy, sxx, sxy, syy, izstatus, x0, y0, az, vzi, vzbot, ztop, zbot, zi, x, y, z, exitFace)
Calculate the particle's local unscaled xyz coordinates after dt.
subroutine calculate_dt(v1, v2, dx, xL, v, dvdx, dt, status, itopbotexit)
Do calculations related to analytical z solution.
@ term_no_exits_sub
terminated in a subcell with no exit face
subroutine, public create_subcell_tri(subcell)
Create a new triangular subcell.
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)
Set coordinates to "canonical" configuration.
subroutine, public traverse_triangle(isolv, tol, texit, alpexit, betexit, itrifaceenter, itrifaceexit, alp1, bet1, alp2, bet2, alpi, beti)
Traverse triangular cell.
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.
Particle tracked by the PRT model.