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

Commit

Permalink
Setting the stage to have the option of using either the flat or tria…
Browse files Browse the repository at this point in the history 
…ngular versions of the two big interaction loops.
  • Loading branch information
@daminton
daminton committed Sep 16, 2021
1 parent e2e1cde commit 25385d2
Show file tree
Hide file tree
Showing 4 changed files with 144 additions and 12 deletions.
57 changes: 53 additions & 4 deletions src/kick/kick.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -13,7 +13,7 @@ module subroutine kick_getacch_int_pl(self)
! Arguments
class(swiftest_pl), intent(inout) :: self !! Swiftest massive body object

call kick_getacch_int_all_pl(self%nbody, self%nplpl, self%k_plpl, self%xh, self%Gmass, self%radius, self%ah)
call kick_getacch_int_all_flat_pl(self%nbody, self%nplpl, self%k_plpl, self%xh, self%Gmass, self%radius, self%ah)

return
end subroutine kick_getacch_int_pl
Expand Down Expand Up @@ -41,10 +41,11 @@ module subroutine kick_getacch_int_tp(self, GMpl, xhp, npl)
end subroutine kick_getacch_int_tp


module subroutine kick_getacch_int_all_pl(npl, nplpl, k_plpl, x, Gmass, radius, acc)
module subroutine kick_getacch_int_all_flat_pl(npl, nplpl, k_plpl, x, Gmass, radius, acc)
!! author: David A. Minton
!!
!! Compute direct cross (third) term heliocentric accelerations for massive bodies, with parallelization
!! Compute direct cross (third) term heliocentric accelerations for massive bodies, with parallelization.
!! This is the flattened (single loop) version that uses the k_plpl interaction pair index array
!!
!! 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
Expand Down Expand Up @@ -88,7 +89,55 @@ module subroutine kick_getacch_int_all_pl(npl, nplpl, k_plpl, x, Gmass, radius,
end do

return
end subroutine kick_getacch_int_all_pl
end subroutine kick_getacch_int_all_flat_pl


module subroutine kick_getacch_int_all_triangular_pl(npl, x, Gmass, radius, acc)
!! author: David A. Minton
!!
!! Compute direct cross (third) term heliocentric accelerations for massive bodies, with parallelization.
!! This is the upper triangular matrix (double loop) version.
!!
!! 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
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
real(DP), dimension(:,:), intent(inout) :: acc !! Acceleration vector array
! Internals
real(DP), dimension(NDIM,npl) :: ahi, ahj
integer(I4B) :: i, j
real(DP) :: rji2, rlim2
real(DP) :: xr, yr, zr

ahi(:,:) = 0.0_DP
ahj(:,:) = 0.0_DP

!$omp parallel do default(private) schedule(static)&
!$omp shared(npl, x, Gmass, radius) &
!$omp lastprivate(rji2, rlim2, xr, yr, zr) &
!$omp reduction(+:ahi) &
!$omp reduction(-:ahj)
do i = 1, npl
do concurrent(j = i+1:npl)
xr = x(1, j) - x(1, i)
yr = x(2, j) - x(2, i)
zr = x(3, j) - x(3, i)
rji2 = xr**2 + yr**2 + zr**2
rlim2 = (radius(i) + radius(j))**2
if (rji2 > rlim2) call kick_getacch_int_one_pl(rji2, xr, yr, zr, Gmass(i), Gmass(j), ahi(1,i), ahi(2,i), ahi(3,i), ahj(1,j), ahj(2,j), ahj(3,j))
end do
end do
!$omp end parallel do

do concurrent(i = 1:npl)
acc(:,i) = acc(:,i) + ahi(:,i) + ahj(:,i)
end do

return
end subroutine kick_getacch_int_all_triangular_pl


module subroutine kick_getacch_int_all_tp(ntp, npl, xtp, xpl, GMpl, lmask, acc)
Expand Down
13 changes: 11 additions & 2 deletions src/modules/swiftest_classes.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -853,7 +853,7 @@ module subroutine kick_getacch_int_tp(self, GMpl, xhp, npl)
integer(I4B), intent(in) :: npl !! Number of active massive bodies
end subroutine kick_getacch_int_tp

module subroutine kick_getacch_int_all_pl(npl, nplpl, k_plpl, x, Gmass, radius, acc)
module subroutine kick_getacch_int_all_flat_pl(npl, nplpl, k_plpl, x, Gmass, radius, acc)
implicit none
integer(I4B), intent(in) :: npl !! Number of massive bodies
integer(I8B), intent(in) :: nplpl !! Number of massive body interactions to compute
Expand All @@ -862,7 +862,16 @@ module subroutine kick_getacch_int_all_pl(npl, nplpl, k_plpl, x, Gmass, radius,
real(DP), dimension(:), intent(in) :: Gmass !! Array of massive body G*mass
real(DP), dimension(:), intent(in) :: radius !! Array of massive body radii
real(DP), dimension(:,:), intent(inout) :: acc !! Acceleration vector array
end subroutine kick_getacch_int_all_pl
end subroutine kick_getacch_int_all_flat_pl

module subroutine kick_getacch_int_all_triangular_pl(npl, x, Gmass, radius, acc)
implicit none
integer(I4B), intent(in) :: npl !! Number of 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
real(DP), dimension(:,:), intent(inout) :: acc !! Acceleration vector array
end subroutine kick_getacch_int_all_triangular_pl

module subroutine kick_getacch_int_all_tp(ntp, npl, xtp, xpl, GMpl, lmask, acc)
implicit none
Expand Down
82 changes: 78 additions & 4 deletions src/symba/symba_encounter_check.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,10 +2,11 @@
use swiftest
contains

subroutine symba_encounter_check_all(nplplm, k_plpl, x, v, rhill, dt, irec, lencounter, loc_lvdotr)
subroutine symba_encounter_check_all_flat(nplplm, k_plpl, x, v, rhill, dt, irec, lencounter, loc_lvdotr)
!! author: David A. Minton
!!
!! Check for encounters between massive bodies. Split off from the main subroutine for performance
!! 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
Expand Down Expand Up @@ -38,7 +39,80 @@ subroutine symba_encounter_check_all(nplplm, k_plpl, x, v, rhill, dt, irec, len
!$omp end parallel do simd

return
end subroutine symba_encounter_check_all
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)
!! 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
! Internals
integer(I4B) :: i, j, nenci, j0, j1
real(DP) :: xr, yr, zr, vxr, vyr, vzr, rhill1, rhill2
logical, dimension(npl) :: lencounteri, loc_lvdotri
integer(I4B), dimension(npl) :: ind_arr
type lenctype
logical, dimension(:), allocatable :: loc_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 shared(npl, nplm, x, v, rhill, dt, irec, lenc) &
!$omp lastprivate(xr, yr, zr, vxr, vyr, vzr, rhill1, rhill2, lencounteri, loc_lvdotri)
do i = 1, nplm
do concurrent(j = i+1:npl)
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)
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))
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))
lenc(i)%index2(:) = pack(ind_arr(i+1:npl), lencounteri(i+1:npl))
end if
end do
!$omp end parallel do simd

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))
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(:)
j0 = j1 + 1
end if
end do
end if

