28 real(dp),
allocatable,
public :: qextl1(:), qextl2(:), qintl(:)
49 method%subcell => subcell
50 method%name => method%subcell%type
51 method%delegates = .false.
57 deallocate (this%name)
65 real(DP),
intent(in) :: tmax
73 select type (subcell => this%subcell)
79 x_origin = subcell%xOrigin
80 y_origin = subcell%yOrigin
81 z_origin = subcell%zOrigin
82 sinrot = subcell%sinrot
83 cosrot = subcell%cosrot
84 call particle%transform(x_origin, y_origin)
85 call this%track_subcell(subcell, particle, tmax)
86 call particle%transform(x_origin, y_origin, invert=.true.)
105 real(DP),
intent(in) :: tmax
107 real(DP) :: dt, dtexit, texit
108 real(DP) :: t, x, y, z
109 real(DP) :: t0, x0, y0, z0
110 integer(I4B) :: i, exit_face, exit_soln
114 x0 = particle%x / subcell%dx
115 y0 = particle%y / subcell%dy
116 z0 = particle%z / subcell%dz
119 call this%find_exits(particle, subcell)
121 exit_x = this%exit_solutions(1)
122 exit_y = this%exit_solutions(2)
123 exit_z = this%exit_solutions(3)
145 exit_soln = this%pick_exit(particle)
146 if (exit_soln == 0)
then
150 exit_face = this%exit_solutions(exit_soln)%iboundary
151 dtexit = this%exit_solutions(exit_soln)%dt
153 texit = particle%ttrack + dtexit
159 call this%tracktimes%advance()
160 if (this%tracktimes%any())
then
161 do i = this%tracktimes%selection(1), this%tracktimes%selection(2)
162 t = this%tracktimes%times(i)
163 if (t < particle%ttrack) cycle
164 if (t >= texit .or. t >= tmax)
exit
166 x =
new_x(exit_x%v, exit_x%dvdx, subcell%vx1, subcell%vx2, &
167 dt, x0, subcell%dx, exit_x%status == 1)
168 y =
new_x(exit_y%v, exit_y%dvdx, subcell%vy1, subcell%vy2, &
169 dt, y0, subcell%dy, exit_y%status == 1)
170 z =
new_x(exit_z%v, exit_z%dvdx, subcell%vz1, subcell%vz2, &
171 dt, z0, subcell%dz, exit_z%status == 1)
172 particle%x = x * subcell%dx
173 particle%y = y * subcell%dy
174 particle%z = z * subcell%dz
177 call this%usertime(particle)
181 if (texit .gt. tmax)
then
187 x =
new_x(exit_x%v, exit_x%dvdx, subcell%vx1, subcell%vx2, &
188 dt, x0, subcell%dx, exit_x%status == 1)
189 y =
new_x(exit_y%v, exit_y%dvdx, subcell%vy1, subcell%vy2, &
190 dt, y0, subcell%dy, exit_y%status == 1)
191 z =
new_x(exit_z%v, exit_z%dvdx, subcell%vz1, subcell%vz2, &
192 dt, z0, subcell%dz, exit_z%status == 1)
195 particle%advancing = .false.
199 particle%x = x * subcell%dx
200 particle%y = y * subcell%dy
201 particle%z = z * subcell%dz
206 call this%timestep(particle)
213 if ((exit_face .eq. 1) .or. (exit_face .eq. 2))
then
215 y =
new_x(exit_y%v, exit_y%dvdx, subcell%vy1, subcell%vy2, &
216 dt, y0, subcell%dy, exit_y%status == 1)
217 z =
new_x(exit_z%v, exit_z%dvdx, subcell%vz1, subcell%vz2, &
218 dt, z0, subcell%dz, exit_z%status == 1)
219 if (exit_face .eq. 2) x =
done
220 else if ((exit_face .eq. 3) .or. (exit_face .eq. 4))
then
221 x =
new_x(exit_x%v, exit_x%dvdx, subcell%vx1, subcell%vx2, dt, &
222 x0, subcell%dx, exit_x%status == 1)
224 z =
new_x(exit_z%v, exit_z%dvdx, subcell%vz1, subcell%vz2, dt, &
225 z0, subcell%dz, exit_z%status == 1)
226 if (exit_face .eq. 4) y =
done
227 else if ((exit_face .eq. 5) .or. (exit_face .eq. 6))
then
228 x =
new_x(exit_x%v, exit_x%dvdx, subcell%vx1, subcell%vx2, &
229 dt, x0, subcell%dx, exit_x%status == 1)
230 y =
new_x(exit_y%v, exit_y%dvdx, subcell%vy1, subcell%vy2, &
231 dt, y0, subcell%dy, exit_y%status == 1)
233 if (exit_face .eq. 6) z =
done
235 print *,
"programmer error, invalid exit face", exit_face
241 particle%x = x * subcell%dx
242 particle%y = y * subcell%dy
243 particle%z = z * subcell%dz
248 call this%subcellexit(particle)
257 integer(I4B) :: exit_soln
264 if (this%exit_solutions(1)%status < 2)
then
266 dtmin = this%exit_solutions(1)%dt
268 if (this%exit_solutions(2)%status < 2 .and. &
269 this%exit_solutions(2)%dt < dtmin)
then
271 dtmin = this%exit_solutions(2)%dt
273 if (this%exit_solutions(3)%status < 2 .and. &
274 this%exit_solutions(3)%dt < dtmin)
then
286 real(DP) :: x0, y0, z0
291 x0 = particle%x / domain%dx
292 y0 = particle%y / domain%dy
293 z0 = particle%z / domain%dz
296 this%exit_solutions = [ &
297 find_exit(domain%vx1, domain%vx2, domain%dx, x0), &
298 find_exit(domain%vy1, domain%vy2, domain%dy, y0), &
299 find_exit(domain%vz1, domain%vz2, domain%dz, z0) &
303 if (this%exit_solutions(1)%v <
dzero)
then
304 this%exit_solutions(1)%iboundary = 1
305 else if (this%exit_solutions(1)%v >
dzero)
then
306 this%exit_solutions(1)%iboundary = 2
308 if (this%exit_solutions(2)%v <
dzero)
then
309 this%exit_solutions(2)%iboundary = 3
310 else if (this%exit_solutions(2)%v >
dzero)
then
311 this%exit_solutions(2)%iboundary = 4
313 if (this%exit_solutions(3)%v <
dzero)
then
314 this%exit_solutions(3)%iboundary = 5
315 else if (this%exit_solutions(3)%v >
dzero)
then
316 this%exit_solutions(3)%iboundary = 6
324 real(dp),
intent(in) :: v1
325 real(dp),
intent(in) :: v2
326 real(dp),
intent(in) :: dx
327 real(dp),
intent(in) :: xl
332 solution%v, solution%dvdx, solution%dt)
352 integer(I4B) :: status
368 logical(LGP) :: nooutflow
374 if (v2a .lt.
dzero) v2a = -v2a
376 if (v1a .lt.
dzero) v1a = -v1a
379 if (dva .lt.
dzero) dva = -dva
384 if ((v2a .lt. tol) .and. (v1a .lt. tol))
then
395 if (v2a .gt. vv) vv = v2a
397 if (vvv .lt. 1.0d-4)
then
402 if (v1 .gt. zro) dt = (dx - x) / v1
403 if (v1 .lt. zrom) dt = -x / v1
412 v = (
done - xl) * v1 + xl * v2
417 if (v1 .lt.
dzero) nooutflow = .false.
418 if (v2 .gt.
dzero) nooutflow = .false.
428 if ((v1 .le.
dzero) .and. (v2 .ge.
dzero))
then
429 if (abs(v) .le.
dzero)
then
431 if (v2 .le.
dzero) v = -v
441 if (vr .le.
dzero)
then
449 if (v1v2 .gt.
dzero)
then
450 if (v .gt.
dzero) vr = vr2
451 if (v .lt.
dzero) vr = vr1
456 if (dabs(vr) .lt. 1.0d-10)
then
476 pure function new_x(v, dvdx, v1, v2, dt, x, dx, velocity_profile)
result(newx)
478 real(dp),
intent(in) :: v
479 real(dp),
intent(in) :: dvdx
480 real(dp),
intent(in) :: v1
481 real(dp),
intent(in) :: v2
482 real(dp),
intent(in) :: dt
483 real(dp),
intent(in) :: x
484 real(dp),
intent(in) :: dx
485 logical(LGP),
intent(in),
optional :: velocity_profile
488 logical(LGP) :: lprofile
491 if (
present(velocity_profile))
then
492 lprofile = velocity_profile
500 newx = newx + (v1 * dt / dx)
501 else if (v .ne.
dzero)
then
502 newx = newx + (v * (exp(dvdx * dt) -
done) / dvdx / dx)
This module contains simulation constants.
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.
@ ok_exit_constant
exit found, constant velocity
@ ok_exit
exit found using velocity interpolation
@ no_exit_stationary
no exit, zero velocity
@ no_exit_no_outflow
no exit, no outflow
This module defines variable data types.
Particle tracking strategies.
@, public level_subfeature
pure real(dp) function new_x(v, dvdx, v1, v2, dt, x, dx, velocity_profile)
Update a cell-local coordinate based on a time increment.
integer(i4b) function pick_exit(this, particle)
Pick the exit solution with the shortest travel time.
integer(i4b) function, public calculate_dt(v1, v2, dx, xL, v, dvdx, dt)
Calculate particle travel time to exit and exit status.
type(linearexitsolutiontype) function find_exit(v1, v2, dx, xL)
Find an exit solution for one dimension.
subroutine track_subcell(this, subcell, particle, tmax)
Track a particle across a rectangular subcell using Pollock's method.
subroutine, public create_method_subcell_pollock(method)
Create a new Pollock's subcell method.
subroutine find_exits(this, particle, domain)
Compute candidate exit solutions.
subroutine apply_msp(this, particle, tmax)
Apply Pollock's method to a rectangular subcell.
@, public featexit
particle exited a grid feature
@, public timestep
time step ended
@ term_timeout
terminated at stop time or end of simulation
@ term_no_exits_sub
terminated in a subcell with no exit face
subroutine, public create_subcell_rect(subcell)
Create a new rectangular subcell.
logical(lgp), pointer, public endofsimulation
flag indicating end of simulation
Base type for grid cells of a concrete type. Contains a cell-definition which is information shared b...
Base type for exit solutions.
Linear velocity interpolation exit solution.
A generic heterogeneous doubly-linked list.
Abstract base type for subcell tracking methods.
Rectangular subcell tracking method.
Particle tracked by the PRT model.