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

Commit

Permalink
Updated OpenMP directrives and restructured some loops
Browse files Browse the repository at this point in the history
  • Loading branch information
@daminton
daminton committed Nov 16, 2021
1 parent 677060f commit 2e13e18
Show file tree
Hide file tree
Showing 3 changed files with 69 additions and 44 deletions.
75 changes: 51 additions & 24 deletions src/encounter/encounter_check.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1074,23 +1074,43 @@ pure subroutine sweep_dl(i, n1, n2, ntot, ext_ind, ibegi, iendi, ibeg, iend, len
! Internals
integer(I4B) :: j, jbox, dim, jlo, jhi
integer(I4B), dimension(ibegi(1):iendi(1)) :: box
logical, dimension(ibegi(1):iendi(1)) :: lencounterj

lencounteri(:) = .false.
jlo = ibegi(1)
jhi = iendi(1)
box(:) = ext_ind(jlo:jhi)
where(box(:) > ntot)
box(:) = box(:) - ntot

lencounteri(:) = .false.
lencounterj(jlo:jhi) = .false.
box(jlo:jhi) = ext_ind(jlo:jhi)

where(box(jlo:jhi) > ntot)
box(jlo:jhi) = box(jlo:jhi) - ntot
endwhere

! only pairs from the two different lists allowed
where (box(jlo:jhi) > n1)
lencounterj(jlo:jhi) = (i <= n1)
elsewhere
lencounterj(jlo:jhi) = (i > n1)
end where

where (lencounterj(jlo:jhi) .and. (iend(2,box(jlo:jhi)) > ibegi(2)) .and. (ibeg(2,box(jlo:jhi)) < iendi(2)) )
lencounterj(jlo:jhi) = (iend(3,box(jlo:jhi)) > ibegi(3)) .and. (ibeg(3,box(jlo:jhi)) < iendi(3))
end where

do concurrent(jbox = jlo:jhi)
lencounteri(box(jbox)) = lencounterj(jbox)
end do

do concurrent (jbox = jlo:jhi)
j = box(jbox)
if (((i <= n1) .and. (j <= n1)) .or. ((i > n1) .and. (j > n1))) cycle ! only pairs from the two different lists allowed
!do concurrent (jbox = jlo:jhi, ( ( (i <= n1).and.(box(jbox) > n1) ) .or. ( (i > n1).and.(box(jbox) <= n1) ) ) .and. &
! ( ( iend(2,box(j)) >= ibegi(2) ) .and. ( ibeg(2,box(j)) <= iendi(2) ) ) )
!j = box(jbox)
!if (((i <= n1) .and. (j <= n1)) .or. ((i > n1) .and. (j > n1))) cycle ! only pairs from the two different lists allowed
! Check the other dimensions
!lencounteri(j) = all((iend(2:SWEEPDIM,j) > ibegi(2:SWEEPDIM)) .and. (ibeg(2:SWEEPDIM,j) < iendi(2:SWEEPDIM)))
if ((iend(2,j) < ibegi(2)) .or. ibeg(2,j) > iendi(2)) cycle
lencounteri(j) = (iend(SWEEPDIM,j) > ibegi(SWEEPDIM)) .and. (ibeg(SWEEPDIM,j) < iendi(SWEEPDIM))
end do
!if ((iend(2,j) < ibegi(2)) .or. ibeg(2,j) > iendi(2)) cycle
! lencounteri(box(jbox)) = (iend(SWEEPDIM,box(jbox)) > ibegi(SWEEPDIM)) .and. (ibeg(SWEEPDIM,box(jbox)) < iendi(SWEEPDIM))
!end do

return
end subroutine sweep_dl
Expand All @@ -1115,8 +1135,7 @@ subroutine encounter_check_sweep_aabb_all_single_list(n, ext_ind, ibeg, iend, in

!$omp parallel do default(private) schedule(dynamic)&
!$omp shared(ext_ind, ibeg, iend, ind_arr, lenc) &
!$omp firstprivate(n) &
!$omp lastprivate(ibegi, iendi, lencounteri)
!$omp firstprivate(n)
do i = 1, n
ibegi(1) = ibeg(1,i) + 1
iendi(1) = iend(1,i) - 1
Expand Down Expand Up @@ -1148,26 +1167,34 @@ pure subroutine sweep_sl(n, ext_ind, ibegi, iendi, ibeg, iend, lencounteri)
! Internals
integer(I4B) :: j, jbox, dim, jlo, jhi
integer(I4B), dimension(ibegi(1):iendi(1)) :: box

lencounteri(:) = .false.
logical, dimension(ibegi(1):iendi(1)) :: lencounterj

jlo = ibegi(1)
jhi = iendi(1)
lencounteri(:) = .false.
lencounterj(jlo:jhi) = .false.
box(jlo:jhi) = ext_ind(jlo:jhi)

box(:) = ext_ind(jlo:jhi)

where(box(:) > n)
box(:) = box(:) - n
where(box(jlo:jhi) > n)
box(jlo:jhi) = box(jlo:jhi) - n
endwhere

do concurrent(jbox = jlo:jhi) ! Sweep forward until the end of the interval
j = box(jbox)
! Check the other dimensions
!lencounteri(j) = all((iend(2:SWEEPDIM,j) > ibegi(2:SWEEPDIM)) .and. (ibeg(2:SWEEPDIM,j) < iendi(2:SWEEPDIM)))
if ((iend(2,j) < ibegi(2)) .or. ibeg(2,j) > iendi(2)) cycle
lencounteri(j) = (iend(SWEEPDIM,j) > ibegi(SWEEPDIM)) .and. (ibeg(SWEEPDIM,j) < iendi(SWEEPDIM))
where ((iend(2,box(jlo:jhi)) > ibegi(2) ) .and. (ibeg(2,box(jlo:jhi)) < iendi(2)) )
lencounterj(jlo:jhi) = (iend(3,box(jlo:jhi)) > ibegi(3)) .and. (ibeg(3,box(jlo:jhi)) < iendi(3))
end where

do concurrent(jbox = jlo:jhi)
lencounteri(box(jbox)) = lencounterj(jbox)
end do

! do concurrent(jbox = jlo:jhi) ! Sweep forward until the end of the interval
! j = box(jbox)
! ! Check the other dimensions
! !lencounteri(j) = all((iend(2:SWEEPDIM,j) > ibegi(2:SWEEPDIM)) .and. (ibeg(2:SWEEPDIM,j) < iendi(2:SWEEPDIM)))
! if ((iend(2,j) < ibegi(2)) .or. ibeg(2,j) > iendi(2)) cycle
! lencounteri(j) = (iend(SWEEPDIM,j) > ibegi(SWEEPDIM)) .and. (ibeg(SWEEPDIM,j) < iendi(SWEEPDIM))
! end do

return
end subroutine sweep_sl
end subroutine encounter_check_sweep_aabb_all_single_list
Expand Down
18 changes: 3 additions & 15 deletions src/symba/symba_util.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -403,29 +403,17 @@ module subroutine symba_util_flatten_eucl_plpl(self, param)
integer(I8B) :: k, npl, nplm
integer(I4B) :: i, j, err

associate(pl => self, nplpl => self%nplpl, nplplm => self%nplplm)
associate(pl => self, nplplm => self%nplplm)
npl = int(self%nbody, kind=I8B)
select type(param)
class is (symba_parameters)
pl%lmtiny(1:npl) = pl%Gmass(1:npl) < param%GMTINY
end select
nplm = count(.not. pl%lmtiny(1:npl))
pl%nplm = int(nplm, kind=I4B)
nplpl = (npl * (npl - 1_I8B)) / 2_I8B ! number of entries in a strict lower triangle, npl x npl, minus first column
nplplm = nplm * npl - nplm * (nplm + 1_I8B) / 2_I8B ! number of entries in a strict lower triangle, npl x npl, minus first column including only mutually interacting bodies
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, nplpl), stat=err)
if (err /= 0) then ! An error occurred trying to allocate this big array. This probably means it's too big to fit in memory, and so we will force the run back into triangular mode
param%lflatten_interactions = .false.
else
do concurrent (i=1:npl, j=1:npl, j>i)
call util_flatten_eucl_ij_to_k(self%nbody, i, j, k)
self%k_plpl(1, k) = i
self%k_plpl(2, k) = j
end do
end if
end if

