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

Commit

Permalink
Consolidated interaction step into a single parallelized subroutine i…
Browse files Browse the repository at this point in the history 
…nstead of having a different one for SyMBA.
  • Loading branch information
@daminton
daminton committed Aug 23, 2021
1 parent ce023d3 commit e43ae9b
Show file tree
Hide file tree
Showing 4 changed files with 124 additions and 117 deletions.
6 changes: 3 additions & 3 deletions src/io/io.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1054,7 +1054,7 @@ module function io_read_frame_body(self, iu, param) result(ierr)
read(iu, err = 667, iomsg = errmsg) pl%Gmass(:)
pl%mass(:) = pl%Gmass(:) / param%GU
if (param%lrhill_present) read(iu, err = 667, iomsg = errmsg) pl%rhill(:)
if (param%lclose) read(iu, err = 667, iomsg = errmsg) pl%radius(:)
read(iu, err = 667, iomsg = errmsg) pl%radius(:)
if (param%lrotation) then
read(iu, err = 667, iomsg = errmsg) pl%Ip(1, :)
read(iu, err = 667, iomsg = errmsg) pl%Ip(2, :)
Expand All @@ -1081,7 +1081,7 @@ module function io_read_frame_body(self, iu, param) result(ierr)
end if
self%Gmass(i) = real(val, kind=DP)
self%mass(i) = real(val / param%GU, kind=DP)
if (param%lclose) read(iu, *, err = 667, iomsg = errmsg) self%radius(i)
read(iu, *, err = 667, iomsg = errmsg) self%radius(i)
class is (swiftest_tp)
read(iu, *, err = 667, iomsg = errmsg) self%id(i)
end select
Expand Down Expand Up @@ -1508,7 +1508,7 @@ module subroutine io_write_frame_body(self, iu, param)
class is (swiftest_pl) ! Additional output if the passed polymorphic object is a massive body
write(iu, err = 667, iomsg = errmsg) pl%Gmass(1:n)
if (param%lrhill_present) write(iu, err = 667, iomsg = errmsg) pl%rhill(1:n)
if (param%lclose) write(iu, err = 667, iomsg = errmsg) pl%radius(1:n)
write(iu, err = 667, iomsg = errmsg) pl%radius(1:n)
if (param%lrotation) then
write(iu, err = 667, iomsg = errmsg) pl%Ip(1, 1:n)
write(iu, err = 667, iomsg = errmsg) pl%Ip(2, 1:n)
Expand Down
148 changes: 87 additions & 61 deletions src/kick/kick.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -12,37 +12,8 @@ module subroutine kick_getacch_int_pl(self)
implicit none
! Arguments
class(swiftest_pl), intent(inout) :: self
! Internals
integer(I8B) :: k, nplpl
real(DP) :: rji2
real(DP) :: dx, dy, dz
integer(I4B) :: i, j
real(DP), dimension(:,:), pointer :: ah, xh
real(DP), dimension(NDIM,self%nbody) :: ahi, ahj
integer(I4B), dimension(:,:), pointer :: k_plpl
logical, dimension(:), pointer :: lmask
real(DP), dimension(:), pointer :: Gmass

associate(ah => self%ah, xh => self%xh, k_plpl => self%k_plpl, lmask => self%lmask, Gmass => self%Gmass)
nplpl = self%nplpl
ahi(:,:) = 0.0_DP
ahj(:,:) = 0.0_DP
!$omp parallel do default(shared)&
!$omp private(k, i, j, dx, dy, dz, rji2) &
!$omp reduction(+:ahi) &
!$omp reduction(-:ahj)
do k = 1_I8B, nplpl
i = k_plpl(1,k)
j = k_plpl(2,k)
dx = xh(1, j) - xh(1, i)
dy = xh(2, j) - xh(2, i)
dz = xh(3, j) - xh(3, i)
rji2 = dx**2 + dy**2 + dz**2
if (lmask(i) .and. lmask(j)) call kick_getacch_int_one_pl(rji2, dx, dy, dz, Gmass(i), Gmass(j), ahi(1,i), ahi(2,i), ahi(3,i), ahj(1,j), ahj(2,j), ahj(3,j))
end do
!$omp end parallel do
ah(:,1:self%nbody) = ah(:,1:self%nbody) + ahi(:,1:self%nbody) + ahj(:,1:self%nbody)
end associate

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

