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

Commit

Permalink
Merge branch 'debug'
Browse files Browse the repository at this point in the history
  • Loading branch information
@daminton
daminton committed Oct 22, 2021
2 parents a36d7e0 + b09010f commit 3bfff4f
Show file tree
Hide file tree
Showing 2 changed files with 156 additions and 78 deletions.
149 changes: 84 additions & 65 deletions src/encounter/encounter_check.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -447,12 +447,12 @@ subroutine encounter_check_all_sort_and_sweep_plplm(nplm, nplt, xplm, vplm, xplt
ntot = nplm + nplt
n = 2 * ntot
if (ntot /= ntot_last) then

call boundingbox%setup(ntot, ntot_last)

ntot_last = ntot
end if

call timer%reset()
call timer%start()
!$omp parallel do default(private) schedule(static) &
!$omp shared(xplm, xplt, vplm, vplt, rencm, renct, boundingbox) &
!$omp firstprivate(dt, nplm, nplt, ntot)
Expand Down Expand Up @@ -929,16 +929,22 @@ module subroutine encounter_check_sweep_aabb_double_list(self, n1, n2, nenc, ind
integer(I4B), dimension(:), allocatable, intent(out) :: index1 !! List of indices for body 1 in each encounter candidate pair
integer(I4B), dimension(:), allocatable, intent(out) :: index2 !! List of indices for body 2 in each encounter candidate pair
!Internals
integer(I4B) :: i, k, ntot
integer(I4B) :: i, k, ntot, dim
type(encounter_list), dimension(n1+n2) :: lenc !! Array of encounter lists (one encounter list per body)
integer(I4B), dimension(:), allocatable, save :: ind_arr
integer(I4B), dimension(:), allocatable :: ibeg, iend

ntot = n1 + n2
call util_index_array(ind_arr, ntot)
allocate(ibeg(SWEEPDIM * ntot))
allocate(iend(SWEEPDIM * ntot))
do dim = 1, SWEEPDIM
ibeg((dim - 1) * ntot + 1:dim * ntot) = self%aabb(dim)%ibeg(:)
iend((dim - 1) * ntot + 1:dim * ntot) = self%aabb(dim)%iend(:)
end do
! 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

call encounter_check_sweep_aabb_all_double_list(n1, n2, self%aabb(1)%ind(:), self%aabb(1)%ibeg(:), self%aabb(1)%iend(:), self%aabb(2)%ibeg(:), self%aabb(2)%iend(:), ind_arr(:), lenc(:))
call encounter_check_sweep_aabb_all_double_list(n1, n2, self%aabb(1)%ind(:), reshape(ibeg(:), [SWEEPDIM, ntot]), reshape(iend(:), [SWEEPDIM, ntot]), ind_arr(:), lenc(:))

call encounter_check_collapse_ragged_list(lenc, ntot, nenc, index1, index2)

Expand Down Expand Up @@ -966,16 +972,23 @@ module subroutine encounter_check_sweep_aabb_single_list(self, n, nenc, index1,
integer(I4B), dimension(:), allocatable, intent(out) :: index1 !! List of indices for body 1 in each encounter candidate pair
integer(I4B), dimension(:), allocatable, intent(out) :: index2 !! List of indices for body 2 in each encounter candidate pair
!Internals
Integer(I4B) :: i, k
Integer(I4B) :: i, k, dim
type(encounter_list), dimension(n) :: lenc !! Array of encounter lists (one encounter list per body)
integer(I4B), dimension(:), allocatable, save :: ind_arr
type(walltimer) :: timer
integer(I4B), dimension(:), allocatable :: ibeg, iend

call util_index_array(ind_arr, n)
allocate(ibeg(SWEEPDIM * n))
allocate(iend(SWEEPDIM * n))
do dim = 1, SWEEPDIM
ibeg((dim - 1) * n + 1:dim * n) = self%aabb(dim)%ibeg(:)
iend((dim - 1) * n + 1:dim * n) = self%aabb(dim)%iend(:)
end do

! 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
call encounter_check_sweep_aabb_all_single_list(n, self%aabb(1)%ind(:), self%aabb(1)%ibeg(:), self%aabb(1)%iend(:), self%aabb(2)%ibeg(:), self%aabb(2)%iend(:), ind_arr(:), lenc(:))
call encounter_check_sweep_aabb_all_single_list(n, self%aabb(1)%ind(:), reshape(ibeg(:), [SWEEPDIM, n]), reshape(iend(:), [SWEEPDIM, n]), ind_arr(:), lenc(:))

call encounter_check_collapse_ragged_list(lenc, n, nenc, index1, index2)

Expand All @@ -990,34 +1003,34 @@ module subroutine encounter_check_sweep_aabb_single_list(self, n, nenc, index1,
end subroutine encounter_check_sweep_aabb_single_list


subroutine encounter_check_sweep_aabb_all_double_list(n1, n2, ext_ind, ibegx, iendx, ibegy, iendy, ind_arr, lenc)
subroutine encounter_check_sweep_aabb_all_double_list(n1, n2, ext_ind, ibeg, iend, ind_arr, lenc)
!! author: David A. Minton
!!
!! Performs the loop part of the sweep operation. Double list version (e.g. pl-tp or plm-plt)
implicit none
! Arguments
integer(I4B), intent(in) :: n1, n2 !! Number of bodies
integer(I4B), dimension(:), intent(in) :: ext_ind !! Sorted index array of extents
integer(I4B), dimension(:), intent(in) :: ibegx, iendx !! Beginning and ending index lists in the x-dimension
integer(I4B), dimension(:), intent(in) :: ibegy, iendy !! Beginning and ending index lists in the y-dimension
integer(I4B), dimension(:), intent(in) :: ind_arr !! index array for mapping body 2 indexes
type(encounter_list), dimension(:), intent(inout) :: lenc !! Encounter list for the ith body
integer(I4B), intent(in) :: n1, n2 !! Number of bodies
integer(I4B), dimension(:), intent(in) :: ext_ind !! Sorted index array of extents
integer(I4B), dimension(:,:), intent(in) :: ibeg, iend !! Beginning and ending index lists in the n-dimensions
integer(I4B), dimension(:), intent(in) :: ind_arr !! index array for mapping body 2 indexes
type(encounter_list), dimension(:), intent(inout) :: lenc !! Encounter list for the ith body
! Internals
integer(I4B) :: i, ntot
logical, dimension(n1+n2) :: lencounteri
integer(I4B) :: ibegxi, iendxi, ibegyi, iendyi
integer(I4B), dimension(SWEEPDIM) :: ibegi, iendi

ntot = n1 + n2
!$omp parallel do default(private) schedule(guided)&
!$omp shared(ext_ind, ibegx, iendx, ibegy, iendy, ind_arr, lenc) &
!$omp firstprivate(ntot, n1, n2)
!$omp shared(ext_ind, ibeg, iend, ind_arr, lenc) &
!$omp firstprivate(ntot, n1, n2) &
!$omp lastprivate(ibegi, iendi)
do i = 1, ntot
ibegxi = ibegx(i) + 1
iendxi = iendx(i) - 1
if (iendxi >= ibegxi) then
ibegyi = ibegy(i)
iendyi = iendy(i)
call encounter_check_sweep_aabb_one_double_list(i, n1, n2, ntot, ext_ind(:), ibegxi, iendxi, ibegyi, iendyi, ibegx(:), iendx(:), ibegy(:), iendy(:), lencounteri(:))
ibegi(1) = ibeg(1,i) + 1
iendi(1) = iend(1,i) - 1
if (iendi(1) >= ibegi(1)) then
ibegi(2:SWEEPDIM) = ibeg(2:SWEEPDIM,i)
iendi(2:SWEEPDIM) = iend(2:SWEEPDIM,i)
call encounter_check_sweep_aabb_one_double_list(i, n1, n2, ntot, ext_ind(:), ibegi(:), iendi(:), ibeg(:,:), iend(:,:), lencounteri(:))
call encounter_check_make_ragged_list(lencounteri(:), ind_arr(:), lenc(i)%nenc, lenc(i)%index2)
else
lenc(i)%nenc = 0
Expand All @@ -1029,33 +1042,33 @@ subroutine encounter_check_sweep_aabb_all_double_list(n1, n2, ext_ind, ibegx, ie
end subroutine encounter_check_sweep_aabb_all_double_list


subroutine encounter_check_sweep_aabb_all_single_list(n, ext_ind, ibegx, iendx, ibegy, iendy, ind_arr, lenc)
subroutine encounter_check_sweep_aabb_all_single_list(n, ext_ind, ibeg, iend, ind_arr, lenc)
!! author: David A. Minton
!!
!! Performs the loop part of the sweep operation. Single list version (e.g. pl-pl)
implicit none
! Arguments
integer(I4B), intent(in) :: n !! Number of bodies
integer(I4B), dimension(:), intent(in) :: ext_ind !! Sorted index array of extents
integer(I4B), dimension(:), intent(in) :: ibegx, iendx !! Beginning and ending index lists in the x-dimension
integer(I4B), dimension(:), intent(in) :: ibegy, iendy !! Beginning and ending index lists in the y-dimension
integer(I4B), dimension(:), intent(in) :: ind_arr !! index array for mapping body 2 indexes
type(encounter_list), dimension(:), intent(inout) :: lenc !! Encounter list for the ith body
integer(I4B), intent(in) :: n !! Number of bodies
integer(I4B), dimension(:), intent(in) :: ext_ind !! Sorted index array of extents
integer(I4B), dimension(:,:), intent(in) :: ibeg, iend !! Beginning and ending index lists in the n-dimension
integer(I4B), dimension(:), intent(in) :: ind_arr !! index array for mapping body 2 indexes
type(encounter_list), dimension(:), intent(inout) :: lenc !! Encounter list for the ith body
! Internals
integer(I4B) :: i
logical, dimension(n) :: lencounteri
integer(I4B) :: ibegxi, iendxi, ibegyi, iendyi
integer(I4B), dimension(SWEEPDIM) :: ibegi, iendi

!$omp parallel do default(private) schedule(guided)&
!$omp shared(ext_ind, ibegx, iendx, ibegy, iendy, ind_arr, lenc) &
!$omp firstprivate(n)
!$omp shared(ext_ind, ibeg, iend, ind_arr, lenc) &
!$omp firstprivate(n) &
!$omp lastprivate(ibegi, iendi)
do i = 1, n
ibegxi = ibegx(i) + 1
iendxi = iendx(i) - 1
if (iendxi >= ibegxi) then
ibegyi = ibegy(i)
iendyi = iendy(i)
call encounter_check_sweep_aabb_one_single_list(n, ext_ind(:), ibegxi, iendxi, ibegyi, iendyi, ibegx(:), iendx(:), ibegy(:), iendy(:), lencounteri(:))
ibegi(1) = ibeg(1,i) + 1
iendi(1) = iend(1,i) - 1
if (iendi(1) >= ibegi(1)) then
ibegi(2:SWEEPDIM) = ibeg(2:SWEEPDIM,i)
iendi(2:SWEEPDIM) = iend(2:SWEEPDIM,i)
call encounter_check_sweep_aabb_one_single_list(n, ext_ind(:), ibegi(:), iendi(:), ibeg(:,:), iend(:,:), lencounteri(:))
call encounter_check_make_ragged_list(lencounteri(:), ind_arr(:), lenc(i)%nenc, lenc(i)%index2)
else
lenc(i)%nenc = 0
Expand All @@ -1067,60 +1080,66 @@ subroutine encounter_check_sweep_aabb_all_single_list(n, ext_ind, ibegx, iendx,
end subroutine encounter_check_sweep_aabb_all_single_list


pure subroutine encounter_check_sweep_aabb_one_double_list(i, n1, n2, ntot, ext_ind, ibegxi, iendxi, ibegyi, iendyi, ibegx, iendx, ibegy, iendy, lencounteri)
pure subroutine encounter_check_sweep_aabb_one_double_list(i, n1, n2, ntot, ext_ind, ibegi, iendi, ibeg, iend, lencounteri)
!! author: David A. Minton
!!
!! Performs a sweep operation on a single body. Encounters from the same lists not allowed (e.g. pl-tp encounters only)
implicit none
! Arguments
integer(I4B), intent(in) :: i !! The current index of the ith body
integer(I4B), intent(in) :: n1 !! Number of bodies 1
integer(I4B), intent(in) :: n2 !! Number of bodies 2
integer(I4B), intent(in) :: ntot !! n1 + n2
integer(I4B), dimension(:), intent(in) :: ext_ind !! Sorted index array of extents
integer(I4B), intent(in) :: ibegxi, iendxi !! The beginning and ending indices of the ith bounding box in the x-dimension
integer(I4B), intent(in) :: ibegyi, iendyi !! The beginning and ending indices of the ith bounding box in the y-dimension
integer(I4B), dimension(:), intent(in) :: ibegx, iendx !! Beginning and ending index lists in the x-dimension
integer(I4B), dimension(:), intent(in) :: ibegy, iendy !! Beginning and ending index lists in the y-dimensio
logical, dimension(:), intent(out) :: lencounteri !! Encounter list for the ith body
integer(I4B), intent(in) :: i !! The current index of the ith body
integer(I4B), intent(in) :: n1 !! Number of bodies 1
integer(I4B), intent(in) :: n2 !! Number of bodies 2
integer(I4B), intent(in) :: ntot !! n1 + n2
integer(I4B), dimension(:), intent(in) :: ext_ind !! Sorted index array of extents
integer(I4B), dimension(:), intent(in) :: ibegi, iendi !! The beginning and ending indices of the ith bounding box in the n-dimensions
integer(I4B), dimension(:,:), intent(in) :: ibeg, iend !! Beginning and ending index lists in the n-dimensions
logical, dimension(:), intent(out) :: lencounteri !! Encounter list for the ith body
! Internals
integer(I4B) :: j, jbox
integer(I4B) :: j, jbox, dim
logical, dimension(SWEEPDIM) :: lenc

lencounteri(:) = .false.
do concurrent(jbox = ibegxi:iendxi) ! Sweep forward until the end of the interval
lenc(1) = .true.
do concurrent(jbox = ibegi(1):iendi(1)) ! Sweep forward until the end of the interval
j = ext_ind(jbox)
if (j > ntot) j = j - ntot ! If this is an endpoint index, shift it to the correct range
if (((i <= n1) .and. (j <= n1)) .or. ((i > n1) .and. (j > n1))) cycle ! only pairs from the two different lists allowed
! Check the y-dimension
lencounteri(j) = (iendy(j) > ibegyi) .and. (ibegy(j) < iendyi)
do dim = 2, SWEEPDIM
lenc(dim) = (iend(dim,j) > ibegi(dim)) .and. (ibeg(dim,j) < iendi(dim))
end do
lencounteri(j) = all(lenc(:))
end do

return
end subroutine encounter_check_sweep_aabb_one_double_list


pure subroutine encounter_check_sweep_aabb_one_single_list(n, ext_ind, ibegxi, iendxi, ibegyi, iendyi, ibegx, iendx, ibegy, iendy, lencounteri)
pure subroutine encounter_check_sweep_aabb_one_single_list(n, ext_ind, ibegi, iendi, ibeg, iend, lencounteri)
!! author: David A. Minton
!!
!! Performs a sweep operation on a single body. Mutual encounters allowed (e.g. pl-pl)
implicit none
! Arguments
integer(I4B), intent(in) :: n !! Number of bodies
integer(I4B), dimension(:), intent(in) :: ext_ind !! Sorted index array of extents
integer(I4B), intent(in) :: ibegxi, iendxi !! The beginning and ending indices of the ith bounding box in the x-dimension
integer(I4B), intent(in) :: ibegyi, iendyi !! The beginning and ending indices of the ith bounding box in the y-dimension
integer(I4B), dimension(:), intent(in) :: ibegx, iendx !! Beginning and ending index lists in the x-dimension
integer(I4B), dimension(:), intent(in) :: ibegy, iendy !! Beginning and ending index lists in the y-dimension
logical, dimension(:), intent(out) :: lencounteri !! Encounter list for the ith body
integer(I4B), intent(in) :: n !! Number of bodies
integer(I4B), dimension(:), intent(in) :: ext_ind !! Sorted index array of extents
integer(I4B), dimension(:), intent(in) :: ibegi, iendi !! The beginning and ending indices of the ith bounding box in the n-dimensions
integer(I4B), dimension(:,:), intent(in) :: ibeg, iend !! Beginning and ending index lists in the x-dimension
logical, dimension(:), intent(out) :: lencounteri !! Encounter list for the ith body
! Internals
integer(I4B) :: j, jbox
integer(I4B) :: j, jbox, dim
logical, dimension(SWEEPDIM) :: lenc

lencounteri(:) = .false.
do concurrent(jbox = ibegxi:iendxi) ! Sweep forward until the end of the interval
lenc(1) = .true.
do concurrent(jbox = ibegi(1):iendi(1)) ! Sweep forward until the end of the interval
j = ext_ind(jbox)
if (j > n) j = j - n ! If this is an endpoint index, shift it to the correct range
! Check the y-dimension
lencounteri(j) = (iendy(j) > ibegyi) .and. (ibegy(j) < iendyi)
! Check the other dimensions
do dim = 2, SWEEPDIM
lenc(dim) = (iend(dim,j) > ibegi(dim)) .and. (ibeg(dim,j) < iendi(dim))
end do
lencounteri(j) = all(lenc(:))
end do

return
Expand Down
85 changes: 72 additions & 13 deletions src/util/util_sort.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -60,7 +60,7 @@ end subroutine util_sort_body
module pure subroutine util_sort_dp(arr)
!! author: David A. Minton
!!
!! Sort input double precision array in place into ascending numerical order using insertion sort.
!! Sort input DP precision array in place into ascending numerical order using insertion sort.
!! This algorithm works well for partially sorted arrays (which is usually the case here)
!!
implicit none
Expand Down Expand Up @@ -89,7 +89,7 @@ end subroutine util_sort_dp
module pure subroutine util_sort_index_dp(arr, ind)
!! author: David A. Minton
!!
!! Sort input double precision array by index in ascending numerical order using insertion sort.
!! Sort input DP precision array by index in ascending numerical order using insertion sort.
!! This algorithm works well for partially sorted arrays (which is usually the case here).
!! If ind is supplied already allocated, we assume it is an existing index array (e.g. a previously
!! sorted array). If it is not allocated, this subroutine allocates it.
Expand All @@ -99,27 +99,86 @@ module pure subroutine util_sort_index_dp(arr, ind)
real(DP), dimension(:), intent(in) :: arr
integer(I4B), dimension(:), allocatable, intent(inout) :: ind
! Internals
integer(I4B) :: n, i, j, jp, itmp
integer(I4B) :: n, i, j, itmp
real(DP), dimension(:), allocatable :: tmparr

n = size(arr)
if (.not.allocated(ind)) then
allocate(ind(n))
ind = [(i, i=1, n)]
end if
do i = 2, n
itmp = ind(i)
j = i - 1
do while (j >= 1)
if (arr(ind(j)) <= arr(itmp)) exit
ind(j + 1) = ind(j)
allocate(tmparr, source=arr)
call qsort_index_DP(tmparr, ind)

return
end subroutine util_sort_index_dp


recursive pure subroutine qsort_index_DP(arr, ind)
!! author: David A. Minton
!!
!! Sort input DP precision array by index in ascending numerical order using quicksort sort.
implicit none
real(DP), dimension(:), intent(inout) :: arr
integer(I4B),dimension(:),intent(out) :: ind
integer :: iq

if (size(arr) > 1) then
call partition_DP(arr, iq, ind)
call qsort_index_DP(arr(:iq-1),ind(:iq-1))
call qsort_index_DP(arr(iq:), ind(iq:))
end if

return
end subroutine qsort_index_DP


pure subroutine partition_DP(arr, marker, ind)
!! author: David A. Minton
!!
!! Partition function for quicksort on DP type
real(DP), intent(inout), dimension(:) :: arr
integer(I4B), intent(inout), dimension(:), optional :: ind
integer(I4B), intent(out) :: marker
integer(I4B) :: i, j, itmp
real(DP) :: temp
real(DP) :: x ! pivot point

x = arr(1)
i = 0
j = size(arr) + 1

do
j = j - 1
do
if (arr(j) <= x) exit
j = j - 1
end do
ind(j + 1) = itmp
i = i + 1
do
if (arr(i) >= x) exit
i = i + 1
end do
if (i < j) then
! exchange A(i) and A(j)
temp = arr(i)
arr(i) = arr(j)
arr(j) = temp
itmp = ind(i)
ind(i) = ind(j)
ind(j) = itmp
else if (i == j) then
marker = i + 1
return
else
marker = i
return
endif
end do

return
end subroutine util_sort_index_dp

end subroutine Partition_DP

module pure subroutine util_sort_i4b(arr)
!! author: David A. Minton
Expand Down

0 comments on commit 3bfff4f

Please sign in to comment.