call util_flatten_eucl_plpl(pl, param)
end associate

return
Expand Down
20 changes: 15 additions & 5 deletions src/util/util_get_energy_momentum.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -144,7 +144,8 @@ subroutine util_get_energy_potential_flat(npl, nplpl, k_plpl, lmask, GMcb, Gmass
end do

!$omp parallel do default(private) schedule(static)&
!$omp shared(nplpl, k_plpl, xb, mass, Gmass, pepl, lstatpl, lmask)
!$omp shared(k_plpl, xb, mass, Gmass, pepl, lstatpl, lmask) &
!$omp firstprivate(nplpl)
do k = 1, nplpl
i = k_plpl(1,k)
j = k_plpl(2,k)
Expand Down Expand Up @@ -178,7 +179,7 @@ subroutine util_get_energy_potential_triangular(npl, lmask, GMcb, Gmass, mass, x
real(DP), intent(out) :: pe
! Internals
integer(I4B) :: i, j
real(DP), dimension(npl) :: pecb
real(DP), dimension(npl) :: pecb, pepl

! Do the central body potential energy component first
where(.not. lmask(1:npl))
Expand All @@ -189,10 +190,19 @@ subroutine util_get_energy_potential_triangular(npl, lmask, GMcb, Gmass, mass, x
pecb(i) = -GMcb * mass(i) / norm2(xb(:,i))
end do

pe = sum(pecb(1:npl), lmask(1:npl))
do concurrent(i = 1:npl, j = 1:npl, (j > i) .and. lmask(i) .and. lmask(j))
pe = pe - (Gmass(i) * mass(j)) / norm2(xb(:, i) - xb(:, j))
pe = 0.0_DP
!$omp parallel do default(private) schedule(static)&
!$omp shared(lmask, Gmass, mass, xb) &
!$omp firstprivate(npl) &
!$omp reduction(+:pe)
do i = 1, npl
do concurrent(j = i+1:npl, lmask(i) .and. lmask(j))
pepl(j) = - (Gmass(i) * mass(j)) / norm2(xb(:, i) - xb(:, j))
end do
pe = pe + sum(pepl(i+1:npl), lmask(i) .and. lmask(j))
end do
!$omp end parallel do
pe = pe + sum(pecb(1:npl), lmask(1:npl))

return
end subroutine util_get_energy_potential_triangular
Expand Down

0 comments on commit 2e13e18

Please sign in to comment.