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

Commit

Permalink
Fixed adaptive encounter check so that it properly skips at least onc…
Browse files Browse the repository at this point in the history 
…e before timing
  • Loading branch information
@daminton
daminton committed Oct 7, 2021
1 parent 1d41b68 commit bdfc400
Show file tree
Hide file tree
Showing 2 changed files with 91 additions and 51 deletions.
136 changes: 85 additions & 51 deletions src/encounter/encounter_check.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -65,25 +65,19 @@ module subroutine encounter_check_all_plpl(param, npl, x, v, renc, dt, lvdotr, i
! Internals
type(interaction_timer), save :: itimer
logical, save :: lfirst = .true.
logical, save :: lsecond = .false.
logical, save :: skipit = .false. ! This will be used to ensure that the sort & sweep subroutine gets called at least once before timing it so that the extent array is nearly sorted when it is timed
integer(I8B) :: nplpl = 0_I8B

if (param%ladaptive_encounters_plpl) then
if (param%ladaptive_encounters_plpl .and. (.not. skipit)) then
nplpl = (npl * (npl - 1) / 2)
if (nplpl > 0) then
if (lfirst) then
write(itimer%loopname, *) "encounter_check_all_plpl"
write(itimer%looptype, *) "ENCOUNTER_PLPL"
lfirst = .false.
lsecond = .true.
else
if (lsecond) then ! This ensures that the encounter check methods are run at least once prior to timing. Sort and sweep improves on the second pass due to the bounding box extents needing to be nearly sorted
call itimer%time_this_loop(param, nplpl)
lsecond = .false.
else if (itimer%check(param, nplpl)) then
lsecond = .true.
itimer%is_on = .false.
end if
itimer%step_counter = INTERACTION_TIMER_CADENCE
else
if (itimer%check(param, nplpl)) call itimer%time_this_loop(param, nplpl)
end if
else
param%lencounter_sas_plpl = .false.
Expand All @@ -96,8 +90,15 @@ module subroutine encounter_check_all_plpl(param, npl, x, v, renc, dt, lvdotr, i
call encounter_check_all_triangular_plpl(npl, x, v, renc, dt, lvdotr, index1, index2, nenc)
end if

if (.not.lfirst .and. param%ladaptive_encounters_plpl .and. nplpl > 0) then
if (itimer%is_on) call itimer%adapt(param, nplpl)
if (skipit) then
skipit = .false.
else
if (param%ladaptive_encounters_plpl .and. nplpl > 0) then
if (itimer%is_on) then
call itimer%adapt(param, nplpl)
skipit = .true.
end if
end if
end if

return
Expand Down Expand Up @@ -128,7 +129,7 @@ module subroutine encounter_check_all_plplm(param, nplm, nplt, xplm, vplm, xplt,
! Internals
type(interaction_timer), save :: itimer
logical, save :: lfirst = .true.
logical, save :: lsecond = .false.
logical, save :: skipit = .false.
integer(I8B) :: nplplm = 0_I8B
integer(I4B) :: npl, i
logical, dimension(:), allocatable :: plmplt_lvdotr !! Logical flag indicating the sign of v .dot. x in the plm-plt group
Expand All @@ -139,34 +140,28 @@ module subroutine encounter_check_all_plplm(param, nplm, nplt, xplm, vplm, xplt,
integer(I4B), dimension(:), allocatable :: itmp, ind
logical, dimension(:), allocatable :: ltmp

! Start with the fully interacting bodies
allocate(tmp_param, source=param)

if (param%ladaptive_encounters_plpl) then
if (param%ladaptive_encounters_plpl .and. (.not. skipit)) then
npl = nplm + nplt
nplplm = nplm * npl - nplm * (nplm + 1) / 2
if (nplplm > 0) then
if (lfirst) then
write(itimer%loopname, *) "encounter_check_all_plpl"
write(itimer%looptype, *) "ENCOUNTER_PLPL"
lfirst = .false.
lsecond = .true.
else
if (lsecond) then ! This ensures that the encounter check methods are run at least once prior to timing. Sort and sweep improves on the second pass due to the bounding box extents needing to be nearly sorted
call itimer%time_this_loop(param, nplplm)
lsecond = .false.
else if (itimer%check(param, nplplm)) then
lsecond = .true.
itimer%is_on = .false.
end if
itimer%step_counter = INTERACTION_TIMER_CADENCE
else
if (itimer%check(param, nplplm)) call itimer%time_this_loop(param, nplplm)
end if
else
param%lencounter_sas_plpl = .false.
end if
! Turn off adaptive encounter checks for the pl-pl group
tmp_param%ladaptive_encounters_plpl = .false.
end if

allocate(tmp_param, source=param)

! Turn off adaptive encounter checks for the pl-pl group
tmp_param%ladaptive_encounters_plpl = .false.

! Start with the pl-pl group
call encounter_check_all_plpl(tmp_param, nplm, xplm, vplm, rencm, dt, lvdotr, index1, index2, nenc)

Expand All @@ -176,8 +171,15 @@ module subroutine encounter_check_all_plplm(param, nplm, nplt, xplm, vplm, xplt,
call encounter_check_all_triangular_plplm(nplm, nplt, xplm, vplm, xplt, vplt, rencm, renct, dt, plmplt_lvdotr, plmplt_index1, plmplt_index2, plmplt_nenc)
end if

if (.not.lfirst .and. param%ladaptive_encounters_plpl .and. nplplm > 0) then
if (itimer%is_on) call itimer%adapt(param, nplplm)
if (skipit) then
skipit = .false.
else
if (param%ladaptive_encounters_plpl .and. nplplm > 0) then
if (itimer%is_on) then
call itimer%adapt(param, nplplm)
skipit = .true.
end if
end if
end if

if (plmplt_nenc > 0) then ! Consolidate the two lists
Expand Down Expand Up @@ -367,16 +369,16 @@ subroutine encounter_check_all_sort_and_sweep_plpl(npl, x, v, renc, dt, lvdotr,
type(walltimer) :: timer

if (npl <= 1) return

! call timer%reset()
! call timer%start()
! If this is the first time through, build the index lists
n = 2 * npl
if (npl_last /= npl) then
call boundingbox%setup(npl, npl_last)
npl_last = npl
end if

! call timer%reset()
! call timer%start()

!$omp parallel do default(private) schedule(static) &
!$omp shared(x, v, renc, boundingbox) &
!$omp firstprivate(dt, npl, n)
Expand All @@ -392,21 +394,34 @@ subroutine encounter_check_all_sort_and_sweep_plpl(npl, x, v, renc, dt, lvdotr,
x(dim,1:npl) + renc(1:npl) + vshift_max(1:npl) * v(dim,1:npl) * dt])
end do
!$omp end parallel do
! call timer%stop()
! call timer%report(nsubsteps=nthreads, message="AABB sort plpl:")


! call timer%reset()
! call timer%start()
call boundingbox%sweep(npl, nenc, index1, index2)
! call timer%stop()
! call timer%report(nsubsteps=1, message="Sweep plpl:")

if (nenc > 0) then
! Now that we have identified potential pairs, use the narrow-phase process to get the final values
allocate(lencounter(nenc))
allocate(lvdotr(nenc))

! call timer%reset()
! call timer%start()
call encounter_check_all(nenc, index1, index2, x, v, x, v, renc, renc, dt, lencounter, lvdotr)
! call timer%stop()
! call timer%report(nsubsteps=1, message="Narrow plpl:")

! call timer%reset()
! call timer%start()
call encounter_check_reduce_broadphase(npl, nenc, index1, index2, lencounter, lvdotr)
! call timer%stop()
! call timer%report(nsubsteps=1, message="Reduce plpl:")
end if

! call timer%stop()
! call timer%report(nsubsteps=1, message="Sort & Sweep plpl:")
return
end subroutine encounter_check_all_sort_and_sweep_plpl

Expand Down Expand Up @@ -444,6 +459,8 @@ subroutine encounter_check_all_sort_and_sweep_plplm(nplm, nplt, xplm, vplm, xplt

! If this is the first time through, build the index lists
if ((nplm == 0) .or. (nplt == 0)) return
! call timer%reset()
! call timer%start()

ntot = nplm + nplt
n = 2 * ntot
Expand All @@ -454,8 +471,7 @@ subroutine encounter_check_all_sort_and_sweep_plplm(nplm, nplt, xplm, vplm, xplt
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 @@ -483,9 +499,13 @@ subroutine encounter_check_all_sort_and_sweep_plplm(nplm, nplt, xplm, vplm, xplt
end do
!$omp end parallel do
! call timer%stop()
! call timer%report(nsubsteps=nthreads, message="AABB sort plplm:")
! call timer%report(nsubsteps=1, message="Sort plplm:")

! call timer%reset()
! call timer%start()
call boundingbox%sweep(nplm, nplt, nenc, index1, index2)
! call timer%stop()
! call timer%report(nsubsteps=1, message="Sweep plplm:")

if (nenc > 0) then
! Shift tiny body indices back into the range of the input position and velocity arrays
Expand All @@ -495,15 +515,25 @@ subroutine encounter_check_all_sort_and_sweep_plplm(nplm, nplt, xplm, vplm, xplt
allocate(lencounter(nenc))
allocate(lvdotr(nenc))

! call timer%reset()
! call timer%start()
call encounter_check_all(nenc, index1, index2, xplm, vplm, xplt, vplt, rencm, renct, dt, lencounter, lvdotr)
! call timer%stop()
! call timer%report(nsubsteps=1, message="Narrow plplm:")

! call timer%reset()
! call timer%start()

! Shift the tiny body indices back to their natural range
index2(:) = index2(:) + nplm

call encounter_check_reduce_broadphase(ntot, nenc, index1, index2, lencounter, lvdotr)
! call timer%stop()
! call timer%report(nsubsteps=1, message="Reduce plplm:")

end if

! call timer%stop()
! call timer%report(nsubsteps=1, message="Sort & Sweep plplm:")
return
end subroutine encounter_check_all_sort_and_sweep_plplm

Expand Down Expand Up @@ -597,7 +627,6 @@ end subroutine encounter_check_all_sort_and_sweep_pltp


subroutine encounter_check_all_triangular_one(i, n, xi, yi, zi, vxi, vyi, vzi, x, y, z, vx, vy, vz, renci, renc, dt, ind_arr, lenci)
!$omp declare simd(encounter_check_all_triangular_one)
implicit none
! Arguments
integer(I4B), intent(in) :: i
Expand Down Expand Up @@ -660,6 +689,10 @@ subroutine encounter_check_all_triangular_plpl(npl, x, v, renc, dt, lvdotr, inde
logical, dimension(npl) :: lencounteri, lvdotri
integer(I4B), dimension(:), allocatable, save :: ind_arr
type(encounter_list), dimension(npl) :: lenc
type(walltimer) :: timer

! call timer%reset()
! call timer%start()

call util_index_array(ind_arr, npl)

Expand All @@ -677,6 +710,9 @@ subroutine encounter_check_all_triangular_plpl(npl, x, v, renc, dt, lvdotr, inde

call encounter_check_collapse_ragged_list(lenc, npl, nenc, index1, index2, lvdotr)

! call timer%stop()
! call timer%report(nsubsteps=1, message="Triangular plpl:")

return
end subroutine encounter_check_all_triangular_plpl

Expand Down Expand Up @@ -707,6 +743,10 @@ subroutine encounter_check_all_triangular_plplm(nplm, nplt, xplm, vplm, xplt, vp
logical, dimension(nplt) :: lencounteri, lvdotri
integer(I4B), dimension(:), allocatable, save :: ind_arr
type(encounter_list), dimension(nplm) :: lenc
type(walltimer) :: timer

! call timer%reset()
! call timer%start()

call util_index_array(ind_arr, nplt)

Expand All @@ -724,6 +764,9 @@ subroutine encounter_check_all_triangular_plplm(nplm, nplt, xplm, vplm, xplt, vp

call encounter_check_collapse_ragged_list(lenc, nplm, nenc, index1, index2, lvdotr)

! call timer%stop()
! call timer%report(nsubsteps=1, message="Triangular plplm:")

return
end subroutine encounter_check_all_triangular_plplm

Expand Down Expand Up @@ -909,23 +952,19 @@ module subroutine encounter_check_sweep_aabb_double_list(self, n1, n2, nenc, ind
integer(I4B) :: i, k, ntot
type(encounter_list), dimension(n1+n2) :: lenc !! Array of encounter lists (one encounter list per body)
integer(I4B), dimension(:), allocatable, save :: ind_arr
type(walltimer) :: timer

ntot = n1 + n2
call util_index_array(ind_arr, ntot)
! 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 timer%reset()
! call timer%start()

!$omp parallel do default(private) schedule(static)&
!$omp shared(self, lenc, ind_arr) &
!$omp firstprivate(ntot, n1, n2)
do i = 1, ntot
call encounter_check_sweep_aabb_one_double_list(i, 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(i))
end do
!$omp end parallel do
! call timer%stop()
! call timer%report(nsubsteps=nthreads, message="AABB sweep double:")

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

Expand Down Expand Up @@ -955,24 +994,19 @@ module subroutine encounter_check_sweep_aabb_single_list(self, n, nenc, index1,
!Internals
Integer(I4B) :: i, k
type(encounter_list), dimension(n) :: lenc !! Array of encounter lists (one encounter list per body)
type(walltimer) :: timer
integer(I4B), dimension(:), allocatable, save :: ind_arr

call util_index_array(ind_arr, n)

! 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 timer%reset()
! call timer%start()
!$omp parallel do default(private) schedule(static)&
!$omp shared(self, lenc, ind_arr) &
!$omp firstprivate(n)
do i = 1, n
call encounter_check_sweep_aabb_one_single_list(i, n, self%aabb(1)%ind(:), self%aabb(1)%ibeg(:), self%aabb(1)%iend(:), self%aabb(2)%ibeg(:), self%aabb(2)%iend(:), ind_arr(:), lenc(i))
end do
!$omp end parallel do
! call timer%stop()
! call timer%report(nsubsteps=nthreads, message="AABB sweep single:")

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

Expand Down
6 changes: 6 additions & 0 deletions src/symba/symba_encounter_check.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -24,6 +24,7 @@ module function symba_encounter_check_pl(self, param, system, dt, irec) result(l
integer(I4B), dimension(:,:), allocatable :: k_plpl_enc
type(interaction_timer), save :: itimer
logical, save :: lfirst = .true.
type(walltimer) :: timer

if (self%nbody == 0) return

Expand All @@ -32,11 +33,16 @@ module function symba_encounter_check_pl(self, param, system, dt, irec) result(l
npl = pl%nbody
nplm = pl%nplm
nplt = npl - nplm

! call timer%reset()
! call timer%start()
if (nplt == 0) then
call encounter_check_all_plpl(param, npl, pl%xh, pl%vh, pl%renc, dt, lvdotr, index1, index2, nenc)
else
call encounter_check_all_plplm(param, nplm, nplt, pl%xh(:,1:nplm), pl%vh(:,1:nplm), pl%xh(:,nplm+1:npl), pl%vh(:,nplm+1:npl), pl%renc(1:nplm), pl%renc(nplm+1:npl), dt, lvdotr, index1, index2, nenc)
end if
! call timer%stop()
! call timer%report(nsubsteps=nthreads, message="Encounter Check Total:")
lany_encounter = nenc > 0
if (lany_encounter) then
call plplenc_list%resize(nenc)
Expand Down

0 comments on commit bdfc400

Please sign in to comment.