diff --git a/src/encounter/encounter_check.f90 b/src/encounter/encounter_check.f90 index 4b0e91e13..516fb84a9 100644 --- a/src/encounter/encounter_check.f90 +++ b/src/encounter/encounter_check.f90 @@ -360,25 +360,24 @@ subroutine encounter_check_all_sort_and_sweep_plpl(npl, x, v, renc, dt, lvdotr, integer(I4B), intent(out) :: nenc !! Total number of encounters ! Internals integer(I4B) :: i, dim, n - integer(I4B), dimension(:), allocatable, save :: ind_arr integer(I4B), save :: npl_last = 0 type(encounter_bounding_box), save :: boundingbox logical, dimension(:), allocatable :: lencounter integer(I2B), dimension(npl) :: vshift_min, vshift_max + type(walltimer) :: timer if (npl <= 1) return ! If this is the first time through, build the index lists n = 2 * npl if (npl_last /= npl) then - if (allocated(ind_arr)) deallocate(ind_arr) - allocate(ind_arr(npl)) - ind_arr(:) = [(i, i = 1, npl)] - call boundingbox%setup(npl, npl_last) + npl_last = npl end if - !$omp taskloop default(private) num_tasks(SWEEPDIM) & + ! call timer%reset() + ! call timer%start() + !$omp parallel do default(private) schedule(static) & !$omp shared(x, v, renc, boundingbox) & !$omp firstprivate(dt, npl, n) do dim = 1, SWEEPDIM @@ -392,9 +391,11 @@ subroutine encounter_check_all_sort_and_sweep_plpl(npl, x, v, renc, dt, lvdotr, 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]) end do - !$omp end taskloop + !$omp end parallel do + ! call timer%stop() + ! call timer%report(nsubsteps=nthreads, message="AABB sort plpl:") - call boundingbox%sweep(npl, ind_arr, nenc, index1, index2) + call boundingbox%sweep(npl, nenc, index1, index2) if (nenc > 0) then ! Now that we have identified potential pairs, use the narrow-phase process to get the final values @@ -435,11 +436,11 @@ subroutine encounter_check_all_sort_and_sweep_plplm(nplm, nplt, xplm, vplm, xplt ! 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 + type(walltimer) :: timer ! If this is the first time through, build the index lists if ((nplm == 0) .or. (nplt == 0)) return @@ -447,16 +448,15 @@ subroutine encounter_check_all_sort_and_sweep_plplm(nplm, nplt, xplm, vplm, xplt 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) & + + ! call timer%reset() + ! call timer%start() + !$omp parallel do default(private) schedule(static) & !$omp shared(xplm, xplt, vplm, vplt, rencm, renct, boundingbox) & !$omp firstprivate(dt, nplm, nplt, ntot) do dim = 1, SWEEPDIM @@ -481,9 +481,11 @@ subroutine encounter_check_all_sort_and_sweep_plplm(nplm, nplt, xplm, vplm, xplt 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 + !$omp end parallel do + ! call timer%stop() + ! call timer%report(nsubsteps=nthreads, message="AABB sort plplm:") - call boundingbox%sweep(nplm, nplt, ind_arr, nenc, index1, index2) + call boundingbox%sweep(nplm, nplt, nenc, index1, index2) if (nenc > 0) then ! Shift tiny body indices back into the range of the input position and velocity arrays @@ -530,7 +532,6 @@ subroutine encounter_check_all_sort_and_sweep_pltp(npl, ntp, xpl, vpl, xtp, vtp, ! 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(npl) :: vplshift_min, vplshift_max @@ -543,16 +544,11 @@ subroutine encounter_check_all_sort_and_sweep_pltp(npl, ntp, xpl, vpl, xtp, vtp, ntot = npl + ntp 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 parallel do default(private) schedule(static) & !$omp shared(xpl, xtp, vpl, vtp, renc, boundingbox) & !$omp firstprivate(dt, npl, ntp, ntot) do dim = 1, SWEEPDIM @@ -577,9 +573,9 @@ subroutine encounter_check_all_sort_and_sweep_pltp(npl, ntp, xpl, vpl, xtp, vtp, xpl(dim,1:npl) + renc(1:npl) + vplshift_max(1:npl) * vpl(dim,1:npl) * dt, & xtp(dim,1:ntp) + vtpshift_max(1:ntp) * vtp(dim,1:ntp) * dt]) end do - !$omp end taskloop + !$omp end parallel do - call boundingbox%sweep(npl, ntp, ind_arr, nenc, index1, index2) + call boundingbox%sweep(npl, ntp, nenc, index1, index2) if (nenc > 0) then ! Shift test particle indices back into the proper range @@ -626,7 +622,7 @@ subroutine encounter_check_all_triangular_plpl(npl, x, v, renc, dt, lvdotr, inde ind_arr(:) = [(i, i = 1, npl)] !$omp parallel do default(private) schedule(static)& !$omp shared(npl, x, v, renc, dt, lenc, ind_arr) & - !$omp lastprivate(xr, yr, zr, vxr, vyr, vzr, renc12, lencounteri, lvdotri) + !$omp lastprivate(xr, yr, zr, vxr, vyr, vzr, renc12, nenci, lencounteri, lvdotri) do i = 1, npl do concurrent(j = i+1:npl) xr = x(1, j) - x(1, i) @@ -701,7 +697,7 @@ subroutine encounter_check_all_triangular_plplm(nplm, nplt, xplm, vplm, xplt, vp 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) + !$omp lastprivate(xr, yr, zr, vxr, vyr, vzr, renc12, nenci, lencounteri, lvdotri) do i = 1, nplm do concurrent(j = 1:nplt) xr = xplt(1, j) - xplm(1, i) @@ -775,7 +771,7 @@ subroutine encounter_check_all_triangular_pltp(npl, ntp, xpl, vpl, xtp, vtp, ren ind_arr(:) = [(i, i = 1, ntp)] !$omp parallel do default(private) schedule(static)& !$omp shared(npl, ntp, xpl, vpl, xtp, vtp, renc, dt, lenc, ind_arr) & - !$omp lastprivate(xr, yr, zr, vxr, vyr, vzr, lencounteri, lvdotri) + !$omp lastprivate(xr, yr, zr, vxr, vyr, vzr, nenci, lencounteri, lvdotri) do i = 1, npl do concurrent(j = 1:ntp) xr = xtp(1, j) - xpl(1, i) @@ -910,35 +906,26 @@ module subroutine encounter_check_sort_aabb_1D(self, n, extent_arr) integer(I4B), intent(in) :: n !! Number of bodies with extents real(DP), dimension(:), intent(in) :: extent_arr !! Array of extents of size 2*n ! Internals - integer(I4B) :: i, j, ibox, jbox + integer(I4B) :: i, j, k, ibox, jbox + integer(I4B), dimension(2) :: extent logical, dimension(:), allocatable :: lfresh call util_sort(extent_arr, self%ind) - allocate(lfresh(n)) - - ! Determine the interval starting points and sizes - lfresh(:) = .true. ! This will prevent double-counting of pairs - do ibox = 1, 2*n - i = self%ind(ibox) - if (i > n) i = i - n ! If this is an endpoint index, shift it to the correct range - if (.not.lfresh(i)) cycle - do jbox = ibox + 1, 2*n - j = self%ind(jbox) - if (j > n) j = j - n ! If this is an endpoint index, shift it to the correct range - if (j == i) then - lfresh(i) = .false. - self%ibeg(i) = ibox - self%iend(i) = jbox - exit ! We've reached the end of this interval - end if - end do + + do concurrent(k = 1:2*n) + i = self%ind(k) + if (i <= n) then + self%ibeg(i) = k + else + self%iend(i - n) = k + end if end do return end subroutine encounter_check_sort_aabb_1D - module subroutine encounter_check_sweep_aabb_double_list(self, n1, n2, ind_arr, nenc, index1, index2) + module subroutine encounter_check_sweep_aabb_double_list(self, n1, n2, nenc, index1, index2) !! author: David A. Minton !! !! Sweeps the sorted bounding box extents and returns the encounter candidates @@ -947,24 +934,28 @@ module subroutine encounter_check_sweep_aabb_double_list(self, n1, n2, ind_arr, class(encounter_bounding_box), intent(inout) :: self !! Multi-dimensional bounding box structure integer(I4B), intent(in) :: n1 !! Number of bodies 1 integer(I4B), intent(in) :: n2 !! Number of bodies 2 - integer(I4B), dimension(:), intent(in) :: ind_arr !! index array for mapping the body indices (body 2 indexes should be shifted by +n1 in this list) integer(I4B), intent(out) :: nenc !! Total number of encounter candidates integer(I4B), dimension(:), allocatable, intent(out) :: index1 !! List of indices for body 1 in each encounter candidate pair integer(I4B), dimension(:), allocatable, intent(out) :: index2 !! List of indices for body 2 in each encounter candidate pair !Internals Integer(I4B) :: i, k, ntot type(encounter_list), dimension(n1+n2) :: lenc !! Array of encounter lists (one encounter list per body) + type(walltimer) :: timer ntot = n1 + n2 ! 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 + ! call timer%reset() + ! call timer%start() !$omp parallel do default(private) schedule(static)& - !$omp shared(self, lenc, ind_arr) & + !$omp shared(self, lenc) & !$omp firstprivate(ntot, n1, n2) do i = 1, ntot - call encounter_check_sweep_aabb_one_double_list(i, n1, n2, self%aabb(1)%ind(:), self%aabb(1)%ibeg(:), self%aabb(1)%iend(:), self%aabb(2)%ibeg(:), self%aabb(2)%iend(:), ind_arr, lenc(i)) + call encounter_check_sweep_aabb_one_double_list(i, n1, n2, self%aabb(1)%ind(:), self%aabb(1)%ibeg(:), self%aabb(1)%iend(:), self%aabb(2)%ibeg(:), self%aabb(2)%iend(:), self%ind_arr(:), lenc(i)) end do !$omp end parallel do + ! call timer%stop() + ! call timer%report(nsubsteps=nthreads, message="AABB sweep double:") call encounter_check_collapse_ragged_list(lenc, ntot, nenc, index1, index2) @@ -979,7 +970,7 @@ module subroutine encounter_check_sweep_aabb_double_list(self, n1, n2, ind_arr, end subroutine encounter_check_sweep_aabb_double_list - module subroutine encounter_check_sweep_aabb_single_list(self, n, ind_arr, nenc, index1, index2) + module subroutine encounter_check_sweep_aabb_single_list(self, n, nenc, index1, index2) !! author: David A. Minton !! !! Sweeps the sorted bounding box extents and returns the encounter candidates. Mutual encounters @@ -988,23 +979,27 @@ module subroutine encounter_check_sweep_aabb_single_list(self, n, ind_arr, nenc, ! Arguments class(encounter_bounding_box), intent(inout) :: self !! Multi-dimensional bounding box structure integer(I4B), intent(in) :: n !! Number of bodies 1 - integer(I4B), dimension(:), intent(in) :: ind_arr !! index array for mapping body 2 indexes integer(I4B), intent(out) :: nenc !! Total number of encounter candidates integer(I4B), dimension(:), allocatable, intent(out) :: index1 !! List of indices for body 1 in each encounter candidate pair integer(I4B), dimension(:), allocatable, intent(out) :: index2 !! List of indices for body 2 in each encounter candidate pair !Internals Integer(I4B) :: i, k type(encounter_list), dimension(n) :: lenc !! Array of encounter lists (one encounter list per body) + type(walltimer) :: timer ! 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 + ! call timer%reset() + ! call timer%start() !$omp parallel do default(private) schedule(static)& - !$omp shared(self, lenc, ind_arr) & + !$omp shared(self, lenc) & !$omp firstprivate(n) do i = 1, n - call encounter_check_sweep_aabb_one_single_list(i, n, self%aabb(1)%ind(:), self%aabb(1)%ibeg(:), self%aabb(1)%iend(:), self%aabb(2)%ibeg(:), self%aabb(2)%iend(:), ind_arr, lenc(i)) + call encounter_check_sweep_aabb_one_single_list(i, n, self%aabb(1)%ind(:), self%aabb(1)%ibeg(:), self%aabb(1)%iend(:), self%aabb(2)%ibeg(:), self%aabb(2)%iend(:), self%ind_arr(:), lenc(i)) end do !$omp end parallel do + ! call timer%stop() + ! call timer%report(nsubsteps=nthreads, message="AABB sweep single:") call encounter_check_collapse_ragged_list(lenc, n, nenc, index1, index2) @@ -1019,7 +1014,7 @@ module subroutine encounter_check_sweep_aabb_single_list(self, n, ind_arr, nenc, end subroutine encounter_check_sweep_aabb_single_list - subroutine encounter_check_sweep_aabb_one_double_list(i, n1, n2, ext_ind, ibegx, iendx, ibegy, iendy, ind_arr2, lenc) + subroutine encounter_check_sweep_aabb_one_double_list(i, n1, n2, ext_ind, ibegx, iendx, ibegy, iendy, ind_arr, lenc) !! author: David A. Minton !! !! Performs a sweep operation on a single body. Encounters from the same lists not allowed (e.g. pl-tp encounters only) @@ -1031,7 +1026,7 @@ subroutine encounter_check_sweep_aabb_one_double_list(i, n1, n2, ext_ind, ibegx, integer(I4B), dimension(:), intent(in) :: ext_ind !! Sorted index array of extents integer(I4B), dimension(:), intent(in) :: ibegx, iendx !! Beginning and ending index lists in the x-dimension integer(I4B), dimension(:), intent(in) :: ibegy, iendy !! Beginning and ending index lists in the y-dimension - integer(I4B), dimension(:), intent(in) :: ind_arr2 !! index array for mapping body 2 indexes + integer(I4B), dimension(:), intent(in) :: ind_arr !! index array for mapping body 2 indexes type(encounter_list), intent(inout) :: lenc !! Encounter list for the ith body ! Internals integer(I4B) :: ibox, jbox, nbox, j, ybegi, yendi, ntot @@ -1054,7 +1049,7 @@ subroutine encounter_check_sweep_aabb_one_double_list(i, n1, n2, ext_ind, ibegx, lenc%nenc = count(lencounteri(:)) if (lenc%nenc > 0) then allocate(lenc%index2(lenc%nenc)) - lenc%index2(:) = pack(ind_arr2(:), lencounteri(:)) + lenc%index2(:) = pack(ind_arr(:), lencounteri(:)) end if return diff --git a/src/encounter/encounter_setup.f90 b/src/encounter/encounter_setup.f90 index 26a7e70c9..9d747c565 100644 --- a/src/encounter/encounter_setup.f90 +++ b/src/encounter/encounter_setup.f90 @@ -19,6 +19,10 @@ module subroutine encounter_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 + allocate(tmp(n)) + if (n_last > 0) tmp(1:n_last) = self%ind_arr(1:n_last) + call move_alloc(tmp, self%ind_arr) + self%ind_arr(n_last+1:n) = [(k, k = n_last+1, n)] do dim = 1, SWEEPDIM allocate(tmp(next)) if (n_last > 0) tmp(1:next_last) = self%aabb(dim)%ind(1:next_last) @@ -26,6 +30,9 @@ module subroutine encounter_setup_aabb(self, n, n_last) self%aabb(dim)%ind(next_last+1:next) = [(k, k = next_last+1, next)] end do else ! The number of bodies has gone down. Resize and chop of the old indices + allocate(tmp(n)) + tmp(1:n) = self%ind_arr(1:n) + call move_alloc(tmp, self%ind_arr) do dim = 1, SWEEPDIM allocate(tmp(next)) tmp(1:next) = pack(self%aabb(dim)%ind(1:next_last), self%aabb(dim)%ind(1:next_last) <= next) diff --git a/src/main/swiftest_driver.f90 b/src/main/swiftest_driver.f90 index 82b874bac..284706d10 100644 --- a/src/main/swiftest_driver.f90 +++ b/src/main/swiftest_driver.f90 @@ -104,7 +104,7 @@ program swiftest_driver write(*, statusfmt) param%t, tfrac, pl%nbody, nbody_system%tp%nbody end select if (param%lenergy) call nbody_system%conservation_report(param, lterminal=.true.) - call integration_timer%report(nsubsteps=istep_dump, message="Integration steps:", param=param) + call integration_timer%report(nsubsteps=istep_dump, message="Integration steps:") call integration_timer%reset() iout = istep_out diff --git a/src/modules/encounter_classes.f90 b/src/modules/encounter_classes.f90 index 81a1c0e12..f87e35553 100644 --- a/src/modules/encounter_classes.f90 +++ b/src/modules/encounter_classes.f90 @@ -42,6 +42,7 @@ module encounter_classes type encounter_bounding_box type(encounter_bounding_box_1D), dimension(SWEEPDIM) :: aabb + integer(I4B), dimension(:), allocatable :: ind_arr !! Index array of bodies contains procedure :: setup => encounter_setup_aabb !! Setup a new axis-aligned bounding box structure procedure :: sweep_single => encounter_check_sweep_aabb_single_list !! Sweeps the sorted bounding box extents and returns the encounter candidates @@ -144,22 +145,20 @@ module subroutine encounter_check_sort_aabb_1D(self, n, extent_arr) real(DP), dimension(:), intent(in) :: extent_arr !! Array of extents of size 2*n end subroutine encounter_check_sort_aabb_1D - module subroutine encounter_check_sweep_aabb_double_list(self, n1, n2, ind_arr, nenc, index1, index2) + module subroutine encounter_check_sweep_aabb_double_list(self, n1, n2, nenc, index1, index2) implicit none class(encounter_bounding_box), intent(inout) :: self !! Multi-dimensional bounding box structure integer(I4B), intent(in) :: n1 !! Number of bodies 1 integer(I4B), intent(in) :: n2 !! Number of bodies 2 - integer(I4B), dimension(:), intent(in) :: ind_arr !! index array for mapping the body indices (body 2 indexes should be shifted by +n1 in this list) integer(I4B), intent(out) :: nenc !! Total number of encounter candidates integer(I4B), dimension(:), allocatable, intent(out) :: index1 !! List of indices for body 1 in each encounter candidate pair integer(I4B), dimension(:), allocatable, intent(out) :: index2 !! List of indices for body 2 in each encounter candidate pair end subroutine encounter_check_sweep_aabb_double_list - module subroutine encounter_check_sweep_aabb_single_list(self, n, ind_arr, nenc, index1, index2) + module subroutine encounter_check_sweep_aabb_single_list(self, n, nenc, index1, index2) implicit none class(encounter_bounding_box), intent(inout) :: self !! Multi-dimensional bounding box structure integer(I4B), intent(in) :: n !! Number of bodies 1 - integer(I4B), dimension(:), intent(in) :: ind_arr !! index array for mapping body 2 indexes integer(I4B), intent(out) :: nenc !! Total number of encounter candidates integer(I4B), dimension(:), allocatable, intent(out) :: index1 !! List of indices for body 1 in each encounter candidate pair integer(I4B), dimension(:), allocatable, intent(out) :: index2 !! List of indices for body 2 in each encounter candidate pair diff --git a/src/modules/walltime_classes.f90 b/src/modules/walltime_classes.f90 index 7f1fcabc2..7e6c971b4 100644 --- a/src/modules/walltime_classes.f90 +++ b/src/modules/walltime_classes.f90 @@ -51,13 +51,11 @@ module walltime_classes end type interaction_timer interface - module subroutine walltime_report(self, nsubsteps, message, param) - use swiftest_classes, only : swiftest_parameters + module subroutine walltime_report(self, nsubsteps, message) implicit none class(walltimer), intent(inout) :: self !! Walltimer object integer(I4B), intent(in) :: nsubsteps !! Number of substeps used to compute the time per step character(len=*), intent(in) :: message !! Message to prepend to the wall time terminal output - class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters end subroutine walltime_report module subroutine walltime_reset(self) diff --git a/src/walltime/walltime.f90 b/src/walltime/walltime.f90 index 118ecd09b..a8d82c2b5 100644 --- a/src/walltime/walltime.f90 +++ b/src/walltime/walltime.f90 @@ -29,7 +29,7 @@ module subroutine walltime_stop(self) end subroutine walltime_stop - module subroutine walltime_report(self, nsubsteps, message, param) + module subroutine walltime_report(self, nsubsteps, message) !! author: David A. Minton !! !! Prints the elapsed time information to the terminal @@ -38,7 +38,6 @@ module subroutine walltime_report(self, nsubsteps, message, param) class(walltimer), intent(inout) :: self !! Walltimer object integer(I4B), intent(in) :: nsubsteps !! Number of substeps used to compute the time per step character(len=*), intent(in) :: message !! Message to prepend to the wall time terminal output - class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters ! Internals character(len=*), parameter :: walltimefmt = '" Total wall time: ", es12.5, "; Interval wall time: ", es12.5, "; Interval wall time/step: ", es12.5' character(len=STRMAX) :: fmt