From a2e1f86804d9a3f14d0818d99be6b8c302de3dfa Mon Sep 17 00:00:00 2001 From: David A Minton Date: Mon, 4 Oct 2021 10:12:52 -0400 Subject: [PATCH] Restructured code to promote encounter subroutines to its own module. Refactored encounter list type --- Makefile | 1 + .../{encounter.f90 => encounter_check.f90} | 335 +++++++++++------- src/encounter/encounter_io.f90 | 90 +++++ src/encounter/encounter_setup.f90 | 61 ++++ src/encounter/encounter_util.f90 | 141 ++++++++ src/io/io.f90 | 86 ----- src/modules/encounter_classes.f90 | 134 +++++++ src/modules/swiftest.f90 | 1 + src/modules/swiftest_classes.f90 | 117 ------ src/modules/symba_classes.f90 | 55 +-- src/rmvs/rmvs_io.f90 | 2 +- src/setup/setup.f90 | 55 --- src/symba/symba_setup.f90 | 6 +- src/symba/symba_util.f90 | 24 +- src/util/util_append.f90 | 32 -- src/util/util_copy.f90 | 28 -- src/util/util_resize.f90 | 38 -- src/util/util_spill.f90 | 38 -- 18 files changed, 672 insertions(+), 572 deletions(-) rename src/encounter/{encounter.f90 => encounter_check.f90} (75%) create mode 100644 src/encounter/encounter_io.f90 create mode 100644 src/encounter/encounter_setup.f90 create mode 100644 src/encounter/encounter_util.f90 create mode 100644 src/modules/encounter_classes.f90 diff --git a/Makefile b/Makefile index 463b5883f..46d620dae 100644 --- a/Makefile +++ b/Makefile @@ -48,6 +48,7 @@ SWIFTEST_MODULES = swiftest_globals.f90 \ swiftest_operators.f90 \ lambda_function.f90\ swiftest_classes.f90 \ + encounter_classes.f90 \ fraggle_classes.f90 \ whm_classes.f90 \ rmvs_classes.f90 \ diff --git a/src/encounter/encounter.f90 b/src/encounter/encounter_check.f90 similarity index 75% rename from src/encounter/encounter.f90 rename to src/encounter/encounter_check.f90 index 89d72d587..2b724731b 100644 --- a/src/encounter/encounter.f90 +++ b/src/encounter/encounter_check.f90 @@ -1,4 +1,4 @@ -submodule (swiftest_classes) s_encounter +submodule (encounter_classes) s_encounter_check use swiftest contains @@ -23,22 +23,25 @@ module subroutine encounter_check_all_plpl(param, npl, nplm, x, v, renc, dt, lvd ! Internals type(interaction_timer), save :: itimer logical, save :: lfirst = .true. - logical, save :: lsecond = .true. + logical, save :: lsecond = .false. integer(I8B) :: nplplm = 0_I8B if (param%ladaptive_encounters) then nplplm = nplm * npl - nplm * (nplm + 1) / 2 if (nplplm > 0) then - if (lfirst) then ! Skip the first pass to allow the bounding box extents to be sorted from their initial state at least once - lfirst = .false. - param%lencounter_sas = .true. - else if (lsecond) then + if (lfirst) then write(itimer%loopname, *) "encounter_check_all_plpl" write(itimer%looptype, *) "ENCOUNTER" - call itimer%time_this_loop(param, nplplm) - lsecond = .false. + lfirst = .false. + lsecond = .true. else - if (itimer%check(param, nplplm)) call itimer%time_this_loop(param, nplplm) + 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, nplplm) + lsecond = .false. + else if (itimer%check(param, nplplm)) then + lsecond = .true. + itimer%is_on = .false. + end if end if else param%lencounter_sas = .false. @@ -128,8 +131,7 @@ subroutine encounter_check_all_sort_and_sweep_plpl(npl, nplm, x, v, renc, dt, lv type(boundingBox), dimension(SWEEPDIM), save :: aabb logical, dimension(:), allocatable :: lenc_final, lvdotr_final integer(I2B), dimension(npl) :: vshift_min, vshift_max - integer :: ybegi, yendi - logical :: lency + integer(I4B) :: ybegi, yendi if (npl <= 1) return @@ -141,6 +143,9 @@ subroutine encounter_check_all_sort_and_sweep_plpl(npl, nplm, x, v, renc, dt, lv if (allocated(ind_arr)) deallocate(ind_arr) allocate(ind_arr(npl)) ind_arr(:) = [(i, i = 1, npl)] + + ! call aabb(dim)%setup(n, n_last, npl) + do dim = 1, SWEEPDIM if (allocated(aabb(dim)%ibeg)) deallocate(aabb(dim)%ibeg) allocate(aabb(dim)%ibeg(npl)) @@ -322,7 +327,7 @@ 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_pltp(npl, ntp, xpl, vpl, xtp, vtp, renc, dt, lvdotr, index1, index2, nenc) + subroutine encounter_check_all_sort_and_sweep_pltp(npl, ntp, xpl, vpl, xtp, vtp, rencpl, renctp, dt, lvdotr, index1, index2, nenc) !! author: David A. Minton !! !! Check for encounters between massive bodies and test particles. @@ -337,7 +342,8 @@ subroutine encounter_check_all_sort_and_sweep_pltp(npl, ntp, xpl, vpl, xtp, vtp, 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) :: renc !! Critical radii of massive bodies that defines an encounter + real(DP), dimension(:), intent(in) :: rencpl !! Critical radii of massive bodies that defines an encounter + real(DP), dimension(:), intent(in) :: renctp !! Critical radii of test particles 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 @@ -363,137 +369,196 @@ subroutine encounter_check_all_sort_and_sweep_pltp(npl, ntp, xpl, vpl, xtp, vtp, integer(I4B), save :: ntp_last = 0 type boundingBox integer(I4B), dimension(:), allocatable :: ind !! Sorted minimum/maximum extent indices + integer(I4B), dimension(:), allocatable :: ibeg !! Beginning index for box + integer(I4B), dimension(:), allocatable :: iend !! Ending index for box end type type(boundingBox), dimension(SWEEPDIM), save :: aabb logical, dimension(:), allocatable :: lenc_final, lvdotr_final integer(I2B), dimension(npl) :: vplshift_min, vplshift_max integer(I2B), dimension(ntp) :: vtpshift_min, vtpshift_max + integer(I4B) :: ybegi, yendi ! If this is the first time through, build the index lists if ((ntp == 0) .or. (npl == 0)) return - ntot = npl + ntp - n = 2 * ntot - if ((ntp_last /= ntp) .or. (npl_last /= npl)) then - if (ntp_last /= ntp) then - if (allocated(tpind_arr)) deallocate(tpind_arr) - allocate(tpind_arr(ntp)) - tpind_arr(1:ntp) = [(i, i = 1, ntp)] - end if - - if (n > n_last) then ! The number of bodies has grown. Resize and append the new bodies - do dim = 1, SWEEPDIM - allocate(tmp(n)) - if (npl_last > 0) tmp(1:n_last) = aabb(dim)%ind(1:n_last) - call move_alloc(tmp, aabb(dim)%ind) - aabb(dim)%ind(n_last+1:n) = [(k, k = n_last+1, n)] - end do - else if (n < n_last) then ! The number of bodies has gone down. Resize and chop of the old indices - do dim = 1, SWEEPDIM - allocate(tmp(n)) - tmp(1:n) = pack(aabb(dim)%ind(1:n_last), aabb(dim)%ind(1:n_last) <= n) - call move_alloc(tmp, aabb(dim)%ind) - end do - end if - - npl_last = npl - ntp_last = ntp - n_last = n - end if - - do i = 1, SWEEPDIM - where(vpl(i,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(i,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 util_sort([xpl(i,1:npl)-renc(1:npl)+vplshift_min(1:npl)*vpl(i,1:npl)*dt, & - xtp(i,1:ntp) +vtpshift_min(1:ntp)*vtp(i,1:ntp)*dt, & - xpl(i,1:npl)+renc(1:npl)+vplshift_max(1:npl)*vpl(i,1:npl)*dt, & - xtp(i,1:ntp) +vtpshift_max(1:ntp)*vtp(i,1:ntp)*dt], aabb(i)%ind) - end do - - ! Sweep the intervals - dim = 1 - lfresh(:) = .true. ! This will allow us to skip end-points of processed massive bodies - do ibox = 1, n - i = aabb(dim)%ind(ibox) - if (i > ntot) i = i - ntot ! If this is an endpoint index, shift it to the correct range - if (i > npl) cycle ! This is a test particle, so move on - if (.not.lfresh(i)) cycle ! This body has already been evaluated, so move on - lencounteri(:) = .false. - do jbox = ibox + 1, n ! Sweep forward until the end of the interval - j = aabb(dim)%ind(jbox) - if (j > ntot) j = j - ntot ! If this is an endpoint index, shift it to the correct range - if (j == i) exit ! We've reached the end of this interval - if (j > npl) then ! this is an unprocessed test particle - j = j - npl ! Shift into the range of the test particles - lencounteri(j) = .true. ! An overlapping tp/pl is found that has not previously been tallied - end if - end do - lfresh(i) = .false. ! This body has now been processed, so it should no longer show up in future encounter lists - nenci = count(lencounteri(:)) - if (nenci > 0) then - allocate(lenc(i)%index2(nenci)) - lenc(i)%nenc = nenci - lenc(i)%index2(:) = pack(tpind_arr(:), lencounteri(:)) - end if - end do + ! ntot = npl + ntp + ! n = 2 * ntot + ! if ((ntp_last /= ntp) .or. (npl_last /= npl)) then + ! if (ntp_last /= ntp) then + ! if (allocated(tpind_arr)) deallocate(tpind_arr) + ! allocate(tpind_arr(ntp)) + ! tpind_arr(1:ntp) = [(i, i = 1, ntp)] + ! end if + + ! do dim = 1, SWEEPDIM + ! if (allocated(aabb(dim)%ibeg)) deallocate(aabb(dim)%ibeg) + ! allocate(aabb(dim)%ibeg(ntot)) + ! if (allocated(aabb(dim)%iend)) deallocate(aabb(dim)%iend) + ! allocate(aabb(dim)%iend(ntot)) + ! end do + + ! do dim = 1, SWEEPDIM + ! if (npl > npl_last) then ! The number of bodies has grown. Resize and append the new bodies + ! allocate(tmp(n)) + ! if (npl_last > 0) tmp(1:n_last) = aabb(dim)%ind(1:n_last) + ! call move_alloc(tmp, aabb(dim)%ind) + ! aabb(dim)%ind(n_last+1:n) = [(k, k = n_last+1, n)] + ! else ! The number of bodies has gone down. Resize and chop of the old indices + ! allocate(tmp(n)) + ! tmp(1:n) = pack(aabb(dim)%ind(1:n_last), aabb(dim)%ind(1:n_last) <= n) + ! call move_alloc(tmp, aabb(dim)%ind) + ! end if + ! end do + + ! npl_last = npl + ! ntp_last = ntp + ! n_last = n + ! end if + + ! !$omp taskloop default(private) num_tasks(SWEEPDIM) & + ! !$omp shared(xpl, xtp, vpl, vtp, rencpl, renctp, aabb) & + ! !$omp firstprivate(dt, npl, n) + ! 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(i,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 util_sort([xpl(dim,1:npl)-rencpl(1:npl)+vplshift_min(1:npl)*vpl(i,1:npl)*dt, & + ! xtp(dim,1:ntp)-renctp(1:ntp)+vtpshift_min(1:ntp)*vtp(i,1:ntp)*dt, & + ! xpl(dim,1:npl)+rencpl(1:npl)+vplshift_max(1:npl)*vpl(i,1:npl)*dt, & + ! xtp(dim,1:ntp)+renctp(1:ntp)+vtpshift_max(1:ntp)*vtp(i,1:ntp)*dt], aabb(i)%ind) + + + ! ! Determine the interval starting points and sizes + ! lfresh(:) = .true. ! This will prevent double-counting of pairs + ! do ibox = 1, n + ! i = aabb(dim)%ind(ibox) + ! if (i > npl) i = i - npl ! If this is an endpoint index, shift it to the correct range + ! if (.not.lfresh(i)) cycle + ! do jbox = ibox + 1, n + ! j = aabb(dim)%ind(jbox) + ! if (j > npl) j = j - npl ! If this is an endpoint index, shift it to the correct range + ! if (j == i) then + ! lfresh(i) = .false. + ! aabb(dim)%ibeg(i) = ibox + ! aabb(dim)%iend(i) = jbox + ! exit ! We've reached the end of this interval + ! end if + ! end do + ! end do + ! end do + ! !$omp end taskloop + + ! !$omp parallel do default(private) schedule(static)& + ! !$omp shared(aabb, lenc, ind_arr) & + ! !$omp firstprivate(npl, ntp, ntot) + ! do i = 1, ntot + ! ibox = aabb(1)%ibeg(i) + 1 + ! nbox = aabb(1)%iend(i) - 1 + ! ybegi = aabb(2)%ibeg(i) + ! yendi = aabb(2)%iend(i) + ! lencounteri(:) = .false. + ! do concurrent(jbox = ibox:nbox) ! Sweep forward until the end of the interval + ! j = aabb(1)%ind(jbox) + ! if (j > ntot) j = j - ntot ! If this is an endpoint index, shift it to the correct range + ! if (((i <= npl) .and. (j <= npl)) .or. ((i > npl) .and. (j > npl))) cycle + ! ! Check the y-dimension + ! lencounteri(j) = (aabb(2)%iend(j) > ybegi) .and. (aabb(2)%ibeg(j) < yendi) + ! end do + + ! lenc(i)%nenc = count(lencounteri(:)) + ! if (lenc(i)%nenc > 0) then + ! allocate(lenc(i)%index2(lenc(i)%nenc)) + ! lenc(i)%index2(:) = pack(ind_arr(:), lencounteri(:)) + ! end if + ! end do + ! !$omp end parallel do + + ! ! Sweep the intervals + ! lfresh(:) = .true. ! This will allow us to skip end-points of processed massive bodies + ! do i = 1, npl + ntp + ! !i = aabb(1)%ind(ibox) + ! !ibox = aabb(1)%ibeg(i) + 1 + ! nbox = aabb(1)%iend(i) - 1 + ! ybegi = aabb(2)%ibeg(i) + ! yendi = aabb(2)%iend(i) + ! if (i > ntot) i = i - ntot ! If this is an endpoint index, shift it to the correct range + ! if (i > npl) cycle ! This is a test particle, so move on + ! if (.not.lfresh(i)) cycle ! This body has already been evaluated, so move on + ! lencounteri(:) = .false. + ! do jbox = ibox + 1, n ! Sweep forward until the end of the interval + ! j = aabb(1)%ind(jbox) + ! if (j > ntot) j = j - ntot ! If this is an endpoint index, shift it to the correct range + ! if (j == i) exit ! We've reached the end of this interval + ! if (j > npl) then ! this is an unprocessed test particle + ! j = j - npl ! Shift into the range of the test particles + ! lencounteri(j) = (aabb(2)%iend(j) > ybegi) .and. (aabb(2)%ibeg(j) < yendi) + ! end if + ! end do + ! lfresh(i) = .false. ! This body has now been processed, so it should no longer show up in future encounter lists + ! nenci = count(lencounteri(:)) + ! if (nenci > 0) then + ! allocate(lenc(i)%index2(nenci)) + ! lenc(i)%nenc = nenci + ! lenc(i)%index2(:) = pack(tpind_arr(:), lencounteri(:)) + ! end if + ! end do - associate(nenc_arr => lenc(:)%nenc) - nenc = sum(nenc_arr(1:npl)) - end associate - - if (nenc > 0) then - allocate(index1(nenc)) - allocate(index2(nenc)) - allocate(lenc_final(nenc)) - allocate(lvdotr_final(nenc)) - j0 = 1 - do i = 1, npl - if (lenc(i)%nenc > 0) then - j1 = j0 + lenc(i)%nenc - 1 - index1(j0:j1) = i - index2(j0:j1) = lenc(i)%index2(:) - j0 = j1 + 1 - end if - end do - ! Now that we have identified potential pairs, use the narrow-phase process to get the final values - lenc_final(:) = .true. - - do k = 1, nenc - i = index1(k) - j = index2(k) - xr = xtp(1, j) - xpl(1, i) - yr = xtp(2, j) - xpl(2, i) - zr = xtp(3, j) - xpl(3, i) - vxr = vtp(1, j) - vpl(1, i) - vyr = vtp(2, j) - vpl(2, i) - vzr = vtp(3, j) - vpl(3, i) - call encounter_check_one(xr, yr, zr, vxr, vyr, vzr, renc(i), dt, lenc_final(k), lvdotr_final(k)) - end do - - nenc = count(lenc_final(:)) ! Count the true number of encounters - allocate(tmp(nenc)) - tmp(:) = pack(index1(:), lenc_final(:)) - call move_alloc(tmp, index1) - allocate(tmp(nenc)) - tmp(:) = pack(index2(:), lenc_final(:)) - call move_alloc(tmp, index2) - allocate(lvdotr(nenc)) - lvdotr(:) = pack(lvdotr_final(:), lenc_final(:)) - end if + ! associate(nenc_arr => lenc(:)%nenc) + ! nenc = sum(nenc_arr(1:npl)) + ! end associate + + ! if (nenc > 0) then + ! allocate(index1(nenc)) + ! allocate(index2(nenc)) + ! allocate(lenc_final(nenc)) + ! allocate(lvdotr_final(nenc)) + ! j0 = 1 + ! do i = 1, npl + ! if (lenc(i)%nenc > 0) then + ! j1 = j0 + lenc(i)%nenc - 1 + ! index1(j0:j1) = i + ! index2(j0:j1) = lenc(i)%index2(:) + ! j0 = j1 + 1 + ! end if + ! end do + ! ! Now that we have identified potential pairs, use the narrow-phase process to get the final values + ! lenc_final(:) = .true. + + ! do k = 1, nenc + ! i = index1(k) + ! j = index2(k) + ! xr = xtp(1, j) - xpl(1, i) + ! yr = xtp(2, j) - xpl(2, i) + ! zr = xtp(3, j) - xpl(3, i) + ! vxr = vtp(1, j) - vpl(1, i) + ! vyr = vtp(2, j) - vpl(2, i) + ! vzr = vtp(3, j) - vpl(3, i) + ! call encounter_check_one(xr, yr, zr, vxr, vyr, vzr, renc(i), dt, lenc_final(k), lvdotr_final(k)) + ! end do + + ! nenc = count(lenc_final(:)) ! Count the true number of encounters + ! allocate(tmp(nenc)) + ! tmp(:) = pack(index1(:), lenc_final(:)) + ! call move_alloc(tmp, index1) + ! allocate(tmp(nenc)) + ! tmp(:) = pack(index2(:), lenc_final(:)) + ! call move_alloc(tmp, index2) + ! allocate(lvdotr(nenc)) + ! lvdotr(:) = pack(lvdotr_final(:), lenc_final(:)) + ! end if return end subroutine encounter_check_all_sort_and_sweep_pltp @@ -696,4 +761,4 @@ module pure subroutine encounter_check_one(xr, yr, zr, vxr, vyr, vzr, renc, dt, return end subroutine encounter_check_one -end submodule s_encounter +end submodule s_encounter_check diff --git a/src/encounter/encounter_io.f90 b/src/encounter/encounter_io.f90 new file mode 100644 index 000000000..60fdf5304 --- /dev/null +++ b/src/encounter/encounter_io.f90 @@ -0,0 +1,90 @@ +submodule (encounter_classes) s_encounter_io + use swiftest +contains + + + module subroutine encounter_io_write_frame(iu, t, id1, id2, Gmass1, Gmass2, radius1, radius2, xh1, xh2, vh1, vh2) + !! author: David A. Minton + !! + !! Write a single frame of close encounter data to output binary files + !! + !! Adapted from David E. Kaufmann's Swifter routine: io_write_encounter.f90 + !! Adapted from Hal Levison's Swift routine io_write_encounter.f + implicit none + ! Arguments + integer(I4B), intent(in) :: iu !! Open file unit number + real(DP), intent(in) :: t !! Time of encounter + integer(I4B), intent(in) :: id1, id2 !! ids of the two encountering bodies + real(DP), intent(in) :: Gmass1, Gmass2 !! G*mass of the two encountering bodies + real(DP), intent(in) :: radius1, radius2 !! Radii of the two encountering bodies + real(DP), dimension(:), intent(in) :: xh1, xh2 !! Heliocentric position vectors of the two encountering bodies + real(DP), dimension(:), intent(in) :: vh1, vh2 !! Heliocentric velocity vectors of the two encountering bodies + ! Internals + character(len=STRMAX) :: errmsg + + write(iu, err = 667, iomsg = errmsg) t + write(iu, err = 667, iomsg = errmsg) id1, xh1(1), xh1(2), xh1(3), vh1(1), vh1(2), Gmass1, radius1 + write(iu, err = 667, iomsg = errmsg) id2, xh2(1), xh2(2), xh2(3), vh2(1), vh2(2), Gmass2, radius2 + + return + 667 continue + write(*,*) "Error writing encounter file: " // trim(adjustl(errmsg)) + call util_exit(FAILURE) + end subroutine + + module subroutine encounter_io_write_list(self, pl, encbody, param) + implicit none + ! Arguments + class(encounter_list), intent(in) :: self !! Swiftest encounter list object + class(swiftest_pl), intent(in) :: pl !! Swiftest massive body object + class(swiftest_body), intent(in) :: encbody !! Encountering body - Swiftest generic body object (pl or tp) + class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters + ! Internals + logical , save :: lfirst = .true. + integer(I4B) :: k, ierr + character(len=STRMAX) :: errmsg + + if (param%enc_out == "" .or. self%nenc == 0) return + + open(unit = LUN, file = param%enc_out, status = 'OLD', position = 'APPEND', form = 'UNFORMATTED', iostat = ierr, iomsg = errmsg) + if (ierr /= 0) then + if (lfirst) then + open(unit = LUN, file = param%enc_out, status = 'NEW', form = 'UNFORMATTED', err = 667, iomsg = errmsg) + else + goto 667 + end if + end if + lfirst = .false. + + associate(ind1 => self%index1, ind2 => self%index2) + select type(encbody) + class is (swiftest_pl) + do k = 1, self%nenc + call encounter_io_write_frame(LUN, self%t(k), & + pl%id(ind1(k)), encbody%id(ind2(k)), & + pl%Gmass(ind1(k)), encbody%Gmass(ind2(k)), & + pl%radius(ind1(k)), encbody%radius(ind2(k)), & + self%x1(:,k), self%x2(:,k), & + self%v1(:,k), self%v2(:,k)) + end do + class is (swiftest_tp) + do k = 1, self%nenc + call encounter_io_write_frame(LUN, self%t(k), & + pl%id(ind1(k)), encbody%id(ind2(k)), & + pl%Gmass(ind1(k)), 0.0_DP, & + pl%radius(ind1(k)), 0.0_DP, & + self%x1(:,k), self%x2(:,k), & + self%v1(:,k), self%v2(:,k)) + end do + end select + end associate + + close(unit = LUN, err = 667, iomsg = errmsg) + + return + 667 continue + write(*,*) "Error writing encounter file: " // trim(adjustl(errmsg)) + call util_exit(FAILURE) + end subroutine encounter_io_write_list + +end submodule s_encounter_io \ No newline at end of file diff --git a/src/encounter/encounter_setup.f90 b/src/encounter/encounter_setup.f90 new file mode 100644 index 000000000..91d8d52d6 --- /dev/null +++ b/src/encounter/encounter_setup.f90 @@ -0,0 +1,61 @@ +submodule (encounter_classes) s_encounter_setup + use swiftest +contains + + module subroutine encounter_setup_list(self, n) + !! author: David A. Minton + !! + !! A constructor that sets the number of encounters and allocates and initializes all arrays + !! + implicit none + ! Arguments + class(encounter_list), intent(inout) :: self !! Swiftest encounter structure + integer(I4B), intent(in) :: n !! Number of encounters to allocate space for + + if (n < 0) return + + if (allocated(self%lvdotr)) deallocate(self%lvdotr) + if (allocated(self%status)) deallocate(self%status) + if (allocated(self%index1)) deallocate(self%index1) + if (allocated(self%index2)) deallocate(self%index2) + if (allocated(self%id1)) deallocate(self%id1) + if (allocated(self%id2)) deallocate(self%id2) + if (allocated(self%x1)) deallocate(self%x1) + if (allocated(self%x2)) deallocate(self%x2) + if (allocated(self%v1)) deallocate(self%v1) + if (allocated(self%v2)) deallocate(self%v2) + if (allocated(self%t)) deallocate(self%t) + + self%nenc = n + if (n == 0) return + + allocate(self%lvdotr(n)) + allocate(self%status(n)) + allocate(self%index1(n)) + allocate(self%index2(n)) + allocate(self%id1(n)) + allocate(self%id2(n)) + allocate(self%x1(NDIM,n)) + allocate(self%x2(NDIM,n)) + allocate(self%v1(NDIM,n)) + allocate(self%v2(NDIM,n)) + allocate(self%t(n)) + + self%lvdotr(:) = .false. + self%status(:) = INACTIVE + self%index1(:) = 0 + self%index2(:) = 0 + self%id1(:) = 0 + self%id2(:) = 0 + self%x1(:,:) = 0.0_DP + self%x2(:,:) = 0.0_DP + self%v1(:,:) = 0.0_DP + self%v2(:,:) = 0.0_DP + self%t(:) = 0.0_DP + + return + end subroutine encounter_setup_list + +end submodule s_encounter_setup + + \ No newline at end of file diff --git a/src/encounter/encounter_util.f90 b/src/encounter/encounter_util.f90 new file mode 100644 index 000000000..ae4a6d486 --- /dev/null +++ b/src/encounter/encounter_util.f90 @@ -0,0 +1,141 @@ +submodule (encounter_classes) s_encounter_util + use swiftest +contains + + module subroutine encounter_util_append_list(self, source, lsource_mask) + !! author: David A. Minton + !! + !! Append components from one Swiftest body object to another. + !! This method will automatically resize the destination body if it is too small + implicit none + ! Arguments + class(encounter_list), intent(inout) :: self !! Swiftest encounter list object + class(encounter_list), intent(in) :: source !! Source object to append + logical, dimension(:), intent(in) :: lsource_mask !! Logical mask indicating which elements to append to + ! Internals + integer(I4B) :: nold, nsrc + + nold = self%nenc + nsrc = source%nenc + call util_append(self%lvdotr, source%lvdotr, nold, nsrc, lsource_mask) + call util_append(self%status, source%status, nold, nsrc, lsource_mask) + call util_append(self%index1, source%index1, nold, nsrc, lsource_mask) + call util_append(self%index2, source%index2, nold, nsrc, lsource_mask) + call util_append(self%id1, source%id1, nold, nsrc, lsource_mask) + call util_append(self%id2, source%id2, nold, nsrc, lsource_mask) + call util_append(self%x1, source%x1, nold, nsrc, lsource_mask) + call util_append(self%x2, source%x2, nold, nsrc, lsource_mask) + call util_append(self%v1, source%v1, nold, nsrc, lsource_mask) + call util_append(self%v2, source%v2, nold, nsrc, lsource_mask) + call util_append(self%t, source%t, nold, nsrc, lsource_mask) + self%nenc = nold + count(lsource_mask(1:nsrc)) + + return + end subroutine encounter_util_append_list + + + module subroutine encounter_util_copy_list(self, source) + !! author: David A. Minton + !! + !! Copies elements from the source encounter list into self. + implicit none + ! Arguments + class(encounter_list), intent(inout) :: self !! Encounter list + class(encounter_list), intent(in) :: source !! Source object to copy into + + associate(n => source%nenc) + self%nenc = n + self%lvdotr(1:n) = source%lvdotr(1:n) + self%status(1:n) = source%status(1:n) + self%index1(1:n) = source%index1(1:n) + self%index2(1:n) = source%index2(1:n) + self%id1(1:n) = source%id1(1:n) + self%id2(1:n) = source%id2(1:n) + self%x1(:,1:n) = source%x1(:,1:n) + self%x2(:,1:n) = source%x2(:,1:n) + self%v1(:,1:n) = source%v1(:,1:n) + self%v2(:,1:n) = source%v2(:,1:n) + self%t(1:n) = source%t(1:n) + end associate + + return + end subroutine encounter_util_copy_list + + + module subroutine encounter_util_resize_list(self, nnew) + !! author: David A. Minton + !! + !! Checks the current size of the encounter list against the required size and extends it by a factor of 2 more than requested if it is too small. + !! Note: The reason to extend it by a factor of 2 is for performance. When there are many enounters per step, resizing every time you want to add an + !! encounter takes significant computational effort. Resizing by a factor of 2 is a tradeoff between performance (fewer resize calls) and memory managment + !! Memory usage grows by a factor of 2 each time it fills up, but no more. + implicit none + ! Arguments + class(encounter_list), intent(inout) :: self !! Swiftest encounter list + integer(I4B), intent(in) :: nnew !! New size of list needed + ! Internals + class(encounter_list), allocatable :: enc_temp + integer(I4B) :: nold + logical :: lmalloc + + lmalloc = allocated(self%status) + if (lmalloc) then + nold = size(self%status) + else + nold = 0 + end if + if (nnew > nold) then + if (lmalloc) allocate(enc_temp, source=self) + call self%setup(2 * nnew) + if (lmalloc) then + call self%copy(enc_temp) + deallocate(enc_temp) + end if + else + self%status(nnew+1:nold) = INACTIVE + end if + self%nenc = nnew + + return + end subroutine encounter_util_resize_list + + + module subroutine encounter_util_spill_list(self, discards, lspill_list, ldestructive) + !! author: David A. Minton + !! + !! Move spilled (discarded) Swiftest encounter structure from active list to discard list + implicit none + ! Arguments + class(encounter_list), intent(inout) :: self !! Swiftest encounter list + class(encounter_list), intent(inout) :: discards !! Discarded object + logical, dimension(:), intent(in) :: lspill_list !! Logical array of bodies to spill into the discards + logical, intent(in) :: ldestructive !! Logical flag indicating whether or not this operation should alter body by removing the discard list + ! Internals + integer(I4B) :: nenc_old + + associate(keeps => self) + call util_spill(keeps%lvdotr, discards%lvdotr, lspill_list, ldestructive) + call util_spill(keeps%status, discards%status, lspill_list, ldestructive) + call util_spill(keeps%index1, discards%index1, lspill_list, ldestructive) + call util_spill(keeps%index2, discards%index2, lspill_list, ldestructive) + call util_spill(keeps%id1, discards%id1, lspill_list, ldestructive) + call util_spill(keeps%id2, discards%id2, lspill_list, ldestructive) + call util_spill(keeps%x1, discards%x1, lspill_list, ldestructive) + call util_spill(keeps%x2, discards%x2, lspill_list, ldestructive) + call util_spill(keeps%v1, discards%v1, lspill_list, ldestructive) + call util_spill(keeps%v2, discards%v2, lspill_list, ldestructive) + call util_spill(keeps%t, discards%t, lspill_list, ldestructive) + + nenc_old = keeps%nenc + + ! This is the base class, so will be the last to be called in the cascade. + ! Therefore we need to set the nenc values for both the keeps and discareds + discards%nenc = count(lspill_list(1:nenc_old)) + if (ldestructive) keeps%nenc = nenc_old - discards%nenc + end associate + + return + end subroutine encounter_util_spill_list + + +end submodule s_encounter_util \ No newline at end of file diff --git a/src/io/io.f90 b/src/io/io.f90 index c0f21c4ec..93d86c2f0 100644 --- a/src/io/io.f90 +++ b/src/io/io.f90 @@ -1929,62 +1929,6 @@ module subroutine io_write_discard(self, param) end subroutine io_write_discard - module subroutine io_write_encounter(self, pl, encbody, param) - implicit none - ! Arguments - class(swiftest_encounter), intent(in) :: self !! Swiftest encounter list object - class(swiftest_pl), intent(in) :: pl !! Swiftest massive body object - class(swiftest_body), intent(in) :: encbody !! Encountering body - Swiftest generic body object (pl or tp) - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters - ! Internals - logical , save :: lfirst = .true. - integer(I4B) :: k, ierr - character(len=STRMAX) :: errmsg - - if (param%enc_out == "" .or. self%nenc == 0) return - - open(unit = LUN, file = param%enc_out, status = 'OLD', position = 'APPEND', form = 'UNFORMATTED', iostat = ierr, iomsg = errmsg) - if (ierr /= 0) then - if (lfirst) then - open(unit = LUN, file = param%enc_out, status = 'NEW', form = 'UNFORMATTED', err = 667, iomsg = errmsg) - else - goto 667 - end if - end if - lfirst = .false. - - associate(ind1 => self%index1, ind2 => self%index2) - select type(encbody) - class is (swiftest_pl) - do k = 1, self%nenc - call io_write_frame_encounter(LUN, self%t(k), & - pl%id(ind1(k)), encbody%id(ind2(k)), & - pl%Gmass(ind1(k)), encbody%Gmass(ind2(k)), & - pl%radius(ind1(k)), encbody%radius(ind2(k)), & - self%x1(:,k), self%x2(:,k), & - self%v1(:,k), self%v2(:,k)) - end do - class is (swiftest_tp) - do k = 1, self%nenc - call io_write_frame_encounter(LUN, self%t(k), & - pl%id(ind1(k)), encbody%id(ind2(k)), & - pl%Gmass(ind1(k)), 0.0_DP, & - pl%radius(ind1(k)), 0.0_DP, & - self%x1(:,k), self%x2(:,k), & - self%v1(:,k), self%v2(:,k)) - end do - end select - end associate - - close(unit = LUN, err = 667, iomsg = errmsg) - - return - 667 continue - write(*,*) "Error writing encounter file: " // trim(adjustl(errmsg)) - call util_exit(FAILURE) - end subroutine io_write_encounter - - module subroutine io_write_frame_body(self, iu, param) !! author: David A. Minton !! @@ -2091,36 +2035,6 @@ module subroutine io_write_frame_cb(self, iu, param) end subroutine io_write_frame_cb - module subroutine io_write_frame_encounter(iu, t, id1, id2, Gmass1, Gmass2, radius1, radius2, xh1, xh2, vh1, vh2) - !! author: David A. Minton - !! - !! Write a single frame of close encounter data to output binary files - !! - !! Adapted from David E. Kaufmann's Swifter routine: io_write_encounter.f90 - !! Adapted from Hal Levison's Swift routine io_write_encounter.f - implicit none - ! Arguments - integer(I4B), intent(in) :: iu !! Open file unit number - real(DP), intent(in) :: t !! Time of encounter - integer(I4B), intent(in) :: id1, id2 !! ids of the two encountering bodies - real(DP), intent(in) :: Gmass1, Gmass2 !! G*mass of the two encountering bodies - real(DP), intent(in) :: radius1, radius2 !! Radii of the two encountering bodies - real(DP), dimension(:), intent(in) :: xh1, xh2 !! Heliocentric position vectors of the two encountering bodies - real(DP), dimension(:), intent(in) :: vh1, vh2 !! Heliocentric velocity vectors of the two encountering bodies - ! Internals - character(len=STRMAX) :: errmsg - - write(iu, err = 667, iomsg = errmsg) t - write(iu, err = 667, iomsg = errmsg) id1, xh1(1), xh1(2), xh1(3), vh1(1), vh1(2), Gmass1, radius1 - write(iu, err = 667, iomsg = errmsg) id2, xh2(1), xh2(2), xh2(3), vh2(1), vh2(2), Gmass2, radius2 - - return - 667 continue - write(*,*) "Error writing encounter file: " // trim(adjustl(errmsg)) - call util_exit(FAILURE) - end subroutine - - module subroutine io_write_frame_system(self, param) !! author: The Purdue Swiftest Team - David A. Minton, Carlisle A. Wishard, Jennifer L.L. Pouplin, and Jacob R. Elliott !! diff --git a/src/modules/encounter_classes.f90 b/src/modules/encounter_classes.f90 new file mode 100644 index 000000000..a8cbdf440 --- /dev/null +++ b/src/modules/encounter_classes.f90 @@ -0,0 +1,134 @@ +module encounter_classes + !! author: The Purdue Swiftest Team - David A. Minton, Carlisle A. Wishard, Jennifer L.L. Pouplin, and Jacob R. Elliott + !! + !! Definition of classes and methods used to determine close encounters + use swiftest_globals + use swiftest_classes + implicit none + public + + type :: encounter_list + integer(I4B) :: nenc !! Total number of encounters + logical, dimension(:), allocatable :: lvdotr !! relative vdotr flag + integer(I4B), dimension(:), allocatable :: status !! status of the interaction + integer(I4B), dimension(:), allocatable :: index1 !! position of the first body in the encounter + integer(I4B), dimension(:), allocatable :: index2 !! position of the second body in the encounter + integer(I4B), dimension(:), allocatable :: id1 !! id of the first body in the encounter + integer(I4B), dimension(:), allocatable :: id2 !! id of the second body in the encounter + real(DP), dimension(:,:), allocatable :: x1 !! the position of body 1 in the encounter + real(DP), dimension(:,:), allocatable :: x2 !! the position of body 2 in the encounter + real(DP), dimension(:,:), allocatable :: v1 !! the velocity of body 1 in the encounter + real(DP), dimension(:,:), allocatable :: v2 !! the velocity of body 2 in the encounter + real(DP), dimension(:), allocatable :: t !! Time of encounter + contains + procedure :: setup => encounter_setup_list !! A constructor that sets the number of encounters and allocates and initializes all arrays + procedure :: append => encounter_util_append_list !! Appends elements from one structure to another + procedure :: copy => encounter_util_copy_list !! Copies elements from the source encounter list into self. + procedure :: spill => encounter_util_spill_list !! "Spills" bodies from one object to another depending on the results of a mask (uses the PACK intrinsic) + procedure :: resize => encounter_util_resize_list !! Checks the current size of the encounter list against the required size and extends it by a factor of 2 more than requested if it is too small. + procedure :: write => encounter_io_write_list !! Write close encounter data to output binary file + end type encounter_list + + + interface + module subroutine encounter_check_all_plpl(param, npl, nplm, 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 + 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_plpl + + module subroutine encounter_check_all_pltp(param, npl, ntp, xpl, vpl, xtp, vtp, 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) :: 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) :: 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 + end subroutine encounter_check_all_pltp + + module pure subroutine encounter_check_one(xr, yr, zr, vxr, vyr, vzr, renc, dt, lencounter, lvdotr) + !$omp declare simd(encounter_check_one) + implicit none + real(DP), intent(in) :: xr, yr, zr !! Relative distance vector components + real(DP), intent(in) :: vxr, vyr, vzr !! Relative velocity vector components + real(DP), intent(in) :: renc !! Critical encounter distance + real(DP), intent(in) :: dt !! Step size + logical, intent(out) :: lencounter !! Flag indicating that an encounter has occurred + logical, intent(out) :: lvdotr !! Logical flag indicating the direction of the v .dot. r vector + end subroutine encounter_check_one + + module subroutine encounter_io_write_frame(iu, t, id1, id2, Gmass1, Gmass2, radius1, radius2, xh1, xh2, vh1, vh2) + implicit none + integer(I4B), intent(in) :: iu !! Open file unit number + real(DP), intent(in) :: t !! Time of encounter + integer(I4B), intent(in) :: id1, id2 !! ids of the two encountering bodies + real(DP), intent(in) :: Gmass1, Gmass2 !! G*mass of the two encountering bodies + real(DP), intent(in) :: radius1, radius2 !! Radii of the two encountering bodies + real(DP), dimension(:), intent(in) :: xh1, xh2 !! Swiftestcentric position vectors of the two encountering bodies + real(DP), dimension(:), intent(in) :: vh1, vh2 !! Swiftestcentric velocity vectors of the two encountering bodies + end subroutine encounter_io_write_frame + + module subroutine encounter_io_write_list(self, pl, encbody, param) + implicit none + class(encounter_list), intent(in) :: self !! Swiftest encounter list object + class(swiftest_pl), intent(in) :: pl !! Swiftest massive body object + class(swiftest_body), intent(in) :: encbody !! Encountering body - Swiftest generic body object (pl or tp) + class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters + end subroutine encounter_io_write_list + + module subroutine encounter_setup_list(self, n) + implicit none + class(encounter_list), intent(inout) :: self !! Swiftest encounter structure + integer(I4B), intent(in) :: n !! Number of encounters to allocate space for + end subroutine encounter_setup_list + + module subroutine encounter_util_append_list(self, source, lsource_mask) + implicit none + class(encounter_list), intent(inout) :: self !! Swiftest encounter list object + class(encounter_list), intent(in) :: source !! Source object to append + logical, dimension(:), intent(in) :: lsource_mask !! Logical mask indicating which elements to append to + end subroutine encounter_util_append_list + + module subroutine encounter_util_copy_list(self, source) + implicit none + class(encounter_list), intent(inout) :: self !! Encounter list + class(encounter_list), intent(in) :: source !! Source object to copy into + end subroutine encounter_util_copy_list + + module subroutine encounter_util_resize_list(self, nnew) + implicit none + class(encounter_list), intent(inout) :: self !! Swiftest encounter list + integer(I4B), intent(in) :: nnew !! New size of list needed + end subroutine encounter_util_resize_list + + module subroutine encounter_util_spill_list(self, discards, lspill_list, ldestructive) + implicit none + class(encounter_list), intent(inout) :: self !! Swiftest encounter list + class(encounter_list), intent(inout) :: discards !! Discarded object + logical, dimension(:), intent(in) :: lspill_list !! Logical array of bodies to spill into the discards + logical, intent(in) :: ldestructive !! Logical flag indicating whether or not this operation should alter body by removing the discard list + end subroutine encounter_util_spill_list + + end interface + +end module encounter_classes + diff --git a/src/modules/swiftest.f90 b/src/modules/swiftest.f90 index d3f8996f0..3cab7a6fb 100644 --- a/src/modules/swiftest.f90 +++ b/src/modules/swiftest.f90 @@ -13,6 +13,7 @@ module swiftest use fraggle_classes use lambda_function use walltime_classes + use encounter_classes !use advisor_annotate !$ use omp_lib implicit none diff --git a/src/modules/swiftest_classes.f90 b/src/modules/swiftest_classes.f90 index 9e242b821..5c3f25f21 100644 --- a/src/modules/swiftest_classes.f90 +++ b/src/modules/swiftest_classes.f90 @@ -456,28 +456,6 @@ module swiftest_classes generic :: write_frame => write_frame_bin, write_frame_netcdf !! Generic method call for writing a frame of output data end type swiftest_nbody_system - type :: swiftest_encounter - integer(I4B) :: nenc !! Total number of encounters - logical, dimension(:), allocatable :: lvdotr !! relative vdotr flag - integer(I4B), dimension(:), allocatable :: status !! status of the interaction - integer(I4B), dimension(:), allocatable :: index1 !! position of the first body in the encounter - integer(I4B), dimension(:), allocatable :: index2 !! position of the second body in the encounter - integer(I4B), dimension(:), allocatable :: id1 !! id of the first body in the encounter - integer(I4B), dimension(:), allocatable :: id2 !! id of the second body in the encounter - real(DP), dimension(:,:), allocatable :: x1 !! the position of body 1 in the encounter - real(DP), dimension(:,:), allocatable :: x2 !! the position of body 2 in the encounter - real(DP), dimension(:,:), allocatable :: v1 !! the velocity of body 1 in the encounter - real(DP), dimension(:,:), allocatable :: v2 !! the velocity of body 2 in the encounter - real(DP), dimension(:), allocatable :: t !! Time of encounter - contains - procedure :: setup => setup_encounter !! A constructor that sets the number of encounters and allocates and initializes all arrays - procedure :: append => util_append_encounter !! Appends elements from one structure to another - procedure :: copy => util_copy_encounter !! Copies elements from the source encounter list into self. - procedure :: spill => util_spill_encounter !! "Spills" bodies from one object to another depending on the results of a mask (uses the PACK intrinsic) - procedure :: resize => util_resize_encounter !! Checks the current size of the encounter list against the required size and extends it by a factor of 2 more than requested if it is too small. - procedure :: write => io_write_encounter !! Write close encounter data to output binary file - end type swiftest_encounter - abstract interface subroutine abstract_discard_body(self, system, param) import swiftest_body, swiftest_nbody_system, swiftest_parameters @@ -589,49 +567,6 @@ module pure elemental subroutine drift_one(mu, px, py, pz, vx, vy, vz, dt, iflag integer(I4B), intent(out) :: iflag !! iflag : error status flag for Danby drift (0 = OK, nonzero = ERROR) end subroutine drift_one - module subroutine encounter_check_all_plpl(param, npl, nplm, x, v, renc, dt, lvdotr, index1, index2, nenc) - 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 - 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_plpl - - module subroutine encounter_check_all_pltp(param, npl, ntp, xpl, vpl, xtp, vtp, renc, dt, lvdotr, index1, index2, nenc) - 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) :: 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) :: 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 - end subroutine encounter_check_all_pltp - - module pure subroutine encounter_check_one(xr, yr, zr, vxr, vyr, vzr, renc, dt, lencounter, lvdotr) - !$omp declare simd(encounter_check_one) - implicit none - real(DP), intent(in) :: xr, yr, zr !! Relative distance vector components - real(DP), intent(in) :: vxr, vyr, vzr !! Relative velocity vector components - real(DP), intent(in) :: renc !! Critical encounter distance - real(DP), intent(in) :: dt !! Step size - logical, intent(out) :: lencounter !! Flag indicating that an encounter has occurred - logical, intent(out) :: lvdotr !! Logical flag indicating the direction of the v .dot. r vector - end subroutine encounter_check_one - module pure subroutine gr_kick_getaccb_ns_body(self, system, param) implicit none class(swiftest_body), intent(inout) :: self !! Swiftest generic body object @@ -909,14 +844,6 @@ module subroutine io_toupper(string) character(*), intent(inout) :: string !! String to make upper case end subroutine io_toupper - module subroutine io_write_encounter(self, pl, encbody, param) - implicit none - class(swiftest_encounter), intent(in) :: self !! Swiftest encounter list object - class(swiftest_pl), intent(in) :: pl !! Swiftest massive body object - class(swiftest_body), intent(in) :: encbody !! Encountering body - Swiftest generic body object (pl or tp) - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters - end subroutine io_write_encounter - module subroutine io_write_frame_body(self, iu, param) implicit none class(swiftest_body), intent(in) :: self !! Swiftest body object @@ -931,17 +858,6 @@ module subroutine io_write_frame_cb(self, iu, param) class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters end subroutine io_write_frame_cb - module subroutine io_write_frame_encounter(iu, t, id1, id2, Gmass1, Gmass2, radius1, radius2, xh1, xh2, vh1, vh2) - implicit none - integer(I4B), intent(in) :: iu !! Open file unit number - real(DP), intent(in) :: t !! Time of encounter - integer(I4B), intent(in) :: id1, id2 !! ids of the two encountering bodies - real(DP), intent(in) :: Gmass1, Gmass2 !! G*mass of the two encountering bodies - real(DP), intent(in) :: radius1, radius2 !! Radii of the two encountering bodies - real(DP), dimension(:), intent(in) :: xh1, xh2 !! Swiftestcentric position vectors of the two encountering bodies - real(DP), dimension(:), intent(in) :: vh1, vh2 !! Swiftestcentric velocity vectors of the two encountering bodies - end subroutine io_write_frame_encounter - module subroutine io_write_frame_system(self, param) implicit none class(swiftest_nbody_system), intent(inout) :: self !! Swiftest system object @@ -1186,12 +1102,6 @@ module subroutine setup_construct_system(system, param) class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters end subroutine setup_construct_system - module subroutine setup_encounter(self, n) - implicit none - class(swiftest_encounter), intent(inout) :: self !! Swiftest encounter structure - integer(I4B), intent(in) :: n !! Number of encounters to allocate space for - end subroutine setup_encounter - module subroutine setup_finalize_system(self, param) implicit none class(swiftest_nbody_system), intent(inout) :: self !! Swiftest system object @@ -1306,13 +1216,6 @@ module subroutine util_append_body(self, source, lsource_mask) logical, dimension(:), intent(in) :: lsource_mask !! Logical mask indicating which elements to append to end subroutine util_append_body - module subroutine util_append_encounter(self, source, lsource_mask) - implicit none - class(swiftest_encounter), intent(inout) :: self !! Swiftest encounter list object - class(swiftest_encounter), intent(in) :: source !! Source object to append - logical, dimension(:), intent(in) :: lsource_mask !! Logical mask indicating which elements to append to - end subroutine util_append_encounter - module subroutine util_append_pl(self, source, lsource_mask) implicit none class(swiftest_pl), intent(inout) :: self !! Swiftest massive body object @@ -1387,12 +1290,6 @@ module subroutine util_coord_xh2xb_tp(self, cb) class(swiftest_cb), intent(in) :: cb !! Swiftest central body object end subroutine util_coord_xh2xb_tp - module subroutine util_copy_encounter(self, source) - implicit none - class(swiftest_encounter), intent(inout) :: self !! Encounter list - class(swiftest_encounter), intent(in) :: source !! Source object to copy into - end subroutine util_copy_encounter - module subroutine util_copy_particle_info(self, source) implicit none class(swiftest_particle_info), intent(inout) :: self @@ -1583,12 +1480,6 @@ module subroutine util_resize_body(self, nnew) integer(I4B), intent(in) :: nnew !! New size neded end subroutine util_resize_body - module subroutine util_resize_encounter(self, nnew) - implicit none - class(swiftest_encounter), intent(inout) :: self !! Swiftest encounter list - integer(I4B), intent(in) :: nnew !! New size of list needed - end subroutine util_resize_encounter - module subroutine util_resize_pl(self, nnew) implicit none class(swiftest_pl), intent(inout) :: self !! Swiftest massive body object @@ -1899,14 +1790,6 @@ module subroutine util_spill_body(self, discards, lspill_list, ldestructive) logical, intent(in) :: ldestructive !! Logical flag indicating whether or not this operation should alter the keeps array or not end subroutine util_spill_body - module subroutine util_spill_encounter(self, discards, lspill_list, ldestructive) - implicit none - class(swiftest_encounter), intent(inout) :: self !! Swiftest encounter list - class(swiftest_encounter), intent(inout) :: discards !! Discarded object - logical, dimension(:), intent(in) :: lspill_list !! Logical array of bodies to spill into the discards - logical, intent(in) :: ldestructive !! Logical flag indicating whether or not this operation should alter body by removing the discard list - end subroutine util_spill_encounter - module subroutine util_spill_pl(self, discards, lspill_list, ldestructive) implicit none class(swiftest_pl), intent(inout) :: self !! Swiftest massive body object diff --git a/src/modules/symba_classes.f90 b/src/modules/symba_classes.f90 index 96f72f1f2..1c3605953 100644 --- a/src/modules/symba_classes.f90 +++ b/src/modules/symba_classes.f90 @@ -4,9 +4,10 @@ module symba_classes !! Definition of classes and methods specific to the SyMBA integrator !! Adapted from David E. Kaufmann's Swifter routine: module_symba.f90 use swiftest_globals - use swiftest_classes, only : swiftest_parameters, swiftest_base, swiftest_encounter, swiftest_particle_info, netcdf_parameters - use helio_classes, only : helio_cb, helio_pl, helio_tp, helio_nbody_system - use fraggle_classes, only : fraggle_colliders, fraggle_fragments + use swiftest_classes, only : swiftest_parameters, swiftest_base, swiftest_particle_info, netcdf_parameters + use helio_classes, only : helio_cb, helio_pl, helio_tp, helio_nbody_system + use fraggle_classes, only : fraggle_colliders, fraggle_fragments + use encounter_classes, only : encounter_list implicit none public @@ -119,16 +120,16 @@ module symba_classes ! symba_encounter class definitions and method interfaces !******************************************************************************************************************************* !> SyMBA class for tracking close encounters in a step - type, extends(swiftest_encounter) :: symba_encounter + type, extends(encounter_list) :: symba_encounter integer(I4B), dimension(:), allocatable :: level !! encounter recursion level contains procedure :: collision_check => symba_collision_check_encounter !! Checks if a test particle is going to collide with a massive body procedure :: encounter_check => symba_encounter_check !! Checks if massive bodies are going through close encounters with each other procedure :: kick => symba_kick_encounter !! Kick barycentric velocities of active test particles within SyMBA recursion - procedure :: setup => symba_setup_encounter !! A constructor that sets the number of encounters and allocates and initializes all arrays - procedure :: copy => symba_util_copy_encounter !! Copies elements from the source encounter list into self. - procedure :: spill => symba_util_spill_encounter !! "Spills" bodies from one object to another depending on the results of a mask (uses the PACK intrinsic) - procedure :: append => symba_util_append_encounter !! Appends elements from one structure to another + procedure :: setup => symba_setup_encounter_list !! A constructor that sets the number of encounters and allocates and initializes all arrays + procedure :: copy => symba_util_copy_encounter_list !! Copies elements from the source encounter list into self. + procedure :: spill => symba_util_spill_encounter_list !! "Spills" bodies from one object to another depending on the results of a mask (uses the PACK intrinsic) + procedure :: append => symba_util_append_encounter_list !! Appends elements from one structure to another end type symba_encounter !******************************************************************************************************************************** @@ -176,7 +177,7 @@ module symba_classes module function symba_collision_check_encounter(self, system, param, t, dt, irec) result(lany_collision) use swiftest_classes, only : swiftest_parameters implicit none - class(symba_encounter), intent(inout) :: self !! SyMBA pl-tp encounter list object + class(symba_encounter), intent(inout) :: self !! SyMBA pl-tp encounter list object class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters real(DP), intent(in) :: t !! current time @@ -418,11 +419,11 @@ module subroutine symba_setup_pl(self, n, param) class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters end subroutine symba_setup_pl - module subroutine symba_setup_encounter(self,n) + module subroutine symba_setup_encounter_list(self,n) implicit none class(symba_encounter), intent(inout) :: self !! SyMBA pl-tp encounter structure - integer(I4B), intent(in) :: n !! Number of encounters to allocate space for - end subroutine symba_setup_encounter + integer(I4B), intent(in) :: n !! Number of encounters to allocate space for + end subroutine symba_setup_encounter_list module subroutine symba_setup_tp(self, n, param) use swiftest_classes, only : swiftest_parameters @@ -483,12 +484,12 @@ end subroutine symba_util_append_arr_kin end interface interface - module subroutine symba_util_append_encounter(self, source, lsource_mask) + module subroutine symba_util_append_encounter_list(self, source, lsource_mask) implicit none class(symba_encounter), intent(inout) :: self !! SyMBA encounter list object - class(swiftest_encounter), intent(in) :: source !! Source object to append + class(encounter_list), intent(in) :: source !! Source object to append logical, dimension(:), intent(in) :: lsource_mask !! Logical mask indicating which elements to append to - end subroutine symba_util_append_encounter + end subroutine symba_util_append_encounter_list module subroutine symba_util_append_merger(self, source, lsource_mask) use swiftest_classes, only : swiftest_body @@ -514,12 +515,12 @@ module subroutine symba_util_append_tp(self, source, lsource_mask) logical, dimension(:), intent(in) :: lsource_mask !! Logical mask indicating which elements to append to end subroutine symba_util_append_tp - module subroutine symba_util_copy_encounter(self, source) - use swiftest_classes, only : swiftest_encounter + module subroutine symba_util_copy_encounter_list(self, source) + use encounter_classes, only : encounter_list implicit none class(symba_encounter), intent(inout) :: self !! Encounter list - class(swiftest_encounter), intent(in) :: source !! Source object to copy into - end subroutine symba_util_copy_encounter + class(encounter_list), intent(in) :: source !! Source object to copy into + end subroutine symba_util_copy_encounter_list end interface interface util_fill @@ -642,7 +643,7 @@ module subroutine symba_util_spill_arr_kin(keeps, discards, lspill_list, ldestru implicit none type(symba_kinship), dimension(:), allocatable, intent(inout) :: keeps !! Array of values to keep type(symba_kinship), dimension(:), allocatable, intent(inout) :: discards !! Array of discards - logical, dimension(:), intent(in) :: lspill_list !! Logical array of bodies to spill into the discardss + logical, dimension(:), intent(in) :: lspill_list !! Logical array of bodies to spill into the discardss logical, intent(in) :: ldestructive !! Logical flag indicating whether or not this operation should alter the keeps array or not end subroutine symba_util_spill_arr_kin end interface @@ -657,14 +658,14 @@ module subroutine symba_util_spill_pl(self, discards, lspill_list, ldestructive) logical, intent(in) :: ldestructive !! Logical flag indicating whether or not this operation should alter the keeps array or not end subroutine symba_util_spill_pl - module subroutine symba_util_spill_encounter(self, discards, lspill_list, ldestructive) - use swiftest_classes, only : swiftest_encounter + module subroutine symba_util_spill_encounter_list(self, discards, lspill_list, ldestructive) + use encounter_classes, only : encounter_list implicit none - class(symba_encounter), intent(inout) :: self !! SyMBA pl-tp encounter list - class(swiftest_encounter), intent(inout) :: discards !! Discarded object - logical, dimension(:), intent(in) :: lspill_list !! Logical array of bodies to spill into the discards - logical, intent(in) :: ldestructive !! Logical flag indicating whether or not this operation should alter body by removing the discard list - end subroutine symba_util_spill_encounter + class(symba_encounter), intent(inout) :: self !! SyMBA pl-tp encounter list + class(encounter_list), intent(inout) :: discards !! Discarded object + logical, dimension(:), intent(in) :: lspill_list !! Logical array of bodies to spill into the discards + logical, intent(in) :: ldestructive !! Logical flag indicating whether or not this operation should alter body by removing the discard list + end subroutine symba_util_spill_encounter_list module subroutine symba_util_spill_tp(self, discards, lspill_list, ldestructive) use swiftest_classes, only : swiftest_body diff --git a/src/rmvs/rmvs_io.f90 b/src/rmvs/rmvs_io.f90 index 51f077852..64190a6ea 100644 --- a/src/rmvs/rmvs_io.f90 +++ b/src/rmvs/rmvs_io.f90 @@ -33,7 +33,7 @@ module subroutine rmvs_io_write_encounter(t, id1, id2, Gmass1, Gmass2, radius1, call util_exit(FAILURE) end if lfirst = .false. - call io_write_frame_encounter(LUN, t, id1, id2, Gmass1, Gmass2, radius1, radius2, xh1, xh2, vh1, vh2) + call encounter_io_write_frame(LUN, t, id1, id2, Gmass1, Gmass2, radius1, radius2, xh1, xh2, vh1, vh2) close(unit = LUN, iostat = ierr) if (ierr /= 0) then write(*, *) "Swiftest Error:" diff --git a/src/setup/setup.f90 b/src/setup/setup.f90 index b8a47d782..1bc780ecc 100644 --- a/src/setup/setup.f90 +++ b/src/setup/setup.f90 @@ -71,61 +71,6 @@ module subroutine setup_construct_system(system, param) end subroutine setup_construct_system - module subroutine setup_encounter(self, n) - !! author: David A. Minton - !! - !! A constructor that sets the number of encounters and allocates and initializes all arrays - !! - implicit none - ! Arguments - class(swiftest_encounter), intent(inout) :: self !! Swiftest encounter structure - integer(I4B), intent(in) :: n !! Number of encounters to allocate space for - - if (n < 0) return - - if (allocated(self%lvdotr)) deallocate(self%lvdotr) - if (allocated(self%status)) deallocate(self%status) - if (allocated(self%index1)) deallocate(self%index1) - if (allocated(self%index2)) deallocate(self%index2) - if (allocated(self%id1)) deallocate(self%id1) - if (allocated(self%id2)) deallocate(self%id2) - if (allocated(self%x1)) deallocate(self%x1) - if (allocated(self%x2)) deallocate(self%x2) - if (allocated(self%v1)) deallocate(self%v1) - if (allocated(self%v2)) deallocate(self%v2) - if (allocated(self%t)) deallocate(self%t) - - self%nenc = n - if (n == 0) return - - allocate(self%lvdotr(n)) - allocate(self%status(n)) - allocate(self%index1(n)) - allocate(self%index2(n)) - allocate(self%id1(n)) - allocate(self%id2(n)) - allocate(self%x1(NDIM,n)) - allocate(self%x2(NDIM,n)) - allocate(self%v1(NDIM,n)) - allocate(self%v2(NDIM,n)) - allocate(self%t(n)) - - self%lvdotr(:) = .false. - self%status(:) = INACTIVE - self%index1(:) = 0 - self%index2(:) = 0 - self%id1(:) = 0 - self%id2(:) = 0 - self%x1(:,:) = 0.0_DP - self%x2(:,:) = 0.0_DP - self%v1(:,:) = 0.0_DP - self%v2(:,:) = 0.0_DP - self%t(:) = 0.0_DP - - return - end subroutine setup_encounter - - module subroutine setup_finalize_system(self, param) !! author: David A. Minton !! diff --git a/src/symba/symba_setup.f90 b/src/symba/symba_setup.f90 index 3a0543ee1..1c8879747 100644 --- a/src/symba/symba_setup.f90 +++ b/src/symba/symba_setup.f90 @@ -115,7 +115,7 @@ module subroutine symba_setup_pl(self, n, param) end subroutine symba_setup_pl - module subroutine symba_setup_encounter(self, n) + module subroutine symba_setup_encounter_list(self, n) !! author: David A. Minton !! !! A constructor that sets the number of encounters and allocates and initializes all arrays @@ -125,7 +125,7 @@ module subroutine symba_setup_encounter(self, n) class(symba_encounter), intent(inout) :: self !! SyMBA pl-tp encounter structure integer(I4B), intent(in) :: n !! Number of encounters to allocate space for - call setup_encounter(self, n) + call encounter_setup_list(self, n) if (n < 0) return if (allocated(self%level)) deallocate(self%level) @@ -137,7 +137,7 @@ module subroutine symba_setup_encounter(self, n) self%level(:) = -1 return - end subroutine symba_setup_encounter + end subroutine symba_setup_encounter_list module subroutine symba_setup_tp(self, n, param) diff --git a/src/symba/symba_util.f90 b/src/symba/symba_util.f90 index 4ed92030a..a4ad56eec 100644 --- a/src/symba/symba_util.f90 +++ b/src/symba/symba_util.f90 @@ -31,7 +31,7 @@ module subroutine symba_util_append_arr_kin(arr, source, nold, nsrc, lsource_mas end subroutine symba_util_append_arr_kin - module subroutine symba_util_append_encounter(self, source, lsource_mask) + module subroutine symba_util_append_encounter_list(self, source, lsource_mask) !! author: David A. Minton !! !! Append components from one encounter list (pl-pl or pl-tp) body object to another. @@ -39,7 +39,7 @@ module subroutine symba_util_append_encounter(self, source, lsource_mask) implicit none ! Arguments class(symba_encounter), intent(inout) :: self !! SyMBA encounter list object - class(swiftest_encounter), intent(in) :: source !! Source object to append + class(encounter_list), intent(in) :: source !! Source object to append logical, dimension(:), intent(in) :: lsource_mask !! Logical mask indicating which elements to append to associate(nold => self%nenc, nsrc => source%nenc) @@ -47,11 +47,11 @@ module subroutine symba_util_append_encounter(self, source, lsource_mask) class is (symba_encounter) call util_append(self%level, source%level, nold, nsrc, lsource_mask) end select - call util_append_encounter(self, source, lsource_mask) + call encounter_util_append_list(self, source, lsource_mask) end associate return - end subroutine symba_util_append_encounter + end subroutine symba_util_append_encounter_list module subroutine symba_util_append_pl(self, source, lsource_mask) @@ -159,14 +159,14 @@ module subroutine symba_util_append_tp(self, source, lsource_mask) end subroutine symba_util_append_tp - module subroutine symba_util_copy_encounter(self, source) + module subroutine symba_util_copy_encounter_list(self, source) !! author: David A. Minton !! !! Copies elements from the source encounter list into self. implicit none ! Arguments class(symba_encounter), intent(inout) :: self !! Encounter list - class(swiftest_encounter), intent(in) :: source !! Source object to copy into + class(encounter_list), intent(in) :: source !! Source object to copy into select type(source) class is (symba_encounter) @@ -175,10 +175,10 @@ module subroutine symba_util_copy_encounter(self, source) end associate end select - call util_copy_encounter(self, source) + call encounter_util_copy_list(self, source) return - end subroutine symba_util_copy_encounter + end subroutine symba_util_copy_encounter_list module subroutine symba_util_fill_arr_kin(keeps, inserts, lfill_list) @@ -978,7 +978,7 @@ module subroutine symba_util_spill_pl(self, discards, lspill_list, ldestructive) end subroutine symba_util_spill_pl - module subroutine symba_util_spill_encounter(self, discards, lspill_list, ldestructive) + module subroutine symba_util_spill_encounter_list(self, discards, lspill_list, ldestructive) !! author: David A. Minton !! !! Move spilled (discarded) SyMBA encounter structure from active list to discard list @@ -986,7 +986,7 @@ module subroutine symba_util_spill_encounter(self, discards, lspill_list, ldestr implicit none ! Arguments class(symba_encounter), intent(inout) :: self !! SyMBA pl-tp encounter list - class(swiftest_encounter), intent(inout) :: discards !! Discarded object + class(encounter_list), intent(inout) :: discards !! Discarded object logical, dimension(:), intent(in) :: lspill_list !! Logical array of bodies to spill into the discards logical, intent(in) :: ldestructive !! Logical flag indicating whether or not this operation should alter body by removing the discard list @@ -994,7 +994,7 @@ module subroutine symba_util_spill_encounter(self, discards, lspill_list, ldestr select type(discards) class is (symba_encounter) call util_spill(keeps%level, discards%level, lspill_list, ldestructive) - call util_spill_encounter(keeps, discards, lspill_list, ldestructive) + call encounter_util_spill_list(keeps, discards, lspill_list, ldestructive) class default write(*,*) "Invalid object passed to the spill method. Source must be of class symba_encounter or its descendents!" call util_exit(FAILURE) @@ -1002,7 +1002,7 @@ module subroutine symba_util_spill_encounter(self, discards, lspill_list, ldestr end associate return - end subroutine symba_util_spill_encounter + end subroutine symba_util_spill_encounter_list module subroutine symba_util_spill_tp(self, discards, lspill_list, ldestructive) diff --git a/src/util/util_append.f90 b/src/util/util_append.f90 index c518b3df3..69a771c6e 100644 --- a/src/util/util_append.f90 +++ b/src/util/util_append.f90 @@ -222,38 +222,6 @@ module subroutine util_append_body(self, source, lsource_mask) end subroutine util_append_body - module subroutine util_append_encounter(self, source, lsource_mask) - !! author: David A. Minton - !! - !! Append components from one Swiftest body object to another. - !! This method will automatically resize the destination body if it is too small - implicit none - ! Arguments - class(swiftest_encounter), intent(inout) :: self !! Swiftest encounter list object - class(swiftest_encounter), intent(in) :: source !! Source object to append - logical, dimension(:), intent(in) :: lsource_mask !! Logical mask indicating which elements to append to - ! Internals - integer(I4B) :: nold, nsrc - - nold = self%nenc - nsrc = source%nenc - call util_append(self%lvdotr, source%lvdotr, nold, nsrc, lsource_mask) - call util_append(self%status, source%status, nold, nsrc, lsource_mask) - call util_append(self%index1, source%index1, nold, nsrc, lsource_mask) - call util_append(self%index2, source%index2, nold, nsrc, lsource_mask) - call util_append(self%id1, source%id1, nold, nsrc, lsource_mask) - call util_append(self%id2, source%id2, nold, nsrc, lsource_mask) - call util_append(self%x1, source%x1, nold, nsrc, lsource_mask) - call util_append(self%x2, source%x2, nold, nsrc, lsource_mask) - call util_append(self%v1, source%v1, nold, nsrc, lsource_mask) - call util_append(self%v2, source%v2, nold, nsrc, lsource_mask) - call util_append(self%t, source%t, nold, nsrc, lsource_mask) - self%nenc = nold + count(lsource_mask(1:nsrc)) - - return - end subroutine util_append_encounter - - module subroutine util_append_pl(self, source, lsource_mask) !! author: David A. Minton !! diff --git a/src/util/util_copy.f90 b/src/util/util_copy.f90 index aba6f5b0e..04209f9f7 100644 --- a/src/util/util_copy.f90 +++ b/src/util/util_copy.f90 @@ -2,34 +2,6 @@ use swiftest contains - module subroutine util_copy_encounter(self, source) - !! author: David A. Minton - !! - !! Copies elements from the source encounter list into self. - implicit none - ! Arguments - class(swiftest_encounter), intent(inout) :: self !! Encounter list - class(swiftest_encounter), intent(in) :: source !! Source object to copy into - - associate(n => source%nenc) - self%nenc = n - self%lvdotr(1:n) = source%lvdotr(1:n) - self%status(1:n) = source%status(1:n) - self%index1(1:n) = source%index1(1:n) - self%index2(1:n) = source%index2(1:n) - self%id1(1:n) = source%id1(1:n) - self%id2(1:n) = source%id2(1:n) - self%x1(:,1:n) = source%x1(:,1:n) - self%x2(:,1:n) = source%x2(:,1:n) - self%v1(:,1:n) = source%v1(:,1:n) - self%v2(:,1:n) = source%v2(:,1:n) - self%t(1:n) = source%t(1:n) - end associate - - return - end subroutine util_copy_encounter - - module subroutine util_copy_particle_info(self, source) !! author: David A. Minton !! diff --git a/src/util/util_resize.f90 b/src/util/util_resize.f90 index 18b874607..4cfd52358 100644 --- a/src/util/util_resize.f90 +++ b/src/util/util_resize.f90 @@ -309,44 +309,6 @@ module subroutine util_resize_body(self, nnew) end subroutine util_resize_body - module subroutine util_resize_encounter(self, nnew) - !! author: David A. Minton - !! - !! Checks the current size of the encounter list against the required size and extends it by a factor of 2 more than requested if it is too small. - !! Note: The reason to extend it by a factor of 2 is for performance. When there are many enounters per step, resizing every time you want to add an - !! encounter takes significant computational effort. Resizing by a factor of 2 is a tradeoff between performance (fewer resize calls) and memory managment - !! Memory usage grows by a factor of 2 each time it fills up, but no more. - implicit none - ! Arguments - class(swiftest_encounter), intent(inout) :: self !! Swiftest encounter list - integer(I4B), intent(in) :: nnew !! New size of list needed - ! Internals - class(swiftest_encounter), allocatable :: enc_temp - integer(I4B) :: nold - logical :: lmalloc - - lmalloc = allocated(self%status) - if (lmalloc) then - nold = size(self%status) - else - nold = 0 - end if - if (nnew > nold) then - if (lmalloc) allocate(enc_temp, source=self) - call self%setup(2 * nnew) - if (lmalloc) then - call self%copy(enc_temp) - deallocate(enc_temp) - end if - else - self%status(nnew+1:nold) = INACTIVE - end if - self%nenc = nnew - - return - end subroutine util_resize_encounter - - module subroutine util_resize_pl(self, nnew) !! author: David A. Minton !! diff --git a/src/util/util_spill.f90 b/src/util/util_spill.f90 index fca48d1c1..06f5dfaf0 100644 --- a/src/util/util_spill.f90 +++ b/src/util/util_spill.f90 @@ -357,44 +357,6 @@ module subroutine util_spill_body(self, discards, lspill_list, ldestructive) end subroutine util_spill_body - module subroutine util_spill_encounter(self, discards, lspill_list, ldestructive) - !! author: David A. Minton - !! - !! Move spilled (discarded) Swiftest encounter structure from active list to discard list - implicit none - ! Arguments - class(swiftest_encounter), intent(inout) :: self !! Swiftest encounter list - class(swiftest_encounter), intent(inout) :: discards !! Discarded object - logical, dimension(:), intent(in) :: lspill_list !! Logical array of bodies to spill into the discards - logical, intent(in) :: ldestructive !! Logical flag indicating whether or not this operation should alter body by removing the discard list - ! Internals - integer(I4B) :: nenc_old - - associate(keeps => self) - call util_spill(keeps%lvdotr, discards%lvdotr, lspill_list, ldestructive) - call util_spill(keeps%status, discards%status, lspill_list, ldestructive) - call util_spill(keeps%index1, discards%index1, lspill_list, ldestructive) - call util_spill(keeps%index2, discards%index2, lspill_list, ldestructive) - call util_spill(keeps%id1, discards%id1, lspill_list, ldestructive) - call util_spill(keeps%id2, discards%id2, lspill_list, ldestructive) - call util_spill(keeps%x1, discards%x1, lspill_list, ldestructive) - call util_spill(keeps%x2, discards%x2, lspill_list, ldestructive) - call util_spill(keeps%v1, discards%v1, lspill_list, ldestructive) - call util_spill(keeps%v2, discards%v2, lspill_list, ldestructive) - call util_spill(keeps%t, discards%t, lspill_list, ldestructive) - - nenc_old = keeps%nenc - - ! This is the base class, so will be the last to be called in the cascade. - ! Therefore we need to set the nenc values for both the keeps and discareds - discards%nenc = count(lspill_list(1:nenc_old)) - if (ldestructive) keeps%nenc = nenc_old - discards%nenc - end associate - - return - end subroutine util_spill_encounter - - module subroutine util_spill_pl(self, discards, lspill_list, ldestructive) !! author: David A. Minton !!