diff --git a/src/encounter/encounter_check.f90 b/src/encounter/encounter_check.f90 index 69fc05f67..9366e1cf9 100644 --- a/src/encounter/encounter_check.f90 +++ b/src/encounter/encounter_check.f90 @@ -45,7 +45,7 @@ module subroutine encounter_check_all(nenc, index1, index2, x1, v1, x2, v2, renc end subroutine encounter_check_all - module subroutine encounter_check_all_plpl(param, npl, nplm, x, v, renc, dt, lvdotr, index1, index2, nenc) + module subroutine encounter_check_all_plpl(param, npl, x, v, renc, dt, lvdotr, index1, index2, nenc) !! author: David A. Minton !! !! Check for encounters between massive bodies. Choose between the standard triangular or the Sort & Sweep method based on user inputs @@ -54,7 +54,6 @@ module subroutine encounter_check_all_plpl(param, npl, nplm, x, v, renc, dt, lvd ! Arguments class(swiftest_parameters), intent(inout) :: param !! Current Swiftest run configuration parameter5s 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) :: renc !! Critical radii of massive bodies that defines an encounter @@ -67,14 +66,89 @@ module subroutine encounter_check_all_plpl(param, npl, nplm, x, v, renc, dt, lvd type(interaction_timer), save :: itimer logical, save :: lfirst = .true. logical, save :: lsecond = .false. + integer(I8B) :: nplpl = 0_I8B + + if (param%ladaptive_encounters_plpl) then + nplpl = (npl * (npl - 1) / 2) + if (nplpl > 0) then + if (lfirst) then + write(itimer%loopname, *) "encounter_check_all_plpl" + write(itimer%looptype, *) "ENCOUNTER_PLPL" + lfirst = .false. + lsecond = .true. + else + if (lsecond) then ! This ensures that the encounter check methods are run at least once prior to timing. Sort and sweep improves on the second pass due to the bounding box extents needing to be nearly sorted + call itimer%time_this_loop(param, nplpl) + lsecond = .false. + else if (itimer%check(param, nplpl)) then + lsecond = .true. + itimer%is_on = .false. + end if + end if + else + param%lencounter_sas_plpl = .false. + end if + end if + + if (param%lencounter_sas_plpl) then + call encounter_check_all_sort_and_sweep_plpl(npl, x, v, renc, dt, lvdotr, index1, index2, nenc) + else + call encounter_check_all_triangular_plpl(npl, x, v, renc, dt, lvdotr, index1, index2, nenc) + end if + + if (.not.lfirst .and. param%ladaptive_encounters_plpl .and. nplpl > 0) then + if (itimer%is_on) call itimer%adapt(param, nplpl) + end if + + return + end subroutine encounter_check_all_plpl + + + module subroutine encounter_check_all_plplm(param, nplm, nplt, xplm, vplm, xplt, vplt, rencm, renct, dt, lvdotr, index1, index2, nenc) + !! author: David A. Minton + !! + !! Check for encounters between fully interacting massive bodies partially interacting massive bodies. Choose between the standard triangular or the Sort & Sweep method based on user inputs + !! + implicit none + ! Arguments + class(swiftest_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) :: xplm !! Position vectors of fully interacting massive bodies + real(DP), dimension(:,:), intent(in) :: vplm !! Velocity vectors of fully interacting massive bodies + real(DP), dimension(:,:), intent(in) :: xplt !! Position vectors of partially interacting massive bodies + real(DP), dimension(:,:), intent(in) :: vplt !! Velocity vectors of partially interacting massive bodies + real(DP), dimension(:), intent(in) :: rencm !! Critical radii of fully interacting massive bodies that defines an encounter + real(DP), dimension(:), intent(in) :: renct !! Critical radii of partially interacting massive bodies that defines an encounter + real(DP), intent(in) :: dt !! Step size + 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 encounter + integer(I4B), dimension(:), allocatable, intent(out) :: index2 !! List of indices for body 2 in each encounter + integer(I4B), intent(out) :: nenc !! Total number of encounters + ! Internals + type(interaction_timer), save :: itimer + logical, save :: lfirst = .true. + logical, save :: lsecond = .false. integer(I8B) :: nplplm = 0_I8B + integer(I4B) :: npl + logical, dimension(:), allocatable :: plmplt_lvdotr !! Logical flag indicating the sign of v .dot. x in the plm-plt group + integer(I4B), dimension(:), allocatable :: plmplt_index1 !! List of indices for body 1 in each encounter in the plm-plt group + integer(I4B), dimension(:), allocatable :: plmplt_index2 !! List of indices for body 2 in each encounter in the plm-lt group + integer(I4B) :: plmplt_nenc !! Number of encounters of the plm-plt group + class(swiftest_parameters), allocatable :: tmp_param !! Temporary parameter structure to turn off adaptive timer for the pl-pl phase if necessary + integer(I4B), dimension(:), allocatable :: itmp + logical, dimension(:), allocatable :: ltmp - if (param%ladaptive_encounters) then + ! Start with the fully interacting bodies + allocate(tmp_param, source=param) + + if (param%ladaptive_encounters_plpl) then + npl = nplm + nplt nplplm = nplm * npl - nplm * (nplm + 1) / 2 if (nplplm > 0) then if (lfirst) then write(itimer%loopname, *) "encounter_check_all_plpl" - write(itimer%looptype, *) "ENCOUNTER" + write(itimer%looptype, *) "ENCOUNTER_PLPL" lfirst = .false. lsecond = .true. else @@ -87,22 +161,43 @@ module subroutine encounter_check_all_plpl(param, npl, nplm, x, v, renc, dt, lvd end if end if else - param%lencounter_sas = .false. + param%lencounter_sas_plpl = .false. end if + ! Turn off adaptive encounter checks for the pl-pl group + tmp_param%ladaptive_encounters_plpl = .false. end if - if (param%lencounter_sas) then - call encounter_check_all_sort_and_sweep_plpl(npl, nplm, x, v, renc, dt, lvdotr, index1, index2, nenc) + ! Start with the pl-pl group + call encounter_check_all_plpl(tmp_param, nplm, xplm, vplm, rencm, dt, lvdotr, index1, index2, nenc) + + if (param%lencounter_sas_plpl) then + call encounter_check_all_sort_and_sweep_plplm(nplm, nplt, xplm, vplm, xplt, vplt, rencm, renct, dt, plmplt_lvdotr, plmplt_index1, plmplt_index2, plmplt_nenc) else - call encounter_check_all_triangular_plpl(npl, nplm, x, v, renc, dt, lvdotr, index1, index2, nenc) + call encounter_check_all_triangular_plplm(nplm, nplt, xplm, vplm, xplt, vplt, rencm, renct, dt, plmplt_lvdotr, plmplt_index1, plmplt_index2, plmplt_nenc) end if - if (.not.lfirst .and. param%ladaptive_encounters .and. nplplm > 0) then + if (.not.lfirst .and. param%ladaptive_encounters_plpl .and. nplplm > 0) then if (itimer%is_on) call itimer%adapt(param, nplplm) end if + if (plmplt_nenc > 0) then ! Consolidate the two lists + allocate(itmp(nenc+plmplt_nenc)) + itmp(1:nenc) = index1(1:nenc) + itmp(nenc+1:nenc+plmplt_nenc) = plmplt_index1(1:plmplt_nenc) + call move_alloc(itmp, index1) + allocate(itmp(nenc+plmplt_nenc)) + itmp(1:nenc) = index2(1:nenc) + itmp(nenc+1:nenc+plmplt_nenc) = plmplt_index2(1:plmplt_nenc) + call move_alloc(itmp, index2) + allocate(ltmp(nenc+plmplt_nenc)) + ltmp(1:nenc) = lvdotr(1:nenc) + ltmp(nenc+1:nenc+plmplt_nenc) = plmplt_lvdotr(1:plmplt_nenc) + call move_alloc(ltmp, lvdotr) + nenc = nenc + plmplt_nenc + end if + return - end subroutine encounter_check_all_plpl + end subroutine encounter_check_all_plplm module subroutine encounter_check_all_pltp(param, npl, ntp, xpl, vpl, xtp, vtp, renc, dt, lvdotr, index1, index2, nenc) @@ -117,18 +212,51 @@ module subroutine encounter_check_all_pltp(param, npl, ntp, xpl, vpl, xtp, vtp, integer(I4B), intent(in) :: ntp !! Total number of test particles real(DP), dimension(:,:), intent(in) :: xpl !! 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) :: vtp !! Velocity vectors of massive bodies + real(DP), dimension(:,:), intent(in) :: xtp !! 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 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 encounter integer(I4B), dimension(:), allocatable, intent(out) :: index2 !! List of indices for body 2 in each encounter integer(I4B), intent(out) :: nenc !! Total number of encounters + ! Internals + type(interaction_timer), save :: itimer + logical, save :: lfirst = .true. + logical, save :: lsecond = .false. + integer(I8B) :: npltp = 0_I8B + if (param%ladaptive_encounters_pltp) then + npltp = npl * ntp + if (npltp > 0) then + if (lfirst) then + write(itimer%loopname, *) "encounter_check_all_pltp" + write(itimer%looptype, *) "ENCOUNTER_PLTP" + lfirst = .false. + lsecond = .true. + else + if (lsecond) then ! This ensures that the encounter check methods are run at least once prior to timing. Sort and sweep improves on the second pass due to the bounding box extents needing to be nearly sorted + call itimer%time_this_loop(param, npltp) + lsecond = .false. + else if (itimer%check(param, npltp)) then + lsecond = .true. + itimer%is_on = .false. + end if + end if + else + param%lencounter_sas_pltp = .false. + end if + end if - !call encounter_check_all_triangular_pltp(npl, ntp, xpl, vpl, xtp, vtp, renc, dt, lvdotr, index1, index2, nenc) - call encounter_check_all_sort_and_sweep_pltp(npl, ntp, xpl, vpl, xtp, vtp, renc, dt, lvdotr, index1, index2, nenc) + if (param%lencounter_sas_pltp) then + call encounter_check_all_sort_and_sweep_pltp(npl, ntp, xpl, vpl, xtp, vtp, renc, dt, lvdotr, index1, index2, nenc) + else + call encounter_check_all_triangular_pltp(npl, ntp, xpl, vpl, xtp, vtp, renc, dt, lvdotr, index1, index2, nenc) + end if + + if (.not.lfirst .and. param%ladaptive_encounters_pltp .and. npltp > 0) then + if (itimer%is_on) call itimer%adapt(param, npltp) + end if return end subroutine encounter_check_all_pltp @@ -206,7 +334,7 @@ subroutine encounter_check_reduce_broadphase(n, nenc, index1, index2, lencounter end subroutine encounter_check_reduce_broadphase - subroutine encounter_check_all_sort_and_sweep_plpl(npl, nplm, x, v, renc, dt, lvdotr, index1, index2, nenc) + subroutine encounter_check_all_sort_and_sweep_plpl(npl, x, v, renc, dt, lvdotr, index1, index2, nenc) !! author: David A. Minton !! !! Check for encounters between massive bodies. @@ -216,7 +344,6 @@ subroutine encounter_check_all_sort_and_sweep_plpl(npl, nplm, x, v, renc, dt, lv implicit none ! Arguments 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) :: renc !! Critical radii of massive bodies that defines an encounter @@ -277,6 +404,101 @@ subroutine encounter_check_all_sort_and_sweep_plpl(npl, nplm, x, v, renc, dt, lv end subroutine encounter_check_all_sort_and_sweep_plpl + subroutine encounter_check_all_sort_and_sweep_plplm(nplm, nplt, xplm, vplm, xplt, vplt, rencm, renct, dt, lvdotr, index1, index2, nenc) + !! author: David A. Minton + !! + !! Check for encounters between massive bodies and test particles. + !! This is the sort and sweep version + !! References: Adapted from Christer Ericson's _Real Time Collision Detection_ + !! + implicit none + ! Arguments + 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) :: xplm !! Position vectors of fully interacting massive bodies + real(DP), dimension(:,:), intent(in) :: vplm !! Velocity vectors of fully interacting massive bodies + real(DP), dimension(:,:), intent(in) :: xplt !! Position vectors of partially interacting massive bodies + real(DP), dimension(:,:), intent(in) :: vplt !! Velocity vectors of partially interacting massive bodies + real(DP), dimension(:), intent(in) :: rencm !! Critical radii of fully interacting massive bodies that defines an encounter + real(DP), dimension(:), intent(in) :: renct !! Critical radii of partially interacting massive bodies that defines an encounter + real(DP), intent(in) :: dt !! Step size + 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 encounter + integer(I4B), dimension(:), allocatable, intent(out) :: index2 !! List of indices for body 2 in each encounter + integer(I4B), intent(out) :: nenc !! Total number of encounter + ! Internals + type(encounter_bounding_box), save :: boundingbox + integer(I4B) :: i, dim, n, ntot + integer(I4B), dimension(:), allocatable, save :: ind_arr + integer(I4B), save :: ntot_last = 0 + logical, dimension(:), allocatable :: lencounter + integer(I2B), dimension(nplm) :: vplmshift_min, vplmshift_max + integer(I2B), dimension(nplt) :: vpltshift_min, vpltshift_max + + ! If this is the first time through, build the index lists + if ((nplm == 0) .or. (nplt == 0)) return + + ntot = nplm + nplt + n = 2 * ntot + if (ntot /= ntot_last) then + if (allocated(ind_arr)) deallocate(ind_arr) + allocate(ind_arr(ntot)) + ind_arr(1:ntot) = [(i, i = 1, ntot)] + + call boundingbox%setup(ntot, ntot_last) + + ntot_last = ntot + end if + + !$omp taskloop default(private) num_tasks(SWEEPDIM) & + !$omp shared(xplm, xplt, 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, [xplm(dim,1:nplm) - rencm(1:nplm) + vplmshift_min(1:nplm) * vplm(dim,1:nplm) * dt, & + xplm(dim,1:nplt) - renct(1:nplt) + vpltshift_min(1:nplt) * vplt(dim,1:nplt) * dt, & + xplm(dim,1:nplm) + rencm(1:nplm) + vplmshift_max(1:nplm) * vplm(dim,1:nplm) * dt, & + xplt(dim,1:nplt) + renct(1:nplt) + vpltshift_max(1:nplt) * vplt(dim,1:nplt) * dt]) + end do + !$omp end taskloop + + call boundingbox%sweep(nplm, nplt, ind_arr, nenc, index1, index2) + + if (nenc > 0) then + ! Shift tiny body indices back into the range of the input position and velocity arrays + index2(:) = index2(:) - nplm + + ! Now that we have identified potential pairs, use the narrow-phase process to get the final values + allocate(lencounter(nenc)) + allocate(lvdotr(nenc)) + + call encounter_check_all(nenc, index1, index2, xplm, vplm, xplt, vplt, rencm, renct, dt, lencounter, lvdotr) + + call encounter_check_reduce_broadphase(ntot, nenc, index1, index2, lencounter, lvdotr) + + ! SHift the tiny body indices back to their natural range + index2(:) = index2(:) + nplm + end if + + return + end subroutine encounter_check_all_sort_and_sweep_plplm + + subroutine encounter_check_all_sort_and_sweep_pltp(npl, ntp, xpl, vpl, xtp, vtp, renc, dt, lvdotr, index1, index2, nenc) !! author: David A. Minton !! @@ -371,7 +593,7 @@ subroutine encounter_check_all_sort_and_sweep_pltp(npl, ntp, xpl, vpl, xtp, vtp, end subroutine encounter_check_all_sort_and_sweep_pltp - subroutine encounter_check_all_triangular_plpl(npl, nplm, x, v, renc, dt, lvdotr, index1, index2, nenc) + subroutine encounter_check_all_triangular_plpl(npl, x, v, renc, dt, lvdotr, index1, index2, nenc) !! author: David A. Minton !! !! Check for encounters between massive bodies. Split off from the main subroutine for performance @@ -379,7 +601,6 @@ subroutine encounter_check_all_triangular_plpl(npl, nplm, x, v, renc, dt, lvdotr implicit none ! Arguments 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) :: renc !! Critical radii of massive bodies that defines an encounter @@ -393,13 +614,13 @@ subroutine encounter_check_all_triangular_plpl(npl, nplm, x, v, renc, dt, lvdotr real(DP) :: xr, yr, zr, vxr, vyr, vzr, renc12 logical, dimension(npl) :: lencounteri, lvdotri integer(I4B), dimension(npl) :: ind_arr - type(encounter_list), dimension(nplm) :: lenc + type(encounter_list), dimension(npl) :: lenc ind_arr(:) = [(i, i = 1, npl)] !$omp parallel do default(private) schedule(static)& - !$omp shared(npl, nplm, x, v, renc, dt, lenc, ind_arr) & + !$omp shared(npl, x, v, renc, dt, lenc, ind_arr) & !$omp lastprivate(xr, yr, zr, vxr, vyr, vzr, renc12, lencounteri, lvdotri) - do i = 1, nplm + do i = 1, npl do concurrent(j = i+1:npl) xr = x(1, j) - x(1, i) yr = x(2, j) - x(2, i) @@ -421,14 +642,14 @@ subroutine encounter_check_all_triangular_plpl(npl, nplm, x, v, renc, dt, lvdotr !$omp end parallel do associate(nenc_arr => lenc(:)%nenc) - nenc = sum(nenc_arr(1:nplm)) + nenc = sum(nenc_arr(1:npl)) end associate if (nenc > 0) then allocate(lvdotr(nenc)) allocate(index1(nenc)) allocate(index2(nenc)) j0 = 1 - do i = 1, nplm + do i = 1, npl if (lenc(i)%nenc > 0) then j1 = j0 + lenc(i)%nenc - 1 lvdotr(j0:j1) = lenc(i)%lvdotr(:) @@ -443,6 +664,81 @@ subroutine encounter_check_all_triangular_plpl(npl, nplm, x, v, renc, dt, lvdotr end subroutine encounter_check_all_triangular_plpl + subroutine encounter_check_all_triangular_plplm(nplm, nplt, xplm, vplm, xplt, vplt, rencm, renct, dt, lvdotr, index1, index2, nenc) + !! author: David A. Minton + !! + !! Check for encounters between massive bodies and test particles. Split off from the main subroutine for performance + !! This is the upper triangular (double loop) version. + implicit none + ! Arguments + 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) :: xplm !! Position vectors of fully interacting massive bodies + real(DP), dimension(:,:), intent(in) :: vplm !! Velocity vectors of fully interacting massive bodies + real(DP), dimension(:,:), intent(in) :: xplt !! Position vectors of partially interacting massive bodies + real(DP), dimension(:,:), intent(in) :: vplt !! Velocity vectors of partially interacting massive bodies + real(DP), dimension(:), intent(in) :: rencm !! Critical radii of fully interacting massive bodies that defines an encounter + real(DP), dimension(:), intent(in) :: renct !! Critical radii of partially interacting massive bodies that defines an encounter + real(DP), intent(in) :: dt !! Step size + 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 encounter + integer(I4B), dimension(:), allocatable, intent(out) :: index2 !! List of indices for body 2 in each encounter + integer(I4B), intent(out) :: nenc !! Total number of encounters + ! Internals + integer(I4B) :: i, j, nenci, j0, j1 + real(DP) :: xr, yr, zr, vxr, vyr, vzr, renc12 + logical, dimension(nplt) :: lencounteri, lvdotri + integer(I4B), dimension(nplt) :: ind_arr + type(encounter_list), dimension(nplm) :: lenc + + ind_arr(:) = [(i, i = 1, nplt)] + !$omp parallel do default(private) schedule(static)& + !$omp shared(nplm, nplt, xplm, vplm, xplt, vplt, rencm, renct, dt, lenc, ind_arr) & + !$omp lastprivate(xr, yr, zr, vxr, vyr, vzr, renc12, lencounteri, lvdotri) + do i = 1, nplm + do concurrent(j = 1:nplt) + xr = xplt(1, j) - xplm(1, i) + yr = xplt(2, j) - xplm(2, i) + zr = xplt(3, j) - xplm(3, i) + vxr = vplt(1, j) - vplm(1, i) + vyr = vplt(2, j) - vplm(2, i) + vzr = vplt(3, j) - vplm(3, i) + renc12 = rencm(i) + renct(j) + call encounter_check_one(xr, yr, zr, vxr, vyr, vzr, renc12, dt, lencounteri(j), lvdotri(j)) + end do + nenci = count(lencounteri(1:nplt)) + if (nenci > 0) then + allocate(lenc(i)%lvdotr(nenci), lenc(i)%index2(nenci)) + lenc(i)%nenc = nenci + lenc(i)%lvdotr(:) = pack(lvdotri(1:nplt), lencounteri(1:nplt)) + lenc(i)%index2(:) = pack(ind_arr(1:nplt), lencounteri(1:nplt)) + end if + end do + !$omp end parallel do + + associate(nenc_arr => lenc(:)%nenc) + nenc = sum(nenc_arr(1:nplm)) + end associate + if (nenc > 0) then + 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 + lvdotr(j0:j1) = lenc(i)%lvdotr(:) + index1(j0:j1) = i + index2(j0:j1) = lenc(i)%index2(:) + nplm + j0 = j1 + 1 + end if + end do + end if + + return + end subroutine encounter_check_all_triangular_plplm + + subroutine encounter_check_all_triangular_pltp(npl, ntp, xpl, vpl, xtp, vtp, renc, dt, lvdotr, index1, index2, nenc) !! author: David A. Minton !! diff --git a/src/io/io.f90 b/src/io/io.f90 index 93d86c2f0..f0ff6f55c 100644 --- a/src/io/io.f90 +++ b/src/io/io.f90 @@ -650,9 +650,16 @@ module subroutine io_param_reader(self, unit, iotype, v_list, iostat, iomsg) case ("INTERACTION_LOOPS") call io_toupper(param_value) param%interaction_loops = param_value + case ("ENCOUNTER_CHECK_PLPL") + call io_toupper(param_value) + param%encounter_check_plpl = param_value + case ("ENCOUNTER_CHECK_PLTP") + call io_toupper(param_value) + param%encounter_check_pltp = param_value case ("ENCOUNTER_CHECK") call io_toupper(param_value) - param%encounter_check = param_value + param%encounter_check_plpl = param_value + param%encounter_check_pltp = param_value case ("FIRSTKICK") call io_toupper(param_value) if (param_value == "NO" .or. param_value == 'F') param%lfirstkick = .false. @@ -836,28 +843,50 @@ module subroutine io_param_reader(self, unit, iotype, v_list, iostat, iomsg) call io_log_one_message(INTERACTION_TIMER_LOG_OUT, "Diagnostic values: loop style, time count, nplpl, metric") end select + select case(trim(adjustl(param%encounter_check_plpl))) + case("ADAPTIVE") + param%ladaptive_encounters_plpl = .true. + param%lencounter_sas_plpl = .true. + call io_log_start(param, ENCOUNTER_PLPL_TIMER_LOG_OUT, "Encounter check loop timer logfile") + call io_log_one_message(ENCOUNTER_PLPL_TIMER_LOG_OUT, "Diagnostic values: loop style, time count, nplpl, metric") + case("TRIANGULAR") + param%ladaptive_encounters_plpl = .false. + param%lencounter_sas_plpl = .false. + case("SORTSWEEP") + param%ladaptive_encounters_plpl = .false. + param%lencounter_sas_plpl = .true. + case default + write(*,*) "Unknown value for parameter ENCOUNTER_CHECK_PLPL: -> ",trim(adjustl(param%encounter_check_plpl)) + write(*,*) "Must be one of the following: TRIANGULAR, SORTSWEEP, or ADAPTIVE" + write(*,*) "Using default value of ADAPTIVE" + param%encounter_check_plpl = "ADAPTIVE" + param%ladaptive_encounters_plpl = .true. + param%lencounter_sas_plpl = .true. + call io_log_start(param, ENCOUNTER_PLPL_TIMER_LOG_OUT, "Encounter check loop timer logfile") + call io_log_one_message(ENCOUNTER_PLPL_TIMER_LOG_OUT, "Diagnostic values: loop style, time count, nplpl, metric") + end select - select case(trim(adjustl(param%encounter_check))) + select case(trim(adjustl(param%encounter_check_pltp))) case("ADAPTIVE") - param%ladaptive_encounters = .true. - param%lencounter_sas = .true. - call io_log_start(param, ENCOUNTER_TIMER_LOG_OUT, "Encounter check loop timer logfile") - call io_log_one_message(ENCOUNTER_TIMER_LOG_OUT, "Diagnostic values: loop style, time count, nplpl, metric") + param%ladaptive_encounters_pltp = .true. + param%lencounter_sas_pltp = .true. + call io_log_start(param, ENCOUNTER_PLTP_TIMER_LOG_OUT, "Encounter check loop timer logfile") + call io_log_one_message(ENCOUNTER_PLTP_TIMER_LOG_OUT, "Diagnostic values: loop style, time count, npltp, metric") case("TRIANGULAR") - param%ladaptive_encounters = .false. - param%lencounter_sas = .false. + param%ladaptive_encounters_pltp = .false. + param%lencounter_sas_pltp = .false. case("SORTSWEEP") - param%ladaptive_encounters = .false. - param%lencounter_sas = .true. + param%ladaptive_encounters_pltp = .false. + param%lencounter_sas_pltp = .true. case default - write(*,*) "Unknown value for parameter ENCOUNTER_CHECK: -> ",trim(adjustl(param%encounter_check)) + write(*,*) "Unknown value for parameter ENCOUNTER_CHECK_PLTP: -> ",trim(adjustl(param%encounter_check_pltp)) write(*,*) "Must be one of the following: TRIANGULAR, SORTSWEEP, or ADAPTIVE" write(*,*) "Using default value of ADAPTIVE" - param%encounter_check = "ADAPTIVE" - param%ladaptive_encounters = .true. - param%lencounter_sas = .true. - call io_log_start(param, ENCOUNTER_TIMER_LOG_OUT, "Encounter check loop timer logfile") - call io_log_one_message(ENCOUNTER_TIMER_LOG_OUT, "Diagnostic values: loop style, time count, nplpl, metric") + param%encounter_check_pltp = "ADAPTIVE" + param%ladaptive_encounters_pltp = .true. + param%lencounter_sas_pltp = .true. + call io_log_start(param, ENCOUNTER_PLTP_TIMER_LOG_OUT, "Encounter check loop timer logfile") + call io_log_one_message(ENCOUNTER_PLTP_TIMER_LOG_OUT, "Diagnostic values: loop style, time count, npltp, metric") end select iostat = 0 @@ -950,7 +979,8 @@ module subroutine io_param_writer(self, unit, iotype, v_list, iostat, iomsg) call io_param_writer_one("ROTATION", param%lrotation, unit) call io_param_writer_one("TIDES", param%ltides, unit) call io_param_writer_one("INTERACTION_LOOPS", param%interaction_loops, unit) - call io_param_writer_one("ENCOUNTER_CHECK", param%encounter_check, unit) + call io_param_writer_one("ENCOUNTER_CHECK_PLPL", param%encounter_check_plpl, unit) + call io_param_writer_one("ENCOUNTER_CHECK_PLTP", param%encounter_check_pltp, unit) if (param%lenergy) then call io_param_writer_one("FIRSTENERGY", param%lfirstenergy, unit) diff --git a/src/modules/encounter_classes.f90 b/src/modules/encounter_classes.f90 index 971423e32..81a1c0e12 100644 --- a/src/modules/encounter_classes.f90 +++ b/src/modules/encounter_classes.f90 @@ -64,12 +64,11 @@ module subroutine encounter_check_all(nenc, index1, index2, x1, v1, x2, v2, renc logical, dimension(:), intent(out) :: lvdotr !! Logical array indicating which pairs are approaching end subroutine encounter_check_all - module subroutine encounter_check_all_plpl(param, npl, nplm, x, v, renc, dt, lvdotr, index1, index2, nenc) + module subroutine encounter_check_all_plpl(param, npl, x, v, renc, dt, lvdotr, index1, index2, nenc) import swiftest_parameters implicit none class(swiftest_parameters), intent(inout) :: param !! Current Swiftest run configuration parameter5s 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) :: renc !! Critical radii of massive bodies that defines an encounter @@ -80,6 +79,25 @@ module subroutine encounter_check_all_plpl(param, npl, nplm, x, v, renc, dt, lvd integer(I4B), intent(out) :: nenc !! Total number of encounters end subroutine encounter_check_all_plpl + module subroutine encounter_check_all_plplm(param, nplm, nplt, xplm, vplm, xplt, vplt, rencm, renct, dt, lvdotr, index1, index2, nenc) + import swiftest_parameters + implicit none + class(swiftest_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) :: xplm !! Position vectors of fully interacting massive bodies + real(DP), dimension(:,:), intent(in) :: vplm !! Velocity vectors of fully interacting massive bodies + real(DP), dimension(:,:), intent(in) :: xplt !! Position vectors of partially interacting massive bodies + real(DP), dimension(:,:), intent(in) :: vplt !! Velocity vectors of partially interacting massive bodies + real(DP), dimension(:), intent(in) :: rencm !! Critical radii of fully interacting massive bodies that defines an encounter + real(DP), dimension(:), intent(in) :: renct !! Critical radii of partially interacting massive bodies that defines an encounter + real(DP), intent(in) :: dt !! Step size + 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 encounter + integer(I4B), dimension(:), allocatable, intent(out) :: index2 !! List of indices for body 2 in each encounter + integer(I4B), intent(out) :: nenc !! Total number of encounters + end subroutine encounter_check_all_plplm + module subroutine encounter_check_all_pltp(param, npl, ntp, xpl, vpl, xtp, vtp, renc, dt, lvdotr, index1, index2, nenc) import swiftest_parameters implicit none diff --git a/src/modules/swiftest_classes.f90 b/src/modules/swiftest_classes.f90 index 5c3f25f21..cd4484f8f 100644 --- a/src/modules/swiftest_classes.f90 +++ b/src/modules/swiftest_classes.f90 @@ -126,12 +126,15 @@ module swiftest_classes real(DP) :: inv_c2 = -1.0_DP !! Inverse speed of light squared in the system units character(STRMAX) :: energy_out = "" !! Name of output energy and momentum report file character(NAMELEN) :: interaction_loops = "ADAPTIVE" !! Method used to compute interaction loops. Options are "TRIANGULAR", "FLAT", or "ADAPTIVE" - character(NAMELEN) :: encounter_check = "ADAPTIVE" !! Method used to compute encounter checks. Options are "TRIANGULAR", "SORTSWEEP", or "ADAPTIVE" + character(NAMELEN) :: encounter_check_plpl = "ADAPTIVE" !! Method used to compute pl-pl encounter checks. Options are "TRIANGULAR", "SORTSWEEP", or "ADAPTIVE" + character(NAMELEN) :: encounter_check_pltp = "ADAPTIVE" !! Method used to compute pl-tp encounter checks. Options are "TRIANGULAR", "SORTSWEEP", or "ADAPTIVE" ! The following are used internally, and are not set by the user, but instead are determined by the input value of INTERACTION_LOOPS logical :: lflatten_interactions = .false. !! Use the flattened upper triangular matrix for pl-pl interaction loops logical :: ladaptive_interactions = .false. !! Adaptive interaction loop is turned on (choose between TRIANGULAR and FLAT based on periodic timing tests) - logical :: lencounter_sas = .false. !! Use the Sort and Sweep algorithm to prune the encounter list before checking for close encounters - logical :: ladaptive_encounters = .false. !! Adaptive encounter checking is turned on (choose between TRIANGULAR or SORTSWEEP based on periodic timing tests) + logical :: lencounter_sas_plpl = .false. !! Use the Sort and Sweep algorithm to prune the encounter list before checking for close encounters + logical :: lencounter_sas_pltp = .false. !! Use the Sort and Sweep algorithm to prune the encounter list before checking for close encounters + logical :: ladaptive_encounters_plpl = .false. !! Adaptive encounter checking is turned on (choose between TRIANGULAR or SORTSWEEP based on periodic timing tests) + logical :: ladaptive_encounters_pltp = .false. !! Adaptive encounter checking is turned on (choose between TRIANGULAR or SORTSWEEP based on periodic timing tests) ! Logical flags to turn on or off various features of the code logical :: lrhill_present = .false. !! Hill radii are given as an input rather than calculated by the code (can be used to inflate close encounter regions manually) diff --git a/src/modules/walltime_classes.f90 b/src/modules/walltime_classes.f90 index a92fc1183..7f1fcabc2 100644 --- a/src/modules/walltime_classes.f90 +++ b/src/modules/walltime_classes.f90 @@ -9,7 +9,8 @@ module walltime_classes integer(I4B) :: INTERACTION_TIMER_CADENCE = 1000 !! Minimum number of steps to wait before timing an interaction loop in ADAPTIVE mode character(len=*), parameter :: INTERACTION_TIMER_LOG_OUT = "interaction_timer.log" !! Name of log file for recording results of interaction loop timing - character(len=*), parameter :: ENCOUNTER_TIMER_LOG_OUT = "encounter_check_timer.log" !! Name of log file for recording results of encounter check method timing + character(len=*), parameter :: ENCOUNTER_PLPL_TIMER_LOG_OUT = "encounter_check_plpl_timer.log" !! Name of log file for recording results of encounter check method timing + character(len=*), parameter :: ENCOUNTER_PLTP_TIMER_LOG_OUT = "encounter_check_pltp_timer.log" !! Name of log file for recording results of encounter check method timing type :: walltimer integer(I8B) :: count_rate !! Rate at wich the clock ticks diff --git a/src/symba/symba_encounter_check.f90 b/src/symba/symba_encounter_check.f90 index 484667ae8..84c4fc3f0 100644 --- a/src/symba/symba_encounter_check.f90 +++ b/src/symba/symba_encounter_check.f90 @@ -18,7 +18,7 @@ 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, nplm + integer(I4B) :: i, j, nenc, npl, nplm, nplt logical, dimension(:), allocatable :: lencounter, loc_lvdotr, lvdotr integer(I4B), dimension(:), allocatable :: index1, index2 integer(I4B), dimension(:,:), allocatable :: k_plpl_enc @@ -31,7 +31,12 @@ module function symba_encounter_check_pl(self, param, system, dt, irec) result(l npl = pl%nbody nplm = pl%nplm - call encounter_check_all_plpl(param, npl, nplm, pl%xh, pl%vh, pl%renc, dt, lvdotr, index1, index2, nenc) + nplt = npl - nplm + if (nplt == 0) then + call encounter_check_all_plpl(param, npl, pl%xh, pl%vh, pl%renc, dt, lvdotr, index1, index2, nenc) + else + call encounter_check_all_plplm(param, nplm, nplt, pl%xh(:,1:nplm), pl%vh(:,1:nplm), pl%xh(:,nplm+1:npl), pl%vh(:,nplm+1:npl), pl%renc(1:nplm), pl%renc(nplm+1:npl), dt, lvdotr, index1, index2, nenc) + end if lany_encounter = nenc > 0 if (lany_encounter) then call plplenc_list%resize(nenc) diff --git a/src/symba/symba_util.f90 b/src/symba/symba_util.f90 index a4ad56eec..656aed3f8 100644 --- a/src/symba/symba_util.f90 +++ b/src/symba/symba_util.f90 @@ -290,6 +290,10 @@ module subroutine symba_util_flatten_eucl_plpl(self, param) associate(pl => self, nplpl => self%nplpl, 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) / 2) ! number of entries in a strict lower triangle, npl x npl, minus first column diff --git a/src/walltime/walltime.f90 b/src/walltime/walltime.f90 index 107e938a3..88d791489 100644 --- a/src/walltime/walltime.f90 +++ b/src/walltime/walltime.f90 @@ -155,7 +155,7 @@ module subroutine walltime_interaction_adapt(self, param, ninteractions, pl) case("ENCOUNTER") write(advancedstyle, *) "SORTSWEEP " write(standardstyle, *) "TRIANGULAR" - write(logfile,*) ENCOUNTER_TIMER_LOG_OUT + write(logfile,*) ENCOUNTER_PLPL_TIMER_LOG_OUT case default write(logfile,*) "unknown_looptimer.log" end select @@ -251,8 +251,10 @@ module subroutine walltime_interaction_flip_loop_style(self, param, pl) select case(trim(adjustl(self%looptype))) case("INTERACTION") param%lflatten_interactions = .not. param%lflatten_interactions - case("ENCOUNTER") - param%lencounter_sas = .not. param%lencounter_sas + case("ENCOUNTER_PLPL") + param%lencounter_sas_plpl= .not. param%lencounter_sas_plpl + case("ENCOUNTER_PLTP") + param%lencounter_sas_pltp= .not. param%lencounter_sas_pltp end select if (present(pl)) then @@ -286,7 +288,7 @@ module subroutine walltime_interaction_time_this_loop(self, param, ninteractions case("INTERACTION") write(logfile,*) INTERACTION_TIMER_LOG_OUT case("ENCOUNTER") - write(logfile,*) ENCOUNTER_TIMER_LOG_OUT + write(logfile,*) ENCOUNTER_PLPL_TIMER_LOG_OUT case default write(logfile,*) "unknown_looptimer.log" end select @@ -299,16 +301,20 @@ module subroutine walltime_interaction_time_this_loop(self, param, ninteractions select case(trim(adjustl(self%looptype))) case("INTERACTION") self%stage1_is_advanced = param%lflatten_interactions - case("ENCOUNTER") - self%stage1_is_advanced = param%lencounter_sas + case("ENCOUNTER_PLPL") + self%stage1_is_advanced = param%lencounter_sas_plpl + case("ENCOUNTER_PLTP") + self%stage1_is_advanced = param%lencounter_sas_pltp end select call io_log_one_message(logfile, trim(adjustl(self%loopname)) // ": loop timer turned on at t = " // trim(adjustl(tstr))) case(2) select case(trim(adjustl(self%looptype))) case("INTERACTION") param%lflatten_interactions = self%stage1_is_advanced - case("ENCOUNTER") - param%lencounter_sas = self%stage1_is_advanced + case("ENCOUNTER_PLPL") + param%lencounter_sas_plpl= self%stage1_is_advanced + case("ENCOUNTER_PLTP") + param%lencounter_sas_pltp= self%stage1_is_advanced end select call self%flip(param, pl) case default