return
end subroutine kick_getacch_int_pl
Expand All @@ -63,7 +34,6 @@ module pure subroutine kick_getacch_int_tp(self, GMpl, xhp, npl)
integer(I4B), intent(in) :: npl !! Number of active massive bodies
! Internals
integer(I4B) :: i, j
!real(DP), dimension(:,:), allocatable :: aht

if ((self%nbody == 0) .or. (npl == 0)) return

Expand All @@ -73,60 +43,116 @@ module pure subroutine kick_getacch_int_tp(self, GMpl, xhp, npl)
do j = 1, npl
block
real(DP) :: rji2
real(DP) :: dx, dy, dz
dx = tp%xh(1, i) - xhp(1, j)
dy = tp%xh(2, i) - xhp(1, j)
dz = tp%xh(3, i) - xhp(1, j)
rji2 = dx**2 + dy**2 + dz**2
call kick_getacch_int_one_tp(rji2, dx, dy, dz, GMpl(i), tp%ah(1,i), tp%ah(2,i), tp%ah(3,i))
real(DP) :: xr, yr, zr
xr = tp%xh(1, i) - xhp(1, j)
yr = tp%xh(2, i) - xhp(1, j)
zr = tp%xh(3, i) - xhp(1, j)
rji2 = xr**2 + yr**2 + zr**2
call kick_getacch_int_one_tp(rji2, xr, yr, zr, GMpl(i), tp%ah(1,i), tp%ah(2,i), tp%ah(3,i))
end block
end do
end if
end do
!call move_alloc(aht, tp%ah)
end associate

return
end subroutine kick_getacch_int_tp

module pure elemental subroutine kick_getacch_int_one_pl(rji2, dx, dy, dz, Gmi, Gmj, axi, ayi, azi, axj, ayj, azj)

module subroutine kick_getacch_int_all_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
!!
!! 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
real(DP), intent(in) :: rji2
real(DP), intent(in) :: dx, dy, dz
real(DP), intent(in) :: Gmi
real(DP), intent(in) :: Gmj
real(DP), intent(inout) :: axi, ayi, azi
real(DP), intent(inout) :: axj, ayj, azj
integer(I4B), intent(in) :: npl !! Number of massive bodies
integer(I8B), intent(in) :: nplpl !! Number of massive body interactions to compute
integer(I4B), dimension(:,:), intent(in) :: k_plpl !! Array of interaction pair indices (flattened upper triangular matrix)
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
integer(I8B) :: k
real(DP) :: rji2, rlim2
real(DP) :: xr, yr, zr
integer(I4B) :: i, j
real(DP), dimension(NDIM,npl) :: ahi, ahj

ahi(:,:) = 0.0_DP
ahj(:,:) = 0.0_DP
!$omp parallel do default(private)&
!$omp shared(nplpl, k_plpl, x, Gmass, radius) &
!$omp reduction(+:ahi) &
!$omp reduction(-:ahj)
do k = 1_I8B, nplpl
i = k_plpl(1,k)
j = k_plpl(2,k)
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
!$omp end parallel do
acc(:,1:npl) = acc(:,1:npl) + ahi(:,1:npl) + ahj(:,1:npl)
return
end subroutine kick_getacch_int_all_pl


module pure subroutine kick_getacch_int_one_pl(rji2, xr, yr, zr, Gmi, Gmj, axi, ayi, azi, axj, ayj, azj)
!! author: David A. Minton
!!
!! Compute direct cross (third) term heliocentric accelerations for a single pair of massive bodies
!!
!! 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
real(DP), intent(in) :: rji2 !! Square of distance between the two bodies
real(DP), intent(in) :: xr, yr, zr !! Distances between the two bodies in x, y, and z directions
real(DP), intent(in) :: Gmi !! G*mass of body i
real(DP), intent(in) :: Gmj !! G*mass of body j
real(DP), intent(inout) :: axi, ayi, azi !! Acceleration vector components of body i
real(DP), intent(inout) :: axj, ayj, azj !! Acceleration vector components of body j
! Internals
real(DP) :: faci, facj, irij3

irij3 = 1.0_DP / (rji2 * sqrt(rji2))
faci = Gmi * irij3
facj = Gmj * irij3
axi = axi + facj * dx
ayi = ayi + facj * dy
azi = azi + facj * dz
axj = axj - faci * dx
ayj = ayj - faci * dy
azj = azj - faci * dz
axi = axi + facj * xr
ayi = ayi + facj * yr
azi = azi + facj * zr
axj = axj - faci * xr
ayj = ayj - faci * yr
azj = azj - faci * zr

