Skip to content
This repository was archived by the owner on Aug 28, 2024. It is now read-only.

Commit

Permalink
Switched to using 2D sort and sweep.
Browse files Browse the repository at this point in the history
  • Loading branch information
@daminton
daminton committed Oct 1, 2021
1 parent 73407cd commit 39a5924
Showing 1 changed file with 68 additions and 50 deletions.
118 changes: 68 additions & 50 deletions src/encounter/encounter.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -51,7 +51,7 @@ module subroutine encounter_check_all_plpl(param, npl, nplm, x, v, renc, dt, lvd
call encounter_check_all_triangular_plpl(npl, nplm, x, v, renc, dt, lvdotr, index1, index2, nenc)
end if

if (param%ladaptive_encounters .and. nplplm > 0) then
if (.not.lfirst .and. param%ladaptive_encounters .and. nplplm > 0) then
if (itimer%is_on) call itimer%adapt(param, nplplm)
end if

Expand Down Expand Up @@ -107,6 +107,7 @@ subroutine encounter_check_all_sort_and_sweep_plpl(npl, nplm, x, v, renc, dt, lv
integer(I4B), dimension(:), allocatable, intent(out) :: index2 !! List of indices for body 2 in each encounter
integer(I4B), intent(out) :: nenc !! Total number of encounters
! Internals
integer(I4B), parameter :: SWEEPDIM = 2
integer(I4B) :: i, j, k, m, nenci, i0, i1, j0, j1, dim, ibox, jbox, nbox, n, n_last
real(DP) :: xr, yr, zr, vxr, vyr, vzr, renc12
logical, dimension(npl) :: lencounteri, lfresh
Expand All @@ -124,9 +125,11 @@ subroutine encounter_check_all_sort_and_sweep_plpl(npl, nplm, x, v, renc, dt, lv
integer(I4B), dimension(:), allocatable :: ibeg !! Beginning index for box
integer(I4B), dimension(:), allocatable :: iend !! Ending index for box
end type
type(boundingBox), dimension(NDIM), save :: aabb
type(boundingBox), dimension(SWEEPDIM), save :: aabb
logical, dimension(:), allocatable :: lenc_final, lvdotr_final
integer(I2B), dimension(npl) :: vshift_min, vshift_max
integer :: ybegi, yendi
logical :: lency

if (npl <= 1) return

Expand All @@ -137,66 +140,78 @@ subroutine encounter_check_all_sort_and_sweep_plpl(npl, nplm, x, v, renc, dt, lv
if (npl_last /= npl) then
if (allocated(ind_arr)) deallocate(ind_arr)
allocate(ind_arr(npl))
if (allocated(aabb(dim)%ibeg)) deallocate(aabb(dim)%ibeg)
allocate(aabb(dim)%ibeg(npl))
if (allocated(aabb(dim)%iend)) deallocate(aabb(dim)%iend)
allocate(aabb(dim)%iend(npl))
ind_arr(:) = [(i, i = 1, npl)]
if (npl > npl_last) then ! The number of bodies has grown. Resize and append the new bodies
allocate(tmp(n))
if (npl_last > 0) tmp(1:n_last) = aabb(dim)%ind(1:n_last)
call move_alloc(tmp, aabb(dim)%ind)
aabb(dim)%ind(n_last+1:n) = [(k, k = n_last+1, n)]
else ! The number of bodies has gone down. Resize and chop of the old indices
allocate(tmp(n))
tmp(1:n) = pack(aabb(dim)%ind(1:n_last), aabb(dim)%ind(1:n_last) <= n)
call move_alloc(tmp, aabb(dim)%ind)
end if
do dim = 1, SWEEPDIM
if (allocated(aabb(dim)%ibeg)) deallocate(aabb(dim)%ibeg)
allocate(aabb(dim)%ibeg(npl))
if (allocated(aabb(dim)%iend)) deallocate(aabb(dim)%iend)
allocate(aabb(dim)%iend(npl))
end do
do dim = 1, SWEEPDIM
if (npl > npl_last) then ! The number of bodies has grown. Resize and append the new bodies
allocate(tmp(n))
if (npl_last > 0) tmp(1:n_last) = aabb(dim)%ind(1:n_last)
call move_alloc(tmp, aabb(dim)%ind)
aabb(dim)%ind(n_last+1:n) = [(k, k = n_last+1, n)]
else ! The number of bodies has gone down. Resize and chop of the old indices
allocate(tmp(n))
tmp(1:n) = pack(aabb(dim)%ind(1:n_last), aabb(dim)%ind(1:n_last) <= n)
call move_alloc(tmp, aabb(dim)%ind)
end if
end do
npl_last = npl
end if

where(v(dim,1:npl) < 0.0_DP)
vshift_min(1:npl) = 1
vshift_max(1:npl) = 0
elsewhere
vshift_min(1:npl) = 0
vshift_max(1:npl) = 1
end where
call util_sort([x(dim,1:npl)-renc(1:npl)+vshift_min(1:npl)*v(dim,1:npl)*dt, &
x(dim,1:npl)+renc(1:npl)+vshift_max(1:npl)*v(dim,1:npl)*dt], aabb(dim)%ind)

! Determine the interval starting points and sizes

lfresh(:) = .true. ! This will prevent double-counting of pairs
do ibox = 1, n
i = aabb(dim)%ind(ibox)
if (i > npl) i = i - npl ! If this is an endpoint index, shift it to the correct range
if (.not.lfresh(i)) cycle
do jbox = ibox + 1, n
j = aabb(dim)%ind(jbox)
if (j > npl) j = j - npl ! If this is an endpoint index, shift it to the correct range
if (j == i) then
lfresh(i) = .false.
aabb(dim)%ibeg(i) = ibox
aabb(dim)%iend(i) = jbox
exit ! We've reached the end of this interval
end if
!$omp taskloop default(private) num_tasks(SWEEPDIM) &
!$omp shared(x, v, renc, aabb) &
!$omp firstprivate(dt, npl, n)
do dim = 1, SWEEPDIM
where(v(dim,1:npl) < 0.0_DP)
vshift_min(1:npl) = 1
vshift_max(1:npl) = 0
elsewhere
vshift_min(1:npl) = 0
vshift_max(1:npl) = 1
end where
call util_sort([x(dim,1:npl)-renc(1:npl)+vshift_min(1:npl)*v(dim,1:npl)*dt, &
x(dim,1:npl)+renc(1:npl)+vshift_max(1:npl)*v(dim,1:npl)*dt], aabb(dim)%ind)

! Determine the interval starting points and sizes
lfresh(:) = .true. ! This will prevent double-counting of pairs
do ibox = 1, n
i = aabb(dim)%ind(ibox)
if (i > npl) i = i - npl ! If this is an endpoint index, shift it to the correct range
if (.not.lfresh(i)) cycle
do jbox = ibox + 1, n
j = aabb(dim)%ind(jbox)
if (j > npl) j = j - npl ! If this is an endpoint index, shift it to the correct range
if (j == i) then
lfresh(i) = .false.
aabb(dim)%ibeg(i) = ibox
aabb(dim)%iend(i) = jbox
exit ! We've reached the end of this interval
end if
end do
end do
end do
!$omp end taskloop

! Sweep the intervals for each of the massive bodies along one dimension
! This will build a ragged pair of index lists inside of the lenc data structure
!$omp parallel do default(private) schedule(static)&
!$omp shared(aabb, lenc, ind_arr) &
!$omp firstprivate(npl)
do i = 1, npl
ibox = aabb(1)%ibeg(i)
ibox = aabb(1)%ibeg(i) + 1
nbox = aabb(1)%iend(i) - 1
ybegi = aabb(2)%ibeg(i)
yendi = aabb(2)%iend(i)
lencounteri(:) = .false.

do jbox = ibox + 1, nbox ! Sweep forward until the end of the interval
do concurrent(jbox = ibox:nbox) ! Sweep forward until the end of the interval
j = aabb(1)%ind(jbox)
if (j > npl) j = j - npl ! If this is an endpoint index, shift it to the correct range
lencounteri(j) = .true.
! Check the y-dimension
lencounteri(j) = (aabb(2)%iend(j) > ybegi) .and. (aabb(2)%ibeg(j) < yendi)
end do

lenc(i)%nenc = count(lencounteri(:))
Expand All @@ -217,6 +232,7 @@ subroutine encounter_check_all_sort_and_sweep_plpl(npl, nplm, x, v, renc, dt, lv
allocate(lenc_final(nenc))
allocate(lvdotr_final(nenc))
j0 = 1
! Convert the ragged index lists inside the lenc data structure into two linear arrays
do i = 1, npl
nenci = lenc(i)%nenc
if (nenci > 0) then
Expand All @@ -226,6 +242,7 @@ subroutine encounter_check_all_sort_and_sweep_plpl(npl, nplm, x, v, renc, dt, lv
j0 = j1 + 1
end if
end do

! Now that we have identified potential pairs, use the narrow-phase process to get the final values
lenc_final(:) = .true.

Expand Down Expand Up @@ -327,6 +344,7 @@ subroutine encounter_check_all_sort_and_sweep_pltp(npl, ntp, xpl, vpl, xtp, vtp,
integer(I4B), dimension(:), allocatable, intent(out) :: index2 !! List of indices for body 2 in each encounter
integer(I4B), intent(out) :: nenc !! Total number of encounter
! Internals
integer(I4B), parameter :: SWEEPDIM = 2
integer(I4B) :: i, j, k, nenci, j0, j1, dim, ibox, jbox, n, ntot
real(DP) :: xr, yr, zr, vxr, vyr, vzr, renc12
logical, dimension(ntp) :: lencounteri
Expand All @@ -346,7 +364,7 @@ subroutine encounter_check_all_sort_and_sweep_pltp(npl, ntp, xpl, vpl, xtp, vtp,
type boundingBox
integer(I4B), dimension(:), allocatable :: ind !! Sorted minimum/maximum extent indices
end type
type(boundingBox), dimension(NDIM), save :: aabb
type(boundingBox), dimension(SWEEPDIM), save :: aabb
logical, dimension(:), allocatable :: lenc_final, lvdotr_final
integer(I2B), dimension(npl) :: vplshift_min, vplshift_max
integer(I2B), dimension(ntp) :: vtpshift_min, vtpshift_max
Expand All @@ -364,14 +382,14 @@ subroutine encounter_check_all_sort_and_sweep_pltp(npl, ntp, xpl, vpl, xtp, vtp,
end if

if (n > n_last) then ! The number of bodies has grown. Resize and append the new bodies
do dim = 1, NDIM
do dim = 1, SWEEPDIM
allocate(tmp(n))
if (npl_last > 0) tmp(1:n_last) = aabb(dim)%ind(1:n_last)
call move_alloc(tmp, aabb(dim)%ind)
aabb(dim)%ind(n_last+1:n) = [(k, k = n_last+1, n)]
end do
else if (n < n_last) then ! The number of bodies has gone down. Resize and chop of the old indices
do dim = 1, NDIM
do dim = 1, SWEEPDIM
allocate(tmp(n))
tmp(1:n) = pack(aabb(dim)%ind(1:n_last), aabb(dim)%ind(1:n_last) <= n)
call move_alloc(tmp, aabb(dim)%ind)
Expand All @@ -383,7 +401,7 @@ subroutine encounter_check_all_sort_and_sweep_pltp(npl, ntp, xpl, vpl, xtp, vtp,
n_last = n
end if

do i = 1, NDIM
do i = 1, SWEEPDIM
where(vpl(i,1:npl) < 0.0_DP)
vplshift_min(1:npl) = 1
vplshift_max(1:npl) = 0
Expand Down

0 comments on commit 39a5924

Please sign in to comment.