From d0e790351bbfe57be9edc5179744a7ddf243f034 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Mon, 4 Oct 2021 16:21:56 -0400 Subject: [PATCH] Restructured encounter checking to make it more generalizable and modular --- src/encounter/encounter_check.f90 | 289 ++++++++++-------------------- src/encounter/encounter_setup.f90 | 42 +++++ src/encounter/encounter_util.f90 | 152 +++++++++++++++- src/modules/encounter_classes.f90 | 69 ++++++- 4 files changed, 352 insertions(+), 200 deletions(-) diff --git a/src/encounter/encounter_check.f90 b/src/encounter/encounter_check.f90 index 2b724731b..dba8c43fa 100644 --- a/src/encounter/encounter_check.f90 +++ b/src/encounter/encounter_check.f90 @@ -2,6 +2,49 @@ use swiftest contains + module subroutine encounter_check_all(nenc, index1, index2, x1, v1, x2, v2, renc1, renc2, dt, lencounter, lvdotr) + !! author: David A. Minton + !! + !! Check for encounters between n pairs of bodies. + !! This implementation is general for any type of body. For instance, for massive bodies, you would pass x1 = x2, for test particles renc2 is an array of zeros, etc. + !! + implicit none + ! Arguments + integer(I4B), intent(in) :: nenc !! Number of encounters in the encounter lists + integer(I4B), dimension(:), intent(in) :: index1 !! List of indices for body 1 in each encounter + integer(I4B), dimension(:), intent(in) :: index2 !! List of indices for body 2 in each encounter1 + real(DP), dimension(:,:), intent(in) :: x1, v1 !! Array of indices of bodies 1 + real(DP), dimension(:,:), intent(in) :: x2, v2 !! Array of indices of bodies 2 + real(DP), dimension(:), intent(in) :: renc1 !! Radius of encounter regions of bodies 1 + real(DP), dimension(:), intent(in) :: renc2 !! Radius of encounter regions of bodies 2 + real(DP), intent(in) :: dt !! Step size + logical, dimension(:), intent(out) :: lencounter !! Logical array indicating which pairs are in a close encounter state + logical, dimension(:), intent(out) :: lvdotr !! Logical array indicating which pairs are approaching + ! Internals + integer(I4B) :: i, j, k + real(DP) :: xr, yr, zr, vxr, vyr, vzr, renc12 + + !$omp parallel do simd default(firstprivate) schedule(static)& + !$omp shared(lencounter, lvdotr, index1, index2, x1, v1, x2, v2) & + !$omp lastprivate(i, j, xr, yr, zr, vxr, vyr, vzr, renc12) + do k = 1, nenc + i = index1(k) + j = index2(k) + xr = x2(1, j) - x1(1, i) + yr = x2(2, j) - x1(2, i) + zr = x2(3, j) - x1(3, i) + vxr = v2(1, j) - v1(1, i) + vyr = v2(2, j) - v1(2, i) + vzr = v2(3, j) - v1(3, i) + renc12 = renc1(i) + renc1(j) + call encounter_check_one(xr, yr, zr, vxr, vyr, vzr, renc12, dt, lencounter(k), lvdotr(k)) + end do + !$omp end parallel do simd + + return + end subroutine encounter_check_all + + module subroutine encounter_check_all_plpl(param, npl, nplm, x, v, renc, dt, lvdotr, index1, index2, nenc) !! author: David A. Minton !! @@ -110,25 +153,13 @@ 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 + integer(I4B) :: i, j, k, m, dim, n, i0 + logical, dimension(npl) :: lfresh integer(I4B), dimension(:), allocatable, save :: ind_arr - type lenctype - logical, dimension(:), allocatable :: lvdotr - integer(I4B), dimension(:), allocatable :: index2 - integer(I4B) :: nenc - end type - type(lenctype), dimension(npl) :: lenc + type(encounter_list), dimension(npl) :: lenc integer(I4B), dimension(:), allocatable :: tmp, ind integer(I4B), save :: npl_last = 0 - type boundingBox - integer(I4B), dimension(:), allocatable :: ind !! Sorted minimum/maximum extent indices - integer(I4B), dimension(:), allocatable :: ibeg !! Beginning index for box - integer(I4B), dimension(:), allocatable :: iend !! Ending index for box - end type - type(boundingBox), dimension(SWEEPDIM), save :: aabb + type(encounter_bounding_box), save :: boundingbox logical, dimension(:), allocatable :: lenc_final, lvdotr_final integer(I2B), dimension(npl) :: vshift_min, vshift_max integer(I4B) :: ybegi, yendi @@ -138,37 +169,16 @@ subroutine encounter_check_all_sort_and_sweep_plpl(npl, nplm, x, v, renc, dt, lv dim = 1 ! The dimension to sort and sweep (x works well for just about all cases. z is not recommended if objects are in near planar orbits, like the solar system) ! If this is the first time through, build the index lists n = 2 * npl - n_last = 2 * npl_last if (npl_last /= npl) then if (allocated(ind_arr)) deallocate(ind_arr) allocate(ind_arr(npl)) ind_arr(:) = [(i, i = 1, npl)] - ! call aabb(dim)%setup(n, n_last, npl) - - 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 + call boundingbox%setup(npl, npl_last) end if !$omp taskloop default(private) num_tasks(SWEEPDIM) & - !$omp shared(x, v, renc, aabb) & + !$omp shared(x, v, renc, boundingbox) & !$omp firstprivate(dt, npl, n) do dim = 1, SWEEPDIM where(v(dim,1:npl) < 0.0_DP) @@ -178,99 +188,21 @@ subroutine encounter_check_all_sort_and_sweep_plpl(npl, nplm, x, v, renc, dt, lv 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 + call boundingbox%aabb(dim)%sort(npl, [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]) 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) + 1 - nbox = aabb(1)%iend(i) - 1 - ybegi = aabb(2)%ibeg(i) - yendi = aabb(2)%iend(i) - lencounteri(:) = .false. - 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 - ! Check the y-dimension - lencounteri(j) = (aabb(2)%iend(j) > ybegi) .and. (aabb(2)%ibeg(j) < yendi) - end do - - lenc(i)%nenc = count(lencounteri(:)) - if (lenc(i)%nenc > 0) then - allocate(lenc(i)%index2(lenc(i)%nenc)) - lenc(i)%index2(:) = pack(ind_arr(:), lencounteri(:)) - end if - end do - !$omp end parallel do - - associate(nenc_arr => lenc(:)%nenc) - nenc = sum(nenc_arr(1:npl)) - end associate + call boundingbox%sweep(npl, npl, ind_arr, nenc, index1, index2) if (nenc > 0) then - allocate(index1(nenc)) - allocate(index2(nenc)) 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 - j1 = j0 + nenci - 1 - index1(j0:j1) = i - index2(j0:j1) = lenc(i)%index2(:) - 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. - !$omp parallel do simd default(firstprivate) schedule(static)& - !$omp shared(lenc_final, lvdotr_final, index1, index2, x, v) & - !$omp lastprivate(i, j, xr, yr, zr, vxr, vyr, vzr, renc12) - do k = 1, nenc - i = index1(k) - j = index2(k) - if ((i > nplm) .and. (j > nplm)) then - lenc_final(k) = .false. - else - xr = x(1, j) - x(1, i) - yr = x(2, j) - x(2, i) - zr = x(3, j) - x(3, i) - vxr = v(1, j) - v(1, i) - vyr = v(2, j) - v(2, i) - vzr = v(3, j) - v(3, i) - renc12 = renc(i) + renc(j) - call encounter_check_one(xr, yr, zr, vxr, vyr, vzr, renc12, dt, lenc_final(k), lvdotr_final(k)) - end if - end do - !$omp end parallel do simd + call encounter_check_all(nenc, index1, index2, x, v, x, v, renc, renc, dt, lenc_final, lvdotr_final) nenc = count(lenc_final(:)) ! Count the true number of encounters allocate(tmp(nenc)) @@ -350,29 +282,18 @@ 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 + type(encounter_list), dimension(npl) :: lenc + type(encounter_bounding_box), save :: boundingbox 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 logical, dimension(npl) :: lfresh integer(I4B), dimension(:), allocatable, save :: tpind_arr - type lenctype - logical, dimension(:), allocatable :: lvdotr - integer(I4B), dimension(:), allocatable :: index2 - integer(I4B) :: nenc - end type - type(lenctype), dimension(npl) :: lenc integer(I4B), dimension(:), allocatable :: tmp integer(I4B), save :: ntot_last = 0 integer(I4B), save :: n_last = 0 integer(I4B), save :: npl_last = 0 integer(I4B), save :: ntp_last = 0 - type boundingBox - integer(I4B), dimension(:), allocatable :: ind !! Sorted minimum/maximum extent indices - integer(I4B), dimension(:), allocatable :: ibeg !! Beginning index for box - integer(I4B), dimension(:), allocatable :: iend !! Ending index for box - end type - 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 @@ -381,42 +302,25 @@ subroutine encounter_check_all_sort_and_sweep_pltp(npl, ntp, xpl, vpl, xtp, vtp, ! If this is the first time through, build the index lists if ((ntp == 0) .or. (npl == 0)) return - ! ntot = npl + ntp - ! n = 2 * ntot - ! if ((ntp_last /= ntp) .or. (npl_last /= npl)) then - ! if (ntp_last /= ntp) then - ! if (allocated(tpind_arr)) deallocate(tpind_arr) - ! allocate(tpind_arr(ntp)) - ! tpind_arr(1:ntp) = [(i, i = 1, ntp)] - ! end if + ntot = npl + ntp + n = 2 * ntot + if ((ntp_last /= ntp) .or. (npl_last /= npl)) then + if (ntp_last /= ntp) then + if (allocated(tpind_arr)) deallocate(tpind_arr) + allocate(tpind_arr(ntp)) + tpind_arr(1:ntp) = [(i, i = 1, ntp)] + end if - ! do dim = 1, SWEEPDIM - ! if (allocated(aabb(dim)%ibeg)) deallocate(aabb(dim)%ibeg) - ! allocate(aabb(dim)%ibeg(ntot)) - ! if (allocated(aabb(dim)%iend)) deallocate(aabb(dim)%iend) - ! allocate(aabb(dim)%iend(ntot)) - ! end do + call boundingbox%setup(ntot, ntot_last) - ! 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 - ! ntp_last = ntp - ! n_last = n - ! end if + npl_last = npl + ntp_last = ntp + n_last = n + end if + nenc = 0 ! !$omp taskloop default(private) num_tasks(SWEEPDIM) & - ! !$omp shared(xpl, xtp, vpl, vtp, rencpl, renctp, aabb) & + ! !$omp shared(xpl, xtp, vpl, vtp, rencpl, renctp, boundingbox) & ! !$omp firstprivate(dt, npl, n) ! do dim = 1, SWEEPDIM ! where(vpl(dim,1:npl) < 0.0_DP) @@ -438,22 +342,22 @@ subroutine encounter_check_all_sort_and_sweep_pltp(npl, ntp, xpl, vpl, xtp, vtp, ! call util_sort([xpl(dim,1:npl)-rencpl(1:npl)+vplshift_min(1:npl)*vpl(i,1:npl)*dt, & ! xtp(dim,1:ntp)-renctp(1:ntp)+vtpshift_min(1:ntp)*vtp(i,1:ntp)*dt, & ! xpl(dim,1:npl)+rencpl(1:npl)+vplshift_max(1:npl)*vpl(i,1:npl)*dt, & - ! xtp(dim,1:ntp)+renctp(1:ntp)+vtpshift_max(1:ntp)*vtp(i,1:ntp)*dt], aabb(i)%ind) + ! xtp(dim,1:ntp)+renctp(1:ntp)+vtpshift_max(1:ntp)*vtp(i,1:ntp)*dt], boundingbox(i)%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) + ! i = boundingbox(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) + ! j = boundingbox(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 + ! boundingbox(dim)%ibeg(i) = ibox + ! boundingbox(dim)%iend(i) = jbox ! exit ! We've reached the end of this interval ! end if ! end do @@ -462,20 +366,20 @@ subroutine encounter_check_all_sort_and_sweep_pltp(npl, ntp, xpl, vpl, xtp, vtp, ! !$omp end taskloop ! !$omp parallel do default(private) schedule(static)& - ! !$omp shared(aabb, lenc, ind_arr) & + ! !$omp shared(boundingbox, lenc, ind_arr) & ! !$omp firstprivate(npl, ntp, ntot) ! do i = 1, ntot - ! ibox = aabb(1)%ibeg(i) + 1 - ! nbox = aabb(1)%iend(i) - 1 - ! ybegi = aabb(2)%ibeg(i) - ! yendi = aabb(2)%iend(i) + ! ibox = boundingbox(1)%ibeg(i) + 1 + ! nbox = boundingbox(1)%iend(i) - 1 + ! ybegi = boundingbox(2)%ibeg(i) + ! yendi = boundingbox(2)%iend(i) ! lencounteri(:) = .false. ! do concurrent(jbox = ibox:nbox) ! Sweep forward until the end of the interval - ! j = aabb(1)%ind(jbox) + ! j = boundingbox(1)%ind(jbox) ! if (j > ntot) j = j - ntot ! If this is an endpoint index, shift it to the correct range ! if (((i <= npl) .and. (j <= npl)) .or. ((i > npl) .and. (j > npl))) cycle ! ! Check the y-dimension - ! lencounteri(j) = (aabb(2)%iend(j) > ybegi) .and. (aabb(2)%ibeg(j) < yendi) + ! lencounteri(j) = (boundingbox(2)%iend(j) > ybegi) .and. (boundingbox(2)%ibeg(j) < yendi) ! end do ! lenc(i)%nenc = count(lencounteri(:)) @@ -489,22 +393,22 @@ subroutine encounter_check_all_sort_and_sweep_pltp(npl, ntp, xpl, vpl, xtp, vtp, ! ! Sweep the intervals ! lfresh(:) = .true. ! This will allow us to skip end-points of processed massive bodies ! do i = 1, npl + ntp - ! !i = aabb(1)%ind(ibox) - ! !ibox = aabb(1)%ibeg(i) + 1 - ! nbox = aabb(1)%iend(i) - 1 - ! ybegi = aabb(2)%ibeg(i) - ! yendi = aabb(2)%iend(i) + ! !i = boundingbox(1)%ind(ibox) + ! !ibox = boundingbox(1)%ibeg(i) + 1 + ! nbox = boundingbox(1)%iend(i) - 1 + ! ybegi = boundingbox(2)%ibeg(i) + ! yendi = boundingbox(2)%iend(i) ! if (i > ntot) i = i - ntot ! If this is an endpoint index, shift it to the correct range ! if (i > npl) cycle ! This is a test particle, so move on ! if (.not.lfresh(i)) cycle ! This body has already been evaluated, so move on ! lencounteri(:) = .false. ! do jbox = ibox + 1, n ! Sweep forward until the end of the interval - ! j = aabb(1)%ind(jbox) + ! j = boundingbox(1)%ind(jbox) ! if (j > ntot) j = j - ntot ! If this is an endpoint index, shift it to the correct range ! if (j == i) exit ! We've reached the end of this interval ! if (j > npl) then ! this is an unprocessed test particle ! j = j - npl ! Shift into the range of the test particles - ! lencounteri(j) = (aabb(2)%iend(j) > ybegi) .and. (aabb(2)%ibeg(j) < yendi) + ! lencounteri(j) = (boundingbox(2)%iend(j) > ybegi) .and. (boundingbox(2)%ibeg(j) < yendi) ! end if ! end do ! lfresh(i) = .false. ! This body has now been processed, so it should no longer show up in future encounter lists @@ -586,12 +490,7 @@ subroutine encounter_check_all_triangular_plpl(npl, nplm, x, v, renc, dt, lvdotr real(DP) :: xr, yr, zr, vxr, vyr, vzr, renc12 logical, dimension(npl) :: lencounteri, lvdotri integer(I4B), dimension(npl) :: ind_arr - type lenctype - logical, dimension(:), allocatable :: lvdotr - integer(I4B), dimension(:), allocatable :: index2 - integer(I4B) :: nenc - end type - type(lenctype), dimension(nplm) :: lenc + type(encounter_list), dimension(nplm) :: lenc ind_arr(:) = [(i, i = 1, npl)] !$omp parallel do default(private) schedule(static)& @@ -665,13 +564,7 @@ subroutine encounter_check_all_triangular_pltp(npl, ntp, xpl, vpl, xtp, vtp, ren real(DP) :: xr, yr, zr, vxr, vyr, vzr, renc1, renc2 logical, dimension(ntp) :: lencounteri, lvdotri integer(I4B), dimension(ntp) :: ind_arr - type lenctype - logical, dimension(:), allocatable :: lvdotr - integer(I4B), dimension(:), allocatable :: index2 - integer(I4B) :: nenc - end type - type(lenctype), dimension(npl) :: lenc - + type(encounter_list), dimension(npl) :: lenc ind_arr(:) = [(i, i = 1, ntp)] !$omp parallel do default(private) schedule(static)& diff --git a/src/encounter/encounter_setup.f90 b/src/encounter/encounter_setup.f90 index 91d8d52d6..26a7e70c9 100644 --- a/src/encounter/encounter_setup.f90 +++ b/src/encounter/encounter_setup.f90 @@ -2,6 +2,48 @@ use swiftest contains + module subroutine encounter_setup_aabb(self, n, n_last) + !! author: David A. Minton + !! + !! Sets up or modifies an axis-aligned bounding box structure. + implicit none + ! Arguments + class(encounter_bounding_box), intent(inout) :: self !! Swiftest encounter structure + integer(I4B), intent(in) :: n !! Number of objects with bounding box extents + integer(I4B), intent(in) :: n_last !! Number of objects with bounding box extents the previous time this was called + ! Internals + integer(I4B) :: next, next_last, k, dim + integer(I4B), dimension(:), allocatable :: tmp + + next = 2 * n + next_last = 2 * n_last + + if (n > n_last) then ! The number of bodies has grown. Resize and append the new bodies + do dim = 1, SWEEPDIM + allocate(tmp(next)) + if (n_last > 0) tmp(1:next_last) = self%aabb(dim)%ind(1:next_last) + call move_alloc(tmp, self%aabb(dim)%ind) + self%aabb(dim)%ind(next_last+1:next) = [(k, k = next_last+1, next)] + end do + else ! The number of bodies has gone down. Resize and chop of the old indices + do dim = 1, SWEEPDIM + allocate(tmp(next)) + tmp(1:next) = pack(self%aabb(dim)%ind(1:next_last), self%aabb(dim)%ind(1:next_last) <= next) + call move_alloc(tmp, self%aabb(dim)%ind) + end do + end if + + do dim = 1, SWEEPDIM + if (allocated(self%aabb(dim)%ibeg)) deallocate(self%aabb(dim)%ibeg) + allocate(self%aabb(dim)%ibeg(n)) + if (allocated(self%aabb(dim)%iend)) deallocate(self%aabb(dim)%iend) + allocate(self%aabb(dim)%iend(n)) + end do + + return + end subroutine encounter_setup_aabb + + module subroutine encounter_setup_list(self, n) !! author: David A. Minton !! diff --git a/src/encounter/encounter_util.f90 b/src/encounter/encounter_util.f90 index ae4a6d486..9f00cf67c 100644 --- a/src/encounter/encounter_util.f90 +++ b/src/encounter/encounter_util.f90 @@ -11,7 +11,7 @@ module subroutine encounter_util_append_list(self, source, lsource_mask) ! Arguments class(encounter_list), intent(inout) :: self !! Swiftest encounter list object class(encounter_list), intent(in) :: source !! Source object to append - logical, dimension(:), intent(in) :: lsource_mask !! Logical mask indicating which elements to append to + logical, dimension(:), intent(in) :: lsource_mask !! Logical mask indicating which elements to append to ! Internals integer(I4B) :: nold, nsrc @@ -62,6 +62,45 @@ module subroutine encounter_util_copy_list(self, source) end subroutine encounter_util_copy_list + module subroutine encounter_util_collapse_ragged_list(ragged_list, n1, n2, nenc, index1, index2, lvdotr) + !! author: David A. Minton + !! + !! Collapses a ragged index list (one encounter list per body) into a pair of index arrays and a vdotr logical array (optional) + implicit none + ! Arguments + type(encounter_list), dimension(:), intent(in) :: ragged_list !! The ragged encounter list + integer(I4B), intent(in) :: n1 !! Number of bodies 1 + integer(I4B), intent(in) :: n2 !! Number of bodies 2 + integer(I4B), intent(out) :: nenc !! Total number of encountersj + integer(I4B), dimension(:), allocatable, intent(out) :: index1 !! Array of indices for body 1 + integer(I4B), dimension(:), allocatable, intent(out) :: index2 !! Array of indices for body 1 + integer(I4B), dimension(:), allocatable, intent(out), optional :: lvdotr !! Array indicating which bodies are approaching + ! Internals + integer(I4B) :: i, j0, j1, nenci + + associate(nenc_arr => ragged_list(:)%nenc) + nenc = sum(nenc_arr(1:n1)) + end associate + if (nenc == 0) return + + allocate(index1(nenc)) + allocate(index2(nenc)) + j0 = 1 + do i = 1, n1 + nenci = ragged_list(i)%nenc + if (nenci > 0) then + j1 = j0 + nenci - 1 + index1(j0:j1) = i + index2(j0:j1) = ragged_list(i)%index2(:) + if (present(lvdotr)) lvdotr(j0:j1) = ragged_list(i)%lvdotr(:) + j0 = j1 + 1 + end if + end do + + return + end subroutine encounter_util_collapse_ragged_list + + module subroutine encounter_util_resize_list(self, nnew) !! author: David A. Minton !! @@ -100,6 +139,45 @@ module subroutine encounter_util_resize_list(self, nnew) end subroutine encounter_util_resize_list + module subroutine encounter_util_sort_aabb_1D(self, n, extent_arr) + !! author: David A. Minton + !! + !! Sorts the bounding box extents along a single dimension prior to the sweep phase. + !! This subroutine sets the sorted index array (ind) and the beginning/ending index list (beg & end) + implicit none + ! Arguments + class(encounter_bounding_box_1D), intent(inout) :: self !! Bounding box structure along a single dimension + integer(I4B), intent(in) :: n !! Number of bodies with extents + real(DP), dimension(:), intent(in) :: extent_arr !! Array of extents of size 2*n + ! Internals + integer(I4B) :: i, j, ibox, jbox + logical, dimension(:), allocatable :: lfresh + + call util_sort(extent_arr, self%ind) + allocate(lfresh(n)) + + ! Determine the interval starting points and sizes + lfresh(:) = .true. ! This will prevent double-counting of pairs + do ibox = 1, 2*n + i = self%ind(ibox) + if (i > n) i = i - n ! If this is an endpoint index, shift it to the correct range + if (.not.lfresh(i)) cycle + do jbox = ibox + 1, 2*n + j = self%ind(jbox) + if (j > n) j = j - n ! If this is an endpoint index, shift it to the correct range + if (j == i) then + lfresh(i) = .false. + self%ibeg(i) = ibox + self%iend(i) = jbox + exit ! We've reached the end of this interval + end if + end do + end do + + return + end subroutine encounter_util_sort_aabb_1D + + module subroutine encounter_util_spill_list(self, discards, lspill_list, ldestructive) !! author: David A. Minton !! @@ -136,6 +214,78 @@ module subroutine encounter_util_spill_list(self, discards, lspill_list, ldestru return end subroutine encounter_util_spill_list + + + module subroutine encounter_util_sweep_aabb(self, n1, n2, ind_arr2, nenc, index1, index2) + !! author: David A. Minton + !! + !! Sweeps the sorted bounding box extents and returns the encounter candidates + implicit none + ! Arguments + class(encounter_bounding_box), intent(inout) :: self !! Multi-dimensional bounding box structure + integer(I4B), intent(in) :: n1 !! Number of bodies 1 + integer(I4B), intent(in) :: n2 !! Number of bodies 2 + integer(I4B), dimension(:), intent(in) :: ind_arr2 !! index array for mapping body 2 indexes + integer(I4B), intent(out) :: nenc !! Total number of encounter candidates + 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, n + type(encounter_list), dimension(n1) :: lenc !! Array of encounter lists (one encounter list per body) + + ! 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(self, lenc, ind_arr) & + !$omp firstprivate(n1) + do i = 1, n1 + call encounter_util_sweep_aabb_one(i, n1, self%aabb(1)%ind(:), self%aabb(1)%ibeg(:), self%aabb(1)%iend(:), self%aabb(2)%ibeg(:), self%aabb(2)%iend(:), ind_arr2, lenc(i)) + end do + !$omp end parallel do + + call encounter_util_collapse_ragged_list(lenc, n1, n2, nenc, index1, index2) + + return + end subroutine encounter_util_sweep_aabb + + + subroutine encounter_util_sweep_aabb_one(i, n1, ext_ind, ibegx, iendx, ibegy, iendy, ind_arr2, lenc) + !! author: David A. Minton + !! + !! Performs a sweep operation on a single body + 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), 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_arr2 !! index array for mapping body 2 indexes + type(encounter_list), intent(inout) :: lenc !! Encounter list for the ith body + ! Internals + integer(I4B) :: ibox, jbox, nbox, j, ybegi, yendi + logical, dimension(n1) :: lencounteri + + ibox = ibegx(i) + 1 + nbox = iendx(i) - 1 + ybegi = ibegy(i) + yendi = iendy(i) + lencounteri(:) = .false. + do concurrent(jbox = ibox:nbox) ! Sweep forward until the end of the interval + j = ext_ind(jbox) + if (j > n1) j = j - n1 ! If this is an endpoint index, shift it to the correct range + ! Check the y-dimension + lencounteri(j) = (iendy(j) > ybegi) .and. (ibegy(j) < yendi) + end do + + lenc%nenc = count(lencounteri(:)) + if (lenc%nenc > 0) then + allocate(lenc%index2(lenc%nenc)) + lenc%index2(:) = pack(ind_arr2(:), lencounteri(:)) + end if + + return + end subroutine encounter_util_sweep_aabb_one end submodule s_encounter_util \ No newline at end of file diff --git a/src/modules/encounter_classes.f90 b/src/modules/encounter_classes.f90 index a8cbdf440..f31bf2548 100644 --- a/src/modules/encounter_classes.f90 +++ b/src/modules/encounter_classes.f90 @@ -7,6 +7,8 @@ module encounter_classes implicit none public + integer(I4B), parameter :: SWEEPDIM = 2 + type :: encounter_list integer(I4B) :: nenc !! Total number of encounters logical, dimension(:), allocatable :: lvdotr !! relative vdotr flag @@ -29,8 +31,37 @@ module encounter_classes procedure :: write => encounter_io_write_list !! Write close encounter data to output binary file end type encounter_list + type encounter_bounding_box_1D + integer(I4B) :: n !! Number of bodies with extents + integer(I4B), dimension(:), allocatable :: ind !! Sorted minimum/maximum extent indices + integer(I4B), dimension(:), allocatable :: ibeg !! Beginning index for box + integer(I4B), dimension(:), allocatable :: iend !! Ending index for box + contains + procedure :: sort => encounter_util_sort_aabb_1D !! Sorts the bounding box extents along a single dimension prior to the sweep phase + end type + + type encounter_bounding_box + type(encounter_bounding_box_1D), dimension(SWEEPDIM) :: aabb + contains + procedure :: setup => encounter_setup_aabb !! Setup a new axis-aligned bounding box structure + procedure :: sweep => encounter_util_sweep_aabb !! Sweeps the sorted bounding box extents and returns the encounter candidates + end type interface + module subroutine encounter_check_all(nenc, index1, index2, x1, v1, x2, v2, renc1, renc2, dt, lencounter, lvdotr) + implicit none + integer(I4B), intent(in) :: nenc !! Number of encounters in the encounter lists + integer(I4B), dimension(:), intent(in) :: index1 !! List of indices for body 1 in each encounter + integer(I4B), dimension(:), intent(in) :: index2 !! List of indices for body 2 in each encounter1 + real(DP), dimension(:,:), intent(in) :: x1, v1 !! Array of indices of bodies 1 + real(DP), dimension(:,:), intent(in) :: x2, v2 !! Array of indices of bodies 2 + real(DP), dimension(:), intent(in) :: renc1 !! Radius of encounter regions of bodies 1 + real(DP), dimension(:), intent(in) :: renc2 !! Radius of encounter regions of bodies 2 + real(DP), intent(in) :: dt !! Step size + logical, dimension(:), intent(out) :: lencounter !! Logical array indicating which pairs are in a close encounter state + logical, dimension(:), intent(out) :: lvdotr !! Logical array indicating which pairs are approaching + end subroutine encounter_check_all + module subroutine encounter_check_all_plpl(param, npl, nplm, x, v, renc, dt, lvdotr, index1, index2, nenc) import swiftest_parameters implicit none @@ -95,6 +126,13 @@ module subroutine encounter_io_write_list(self, pl, encbody, param) class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters end subroutine encounter_io_write_list + module subroutine encounter_setup_aabb(self, n, n_last) + implicit none + class(encounter_bounding_box), intent(inout) :: self !! Swiftest encounter structure + integer(I4B), intent(in) :: n !! Number of objects with bounding box extents + integer(I4B), intent(in) :: n_last !! Number of objects with bounding box extents the previous time this was called + end subroutine encounter_setup_aabb + module subroutine encounter_setup_list(self, n) implicit none class(encounter_list), intent(inout) :: self !! Swiftest encounter structure @@ -108,6 +146,17 @@ module subroutine encounter_util_append_list(self, source, lsource_mask) logical, dimension(:), intent(in) :: lsource_mask !! Logical mask indicating which elements to append to end subroutine encounter_util_append_list + module subroutine encounter_util_collapse_ragged_list(ragged_list, n1, n2, nenc, index1, index2, lvdotr) + implicit none + type(encounter_list), dimension(:), intent(in) :: ragged_list !! The ragged encounter list + integer(I4B), intent(in) :: n1 !! Number of bodies 1 + integer(I4B), intent(in) :: n2 !! Number of bodies 2 + integer(I4B), intent(out) :: nenc !! Total number of encountersj + integer(I4B), dimension(:), allocatable, intent(out) :: index1 !! Array of indices for body 1 + integer(I4B), dimension(:), allocatable, intent(out) :: index2 !! Array of indices for body 1 + integer(I4B), dimension(:), allocatable, intent(out), optional :: lvdotr !! Array indicating which bodies are approaching + end subroutine encounter_util_collapse_ragged_list + module subroutine encounter_util_copy_list(self, source) implicit none class(encounter_list), intent(inout) :: self !! Encounter list @@ -117,9 +166,27 @@ end subroutine encounter_util_copy_list module subroutine encounter_util_resize_list(self, nnew) implicit none class(encounter_list), intent(inout) :: self !! Swiftest encounter list - integer(I4B), intent(in) :: nnew !! New size of list needed + integer(I4B), intent(in) :: nnew !! New size of list needed end subroutine encounter_util_resize_list + module subroutine encounter_util_sort_aabb_1D(self, n, extent_arr) + implicit none + class(encounter_bounding_box_1D), intent(inout) :: self !! Bounding box structure along a single dimension + integer(I4B), intent(in) :: n !! Number of bodies with extents + real(DP), dimension(:), intent(in) :: extent_arr !! Array of extents of size 2*n + end subroutine encounter_util_sort_aabb_1D + + module subroutine encounter_util_sweep_aabb(self, n1, n2, ind_arr2, nenc, index1, index2) + implicit none + class(encounter_bounding_box), intent(inout) :: self !! Multi-dimensional bounding box structure + integer(I4B), intent(in) :: n1 !! Number of bodies 1 + integer(I4B), intent(in) :: n2 !! Number of bodies 2 + integer(I4B), dimension(:), intent(in) :: ind_arr2 !! index array for mapping body 2 indexes + integer(I4B), intent(out) :: nenc !! Total number of encounter candidates + 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 + end subroutine encounter_util_sweep_aabb + module subroutine encounter_util_spill_list(self, discards, lspill_list, ldestructive) implicit none class(encounter_list), intent(inout) :: self !! Swiftest encounter list