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

Commit

Permalink
Added ability to switch between flat and triangular versions of doubl…
Browse files Browse the repository at this point in the history 
…e loops
  • Loading branch information
@daminton
daminton committed Sep 16, 2021
1 parent ea1f5f8 commit cc97508
Show file tree
Hide file tree
Showing 7 changed files with 159 additions and 90 deletions.
15 changes: 10 additions & 5 deletions src/kick/kick.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -14,7 +14,11 @@ module subroutine kick_getacch_int_pl(self, param)
class(swiftest_pl), intent(inout) :: self !! Swiftest massive body object
class(swiftest_parameters), intent(in) :: param !! Current swiftest run configuration parameters

call kick_getacch_int_all_flat_pl(self%nbody, self%nplpl, self%k_plpl, self%xh, self%Gmass, self%radius, self%ah)
if (param%lflatten_interactions) then
call kick_getacch_int_all_flat_pl(self%nbody, self%nplpl, self%k_plpl, self%xh, self%Gmass, self%radius, self%ah)
else
call kick_getacch_int_all_triangular_pl(self%nbody, self%nbody, self%xh, self%Gmass, self%radius, self%ah)
end if

return
end subroutine kick_getacch_int_pl
Expand Down Expand Up @@ -94,7 +98,7 @@ module subroutine kick_getacch_int_all_flat_pl(npl, nplpl, k_plpl, x, Gmass, rad
end subroutine kick_getacch_int_all_flat_pl


module subroutine kick_getacch_int_all_triangular_pl(npl, x, Gmass, radius, acc)
module subroutine kick_getacch_int_all_triangular_pl(npl, nplm, x, Gmass, radius, acc)
!! author: David A. Minton
!!
!! Compute direct cross (third) term heliocentric accelerations for massive bodies, with parallelization.
Expand All @@ -103,7 +107,8 @@ module subroutine kick_getacch_int_all_triangular_pl(npl, x, Gmass, radius, acc)
!! Adapted from Hal Levison's Swift routine getacch_ah3.f
!! Adapted from David E. Kaufmann's Swifter routine whm_kick_getacch_ah3.f90 and helio_kick_getacch_int.f9
implicit none
integer(I4B), intent(in) :: npl !! Number of massive bodies
integer(I4B), intent(in) :: npl !! Total number of massive bodies
integer(I4B), intent(in) :: nplm !! Number of fully interacting massive bodies
real(DP), dimension(:,:), intent(in) :: x !! Position vector array
real(DP), dimension(:), intent(in) :: Gmass !! Array of massive body G*mass
real(DP), dimension(:), intent(in) :: radius !! Array of massive body radii
Expand All @@ -118,11 +123,11 @@ module subroutine kick_getacch_int_all_triangular_pl(npl, x, Gmass, radius, acc)
ahj(:,:) = 0.0_DP

!$omp parallel do default(private) schedule(static)&
!$omp shared(npl, x, Gmass, radius) &
!$omp shared(npl, nplm, x, Gmass, radius) &
!$omp lastprivate(rji2, rlim2, xr, yr, zr) &
!$omp reduction(+:ahi) &
!$omp reduction(-:ahj)
do i = 1, npl
do i = 1, nplm
do concurrent(j = i+1:npl)
xr = x(1, j) - x(1, i)
yr = x(2, j) - x(2, i)
Expand Down
5 changes: 3 additions & 2 deletions src/modules/swiftest_classes.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -867,9 +867,10 @@ module subroutine kick_getacch_int_all_flat_pl(npl, nplpl, k_plpl, x, Gmass, rad
real(DP), dimension(:,:), intent(inout) :: acc !! Acceleration vector array
end subroutine kick_getacch_int_all_flat_pl

module subroutine kick_getacch_int_all_triangular_pl(npl, x, Gmass, radius, acc)
module subroutine kick_getacch_int_all_triangular_pl(npl, nplm, x, Gmass, radius, acc)
implicit none
integer(I4B), intent(in) :: npl !! Number of massive bodies
integer(I4B), intent(in) :: npl !! Total number of massive bodies
integer(I4B), intent(in) :: nplm !! Number of fully interacting massive bodies
real(DP), dimension(:,:), intent(in) :: x !! Position vector array
real(DP), dimension(:), intent(in) :: Gmass !! Array of massive body G*mass
real(DP), dimension(:), intent(in) :: radius !! Array of massive body radii
Expand Down
136 changes: 77 additions & 59 deletions src/symba/symba_encounter_check.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -8,13 +8,15 @@ subroutine symba_encounter_check_all_flat(nplplm, k_plpl, x, v, rhill, dt, irec
!! Check for encounters between massive bodies. Split off from the main subroutine for performance.
!! This is the flat (single loop) version.
implicit none
integer(I8B), intent(in) :: nplplm
integer(I4B), dimension(:,:), intent(in) :: k_plpl
real(DP), dimension(:,:), intent(in) :: x, v
real(DP), dimension(:), intent(in) :: rhill
real(DP), intent(in) :: dt
integer(I4B), intent(in) :: irec
logical, dimension(:), intent(out) :: lencounter, loc_lvdotr
integer(I8B), intent(in) :: nplplm !! Total number of plm-pl interactions
integer(I4B), dimension(:,:), intent(in) :: k_plpl !! List of all pl-pl interactions
real(DP), dimension(:,:), intent(in) :: x !! Position vectors of massive bodies
real(DP), dimension(:,:), intent(in) :: v !! Velocity vectors of massive bodies
real(DP), dimension(:), intent(in) :: rhill !! Hill's radii of massive bodies
real(DP), intent(in) :: dt !! Step size
integer(I4B), intent(in) :: irec !! Current recursion depth
logical, dimension(:), intent(out) :: lencounter !! Logical array indicating which pair is in an encounter state
logical, dimension(:), intent(out) :: loc_lvdotr !! Logical array indicating the sign of v .dot. x for each encounter
! Internals
integer(I8B) :: k
integer(I4B) :: i, j
Expand Down Expand Up @@ -42,36 +44,40 @@ subroutine symba_encounter_check_all_flat(nplplm, k_plpl, x, v, rhill, dt, irec
end subroutine symba_encounter_check_all_flat


subroutine symba_encounter_check_all_triangular(npl, nplm, x, v, rhill, dt, irec, loc_lvdotr, k_plplenc, nenc)
subroutine symba_encounter_check_all_triangular(npl, nplm, x, v, rhill, dt, irec, lvdotr, index1, index2, nenc)
!! author: David A. Minton
!!
!! Check for encounters between massive bodies. Split off from the main subroutine for performance
!! This is the upper triangular (double loop) version.
implicit none
integer(I4B), intent(in) :: npl, nplm
real(DP), dimension(:,:), intent(in) :: x, v
real(DP), dimension(:), intent(in) :: rhill
real(DP), intent(in) :: dt
integer(I4B), intent(in) :: irec
logical, dimension(:), allocatable, intent(out) :: loc_lvdotr
integer(I4B), dimension(:,:), allocatable, intent(out) :: k_plplenc
integer(I4B), intent(out) :: nenc
integer(I4B), intent(in) :: npl !! Total number of massive bodies
integer(I4B), intent(in) :: nplm !! Number of fully interacting massive bodies
real(DP), dimension(:,:), intent(in) :: x !! Position vectors of massive bodies
real(DP), dimension(:,:), intent(in) :: v !! Velocity vectors of massive bodies
real(DP), dimension(:), intent(in) :: rhill !! Hill's radii of massive bodies
real(DP), intent(in) :: dt !! Step size
integer(I4B), intent(in) :: irec !! Current recursion depth
logical, dimension(:), allocatable, intent(out) :: lvdotr !! Logical flag indicating the sign of v .dot. x
integer(I4B), dimension(:), allocatable, intent(out) :: index1 !! List of indices for body 1 in each interaction
integer(I4B), dimension(:), allocatable, intent(out) :: index2 !! List of indices for body 2 in each interaction
integer(I4B), intent(out) :: nenc !! Total number of interactions
! Internals
integer(I4B) :: i, j, nenci, j0, j1
real(DP) :: xr, yr, zr, vxr, vyr, vzr, rhill1, rhill2
logical, dimension(npl) :: lencounteri, loc_lvdotri
logical, dimension(npl) :: lencounteri, lvdotri
integer(I4B), dimension(npl) :: ind_arr
type lenctype
logical, dimension(:), allocatable :: loc_lvdotr
logical, dimension(:), allocatable :: lvdotr
integer(I4B), dimension(:), allocatable :: index2
integer(I4B) :: nenc
end type
type(lenctype), dimension(nplm) :: lenc

ind_arr(:) = [(i, i = 1, npl)]
!$omp parallel do simd default(private) schedule(static)&
!$omp parallel do default(private) schedule(static)&
!$omp shared(npl, nplm, x, v, rhill, dt, irec, lenc) &
!$omp lastprivate(xr, yr, zr, vxr, vyr, vzr, rhill1, rhill2, lencounteri, loc_lvdotri)
!$omp firstprivate(ind_arr) &
!$omp lastprivate(xr, yr, zr, vxr, vyr, vzr, rhill1, rhill2, lencounteri, lvdotri)
do i = 1, nplm
do concurrent(j = i+1:npl)
xr = x(1, j) - x(1, i)
Expand All @@ -82,30 +88,31 @@ subroutine symba_encounter_check_all_triangular(npl, nplm, x, v, rhill, dt, ire
vzr = v(3, j) - v(3, i)
rhill1 = rhill(i)
rhill2 = rhill(j)
call symba_encounter_check_one(xr, yr, zr, vxr, vyr, vzr, rhill1, rhill2, dt, irec, lencounteri(j), loc_lvdotri(j))
call symba_encounter_check_one(xr, yr, zr, vxr, vyr, vzr, rhill1, rhill2, dt, irec, lencounteri(j), lvdotri(j))
end do
lenc(i)%nenc = count(lencounteri(i+1:npl))
if (lenc(i)%nenc > 0) then
allocate(lenc(i)%loc_lvdotr(nenci), lenc(i)%index2(nenci))
lenc(i)%loc_lvdotr(:) = pack(loc_lvdotri(i+1:npl), lencounteri(i+1:npl))
allocate(lenc(i)%lvdotr(nenci), lenc(i)%index2(nenci))
lenc(i)%lvdotr(:) = pack(lvdotri(i+1:npl), lencounteri(i+1:npl))
lenc(i)%index2(:) = pack(ind_arr(i+1:npl), lencounteri(i+1:npl))
end if
end do
!$omp end parallel do simd
!$omp end parallel do

associate(nenc_arr => lenc(:)%nenc)
nenc = sum(nenc_arr(1:nplm))
end associate
if (nenc > 0) then
allocate(loc_lvdotr(nenc))
allocate(k_plplenc(2,nenc))
allocate(lvdotr(nenc))
allocate(index1(nenc))
allocate(index2(nenc))
j0 = 1
do i = 1, nplm
if (lenc(i)%nenc > 0) then
j1 = j0 + lenc(i)%nenc - 1
loc_lvdotr(j0:j1) = lenc(i)%loc_lvdotr(:)
k_plplenc(1,j0:j1) = i
k_plplenc(2,j0:j1) = lenc(i)%index2(:)
lvdotr(j0:j1) = lenc(i)%lvdotr(:)
index1(j0:j1) = i
index2(j0:j1) = lenc(i)%index2(:)
j0 = j1 + 1
end if
end do
Expand All @@ -131,25 +138,26 @@ module function symba_encounter_check_pl(self, param, system, dt, irec) result(l
logical :: lany_encounter !! Returns true if there is at least one close encounter
! Internals
integer(I8B) :: k, nplplm, kenc
integer(I4B) :: i, j, nenc, npl
integer(I4B) :: i, j, nenc, npl, nplm
logical, dimension(:), allocatable :: lencounter, loc_lvdotr, lvdotr
integer(I4B), dimension(:), allocatable :: index1, index2
integer(I4B), dimension(:,:), allocatable :: k_plpl_enc

if (self%nbody == 0) return

associate(pl => self)
nplplm = pl%nplplm
associate(pl => self, plplenc_list => system%plplenc_list)
npl = pl%nbody
allocate(lencounter(nplplm))
allocate(loc_lvdotr(nplplm))
if (param%lflatten_interactions) then
nplplm = pl%nplplm
allocate(lencounter(nplplm))
allocate(loc_lvdotr(nplplm))

call symba_encounter_check_all_flat(nplplm, pl%k_plpl, pl%xh, pl%vh, pl%rhill, dt, irec, lencounter, loc_lvdotr)
call symba_encounter_check_all_flat(nplplm, pl%k_plpl, pl%xh, pl%vh, pl%rhill, dt, irec, lencounter, loc_lvdotr)

nenc = count(lencounter(:))
nenc = count(lencounter(:))

lany_encounter = nenc > 0
if (lany_encounter) then
associate(plplenc_list => system%plplenc_list)
lany_encounter = nenc > 0
if (lany_encounter) then
call plplenc_list%resize(nenc)
allocate(lvdotr(nenc))
allocate(index1(nenc))
Expand All @@ -160,28 +168,38 @@ module function symba_encounter_check_pl(self, param, system, dt, irec) result(l
deallocate(lencounter, loc_lvdotr)
call move_alloc(lvdotr, plplenc_list%lvdotr)
call move_alloc(index1, plplenc_list%index1)
call move_alloc(index2, plplenc_list%index2)
do k = 1, nenc
i = plplenc_list%index1(k)
j = plplenc_list%index2(k)
call util_index_eucl_ij_to_k(npl, i, j, kenc)
plplenc_list%kidx(k) = kenc
plplenc_list%id1(k) = pl%id(plplenc_list%index1(k))
plplenc_list%id2(k) = pl%id(plplenc_list%index2(k))
plplenc_list%status(k) = ACTIVE
plplenc_list%level(k) = irec
pl%lencounter(i) = .true.
pl%lencounter(j) = .true.
pl%levelg(i) = irec
pl%levelm(i) = irec
pl%levelg(j) = irec
pl%levelm(j) = irec
pl%nplenc(i) = pl%nplenc(i) + 1
pl%nplenc(j) = pl%nplenc(j) + 1
end do
end associate
call move_alloc(index2, plplenc_list%index2)
end if
else
nplm = pl%nplm
call symba_encounter_check_all_triangular(npl, nplm, pl%xh, pl%vh, pl%rhill, dt, irec, lvdotr, index1, index2, nenc)
call move_alloc(lvdotr, plplenc_list%lvdotr)
call move_alloc(index1, plplenc_list%index1)
call move_alloc(index2, plplenc_list%index2)
end if

if (lany_encounter) then
do k = 1, nenc
i = plplenc_list%index1(k)
j = plplenc_list%index2(k)
call util_index_eucl_ij_to_k(npl, i, j, kenc)
plplenc_list%kidx(k) = kenc
plplenc_list%id1(k) = pl%id(plplenc_list%index1(k))
plplenc_list%id2(k) = pl%id(plplenc_list%index2(k))
plplenc_list%status(k) = ACTIVE
plplenc_list%level(k) = irec
pl%lencounter(i) = .true.
pl%lencounter(j) = .true.
pl%levelg(i) = irec
pl%levelm(i) = irec
pl%levelg(j) = irec
pl%levelm(j) = irec
pl%nplenc(i) = pl%nplenc(i) + 1
pl%nplenc(j) = pl%nplenc(j) + 1
end do
end if
end associate

return
end function symba_encounter_check_pl

Expand Down
11 changes: 7 additions & 4 deletions src/symba/symba_kick.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -14,7 +14,11 @@ module subroutine symba_kick_getacch_int_pl(self, param)
class(symba_pl), intent(inout) :: self !! SyMBA massive body object
class(swiftest_parameters), intent(in) :: param !! Current swiftest run configuration parameter

call kick_getacch_int_all_flat_pl(self%nbody, self%nplplm, self%k_plpl, self%xh, self%Gmass, self%radius, self%ah)
if (param%lflatten_interactions) then
call kick_getacch_int_all_flat_pl(self%nbody, self%nplplm, self%k_plpl, self%xh, self%Gmass, self%radius, self%ah)
else
call kick_getacch_int_all_triangular_pl(self%nbody, self%nplm, self%xh, self%Gmass, self%radius, self%ah)
end if

return
end subroutine symba_kick_getacch_int_pl
Expand Down Expand Up @@ -44,15 +48,15 @@ module subroutine symba_kick_getacch_pl(self, system, param, t, lbeg)
if (self%nbody == 0) return
select type(system)
class is (symba_nbody_system)
associate(pl => self, npl => self%nbody, plplenc_list => system%plplenc_list, radius => self%radius)
associate(pl => self, npl => self%nbody, nplm => self%nplm, plplenc_list => system%plplenc_list, radius => self%radius)
! Apply kicks to all bodies (including those in the encounter list)
call helio_kick_getacch_pl(pl, system, param, t, lbeg)
if (plplenc_list%nenc > 0) then
! Remove kicks from bodies involved currently in the encounter list, as these are dealt with separately.
ah_enc(:,:) = 0.0_DP
nplplenc = int(plplenc_list%nenc, kind=I8B)
allocate(k_plpl_enc(2,nplplenc))
k_plpl_enc(:,1:nplplenc) = pl%k_plpl(:,plplenc_list%kidx(1:nplplenc))
ah_enc(:,:) = 0.0_DP
call kick_getacch_int_all_flat_pl(npl, nplplenc, k_plpl_enc, pl%xh, pl%Gmass, pl%radius, ah_enc)
pl%ah(:,1:npl) = pl%ah(:,1:npl) - ah_enc(:,1:npl)
end if
Expand All @@ -63,7 +67,6 @@ module subroutine symba_kick_getacch_pl(self, system, param, t, lbeg)
return
end subroutine symba_kick_getacch_pl


module subroutine symba_kick_getacch_tp(self, system, param, t, lbeg)
!! author: David A. Minton
!!
Expand Down
19 changes: 10 additions & 9 deletions src/symba/symba_util.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -292,18 +292,19 @@ module subroutine symba_util_index_eucl_plpl(self, param)
npl = int(self%nbody, kind=I8B)
nplm = count(.not. pl%lmtiny(1:npl))
pl%nplm = int(nplm, kind=I4B)

pl%nplpl = (npl * (npl - 1) / 2) ! number of entries in a strict lower triangle, npl x npl, minus first column
pl%nplplm = nplm * npl - nplm * (nplm + 1) / 2 ! number of entries in a strict lower triangle, npl x npl, minus first column including only mutually interacting bodies
if (allocated(self%k_plpl)) deallocate(self%k_plpl) ! Reset the index array if it's been set previously
allocate(self%k_plpl(2, pl%nplpl))
do concurrent (i = 1:npl)
do concurrent (j = i+1:npl)
call util_index_eucl_ij_to_k(npl, i, j, k)
self%k_plpl(1, k) = i
self%k_plpl(2, k) = j
if (param%lflatten_interactions) then
if (allocated(self%k_plpl)) deallocate(self%k_plpl) ! Reset the index array if it's been set previously
allocate(self%k_plpl(2, pl%nplpl))
do concurrent (i = 1:npl)
do concurrent (j = i+1:npl)
call util_index_eucl_ij_to_k(npl, i, j, k)
self%k_plpl(1, k) = i
self%k_plpl(2, k) = j
end do
end do
end do
end if
end associate

return
Expand Down
Loading

0 comments on commit cc97508

Please sign in to comment.