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

Commit

Permalink
Converted the cross term in the WHM accelration to use the flattened …
Browse files Browse the repository at this point in the history
…array. Tests pass
  • Loading branch information
daminton committed Jul 23, 2021
1 parent 5cb07aa commit e47a6f1
Showing 1 changed file with 7 additions and 8 deletions.
15 changes: 7 additions & 8 deletions src/whm/whm_getacch.f90
Original file line number Diff line number Diff line change
Expand Up @@ -187,28 +187,27 @@ pure subroutine whm_getacch_ah3(pl)
implicit none

class(whm_pl), intent(inout) :: pl
integer(I4B) :: i, j
integer(I4B) :: k
real(DP) :: rji2, irij3, faci, facj
real(DP), dimension(NDIM) :: dx
real(DP), dimension(:,:), allocatable :: ah3

associate(npl => pl%nbody)
associate(npl => pl%nbody, nplpl => pl%nplpl)
allocate(ah3, mold=pl%ah)
ah3(:, :) = 0.0_DP

do i = 1, npl - 1
do j = i + 1, npl
do k = 1, nplpl
associate(i => pl%k_eucl(1, k), j => pl%k_eucl(2, k))
dx(:) = pl%xh(:, j) - pl%xh(:, i)
rji2 = dot_product(dx(:), dx(:))
irij3 = 1.0_DP / (rji2 * sqrt(rji2))
faci = pl%Gmass(i) * irij3
facj = pl%Gmass(j) * irij3
ah3(:, i) = ah3(:, i) + facj * dx(:)
ah3(:, j) = ah3(:, j) - faci * dx(:)
end do
end associate
end do
do i = 1, NDIM
pl%ah(i, 1:npl) = pl%ah(i, 1:npl) + ah3(i, 1:npl)
do concurrent (k = 1:npl)
pl%ah(:, k) = pl%ah(:, k) + ah3(:, k)
end do
deallocate(ah3)
end associate
Expand Down

0 comments on commit e47a6f1

Please sign in to comment.