From 5e6eeb0b2bb256ccceefc1213c325ec0b2515f18 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Wed, 8 Mar 2023 08:54:45 -0500 Subject: [PATCH] Simplified sort and sweep to use radial distance and set up for testing --- src/encounter/encounter_check.f90 | 270 ++++++++++++---------------- src/encounter/encounter_module.f90 | 16 +- src/encounter/encounter_util.f90 | 32 ++-- src/rmvs/rmvs_encounter_check.f90 | 2 +- src/swiftest/swiftest_module.f90 | 4 +- src/swiftest/swiftest_orbel.f90 | 4 +- src/symba/symba_encounter_check.f90 | 2 +- src/symba/symba_step.f90 | 2 +- 8 files changed, 141 insertions(+), 191 deletions(-) diff --git a/src/encounter/encounter_check.f90 b/src/encounter/encounter_check.f90 index 917fcae45..aa8da45a5 100644 --- a/src/encounter/encounter_check.f90 +++ b/src/encounter/encounter_check.f90 @@ -11,18 +11,18 @@ use swiftest contains - module subroutine encounter_check_all_plpl(param, npl, x, v, renc, dt, nenc, index1, index2, lvdotr) + module subroutine encounter_check_all_plpl(param, npl, r, v, renc, dt, nenc, index1, index2, lvdotr) !! author: David A. Minton !! !! Check for encounters between massive bodies. Choose between the standard triangular or the Sort & Sweep method based on user inputs !! implicit none ! Arguments - class(base_parameters), intent(inout) :: param !! Current Swiftest run configuration parameter5s + class(base_parameters), intent(inout) :: param !! Current Swiftest run configuration parameter5s integer(I4B), intent(in) :: npl !! Total number of 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) :: renc !! Critical radii of massive bodies that defines an encounter + real(DP), dimension(:,:), intent(in) :: r !! Position vectors of massive bodies + real(DP), dimension(:,:), intent(in) :: v !! Velocity vectors of massive bodies real(DP), intent(in) :: dt !! Step size integer(I8B), intent(out) :: nenc !! Total number of encounters integer(I4B), dimension(:), allocatable, intent(out) :: index1 !! List of indices for body 1 in each encounter @@ -35,7 +35,8 @@ module subroutine encounter_check_all_plpl(param, npl, x, v, renc, dt, nenc, ind integer(I8B) :: nplpl = 0_I8B - call encounter_check_all_triangular_plpl(npl, x, v, renc, dt, nenc, index1, index2, lvdotr) + !call encounter_check_all_triangular_plpl(npl, r, v, renc, dt, nenc, index1, index2, lvdotr) + call encounter_check_all_sort_and_sweep_plpl(npl, r, v, renc, dt, nenc, index1, index2, lvdotr) return end subroutine encounter_check_all_plpl @@ -48,7 +49,7 @@ module subroutine encounter_check_all_plplm(param, nplm, nplt, rplm, vplm, rplt, !! implicit none ! Arguments - class(base_parameters), intent(inout) :: param !! Current Swiftest run configuration parameter5s + class(base_parameters), intent(inout) :: param !! Current Swiftest run configuration parameter5s integer(I4B), intent(in) :: nplm !! Total number of fully interacting massive bodies integer(I4B), intent(in) :: nplt !! Total number of partially interacting masive bodies (GM < GMTINY) real(DP), dimension(:,:), intent(in) :: rplm !! Position vectors of fully interacting massive bodies @@ -113,19 +114,19 @@ module subroutine encounter_check_all_plplm(param, nplm, nplt, rplm, vplm, rplt, end subroutine encounter_check_all_plplm - module subroutine encounter_check_all_pltp(param, npl, ntp, rpl, vpl, xtp, vtp, renc, dt, nenc, index1, index2, lvdotr) + module subroutine encounter_check_all_pltp(param, npl, ntp, rpl, vpl, rtp, vtp, renc, dt, nenc, index1, index2, lvdotr) !! author: David A. Minton !! !! Check for encounters between massive bodies and test particles. Choose between the standard triangular or the Sort & Sweep method based on user inputs !! implicit none ! Arguments - class(base_parameters), intent(inout) :: param !! Current Swiftest run configuration parameter5s + class(base_parameters), intent(inout) :: param !! Current Swiftest run configuration parameter5s integer(I4B), intent(in) :: npl !! Total number of massive bodies integer(I4B), intent(in) :: ntp !! Total number of test particles real(DP), dimension(:,:), intent(in) :: rpl !! Position vectors of massive bodies real(DP), dimension(:,:), intent(in) :: vpl !! Velocity vectors of massive bodies - real(DP), dimension(:,:), intent(in) :: xtp !! Position vectors of test particlse + real(DP), dimension(:,:), intent(in) :: rtp !! Position vectors of test particlse real(DP), dimension(:,:), intent(in) :: vtp !! Velocity vectors of test particles real(DP), dimension(:), intent(in) :: renc !! Critical radii of massive bodies that defines an encounter real(DP), intent(in) :: dt !! Step size @@ -139,13 +140,13 @@ module subroutine encounter_check_all_pltp(param, npl, ntp, rpl, vpl, xtp, vtp, logical, save :: lsecond = .false. integer(I8B) :: npltp = 0_I8B - call encounter_check_all_triangular_pltp(npl, ntp, rpl, vpl, xtp, vtp, renc, dt, nenc, index1, index2, lvdotr) + call encounter_check_all_triangular_pltp(npl, ntp, rpl, vpl, rtp, vtp, renc, dt, nenc, index1, index2, lvdotr) return end subroutine encounter_check_all_pltp - subroutine encounter_check_all_sort_and_sweep_plpl(npl, x, v, renc, dt, nenc, index1, index2, lvdotr) + subroutine encounter_check_all_sort_and_sweep_plpl(npl, r, v, renc, dt, nenc, index1, index2, lvdotr) !! author: David A. Minton !! !! Check for encounters between massive bodies. @@ -155,7 +156,7 @@ subroutine encounter_check_all_sort_and_sweep_plpl(npl, x, v, renc, dt, nenc, in implicit none ! Arguments integer(I4B), intent(in) :: npl !! Total number of massive bodies - real(DP), dimension(:,:), intent(in) :: x !! Position vectors of massive bodies + real(DP), dimension(:,:), intent(in) :: r !! Position vectors of massive bodies real(DP), dimension(:,:), intent(in) :: v !! Velocity vectors of massive bodies real(DP), dimension(:), intent(in) :: renc !! Critical radii of massive bodies that defines an encounter real(DP), intent(in) :: dt !! Step size @@ -164,10 +165,11 @@ subroutine encounter_check_all_sort_and_sweep_plpl(npl, x, v, renc, dt, nenc, in integer(I4B), dimension(:), allocatable, intent(out) :: index2 !! List of indices for body 2 in each encounter logical, dimension(:), allocatable, intent(out) :: lvdotr !! Logical flag indicating the sign of v .dot. x ! Internals - integer(I4B) :: dim, n + integer(I4B) :: i, n integer(I4B), save :: npl_last = 0 type(encounter_bounding_box), save :: boundingbox - integer(I2B), dimension(npl) :: vshift_min, vshift_max + real(DP), dimension(npl) :: rmin,rmax + real(DP) :: rmag if (npl == 0) return @@ -178,23 +180,15 @@ subroutine encounter_check_all_sort_and_sweep_plpl(npl, x, v, renc, dt, nenc, in npl_last = npl end if - !$omp parallel do default(private) schedule(static) & - !$omp shared(x, v, renc, boundingbox) & - !$omp firstprivate(dt, npl, n) - do dim = 1, SWEEPDIM - where(v(dim,1:npl) < 0.0_DP) - vshift_min(1:npl) = 1 - vshift_max(1:npl) = 0 - elsewhere - vshift_min(1:npl) = 0 - vshift_max(1:npl) = 1 - end where - call boundingbox%aabb(dim)%sort(npl, [x(dim,1:npl) - renc(1:npl) + vshift_min(1:npl) * v(dim,1:npl) * dt, & - x(dim,1:npl) + renc(1:npl) + vshift_max(1:npl) * v(dim,1:npl) * dt]) + do concurrent (i = 1:npl) + rmag = .mag.r(:,i) + rmax(i) = rmag + RSWEEP_FACTOR * renc(i) + rmin(i) = rmag - RSWEEP_FACTOR * renc(i) end do - !$omp end parallel do - call boundingbox%sweep(npl, x, v, renc, dt, nenc, index1, index2, lvdotr) + call boundingbox%aabb%sort(npl, [rmin,rmax]) + + call boundingbox%sweep(npl, r, v, renc, dt, nenc, index1, index2, lvdotr) return end subroutine encounter_check_all_sort_and_sweep_plpl @@ -224,10 +218,10 @@ subroutine encounter_check_all_sort_and_sweep_plplm(nplm, nplt, rplm, vplm, rplt logical, dimension(:), allocatable, intent(out) :: lvdotr !! Logical flag indicating the sign of v .dot. x ! Internals type(encounter_bounding_box), save :: boundingbox - integer(I4B) :: dim, n, ntot + integer(I4B) :: i, n, ntot integer(I4B), save :: ntot_last = 0 - integer(I2B), dimension(nplm) :: vplmshift_min, vplmshift_max - integer(I2B), dimension(nplt) :: vpltshift_min, vpltshift_max + real(DP), dimension(nplm+nplt) :: rmin,rmax + real(DP) :: rmag ! If this is the first time through, build the index lists if ((nplm == 0) .or. (nplt == 0)) return @@ -239,32 +233,18 @@ subroutine encounter_check_all_sort_and_sweep_plplm(nplm, nplt, rplm, vplm, rplt ntot_last = ntot end if - !$omp parallel do default(private) schedule(static) & - !$omp shared(rplm, rplt, vplm, vplt, rencm, renct, boundingbox) & - !$omp firstprivate(dt, nplm, nplt, ntot) - do dim = 1, SWEEPDIM - where(vplm(dim,1:nplm) < 0.0_DP) - vplmshift_min(1:nplm) = 1 - vplmshift_max(1:nplm) = 0 - elsewhere - vplmshift_min(1:nplm) = 0 - vplmshift_max(1:nplm) = 1 - end where - - where(vplt(dim,1:nplt) < 0.0_DP) - vpltshift_min(1:nplt) = 1 - vpltshift_max(1:nplt) = 0 - elsewhere - vpltshift_min(1:nplt) = 0 - vpltshift_max(1:nplt) = 1 - end where - - call boundingbox%aabb(dim)%sort(ntot, [rplm(dim,1:nplm) - rencm(1:nplm) + vplmshift_min(1:nplm) * vplm(dim,1:nplm) * dt, & - rplt(dim,1:nplt) - renct(1:nplt) + vpltshift_min(1:nplt) * vplt(dim,1:nplt) * dt, & - rplm(dim,1:nplm) + rencm(1:nplm) + vplmshift_max(1:nplm) * vplm(dim,1:nplm) * dt, & - rplt(dim,1:nplt) + renct(1:nplt) + vpltshift_max(1:nplt) * vplt(dim,1:nplt) * dt]) + do concurrent (i = 1:nplm) + rmag = .mag.rplm(:,i) + rmax(i) = rmag + RSWEEP_FACTOR * rencm(i) + rmin(i) = rmag - RSWEEP_FACTOR * rencm(i) end do - !$omp end parallel do + do concurrent (i = 1:nplt) + rmag = .mag.rplt(:,i) + rmax(nplm+i) = rmag + RSWEEP_FACTOR * renct(i) + rmin(nplm+i) = rmag - RSWEEP_FACTOR * renct(i) + end do + + call boundingbox%aabb%sort(ntot, [rmin, rmax]) call boundingbox%sweep(nplm, nplt, rplm, vplm, rplt, vplt, rencm, renct, dt, nenc, index1, index2, lvdotr) @@ -272,7 +252,7 @@ subroutine encounter_check_all_sort_and_sweep_plplm(nplm, nplt, rplm, vplm, rplt end subroutine encounter_check_all_sort_and_sweep_plplm - subroutine encounter_check_all_sort_and_sweep_pltp(npl, ntp, rpl, vpl, xtp, vtp, rencpl, dt, nenc, index1, index2, lvdotr) + subroutine encounter_check_all_sort_and_sweep_pltp(npl, ntp, rpl, vpl, rtp, vtp, rencpl, dt, nenc, index1, index2, lvdotr) !! author: David A. Minton !! !! Check for encounters between massive bodies and test particles. @@ -285,7 +265,7 @@ subroutine encounter_check_all_sort_and_sweep_pltp(npl, ntp, rpl, vpl, xtp, vtp, integer(I4B), intent(in) :: ntp !! Total number of test particles real(DP), dimension(:,:), intent(in) :: rpl !! Position vectors of massive bodies real(DP), dimension(:,:), intent(in) :: vpl !! Velocity vectors of massive bodies - real(DP), dimension(:,:), intent(in) :: xtp !! Position vectors of massive bodies + real(DP), dimension(:,:), intent(in) :: rtp !! Position vectors of massive bodies real(DP), dimension(:,:), intent(in) :: vtp !! Velocity vectors of massive bodies real(DP), dimension(:), intent(in) :: rencpl !! Critical radii of massive bodies that defines an encounter real(DP), intent(in) :: dt !! Step size @@ -295,11 +275,11 @@ subroutine encounter_check_all_sort_and_sweep_pltp(npl, ntp, rpl, vpl, xtp, vtp, logical, dimension(:), allocatable, intent(out) :: lvdotr !! Logical flag indicating the sign of v .dot. x ! Internals type(encounter_bounding_box), save :: boundingbox - integer(I4B) :: dim, n, ntot + integer(I4B) :: i, n, ntot integer(I4B), save :: ntot_last = 0 - integer(I2B), dimension(npl) :: vplshift_min, vplshift_max - integer(I2B), dimension(ntp) :: vtpshift_min, vtpshift_max - real(DP), dimension(ntp) :: renctp + real(DP), dimension(npl+ntp) :: rmin,rmax + real(DP), dimension(ntp) :: renctp + real(DP) :: rmag ! If this is the first time through, build the index lists if ((ntp == 0) .or. (npl == 0)) return @@ -312,35 +292,21 @@ subroutine encounter_check_all_sort_and_sweep_pltp(npl, ntp, rpl, vpl, xtp, vtp, end if renctp(:) = 0.0_DP - - !$omp parallel do default(private) schedule(static) & - !$omp shared(rpl, xtp, vpl, vtp, rencpl, renctp, boundingbox) & - !$omp firstprivate(dt, npl, ntp, ntot) - do dim = 1, SWEEPDIM - where(vpl(dim,1:npl) < 0.0_DP) - vplshift_min(1:npl) = 1 - vplshift_max(1:npl) = 0 - elsewhere - vplshift_min(1:npl) = 0 - vplshift_max(1:npl) = 1 - end where - - where(vtp(dim,1:ntp) < 0.0_DP) - vtpshift_min(1:ntp) = 1 - vtpshift_max(1:ntp) = 0 - elsewhere - vtpshift_min(1:ntp) = 0 - vtpshift_max(1:ntp) = 1 - end where - - call boundingbox%aabb(dim)%sort(ntot, [rpl(dim,1:npl) - rencpl(1:npl) + vplshift_min(1:npl) * vpl(dim,1:npl) * dt, & - xtp(dim,1:ntp) - renctp(1:ntp) + vtpshift_min(1:ntp) * vtp(dim,1:ntp) * dt, & - rpl(dim,1:npl) + rencpl(1:npl) + vplshift_max(1:npl) * vpl(dim,1:npl) * dt, & - xtp(dim,1:ntp) + renctp(1:ntp) + vtpshift_max(1:ntp) * vtp(dim,1:ntp) * dt]) + + do concurrent (i = 1:npl) + rmag = .mag.rpl(:,i) + rmax(i) = rmag + RSWEEP_FACTOR * rencpl(i) + rmin(i) = rmag - RSWEEP_FACTOR * rencpl(i) + end do + do concurrent (i = 1:ntp) + rmag = .mag.rtp(:,i) + rmax(npl+i) = rmag + RSWEEP_FACTOR * renctp(i) + rmin(npl+i) = rmag - RSWEEP_FACTOR * renctp(i) end do - !$omp end parallel do - call boundingbox%sweep(npl, ntp, rpl, vpl, xtp, vtp, rencpl, renctp, dt, nenc, index1, index2, lvdotr) + call boundingbox%aabb%sort(ntot, [rmin, rmax]) + + call boundingbox%sweep(npl, ntp, rpl, vpl, rtp, vtp, rencpl, renctp, dt, nenc, index1, index2, lvdotr) return end subroutine encounter_check_all_sort_and_sweep_pltp @@ -445,7 +411,7 @@ subroutine encounter_check_all_triangular_one(i, n, xi, yi, zi, vxi, vyi, vzi, x end subroutine encounter_check_all_triangular_one - subroutine encounter_check_all_triangular_plpl(npl, x, v, renc, dt, nenc, index1, index2, lvdotr) + subroutine encounter_check_all_triangular_plpl(npl, r, v, renc, dt, nenc, index1, index2, lvdotr) !! author: David A. Minton !! !! Check for encounters between massive bodies. Split off from the main subroutine for performance @@ -453,9 +419,9 @@ subroutine encounter_check_all_triangular_plpl(npl, x, v, renc, dt, nenc, index1 implicit none ! Arguments integer(I4B), intent(in) :: npl !! Total number of massive bodies - real(DP), dimension(:,:), intent(in) :: x !! Position vectors of massive bodies + real(DP), dimension(:,:), intent(in) :: r !! Position vectors of massive bodies real(DP), dimension(:,:), intent(in) :: v !! Velocity vectors of massive bodies - real(DP), dimension(:), intent(in) :: renc !! Critical radii of massive bodies that defines an encounter + real(DP), dimension(:), intent(in) :: renc !! Critical radii of massive bodies that defines an encounter real(DP), intent(in) :: dt !! Step size integer(I8B), intent(out) :: nenc !! Total number of encounters integer(I4B), dimension(:), allocatable, intent(out) :: index1 !! List of indices for body 1 in each encounter @@ -469,12 +435,12 @@ subroutine encounter_check_all_triangular_plpl(npl, x, v, renc, dt, nenc, index1 call swiftest_util_index_array(ind_arr, npl) !$omp parallel do default(private) schedule(static)& - !$omp shared(x, v, renc, lenc, ind_arr) & + !$omp shared(r, v, renc, lenc, ind_arr) & !$omp firstprivate(npl, dt) do i = 1,npl - call encounter_check_all_triangular_one(i, npl, x(1,i), x(2,i), x(3,i), & + call encounter_check_all_triangular_one(i, npl, r(1,i), r(2,i), r(3,i), & v(1,i), v(2,i), v(3,i), & - x(1,:), x(2,:), x(3,:), & + r(1,:), r(2,:), r(3,:), & v(1,:), v(2,:), v(3,:), & renc(i), renc(:), dt, ind_arr(:), lenc(i)) if (lenc(i)%nenc > 0) lenc(i)%index1(:) = i @@ -483,7 +449,6 @@ subroutine encounter_check_all_triangular_plpl(npl, x, v, renc, dt, nenc, index1 call encounter_check_collapse_ragged_list(lenc, npl, nenc, index1, index2, lvdotr) - nenc = 0 return end subroutine encounter_check_all_triangular_plpl @@ -820,9 +785,8 @@ module subroutine encounter_check_sweep_aabb_double_list(self, n1, n2, r1, v1, r ! Internals integer(I4B) :: ii, i, ntot, nbox, dim logical, dimension(n1+n2) :: loverlap - logical, dimension(SWEEPDIM,n1+n2) :: loverlap_by_dimension - logical, dimension(SWEEPDIM,2*(n1+n2)) :: llist1 - integer(I4B), dimension(SWEEPDIM,2*(n1+n2)) :: ext_ind + logical, dimension(2*(n1+n2)) :: llist1 + integer(I4B), dimension(2*(n1+n2)) :: ext_ind type(collision_list_pltp), dimension(n1+n2) :: lenc !! Array of encounter lists (one encounter list per body) integer(I4B), dimension(:), allocatable, save :: ind_arr integer(I8B) :: ibeg, iend @@ -831,39 +795,33 @@ module subroutine encounter_check_sweep_aabb_double_list(self, n1, n2, r1, v1, r ntot = n1 + n2 call swiftest_util_index_array(ind_arr, ntot) - do concurrent(dim = 1:SWEEPDIM) - loverlap_by_dimension(dim,:) = (self%aabb(dim)%ibeg(:) + 1_I8B) < (self%aabb(dim)%iend(:) - 1_I8B) - where(self%aabb(dim)%ind(:) > ntot) - ext_ind(dim,:) = self%aabb(dim)%ind(:) - ntot - elsewhere - ext_ind(dim,:) = self%aabb(dim)%ind(:) - endwhere - end do - llist1(:,:) = ext_ind(:,:) <= n1 - where(.not.llist1(:,:)) ext_ind(:,:) = ext_ind(:,:) - n1 + loverlap(:) = (self%aabb%ibeg(:) + 1_I8B) < (self%aabb%iend(:) - 1_I8B) + where(self%aabb%ind(:) > ntot) + ext_ind(:) = self%aabb%ind(:) - ntot + elsewhere + ext_ind(:) = self%aabb%ind(:) + endwhere + llist1(:) = ext_ind(:) <= n1 + where(.not.llist1(:)) ext_ind(:) = ext_ind(:) - n1 - loverlap(:) = loverlap_by_dimension(1,:) - do dim = 2, SWEEPDIM - loverlap(:) = loverlap(:) .and. loverlap_by_dimension(dim,:) - end do dim = 1 - where(llist1(dim,:)) - xind(:) = r1(1,ext_ind(dim,:)) - yind(:) = r1(2,ext_ind(dim,:)) - zind(:) = r1(3,ext_ind(dim,:)) - vxind(:) = v1(1,ext_ind(dim,:)) - vyind(:) = v1(2,ext_ind(dim,:)) - vzind(:) = v1(3,ext_ind(dim,:)) - rencind(:) = renc1(ext_ind(dim,:)) + where(llist1(:)) + xind(:) = r1(1,ext_ind(:)) + yind(:) = r1(2,ext_ind(:)) + zind(:) = r1(3,ext_ind(:)) + vxind(:) = v1(1,ext_ind(:)) + vyind(:) = v1(2,ext_ind(:)) + vzind(:) = v1(3,ext_ind(:)) + rencind(:) = renc1(ext_ind(:)) elsewhere - xind(:) = r2(1,ext_ind(dim,:)) - yind(:) = r2(2,ext_ind(dim,:)) - zind(:) = r2(3,ext_ind(dim,:)) - vxind(:) = v2(1,ext_ind(dim,:)) - vyind(:) = v2(2,ext_ind(dim,:)) - vzind(:) = v2(3,ext_ind(dim,:)) - rencind(:) = renc2(ext_ind(dim,:)) + xind(:) = r2(1,ext_ind(:)) + yind(:) = r2(2,ext_ind(:)) + zind(:) = r2(3,ext_ind(:)) + vxind(:) = v2(1,ext_ind(:)) + vyind(:) = v2(2,ext_ind(:)) + vzind(:) = v2(3,ext_ind(:)) + rencind(:) = renc2(ext_ind(:)) endwhere where(.not.loverlap(:)) lenc(:)%nenc = 0 @@ -875,14 +833,14 @@ module subroutine encounter_check_sweep_aabb_double_list(self, n1, n2, r1, v1, r !$omp do schedule(static) do i = 1, n1 if (loverlap(i)) then - ibeg = self%aabb(dim)%ibeg(i) + 1_I8B - iend = self%aabb(dim)%iend(i) - 1_I8B + ibeg = self%aabb%ibeg(i) + 1_I8B + iend = self%aabb%iend(i) - 1_I8B nbox = iend - ibeg + 1 call encounter_check_all_sweep_one(i, nbox, r1(1,i), r1(2,i), r1(3,i), v1(1,i), v1(2,i), v1(3,i), & xind(ibeg:iend), yind(ibeg:iend), zind(ibeg:iend),& vxind(ibeg:iend), vyind(ibeg:iend), vzind(ibeg:iend), & - renc1(i), rencind(ibeg:iend), dt, ext_ind(dim,ibeg:iend), & - .not.llist1(dim,ibeg:iend), lenc(i)%nenc, lenc(i)%index1, lenc(i)%index2, lenc(i)%lvdotr) + renc1(i), rencind(ibeg:iend), dt, ext_ind(ibeg:iend), & + .not.llist1(ibeg:iend), lenc(i)%nenc, lenc(i)%index1, lenc(i)%index2, lenc(i)%lvdotr) end if end do !$omp end do nowait @@ -891,15 +849,15 @@ module subroutine encounter_check_sweep_aabb_double_list(self, n1, n2, r1, v1, r !$omp do schedule(static) do i = n1+1, ntot if (loverlap(i)) then - ibeg = self%aabb(dim)%ibeg(i) + 1_I8B - iend = self%aabb(dim)%iend(i) - 1_I8B + ibeg = self%aabb%ibeg(i) + 1_I8B + iend = self%aabb%iend(i) - 1_I8B nbox = iend - ibeg + 1 ii = i - n1 call encounter_check_all_sweep_one(ii, nbox, r1(1,ii), r1(2,ii), r1(3,ii), v2(1,ii), v2(2,ii), v2(3,ii), & xind(ibeg:iend), yind(ibeg:iend), zind(ibeg:iend),& vxind(ibeg:iend), vyind(ibeg:iend), vzind(ibeg:iend), & - renc2(ii), rencind(ibeg:iend), dt, ext_ind(dim,ibeg:iend), & - llist1(dim,ibeg:iend), lenc(i)%nenc, lenc(i)%index2, lenc(i)%index1, lenc(i)%lvdotr) + renc2(ii), rencind(ibeg:iend), dt, ext_ind(ibeg:iend), & + llist1(ibeg:iend), lenc(i)%nenc, lenc(i)%index2, lenc(i)%index1, lenc(i)%lvdotr) end if end do !$omp end do nowait @@ -914,7 +872,7 @@ module subroutine encounter_check_sweep_aabb_double_list(self, n1, n2, r1, v1, r end subroutine encounter_check_sweep_aabb_double_list - module subroutine encounter_check_sweep_aabb_single_list(self, n, x, v, renc, dt, nenc, index1, index2, lvdotr) + module subroutine encounter_check_sweep_aabb_single_list(self, n, r, v, renc, dt, nenc, index1, index2, lvdotr) !! author: David A. Minton !! !! Sweeps the sorted bounding box extents and returns the true encounters (combines broad and narrow phases) @@ -923,7 +881,7 @@ module subroutine encounter_check_sweep_aabb_single_list(self, n, x, v, renc, dt ! Arguments class(encounter_bounding_box), intent(inout) :: self !! Multi-dimensional bounding box structure integer(I4B), intent(in) :: n !! Number of bodies - real(DP), dimension(:,:), intent(in) :: x, v !! Array of position and velocity vectors + real(DP), dimension(:,:), intent(in) :: r, v !! Array of position and velocity vectors real(DP), dimension(:), intent(in) :: renc !! Radius of encounter regions of bodies 1 real(DP), intent(in) :: dt !! Step size integer(I8B), intent(out) :: nenc !! Total number of encounter candidates @@ -936,7 +894,7 @@ module subroutine encounter_check_sweep_aabb_single_list(self, n, x, v, renc, dt logical, dimension(n) :: loverlap logical, dimension(2*n) :: lencounteri real(DP), dimension(2*n) :: xind, yind, zind, vxind, vyind, vzind, rencind - integer(I4B), dimension(SWEEPDIM,2*n) :: ext_ind + integer(I4B), dimension(2*n) :: ext_ind type(collision_list_plpl), dimension(n) :: lenc !! Array of encounter lists (one encounter list per body) integer(I4B), dimension(:), allocatable, save :: ind_arr integer(I8B) :: ibeg, iend @@ -946,36 +904,36 @@ module subroutine encounter_check_sweep_aabb_single_list(self, n, x, v, renc, dt ! Sweep the intervals for each of the massive bodies along one dimension ! This will build a ragged pair of index lists inside of the lenc data structure - where(self%aabb(dim)%ind(:) > n) - ext_ind(dim,:) = self%aabb(1)%ind(:) - n + where(self%aabb%ind(:) > n) + ext_ind(:) = self%aabb%ind(:) - n elsewhere - ext_ind(dim,:) = self%aabb(1)%ind(:) + ext_ind(:) = self%aabb%ind(:) endwhere - xind(:) = x(1,ext_ind(dim,:)) - yind(:) = x(2,ext_ind(dim,:)) - zind(:) = x(3,ext_ind(dim,:)) - vxind(:) = v(1,ext_ind(dim,:)) - vyind(:) = v(2,ext_ind(dim,:)) - vzind(:) = v(3,ext_ind(dim,:)) - rencind(:) = renc(ext_ind(dim,:)) + xind(:) = r(1,ext_ind(:)) + yind(:) = r(2,ext_ind(:)) + zind(:) = r(3,ext_ind(:)) + vxind(:) = v(1,ext_ind(:)) + vyind(:) = v(2,ext_ind(:)) + vzind(:) = v(3,ext_ind(:)) + rencind(:) = renc(ext_ind(:)) - loverlap(:) = (self%aabb(dim)%ibeg(:) + 1_I8B) < (self%aabb(dim)%iend(:) - 1_I8B) + loverlap(:) = (self%aabb%ibeg(:) + 1_I8B) < (self%aabb%iend(:) - 1_I8B) where(.not.loverlap(:)) lenc(:)%nenc = 0 !$omp parallel do default(private) schedule(static)& - !$omp shared(self, ext_ind, lenc, loverlap, x, v, renc, xind, yind, zind, vxind, vyind, vzind, rencind) & + !$omp shared(self, ext_ind, lenc, loverlap, r, v, renc, xind, yind, zind, vxind, vyind, vzind, rencind) & !$omp firstprivate(n, dt, dim) do i = 1, n if (loverlap(i)) then - ibeg = self%aabb(dim)%ibeg(i) + 1_I8B - iend = self%aabb(dim)%iend(i) - 1_I8B + ibeg = self%aabb%ibeg(i) + 1_I8B + iend = self%aabb%iend(i) - 1_I8B nbox = int(iend - ibeg + 1, kind=I4B) lencounteri(ibeg:iend) = .true. - call encounter_check_all_sweep_one(i, nbox, x(1,i), x(2,i), x(3,i), v(1,i), v(2,i), v(3,i), & + call encounter_check_all_sweep_one(i, nbox, r(1,i), r(2,i), r(3,i), v(1,i), v(2,i), v(3,i), & xind(ibeg:iend), yind(ibeg:iend), zind(ibeg:iend),& vxind(ibeg:iend), vyind(ibeg:iend), vzind(ibeg:iend), & - renc(i), rencind(ibeg:iend), dt, ext_ind(dim,ibeg:iend), & + renc(i), rencind(ibeg:iend), dt, ext_ind(ibeg:iend), & lencounteri(ibeg:iend), lenc(i)%nenc, lenc(i)%index1, lenc(i)%index2, lenc(i)%lvdotr) end if end do diff --git a/src/encounter/encounter_module.f90 b/src/encounter/encounter_module.f90 index 720627485..f3c24c763 100644 --- a/src/encounter/encounter_module.f90 +++ b/src/encounter/encounter_module.f90 @@ -17,8 +17,8 @@ module encounter implicit none public - integer(I4B), parameter :: SWEEPDIM = 3 character(len=*), parameter :: ENCOUNTER_OUTFILE = 'encounters.nc' !! Name of NetCDF output file for encounter information + real(DP), parameter :: RSWEEP_FACTOR = 1.1_DP type, abstract :: encounter_list integer(I8B) :: nenc = 0 !! Total number of encounters @@ -95,7 +95,7 @@ module encounter type encounter_bounding_box - type(encounter_bounding_box_1D), dimension(SWEEPDIM) :: aabb + type(encounter_bounding_box_1D) :: aabb contains procedure :: dealloc => encounter_util_dealloc_bounding_box !! Deallocates all allocatables procedure :: setup => encounter_util_setup_aabb !! Setup a new axis-aligned bounding box structure @@ -107,12 +107,12 @@ module encounter interface - module subroutine encounter_check_all_plpl(param, npl, x, v, renc, dt, nenc, index1, index2, lvdotr) + module subroutine encounter_check_all_plpl(param, npl, r, v, renc, dt, nenc, index1, index2, lvdotr) use base, only: base_parameters implicit none class(base_parameters), intent(inout) :: param !! Current Swiftest run configuration parameter5s integer(I4B), intent(in) :: npl !! Total number of massive bodies - real(DP), dimension(:,:), intent(in) :: x !! Position vectors of massive bodies + real(DP), dimension(:,:), intent(in) :: r !! Position vectors of massive bodies real(DP), dimension(:,:), intent(in) :: v !! Velocity vectors of massive bodies real(DP), dimension(:), intent(in) :: renc !! Critical radii of massive bodies that defines an encounter real(DP), intent(in) :: dt !! Step size @@ -142,7 +142,7 @@ module subroutine encounter_check_all_plplm(param, nplm, nplt, rplm, vplm, rplt, logical, dimension(:), allocatable, intent(out) :: lvdotr !! Logical flag indicating the sign of v .dot. x end subroutine encounter_check_all_plplm - module subroutine encounter_check_all_pltp(param, npl, ntp, rpl, vpl, xtp, vtp, renc, dt, nenc, index1, index2, lvdotr) + module subroutine encounter_check_all_pltp(param, npl, ntp, rpl, vpl, rtp, vtp, renc, dt, nenc, index1, index2, lvdotr) use base, only: base_parameters implicit none class(base_parameters), intent(inout) :: param !! Current Swiftest run configuration parameter5s @@ -150,7 +150,7 @@ module subroutine encounter_check_all_pltp(param, npl, ntp, rpl, vpl, xtp, vtp, integer(I4B), intent(in) :: ntp !! Total number of test particles real(DP), dimension(:,:), intent(in) :: rpl !! Position vectors of massive bodies real(DP), dimension(:,:), intent(in) :: vpl !! Velocity vectors of massive bodies - real(DP), dimension(:,:), intent(in) :: xtp !! Position vectors of massive bodies + real(DP), dimension(:,:), intent(in) :: rtp !! Position vectors of massive bodies real(DP), dimension(:,:), intent(in) :: vtp !! Velocity vectors of massive bodies real(DP), dimension(:), intent(in) :: renc !! Critical radii of massive bodies that defines an encounter real(DP), intent(in) :: dt !! Step size @@ -205,11 +205,11 @@ module subroutine encounter_check_sweep_aabb_double_list(self, n1, n2, r1, v1, r logical, dimension(:), allocatable, intent(out) :: lvdotr !! Logical array indicating which pairs are approaching end subroutine encounter_check_sweep_aabb_double_list - module subroutine encounter_check_sweep_aabb_single_list(self, n, x, v, renc, dt, nenc, index1, index2, lvdotr) + module subroutine encounter_check_sweep_aabb_single_list(self, n, r, v, renc, dt, nenc, index1, index2, lvdotr) implicit none class(encounter_bounding_box), intent(inout) :: self !! Multi-dimensional bounding box structure integer(I4B), intent(in) :: n !! Number of bodies - real(DP), dimension(:,:), intent(in) :: x, v !! Array of position and velocity vectors + real(DP), dimension(:,:), intent(in) :: r, v !! Array of position and velocity vectors real(DP), dimension(:), intent(in) :: renc !! Radius of encounter regions of bodies 1 real(DP), intent(in) :: dt !! Step size integer(I8B), intent(out) :: nenc !! Total number of encounter candidates diff --git a/src/encounter/encounter_util.f90 b/src/encounter/encounter_util.f90 index 941ed6e16..01848d571 100644 --- a/src/encounter/encounter_util.f90 +++ b/src/encounter/encounter_util.f90 @@ -104,9 +104,7 @@ module subroutine encounter_util_dealloc_bounding_box(self) ! Internals integer(I4B) :: i - do i = 1,NDIM - call self%aabb(i)%dealloc() - end do + call self%aabb%dealloc() return end subroutine encounter_util_dealloc_bounding_box @@ -351,26 +349,20 @@ module subroutine encounter_util_setup_aabb(self, n, n_last) next_last = 2 * n_last if (n > n_last) then ! The number of bodies has grown. Resize and append the new bodies - do dim = 1, SWEEPDIM - allocate(itmp(next)) - if (n_last > 0) itmp(1:next_last) = self%aabb(dim)%ind(1:next_last) - call move_alloc(itmp, self%aabb(dim)%ind) - self%aabb(dim)%ind(next_last+1:next) = [(k, k = next_last+1, next)] - end do + allocate(itmp(next)) + if (n_last > 0) itmp(1:next_last) = self%aabb%ind(1:next_last) + call move_alloc(itmp, self%aabb%ind) + self%aabb%ind(next_last+1:next) = [(k, k = next_last+1, next)] else ! The number of bodies has gone down. Resize and chop of the old indices - do dim = 1, SWEEPDIM - allocate(itmp(next)) - itmp(1:next) = pack(self%aabb(dim)%ind(1:next_last), self%aabb(dim)%ind(1:next_last) <= next) - call move_alloc(itmp, self%aabb(dim)%ind) - end do + allocate(itmp(next)) + itmp(1:next) = pack(self%aabb%ind(1:next_last), self%aabb%ind(1:next_last) <= next) + call move_alloc(itmp, self%aabb%ind) end if - do dim = 1, SWEEPDIM - if (allocated(self%aabb(dim)%ibeg)) deallocate(self%aabb(dim)%ibeg) - allocate(self%aabb(dim)%ibeg(n)) - if (allocated(self%aabb(dim)%iend)) deallocate(self%aabb(dim)%iend) - allocate(self%aabb(dim)%iend(n)) - end do + if (allocated(self%aabb%ibeg)) deallocate(self%aabb%ibeg) + allocate(self%aabb%ibeg(n)) + if (allocated(self%aabb%iend)) deallocate(self%aabb%iend) + allocate(self%aabb%iend(n)) return end subroutine encounter_util_setup_aabb diff --git a/src/rmvs/rmvs_encounter_check.f90 b/src/rmvs/rmvs_encounter_check.f90 index 00aafd1fb..0880cb136 100644 --- a/src/rmvs/rmvs_encounter_check.f90 +++ b/src/rmvs/rmvs_encounter_check.f90 @@ -40,7 +40,7 @@ module function rmvs_encounter_check_tp(self, param, nbody_system, dt) result(le select type(pl => nbody_system%pl) class is (rmvs_pl) - associate(tp => self, ntp => self%nbody, npl => pl%nbody) + associate(tp => self, ntp => self%nbody, npl => pl%nbody, cb => nbody_system%cb) tp%plencP(1:ntp) = 0 call encounter_check_all_pltp(param, npl, ntp, pl%rbeg, pl%vbeg, tp%rh, tp%vh, pl%renc, dt, & nenc, index1, index2, lvdotr) diff --git a/src/swiftest/swiftest_module.f90 b/src/swiftest/swiftest_module.f90 index d360db56b..7913fd1ae 100644 --- a/src/swiftest/swiftest_module.f90 +++ b/src/swiftest/swiftest_module.f90 @@ -140,8 +140,8 @@ module swiftest procedure :: write_info => swiftest_io_netcdf_write_info_body !! Dump contents of particle information metadata to file procedure :: accel_obl => swiftest_obl_acc_body !! Compute the barycentric accelerations of bodies due to the oblateness of the central body procedure :: el2xv => swiftest_orbel_el2xv_vec !! Convert orbital elements to position and velocity vectors - procedure :: xv2el => swiftest_orbel_xv2el_vec !! Convert position and velocity vectors to orbital elements - procedure :: setup => swiftest_util_setup_body !! A constructor that sets the number of bodies and allocates all allocatable arrays + procedure :: xv2el => swiftest_orbel_xv2el_vec !! Convert position and velocity vectors to orbital elements + procedure :: setup => swiftest_util_setup_body !! A constructor that sets the number of bodies and allocates all allocatable arrays procedure :: accel_user => swiftest_user_kick_getacch_body !! Add user-supplied heliocentric accelerations to planets procedure :: append => swiftest_util_append_body !! Appends elements from one structure to another procedure :: dealloc => swiftest_util_dealloc_body !! Deallocates all allocatable arrays diff --git a/src/swiftest/swiftest_orbel.f90 b/src/swiftest/swiftest_orbel.f90 index bf72c91b2..8e351e911 100644 --- a/src/swiftest/swiftest_orbel.f90 +++ b/src/swiftest/swiftest_orbel.f90 @@ -680,7 +680,7 @@ real(DP) pure function swiftest_orbel_fhybrid(e,n) abn = n if(n < 0._DP) abn = -abn - if(abn < 0.636_DP * e -0.6_DP) then + if(abn < 0.636_DP * e - 0.6_DP) then swiftest_orbel_fhybrid = swiftest_orbel_flon(e,n) else swiftest_orbel_fhybrid = swiftest_orbel_fget(e,n) @@ -715,7 +715,7 @@ pure elemental module subroutine swiftest_orbel_xv2aeq(mu, px, py, pz, vx, vy, v q = 0.0_DP x = [px, py, pz] v = [vx, vy, vz] - r = sqrt(dot_product(x(:), x(:))) + r = .mag.x(:) v2 = dot_product(v(:), v(:)) hvec(:) = x(:) .cross. v(:) h2 = dot_product(hvec(:), hvec(:)) diff --git a/src/symba/symba_encounter_check.f90 b/src/symba/symba_encounter_check.f90 index ca7761ff0..895a9f34d 100644 --- a/src/symba/symba_encounter_check.f90 +++ b/src/symba/symba_encounter_check.f90 @@ -249,7 +249,7 @@ module function symba_encounter_check_tp(self, param, nbody_system, dt, irec) re lany_encounter = .false. if (self%nbody == 0) return - associate(tp => self, ntp => self%nbody, pl => nbody_system%pl, npl => nbody_system%pl%nbody) + associate(tp => self, ntp => self%nbody, pl => nbody_system%pl, npl => nbody_system%pl%nbody, cb => nbody_system%cb) call pl%set_renc(irec) call encounter_check_all_pltp(param, npl, ntp, pl%rh, pl%vb, tp%rh, tp%vb, pl%renc, dt, nenc, index1, index2, lvdotr) diff --git a/src/symba/symba_step.f90 b/src/symba/symba_step.f90 index c811b1968..6e72c3e95 100644 --- a/src/symba/symba_step.f90 +++ b/src/symba/symba_step.f90 @@ -27,6 +27,7 @@ module subroutine symba_step_system(self, param, t, dt) real(DP), intent(in) :: dt !! Current stepsize ! Internals logical :: lencounter + type(walltimer) :: timer1,timer2,timer3 !! Object used for computing elapsed wall time select type(pl => self%pl) class is (symba_pl) @@ -48,7 +49,6 @@ module subroutine symba_step_system(self, param, t, dt) end select end select end select - return end subroutine symba_step_system