return
end subroutine kick_getacch_int_one_pl


!module pure elemental subroutine kick_getacch_int_one_tp(rji2, dx, dy, dz, GMpl, ax, ay, az)
module pure subroutine kick_getacch_int_one_tp(rji2, dx, dy, dz, GMpl, ax, ay, az)
module pure subroutine kick_getacch_int_one_tp(rji2, xr, yr, zr, GMpl, ax, ay, az)
!! author: David A. Minton
!!
!! Compute direct cross (third) term heliocentric accelerations of a single test particle massive body pair.
!!
!! Adapted from Hal Levison's Swift routine getacch_ah3_tp.f
!! Adapted from David E. Kaufmann's Swifter routine whm_kick_getacch_ah3.f90 and helio_kick_getacch_int_tp.f90
implicit none
real(DP), intent(in) :: rji2
real(DP), intent(in) :: dx, dy, dz
real(DP), intent(in) :: GMpl
real(DP), intent(inout) :: ax, ay, az
real(DP), intent(in) :: rji2 !! Square of distance between the test particle and massive body
real(DP), intent(in) :: xr, yr, zr !! Distances between the two bodies in x, y, and z directions
real(DP), intent(in) :: Gmpl !! G*mass of massive body
real(DP), intent(inout) :: ax, ay, az !! Acceleration vector components of test particle
! Internals
real(DP) :: fac

fac = GMpl / (rji2 * sqrt(rji2))
ax = ax - fac * dx
ay = ay - fac * dy
az = az - fac * dz
ax = ax - fac * xr
ay = ay - fac * yr
az = az - fac * zr
return
end subroutine kick_getacch_int_one_tp

Expand Down
41 changes: 26 additions & 15 deletions src/modules/swiftest_classes.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -709,29 +709,40 @@ module subroutine io_write_frame_system(self, iu, param)
class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters
end subroutine io_write_frame_system

module pure elemental subroutine kick_getacch_int_one_pl(rji2, dx, dy, dz, Gmi, Gmj, axi, ayi, azi, axj, ayj, azj)
implicit none
real(DP), intent(in) :: rji2
real(DP), intent(in) :: dx, dy, dz
real(DP), intent(in) :: Gmi
real(DP), intent(in) :: Gmj
real(DP), intent(inout) :: axi, ayi, azi
real(DP), intent(inout) :: axj, ayj, azj
module subroutine kick_getacch_int_all_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
integer(I4B), dimension(:,:), intent(in) :: k_plpl !! Array of interaction pair indices (flattened upper triangular matrix)
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_pl

module pure subroutine kick_getacch_int_one_pl(rji2, xr, yr, zr, Gmi, Gmj, axi, ayi, azi, axj, ayj, azj)
!$omp declare simd(kick_getacch_int_one_pl)
implicit none
real(DP), intent(in) :: rji2 !! Square of distance between the two bodies
real(DP), intent(in) :: xr, yr, zr !! Distances between the two bodies in x, y, and z directions
real(DP), intent(in) :: Gmi !! G*mass of body i
real(DP), intent(in) :: Gmj !! G*mass of body j
real(DP), intent(inout) :: axi, ayi, azi !! Acceleration vector components of body i
real(DP), intent(inout) :: axj, ayj, azj !! Acceleration vector components of body j
end subroutine kick_getacch_int_one_pl

!module pure elemental subroutine kick_getacch_int_one_tp(rji2, dx, dy, dz, Gmpl, ax, ay, az)
module pure subroutine kick_getacch_int_one_tp(rji2, dx, dy, dz, Gmpl, ax, ay, az)
module pure subroutine kick_getacch_int_one_tp(rji2, xr, yr, zr, Gmpl, ax, ay, az)
!$omp declare simd(kick_getacch_int_one_tp)
implicit none
real(DP), intent(in) :: rji2
real(DP), intent(in) :: dx, dy, dz
real(DP), intent(in) :: Gmpl
real(DP), intent(inout) :: ax, ay, az
real(DP), intent(in) :: rji2 !! Square of distance between the test particle and massive body
real(DP), intent(in) :: xr, yr, zr !! Distances between the two bodies in x, y, and z directions
real(DP), intent(in) :: Gmpl !! G*mass of massive body
real(DP), intent(inout) :: ax, ay, az !! Acceleration vector components of test particle
end subroutine kick_getacch_int_one_tp

