diff --git a/src/kick/kick.f90 b/src/kick/kick.f90 index de99a9bb0..ea21c42e8 100644 --- a/src/kick/kick.f90 +++ b/src/kick/kick.f90 @@ -14,7 +14,11 @@ module subroutine kick_getacch_int_pl(self, param) class(swiftest_pl), intent(inout) :: self !! Swiftest massive body object class(swiftest_parameters), intent(in) :: param !! Current swiftest run configuration parameters - call kick_getacch_int_all_flat_pl(self%nbody, self%nplpl, self%k_plpl, self%xh, self%Gmass, self%radius, self%ah) + if (param%lflatten_interactions) then + call kick_getacch_int_all_flat_pl(self%nbody, self%nplpl, self%k_plpl, self%xh, self%Gmass, self%radius, self%ah) + else + call kick_getacch_int_all_triangular_pl(self%nbody, self%nbody, self%xh, self%Gmass, self%radius, self%ah) + end if return end subroutine kick_getacch_int_pl @@ -94,7 +98,7 @@ module subroutine kick_getacch_int_all_flat_pl(npl, nplpl, k_plpl, x, Gmass, rad end subroutine kick_getacch_int_all_flat_pl - module subroutine kick_getacch_int_all_triangular_pl(npl, x, Gmass, radius, acc) + module subroutine kick_getacch_int_all_triangular_pl(npl, nplm, x, Gmass, radius, acc) !! author: David A. Minton !! !! Compute direct cross (third) term heliocentric accelerations for massive bodies, with parallelization. @@ -103,7 +107,8 @@ module subroutine kick_getacch_int_all_triangular_pl(npl, x, Gmass, radius, acc) !! 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 + integer(I4B), intent(in) :: npl !! Total number of massive bodies + integer(I4B), intent(in) :: nplm !! Number of fully interacting 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 @@ -118,11 +123,11 @@ module subroutine kick_getacch_int_all_triangular_pl(npl, x, Gmass, radius, acc) ahj(:,:) = 0.0_DP !$omp parallel do default(private) schedule(static)& - !$omp shared(npl, x, Gmass, radius) & + !$omp shared(npl, nplm, x, Gmass, radius) & !$omp lastprivate(rji2, rlim2, xr, yr, zr) & !$omp reduction(+:ahi) & !$omp reduction(-:ahj) - do i = 1, npl + do i = 1, nplm do concurrent(j = i+1:npl) xr = x(1, j) - x(1, i) yr = x(2, j) - x(2, i) diff --git a/src/modules/swiftest_classes.f90 b/src/modules/swiftest_classes.f90 index a88fb7d7a..f4f9b2698 100644 --- a/src/modules/swiftest_classes.f90 +++ b/src/modules/swiftest_classes.f90 @@ -867,9 +867,10 @@ module subroutine kick_getacch_int_all_flat_pl(npl, nplpl, k_plpl, x, Gmass, rad real(DP), dimension(:,:), intent(inout) :: acc !! Acceleration vector array end subroutine kick_getacch_int_all_flat_pl - module subroutine kick_getacch_int_all_triangular_pl(npl, x, Gmass, radius, acc) + module subroutine kick_getacch_int_all_triangular_pl(npl, nplm, x, Gmass, radius, acc) implicit none - integer(I4B), intent(in) :: npl !! Number of massive bodies + integer(I4B), intent(in) :: npl !! Total number of massive bodies + integer(I4B), intent(in) :: nplm !! Number of fully interacting 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 diff --git a/src/symba/symba_encounter_check.f90 b/src/symba/symba_encounter_check.f90 index 82c4c1890..1df9ca03c 100644 --- a/src/symba/symba_encounter_check.f90 +++ b/src/symba/symba_encounter_check.f90 @@ -8,13 +8,15 @@ subroutine symba_encounter_check_all_flat(nplplm, k_plpl, x, v, rhill, dt, irec !! 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 - 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(:), intent(out) :: lencounter, loc_lvdotr + integer(I8B), intent(in) :: nplplm !! Total number of plm-pl interactions + integer(I4B), dimension(:,:), intent(in) :: k_plpl !! List of all pl-pl interactions + real(DP), dimension(:,:), intent(in) :: x !! Position vectors of massive bodies + real(DP), dimension(:,:), intent(in) :: v !! Velocity vectors of massive bodies + real(DP), dimension(:), intent(in) :: rhill !! Hill's radii of massive bodies + real(DP), intent(in) :: dt !! Step size + integer(I4B), intent(in) :: irec !! Current recursion depth + logical, dimension(:), intent(out) :: lencounter !! Logical array indicating which pair is in an encounter state + logical, dimension(:), intent(out) :: loc_lvdotr !! Logical array indicating the sign of v .dot. x for each encounter ! Internals integer(I8B) :: k integer(I4B) :: i, j @@ -42,36 +44,40 @@ subroutine symba_encounter_check_all_flat(nplplm, k_plpl, x, v, rhill, dt, irec 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) + subroutine symba_encounter_check_all_triangular(npl, nplm, x, v, rhill, dt, irec, lvdotr, index1, index2, 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 + integer(I4B), intent(in) :: npl !! Total number of massive bodies + integer(I4B), intent(in) :: nplm !! Number of fully interacting massive bodies + real(DP), dimension(:,:), intent(in) :: x !! Position vectors of massive bodies + real(DP), dimension(:,:), intent(in) :: v !! Velocity vectors of massive bodies + real(DP), dimension(:), intent(in) :: rhill !! Hill's radii of massive bodies + real(DP), intent(in) :: dt !! Step size + integer(I4B), intent(in) :: irec !! Current recursion depth + logical, dimension(:), allocatable, intent(out) :: lvdotr !! Logical flag indicating the sign of v .dot. x + integer(I4B), dimension(:), allocatable, intent(out) :: index1 !! List of indices for body 1 in each interaction + integer(I4B), dimension(:), allocatable, intent(out) :: index2 !! List of indices for body 2 in each interaction + integer(I4B), intent(out) :: nenc !! Total number of interactions ! Internals integer(I4B) :: i, j, nenci, j0, j1 real(DP) :: xr, yr, zr, vxr, vyr, vzr, rhill1, rhill2 - logical, dimension(npl) :: lencounteri, loc_lvdotri + logical, dimension(npl) :: lencounteri, lvdotri integer(I4B), dimension(npl) :: ind_arr type lenctype - logical, dimension(:), allocatable :: loc_lvdotr + logical, dimension(:), allocatable :: 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 parallel do 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) + !$omp firstprivate(ind_arr) & + !$omp lastprivate(xr, yr, zr, vxr, vyr, vzr, rhill1, rhill2, lencounteri, lvdotri) do i = 1, nplm do concurrent(j = i+1:npl) xr = x(1, j) - x(1, i) @@ -82,30 +88,31 @@ subroutine symba_encounter_check_all_triangular(npl, nplm, x, v, rhill, dt, ire 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)) + call symba_encounter_check_one(xr, yr, zr, vxr, vyr, vzr, rhill1, rhill2, dt, irec, lencounteri(j), 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)) + allocate(lenc(i)%lvdotr(nenci), lenc(i)%index2(nenci)) + lenc(i)%lvdotr(:) = pack(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 + !$omp end parallel do 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)) + allocate(lvdotr(nenc)) + allocate(index1(nenc)) + allocate(index2(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(:) + lvdotr(j0:j1) = lenc(i)%lvdotr(:) + index1(j0:j1) = i + index2(j0:j1) = lenc(i)%index2(:) j0 = j1 + 1 end if end do @@ -131,25 +138,26 @@ module function symba_encounter_check_pl(self, param, system, dt, irec) result(l logical :: lany_encounter !! Returns true if there is at least one close encounter ! Internals integer(I8B) :: k, nplplm, kenc - integer(I4B) :: i, j, nenc, npl + integer(I4B) :: i, j, nenc, npl, nplm logical, dimension(:), allocatable :: lencounter, loc_lvdotr, lvdotr integer(I4B), dimension(:), allocatable :: index1, index2 + integer(I4B), dimension(:,:), allocatable :: k_plpl_enc if (self%nbody == 0) return - associate(pl => self) - nplplm = pl%nplplm + associate(pl => self, plplenc_list => system%plplenc_list) npl = pl%nbody - allocate(lencounter(nplplm)) - allocate(loc_lvdotr(nplplm)) + if (param%lflatten_interactions) then + nplplm = pl%nplplm + allocate(lencounter(nplplm)) + allocate(loc_lvdotr(nplplm)) - call symba_encounter_check_all_flat(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(:)) + nenc = count(lencounter(:)) - lany_encounter = nenc > 0 - if (lany_encounter) then - associate(plplenc_list => system%plplenc_list) + lany_encounter = nenc > 0 + if (lany_encounter) then call plplenc_list%resize(nenc) allocate(lvdotr(nenc)) allocate(index1(nenc)) @@ -160,28 +168,38 @@ module function symba_encounter_check_pl(self, param, system, dt, irec) result(l deallocate(lencounter, loc_lvdotr) call move_alloc(lvdotr, plplenc_list%lvdotr) call move_alloc(index1, plplenc_list%index1) - call move_alloc(index2, plplenc_list%index2) - do k = 1, nenc - i = plplenc_list%index1(k) - j = plplenc_list%index2(k) - call util_index_eucl_ij_to_k(npl, i, j, kenc) - plplenc_list%kidx(k) = kenc - plplenc_list%id1(k) = pl%id(plplenc_list%index1(k)) - plplenc_list%id2(k) = pl%id(plplenc_list%index2(k)) - plplenc_list%status(k) = ACTIVE - plplenc_list%level(k) = irec - pl%lencounter(i) = .true. - pl%lencounter(j) = .true. - pl%levelg(i) = irec - pl%levelm(i) = irec - pl%levelg(j) = irec - pl%levelm(j) = irec - pl%nplenc(i) = pl%nplenc(i) + 1 - pl%nplenc(j) = pl%nplenc(j) + 1 - end do - end associate + call move_alloc(index2, plplenc_list%index2) + end if + else + nplm = pl%nplm + call symba_encounter_check_all_triangular(npl, nplm, pl%xh, pl%vh, pl%rhill, dt, irec, lvdotr, index1, index2, nenc) + call move_alloc(lvdotr, plplenc_list%lvdotr) + call move_alloc(index1, plplenc_list%index1) + call move_alloc(index2, plplenc_list%index2) + end if + + if (lany_encounter) then + do k = 1, nenc + i = plplenc_list%index1(k) + j = plplenc_list%index2(k) + call util_index_eucl_ij_to_k(npl, i, j, kenc) + plplenc_list%kidx(k) = kenc + plplenc_list%id1(k) = pl%id(plplenc_list%index1(k)) + plplenc_list%id2(k) = pl%id(plplenc_list%index2(k)) + plplenc_list%status(k) = ACTIVE + plplenc_list%level(k) = irec + pl%lencounter(i) = .true. + pl%lencounter(j) = .true. + pl%levelg(i) = irec + pl%levelm(i) = irec + pl%levelg(j) = irec + pl%levelm(j) = irec + pl%nplenc(i) = pl%nplenc(i) + 1 + pl%nplenc(j) = pl%nplenc(j) + 1 + end do end if end associate + return end function symba_encounter_check_pl diff --git a/src/symba/symba_kick.f90 b/src/symba/symba_kick.f90 index 2011ee0c2..414c89805 100644 --- a/src/symba/symba_kick.f90 +++ b/src/symba/symba_kick.f90 @@ -14,7 +14,11 @@ module subroutine symba_kick_getacch_int_pl(self, param) class(symba_pl), intent(inout) :: self !! SyMBA massive body object class(swiftest_parameters), intent(in) :: param !! Current swiftest run configuration parameter - call kick_getacch_int_all_flat_pl(self%nbody, self%nplplm, self%k_plpl, self%xh, self%Gmass, self%radius, self%ah) + if (param%lflatten_interactions) then + call kick_getacch_int_all_flat_pl(self%nbody, self%nplplm, self%k_plpl, self%xh, self%Gmass, self%radius, self%ah) + else + call kick_getacch_int_all_triangular_pl(self%nbody, self%nplm, self%xh, self%Gmass, self%radius, self%ah) + end if return end subroutine symba_kick_getacch_int_pl @@ -44,15 +48,15 @@ module subroutine symba_kick_getacch_pl(self, system, param, t, lbeg) if (self%nbody == 0) return select type(system) class is (symba_nbody_system) - associate(pl => self, npl => self%nbody, plplenc_list => system%plplenc_list, radius => self%radius) + associate(pl => self, npl => self%nbody, nplm => self%nplm, plplenc_list => system%plplenc_list, radius => self%radius) ! Apply kicks to all bodies (including those in the encounter list) call helio_kick_getacch_pl(pl, system, param, t, lbeg) if (plplenc_list%nenc > 0) then ! Remove kicks from bodies involved currently in the encounter list, as these are dealt with separately. + ah_enc(:,:) = 0.0_DP nplplenc = int(plplenc_list%nenc, kind=I8B) 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_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 @@ -63,7 +67,6 @@ module subroutine symba_kick_getacch_pl(self, system, param, t, lbeg) return end subroutine symba_kick_getacch_pl - module subroutine symba_kick_getacch_tp(self, system, param, t, lbeg) !! author: David A. Minton !! diff --git a/src/symba/symba_util.f90 b/src/symba/symba_util.f90 index 752028cf7..0d062894f 100644 --- a/src/symba/symba_util.f90 +++ b/src/symba/symba_util.f90 @@ -292,18 +292,19 @@ module subroutine symba_util_index_eucl_plpl(self, param) npl = int(self%nbody, kind=I8B) nplm = count(.not. pl%lmtiny(1:npl)) pl%nplm = int(nplm, kind=I4B) - pl%nplpl = (npl * (npl - 1) / 2) ! number of entries in a strict lower triangle, npl x npl, minus first column pl%nplplm = nplm * npl - nplm * (nplm + 1) / 2 ! number of entries in a strict lower triangle, npl x npl, minus first column including only mutually interacting bodies - if (allocated(self%k_plpl)) deallocate(self%k_plpl) ! Reset the index array if it's been set previously - allocate(self%k_plpl(2, pl%nplpl)) - do concurrent (i = 1:npl) - do concurrent (j = i+1:npl) - call util_index_eucl_ij_to_k(npl, i, j, k) - self%k_plpl(1, k) = i - self%k_plpl(2, k) = j + 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, pl%nplpl)) + do concurrent (i = 1:npl) + do concurrent (j = i+1:npl) + call util_index_eucl_ij_to_k(npl, i, j, k) + self%k_plpl(1, k) = i + self%k_plpl(2, k) = j + end do end do - end do + end if end associate return diff --git a/src/util/util_get_energy_momentum.f90 b/src/util/util_get_energy_momentum.f90 index 677de5646..3474e2921 100644 --- a/src/util/util_get_energy_momentum.f90 +++ b/src/util/util_get_energy_momentum.f90 @@ -78,8 +78,12 @@ module subroutine util_get_energy_momentum_system(self, param) kespincb = 0.0_DP kespinpl(:) = 0.0_DP end if - - call util_get_energy_potential(npl, nplpl, pl%k_plpl, pl%lmask, cb%Gmass, pl%Gmass, pl%mass, pl%xb, system%pe) + + if (param%lflatten_interactions) then + call util_get_energy_potential_flat(npl, nplpl, pl%k_plpl, pl%lmask, cb%Gmass, pl%Gmass, pl%mass, pl%xb, system%pe) + else + call util_get_energy_potential_triangular(npl, pl%lmask, cb%Gmass, pl%Gmass, pl%mass, pl%xb, system%pe) + end if ! Potential energy from the oblateness term if (param%loblatecb) then @@ -108,7 +112,7 @@ module subroutine util_get_energy_momentum_system(self, param) end subroutine util_get_energy_momentum_system - subroutine util_get_energy_potential(npl, nplpl, k_plpl, lmask, GMcb, Gmass, mass, xb, pe) + subroutine util_get_energy_potential_flat(npl, nplpl, k_plpl, lmask, GMcb, Gmass, mass, xb, pe) !! author: David A. Minton !! !! Compute total system potential energy @@ -156,6 +160,41 @@ subroutine util_get_energy_potential(npl, nplpl, k_plpl, lmask, GMcb, Gmass, mas pe = sum(pepl(:), lstatpl(:)) + sum(pecb(1:npl), lmask(1:npl)) return - end subroutine util_get_energy_potential + end subroutine util_get_energy_potential_flat + + + subroutine util_get_energy_potential_triangular(npl, lmask, GMcb, Gmass, mass, xb, pe) + !! author: David A. Minton + !! + !! Compute total system potential energy + implicit none + ! Arguments + integer(I4B), intent(in) :: npl + logical, dimension(:), intent(in) :: lmask + real(DP), intent(in) :: GMcb + real(DP), dimension(:), intent(in) :: Gmass + real(DP), dimension(:), intent(in) :: mass + real(DP), dimension(:,:), intent(in) :: xb + real(DP), intent(out) :: pe + ! Internals + integer(I4B) :: i, j + real(DP), dimension(npl) :: pecb + + ! Do the central body potential energy component first + where(.not. lmask(1:npl)) + pecb(1:npl) = 0.0_DP + end where + + do concurrent(i = 1:npl, lmask(i)) + 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)) + end do + + return + end subroutine util_get_energy_potential_triangular end submodule s_util_get_energy_momentum diff --git a/src/util/util_index.f90 b/src/util/util_index.f90 index 1ee60e400..d0dd83896 100644 --- a/src/util/util_index.f90 +++ b/src/util/util_index.f90 @@ -82,13 +82,15 @@ module subroutine util_index_eucl_plpl(self, param) npl = int(self%nbody, kind=I8B) associate(nplpl => self%nplpl) nplpl = (npl * (npl - 1) / 2) ! number of entries in a strict lower triangle, npl x npl - 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)) - do concurrent (i=1:npl, j=1:npl, j>i) - call util_index_eucl_ij_to_k(npl, i, j, k) - self%k_plpl(1, k) = i - self%k_plpl(2, k) = j - end do + 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)) + do concurrent (i=1:npl, j=1:npl, j>i) + call util_index_eucl_ij_to_k(npl, i, j, k) + self%k_plpl(1, k) = i + self%k_plpl(2, k) = j + end do + end if end associate return