diff --git a/src/encounter/encounter_check.f90 b/src/encounter/encounter_check.f90 index 1f7550125..a29aab6eb 100644 --- a/src/encounter/encounter_check.f90 +++ b/src/encounter/encounter_check.f90 @@ -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 @@ -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 @@ -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 diff --git a/src/symba/symba_util.f90 b/src/symba/symba_util.f90 index 5a985faea..fcb576d0f 100644 --- a/src/symba/symba_util.f90 +++ b/src/symba/symba_util.f90 @@ -403,7 +403,7 @@ 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) @@ -411,21 +411,9 @@ module subroutine symba_util_flatten_eucl_plpl(self, param) 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 diff --git a/src/util/util_get_energy_momentum.f90 b/src/util/util_get_energy_momentum.f90 index 3474e2921..b935a680c 100644 --- a/src/util/util_get_energy_momentum.f90 +++ b/src/util/util_get_energy_momentum.f90 @@ -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) @@ -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)) @@ -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