module subroutine kick_getacch_int_pl(self)
implicit none
class(swiftest_pl), intent(inout) :: self
class(swiftest_pl), intent(inout) :: self !! Swiftest massive body object
end subroutine kick_getacch_int_pl

module pure subroutine kick_getacch_int_tp(self, GMpl, xhp, npl)
Expand Down
46 changes: 8 additions & 38 deletions src/symba/symba_kick.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -12,38 +12,8 @@ module subroutine symba_kick_getacch_int_pl(self)
implicit none
! Arguments
class(symba_pl), intent(inout) :: self
! Internals
integer(I8B) :: k, nplplm
real(DP) :: rji2, rlim2
real(DP) :: dx, dy, dz
integer(I4B) :: i, j
real(DP), dimension(:,:), pointer :: ah, xh
real(DP), dimension(NDIM,self%nbody) :: ahi, ahj
integer(I4B), dimension(:,:), pointer :: k_plpl
logical, dimension(:), pointer :: lmask
real(DP), dimension(:), pointer :: Gmass, radius

associate(ah => self%ah, xh => self%xh, k_plpl => self%k_plpl, lmask => self%lmask, Gmass => self%Gmass, radius => self%radius)
nplplm = self%nplplm
ahi(:,:) = 0.0_DP
ahj(:,:) = 0.0_DP
!$omp parallel do default(shared)&
!$omp private(k, i, j, dx, dy, dz, rji2, rlim2) &
!$omp reduction(+:ahi) &
!$omp reduction(-:ahj)
do k = 1_I8B, nplplm
i = k_plpl(1,k)
j = k_plpl(2,k)
dx = xh(1, j) - xh(1, i)
dy = xh(2, j) - xh(2, i)
dz = xh(3, j) - xh(3, i)
rji2 = dx**2 + dy**2 + dz**2
rlim2 = (radius(i)+radius(j))**2
if ((rji2 > rlim2) .and. lmask(i) .and. lmask(j)) call kick_getacch_int_one_pl(rji2, dx, dy, dz, Gmass(i), Gmass(j), ahi(1,i), ahi(2,i), ahi(3,i), ahj(1,j), ahj(2,j), ahj(3,j))
end do
!$omp end parallel do
ah(:,1:self%nbody) = ah(:,1:self%nbody) + ahi(:,1:self%nbody) + ahj(:,1:self%nbody)
end associate
call kick_getacch_int_all_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 All @@ -66,7 +36,7 @@ module subroutine symba_kick_getacch_pl(self, system, param, t, lbeg)
! Internals
integer(I4B) :: i, j
integer(I8B) :: k, nplplenc
real(DP) :: rji2, rlim2, dx, dy, dz
real(DP) :: rji2, rlim2, xr, yr, zr
real(DP), dimension(NDIM,self%nbody) :: ahi, ahj
class(symba_plplenc), pointer :: plplenc_list
real(DP), dimension(:), pointer :: Gmass, radius
Expand All @@ -82,18 +52,18 @@ module subroutine symba_kick_getacch_pl(self, system, param, t, lbeg)
ahi(:,:) = 0.0_DP
ahj(:,:) = 0.0_DP
!$omp parallel do default(shared)&
!$omp private(k, i, j, dx, dy, dz, rji2, rlim2) &
!$omp private(k, i, j, xr, yr, zr, rji2, rlim2) &
!$omp reduction(+:ahi) &
!$omp reduction(-:ahj)
do k = 1_I8B, nplplenc
i = plplenc_list%index1(k)
j = plplenc_list%index2(k)
dx = xh(1, j) - xh(1, i)
dy = xh(2, j) - xh(2, i)
dz = xh(3, j) - xh(3, i)
rji2 = dx**2 + dy**2 + dz**2
xr = xh(1, j) - xh(1, i)
yr = xh(2, j) - xh(2, i)
zr = xh(3, j) - xh(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, dx, dy, dz, Gmass(i), Gmass(j), ahi(1,i), ahi(2,i), ahi(3,i), ahj(1,j), ahj(2,j), ahj(3,j))
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
!$omp end parallel do
ah(:,1:self%nbody) = ah(:,1:self%nbody) - ahi(:,1:self%nbody) - ahj(:,1:self%nbody)
Expand Down

0 comments on commit e43ae9b

Please sign in to comment.