return
end subroutine symba_encounter_check_all_triangular


module function symba_encounter_check_pl(self, system, dt, irec) result(lany_encounter)
Expand Down Expand Up @@ -68,7 +142,7 @@ module function symba_encounter_check_pl(self, system, dt, irec) result(lany_enc
allocate(lencounter(nplplm))
allocate(loc_lvdotr(nplplm))

call symba_encounter_check_all(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(:))

Expand Down
4 changes: 2 additions & 2 deletions src/symba/symba_kick.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -13,7 +13,7 @@ module subroutine symba_kick_getacch_int_pl(self)
! Arguments
class(symba_pl), intent(inout) :: self

call kick_getacch_int_all_pl(self%nbody, self%nplplm, self%k_plpl, self%xh, self%Gmass, self%radius, self%ah)
call kick_getacch_int_all_flat_pl(self%nbody, self%nplplm, self%k_plpl, self%xh, self%Gmass, self%radius, self%ah)

return
end subroutine symba_kick_getacch_int_pl
Expand Down Expand Up @@ -52,7 +52,7 @@ module subroutine symba_kick_getacch_pl(self, system, param, t, lbeg)
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_pl(npl, nplplenc, k_plpl_enc, pl%xh, pl%Gmass, pl%radius, ah_enc)
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 Down

0 comments on commit 25385d2

Please sign in to comment.