From 8af4da899d43ffe810e0bbb3a9aa392148f1fe26 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Tue, 13 Dec 2022 14:08:43 -0500 Subject: [PATCH] Refactored x1,2 -> r1,2 --- src/encounter/encounter_check.f90 | 24 ++-- src/encounter/encounter_setup.f90 | 8 +- src/encounter/encounter_util.f90 | 209 +++++++++++++++++----------- src/modules/encounter_classes.f90 | 10 +- src/symba/symba_collision.f90 | 30 ++-- src/symba/symba_encounter_check.f90 | 4 +- src/symba/symba_util.f90 | 12 +- 7 files changed, 176 insertions(+), 121 deletions(-) diff --git a/src/encounter/encounter_check.f90 b/src/encounter/encounter_check.f90 index 6d866fb50..4e60ecf4f 100644 --- a/src/encounter/encounter_check.f90 +++ b/src/encounter/encounter_check.f90 @@ -891,7 +891,7 @@ pure module subroutine encounter_check_sort_aabb_1D(self, n, extent_arr) end subroutine encounter_check_sort_aabb_1D - module subroutine encounter_check_sweep_aabb_double_list(self, n1, n2, x1, v1, x2, v2, renc1, renc2, dt, & + module subroutine encounter_check_sweep_aabb_double_list(self, n1, n2, r1, v1, r2, v2, renc1, renc2, dt, & nenc, index1, index2, lvdotr) !! author: David A. Minton !! @@ -902,8 +902,8 @@ module subroutine encounter_check_sweep_aabb_double_list(self, n1, n2, x1, v1, x 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 - real(DP), dimension(:,:), intent(in) :: x1, v1 !! Array of position and velocity vectorrs for bodies 1 - real(DP), dimension(:,:), intent(in) :: x2, v2 !! Array of position and velocity vectorrs for bodies 2 + real(DP), dimension(:,:), intent(in) :: r1, v1 !! Array of position and velocity vectorrs for bodies 1 + real(DP), dimension(:,:), intent(in) :: r2, v2 !! Array of position and velocity vectorrs for bodies 2 real(DP), dimension(:), intent(in) :: renc1 !! Radius of encounter regions of bodies 1 real(DP), dimension(:), intent(in) :: renc2 !! Radius of encounter regions of bodies 2 real(DP), intent(in) :: dt !! Step size @@ -943,17 +943,17 @@ module subroutine encounter_check_sweep_aabb_double_list(self, n1, n2, x1, v1, x dim = 1 where(llist1(dim,:)) - xind(:) = x1(1,ext_ind(dim,:)) - yind(:) = x1(2,ext_ind(dim,:)) - zind(:) = x1(3,ext_ind(dim,:)) + xind(:) = r1(1,ext_ind(dim,:)) + yind(:) = r1(2,ext_ind(dim,:)) + zind(:) = r1(3,ext_ind(dim,:)) vxind(:) = v1(1,ext_ind(dim,:)) vyind(:) = v1(2,ext_ind(dim,:)) vzind(:) = v1(3,ext_ind(dim,:)) rencind(:) = renc1(ext_ind(dim,:)) elsewhere - xind(:) = x2(1,ext_ind(dim,:)) - yind(:) = x2(2,ext_ind(dim,:)) - zind(:) = x2(3,ext_ind(dim,:)) + xind(:) = r2(1,ext_ind(dim,:)) + yind(:) = r2(2,ext_ind(dim,:)) + zind(:) = r2(3,ext_ind(dim,:)) vxind(:) = v2(1,ext_ind(dim,:)) vyind(:) = v2(2,ext_ind(dim,:)) vzind(:) = v2(3,ext_ind(dim,:)) @@ -962,7 +962,7 @@ module subroutine encounter_check_sweep_aabb_double_list(self, n1, n2, x1, v1, x where(.not.loverlap(:)) lenc(:)%nenc = 0 !$omp parallel default(private) & - !$omp shared(self, ext_ind, lenc, loverlap, x1, v1, x2, v2, renc1, renc2, xind, yind, zind, vxind, vyind, vzind, rencind, llist1) & + !$omp shared(self, ext_ind, lenc, loverlap, r1, v1, r2, v2, renc1, renc2, xind, yind, zind, vxind, vyind, vzind, rencind, llist1) & !$omp firstprivate(ntot, n1, n2, dt, dim) ! Do the first group of bodies (i is in list 1, all the others are from list 2) @@ -972,7 +972,7 @@ module subroutine encounter_check_sweep_aabb_double_list(self, n1, n2, x1, v1, x ibeg = self%aabb(dim)%ibeg(i) + 1_I8B iend = self%aabb(dim)%iend(i) - 1_I8B nbox = iend - ibeg + 1 - call encounter_check_all_sweep_one(i, nbox, x1(1,i), x1(2,i), x1(3,i), v1(1,i), v1(2,i), v1(3,i), & + call encounter_check_all_sweep_one(i, nbox, r1(1,i), r1(2,i), r1(3,i), v1(1,i), v1(2,i), v1(3,i), & xind(ibeg:iend), yind(ibeg:iend), zind(ibeg:iend),& vxind(ibeg:iend), vyind(ibeg:iend), vzind(ibeg:iend), & renc1(i), rencind(ibeg:iend), dt, ext_ind(dim,ibeg:iend), & @@ -989,7 +989,7 @@ module subroutine encounter_check_sweep_aabb_double_list(self, n1, n2, x1, v1, x iend = self%aabb(dim)%iend(i) - 1_I8B nbox = iend - ibeg + 1 ii = i - n1 - call encounter_check_all_sweep_one(ii, nbox, x2(1,ii), x2(2,ii), x2(3,ii), v2(1,ii), v2(2,ii), v2(3,ii), & + call encounter_check_all_sweep_one(ii, nbox, r1(1,ii), r1(2,ii), r1(3,ii), v2(1,ii), v2(2,ii), v2(3,ii), & xind(ibeg:iend), yind(ibeg:iend), zind(ibeg:iend),& vxind(ibeg:iend), vyind(ibeg:iend), vzind(ibeg:iend), & renc2(ii), rencind(ibeg:iend), dt, ext_ind(dim,ibeg:iend), & diff --git a/src/encounter/encounter_setup.f90 b/src/encounter/encounter_setup.f90 index 4423952cd..18b60d229 100644 --- a/src/encounter/encounter_setup.f90 +++ b/src/encounter/encounter_setup.f90 @@ -79,8 +79,8 @@ module subroutine encounter_setup_list(self, 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%r1(NDIM,n)) + allocate(self%r2(NDIM,n)) allocate(self%v1(NDIM,n)) allocate(self%v2(NDIM,n)) @@ -91,8 +91,8 @@ module subroutine encounter_setup_list(self, n) self%index2(:) = 0 self%id1(:) = 0 self%id2(:) = 0 - self%x1(:,:) = 0.0_DP - self%x2(:,:) = 0.0_DP + self%r1(:,:) = 0.0_DP + self%r2(:,:) = 0.0_DP self%v1(:,:) = 0.0_DP self%v2(:,:) = 0.0_DP diff --git a/src/encounter/encounter_util.f90 b/src/encounter/encounter_util.f90 index 4386ec530..5d9e93a2e 100644 --- a/src/encounter/encounter_util.f90 +++ b/src/encounter/encounter_util.f90 @@ -33,8 +33,8 @@ module subroutine encounter_util_append_list(self, source, 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%r1, source%r1, nold, nsrc, lsource_mask) + call util_append(self%r2, source%r2, 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) self%nenc = nold + count(lsource_mask(1:nsrc)) @@ -62,8 +62,8 @@ module subroutine encounter_util_copy_list(self, source) 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%r1(:,1:n) = source%r1(:,1:n) + self%r2(:,1:n) = source%r2(:,1:n) self%v1(:,1:n) = source%v1(:,1:n) self%v2(:,1:n) = source%v2(:,1:n) end associate @@ -103,8 +103,8 @@ module subroutine encounter_util_dealloc_list(self) 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%r1)) deallocate(self%r1) + if (allocated(self%r2)) deallocate(self%r2) if (allocated(self%v1)) deallocate(self%v1) if (allocated(self%v2)) deallocate(self%v2) @@ -385,8 +385,8 @@ module subroutine encounter_util_spill_list(self, discards, lspill_list, ldestru 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%r1, discards%r1, lspill_list, ldestructive) + call util_spill(keeps%r2, discards%r2, lspill_list, ldestructive) call util_spill(keeps%v1, discards%v1, lspill_list, ldestructive) call util_spill(keeps%v2, discards%v2, lspill_list, ldestructive) @@ -570,7 +570,10 @@ module subroutine encounter_util_snapshot_encounter(self, param, system, t, arg) character(*), intent(in), optional :: arg !! Optional argument (needed for extended storage type used in collision snapshots) ! Arguments class(encounter_snapshot), allocatable :: snapshot - integer(I4B) :: i, npl_snap, ntp_snap + integer(I4B) :: i, j, k, npl_snap, ntp_snap + real(DP), dimension(NDIM) :: rrel, vrel, rcom, vcom + real(DP) :: mu, a, q, capm, tperi + real(DP), dimension(NDIM,2) :: rb,vb if (.not.present(t)) then write(*,*) "encounter_util_snapshot_encounter requires `t` to be passed" @@ -582,91 +585,139 @@ module subroutine encounter_util_snapshot_encounter(self, param, system, t, arg) return end if - select type (system) - class is (symba_nbody_system) - select case(arg) - case("trajectory") + select type(param) + class is (symba_parameters) + select type (system) + class is (symba_nbody_system) select type(pl => system%pl) class is (symba_pl) select type (tp => system%tp) class is (symba_tp) associate(npl => pl%nbody, ntp => tp%nbody) - allocate(encounter_snapshot :: snapshot) - snapshot%t = t - snapshot%iloop = param%iloop - if (npl + ntp == 0) return - npl_snap = npl - ntp_snap = ntp - + allocate(encounter_snapshot :: snapshot) allocate(snapshot%pl, mold=pl) allocate(snapshot%tp, mold=tp) - select type(pl_snap => snapshot%pl) - class is (symba_pl) - if (npl > 0) then - pl%lmask(1:npl) = pl%status(1:npl) /= INACTIVE .and. pl%levelg(1:npl) == system%irec - npl_snap = count(pl%lmask(1:npl)) - end if - if (ntp > 0) then - tp%lmask(1:ntp) = tp%status(1:ntp) /= INACTIVE .and. tp%levelg(1:ntp) == system%irec - ntp_snap = count(tp%lmask(1:ntp)) - end if - pl_snap%nbody = npl_snap - end select + snapshot%iloop = param%iloop select type(pl_snap => snapshot%pl) class is (symba_pl) - ! Take snapshot of the currently encountering massive bodies - if (npl_snap > 0) then - call pl_snap%setup(npl_snap, param) - pl_snap%levelg(:) = pack(pl%levelg(1:npl), pl%lmask(1:npl)) - pl_snap%id(:) = pack(pl%id(1:npl), pl%lmask(1:npl)) - pl_snap%info(:) = pack(pl%info(1:npl), pl%lmask(1:npl)) - pl_snap%Gmass(:) = pack(pl%Gmass(1:npl), pl%lmask(1:npl)) - do i = 1, NDIM - pl_snap%rh(i,:) = pack(pl%rh(i,1:npl), pl%lmask(1:npl)) - pl_snap%vh(i,:) = pack(pl%vb(i,1:npl), pl%lmask(1:npl)) - end do - if (param%lclose) then - pl_snap%radius(:) = pack(pl%radius(1:npl), pl%lmask(1:npl)) - end if - - if (param%lrotation) then - do i = 1, NDIM - pl_snap%Ip(i,:) = pack(pl%Ip(i,1:npl), pl%lmask(1:npl)) - pl_snap%rot(i,:) = pack(pl%rot(i,1:npl), pl%lmask(1:npl)) - end do - end if - call pl_snap%sort("id", ascending=.true.) - end if - end select - - select type(tp_snap => snapshot%tp) - class is (symba_tp) - ! Take snapshot of the currently encountering test particles - tp_snap%nbody = ntp_snap - if (ntp_snap > 0) then - call tp_snap%setup(ntp_snap, param) - tp_snap%id(:) = pack(tp%id(1:ntp), tp%lmask(1:ntp)) - tp_snap%info(:) = pack(tp%info(1:ntp), tp%lmask(1:ntp)) - do i = 1, NDIM - tp_snap%rh(i,:) = pack(tp%rh(i,1:ntp), tp%lmask(1:ntp)) - tp_snap%vh(i,:) = pack(tp%vh(i,1:ntp), tp%lmask(1:ntp)) - end do - end if + select type(tp_snap => snapshot%tp) + class is (symba_tp) + + select case(arg) + case("trajectory") + snapshot%t = t + + npl_snap = npl + ntp_snap = ntp + + if (npl > 0) then + pl%lmask(1:npl) = pl%status(1:npl) /= INACTIVE .and. pl%levelg(1:npl) == system%irec + npl_snap = count(pl%lmask(1:npl)) + end if + if (ntp > 0) then + tp%lmask(1:ntp) = tp%status(1:ntp) /= INACTIVE .and. tp%levelg(1:ntp) == system%irec + ntp_snap = count(tp%lmask(1:ntp)) + end if + + pl_snap%nbody = npl_snap + + ! Take snapshot of the currently encountering massive bodies + if (npl_snap > 0) then + call pl_snap%setup(npl_snap, param) + pl_snap%levelg(:) = pack(pl%levelg(1:npl), pl%lmask(1:npl)) + pl_snap%id(:) = pack(pl%id(1:npl), pl%lmask(1:npl)) + pl_snap%info(:) = pack(pl%info(1:npl), pl%lmask(1:npl)) + pl_snap%Gmass(:) = pack(pl%Gmass(1:npl), pl%lmask(1:npl)) + do i = 1, NDIM + pl_snap%rh(i,:) = pack(pl%rh(i,1:npl), pl%lmask(1:npl)) + pl_snap%vh(i,:) = pack(pl%vb(i,1:npl), pl%lmask(1:npl)) + end do + if (param%lclose) then + pl_snap%radius(:) = pack(pl%radius(1:npl), pl%lmask(1:npl)) + end if + + if (param%lrotation) then + do i = 1, NDIM + pl_snap%Ip(i,:) = pack(pl%Ip(i,1:npl), pl%lmask(1:npl)) + pl_snap%rot(i,:) = pack(pl%rot(i,1:npl), pl%lmask(1:npl)) + end do + end if + call pl_snap%sort("id", ascending=.true.) + end if + + ! Take snapshot of the currently encountering test particles + tp_snap%nbody = ntp_snap + if (ntp_snap > 0) then + call tp_snap%setup(ntp_snap, param) + tp_snap%id(:) = pack(tp%id(1:ntp), tp%lmask(1:ntp)) + tp_snap%info(:) = pack(tp%info(1:ntp), tp%lmask(1:ntp)) + do i = 1, NDIM + tp_snap%rh(i,:) = pack(tp%rh(i,1:ntp), tp%lmask(1:ntp)) + tp_snap%vh(i,:) = pack(tp%vh(i,1:ntp), tp%lmask(1:ntp)) + end do + end if + + ! Save the snapshot + param%encounter_history%nid = param%encounter_history%nid + ntp_snap + npl_snap + call encounter_util_save_encounter(param%encounter_history,snapshot,t) + case("closest") + associate(plplenc_list => system%plplenc_list, pltpenc_list => system%pltpenc_list) + if (any(plplenc_list%lclosest(:))) then + call pl_snap%setup(2, param) + do k = 1, plplenc_list%nenc + if (plplenc_list%lclosest(k)) then + i = plplenc_list%index1(k) + j = plplenc_list%index2(k) + pl_snap%levelg(:) = pl%levelg([i,j]) + pl_snap%id(:) = pl%id([i,j]) + pl_snap%info(:) = pl%info([i,j]) + pl_snap%Gmass(:) = pl%Gmass([i,j]) + mu = sum(pl_snap%Gmass(:)) + if (param%lclose) pl_snap%radius(:) = pl%radius([i,j]) + if (param%lrotation) then + do i = 1, NDIM + pl_snap%Ip(i,:) = pl%Ip(i,[i,j]) + pl_snap%rot(i,:) = pl%rot(i,[i,j]) + end do + end if + + ! Compute pericenter passage time to get the closest approach parameters + rrel(:) = plplenc_list%r2(:,k) - plplenc_list%r1(:,k) + vrel(:) = plplenc_list%v2(:,k) - plplenc_list%v1(:,k) + call orbel_xv2aqt(mu, rrel(1), rrel(2), rrel(3), vrel(1), vrel(2), vrel(3), a, q, capm, tperi) + snapshot%t = t + tperi + + ! Computer the center mass of the pair + + ! do i = 1, NDIM + ! pl_snap%rh(i,:) = pl%rh(i,[i,j]) + ! pl_snap%vh(i,:) = pl%vb(i,1:npl), pl%lmask(1:npl)) + ! end do + + call pl_snap%sort("id", ascending=.true.) + call encounter_util_save_encounter(param%encounter_history,snapshot,snapshot%t) + end if + end do + + plplenc_list%lclosest(:) = .false. + end if + + if (any(pltpenc_list%lclosest(:))) then + do k = 1, pltpenc_list%nenc + end do + pltpenc_list%lclosest(:) = .false. + end if + end associate + case default + write(*,*) "encounter_util_snapshot_encounter requires `arg` to be either `trajectory` or `closest`" + end select + end select end select end associate - ! Save the snapshot - select type(param) - class is (symba_parameters) - param%encounter_history%nid = param%encounter_history%nid + ntp_snap + npl_snap - call encounter_util_save_encounter(param%encounter_history,snapshot,t) - end select end select end select - case("closest") - case default - write(*,*) "encounter_util_snapshot_encounter requires `arg` to be either `trajectory` or `closest`" end select end select diff --git a/src/modules/encounter_classes.f90 b/src/modules/encounter_classes.f90 index f566a8070..2b313eeb5 100644 --- a/src/modules/encounter_classes.f90 +++ b/src/modules/encounter_classes.f90 @@ -29,8 +29,8 @@ module encounter_classes 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 :: r1 !! the position of body 1 in the encounter + real(DP), dimension(:,:), allocatable :: r2 !! 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 contains @@ -186,14 +186,14 @@ pure 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, x1, v1, x2, v2, renc1, renc2, dt, & + module subroutine encounter_check_sweep_aabb_double_list(self, n1, n2, r1, v1, r2, v2, renc1, renc2, dt, & nenc, index1, index2, lvdotr) 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 - real(DP), dimension(:,:), intent(in) :: x1, v1 !! Array of indices of bodies 1 - real(DP), dimension(:,:), intent(in) :: x2, v2 !! Array of indices of bodies 2 + real(DP), dimension(:,:), intent(in) :: r1, v1 !! Array of indices of bodies 1 + real(DP), dimension(:,:), intent(in) :: r2, v2 !! Array of indices of bodies 2 real(DP), dimension(:), intent(in) :: renc1 !! Radius of encounter regions of bodies 1 real(DP), dimension(:), intent(in) :: renc2 !! Radius of encounter regions of bodies 2 real(DP), intent(in) :: dt !! Step size diff --git a/src/symba/symba_collision.f90 b/src/symba/symba_collision.f90 index 9015d3403..abfec217d 100644 --- a/src/symba/symba_collision.f90 +++ b/src/symba/symba_collision.f90 @@ -280,11 +280,11 @@ module subroutine symba_collision_check_encounter(self, system, param, t, dt, ir integer(I4B), intent(in) :: irec !! Current recursion level logical, intent(out) :: lany_collision !! Returns true if cany pair of encounters resulted in a collision ! Internals - logical, dimension(:), allocatable :: lcollision, lclosest, lmask + logical, dimension(:), allocatable :: lcollision, lmask real(DP), dimension(NDIM) :: xr, vr integer(I4B) :: i, j, k, nenc real(DP) :: rlim, Gmtot - logical :: isplpl + logical :: isplpl, lany_closest character(len=STRMAX) :: timestr, idstri, idstrj, message class(symba_encounter), allocatable :: tmp @@ -316,8 +316,7 @@ module subroutine symba_collision_check_encounter(self, system, param, t, dt, ir allocate(lcollision(nenc)) lcollision(:) = .false. - allocate(lclosest(nenc)) - lclosest(:) = .false. + self%lclosest(:) = .false. if (isplpl) then do concurrent(k = 1:nenc, lmask(k)) @@ -327,7 +326,7 @@ module subroutine symba_collision_check_encounter(self, system, param, t, dt, ir vr(:) = pl%vb(:, i) - pl%vb(:, j) rlim = pl%radius(i) + pl%radius(j) Gmtot = pl%Gmass(i) + pl%Gmass(j) - call symba_collision_check_one(xr(1), xr(2), xr(3), vr(1), vr(2), vr(3), Gmtot, rlim, dt, self%lvdotr(k), lcollision(k), lclosest(k)) + call symba_collision_check_one(xr(1), xr(2), xr(3), vr(1), vr(2), vr(3), Gmtot, rlim, dt, self%lvdotr(k), lcollision(k), self%lclosest(k)) end do else do concurrent(k = 1:nenc, lmask(k)) @@ -335,23 +334,28 @@ module subroutine symba_collision_check_encounter(self, system, param, t, dt, ir j = self%index2(k) xr(:) = pl%rh(:, i) - tp%rh(:, j) vr(:) = pl%vb(:, i) - tp%vb(:, j) - call symba_collision_check_one(xr(1), xr(2), xr(3), vr(1), vr(2), vr(3), pl%Gmass(i), pl%radius(i), dt, self%lvdotr(k), lcollision(k), lclosest(k)) + call symba_collision_check_one(xr(1), xr(2), xr(3), vr(1), vr(2), vr(3), pl%Gmass(i), pl%radius(i), dt, self%lvdotr(k), lcollision(k), self%lclosest(k)) end do end if lany_collision = any(lcollision(:)) + lany_closest = (param%lenc_save_closest .and. any(self%lclosest(:))) - if (lany_collision) then + + if (lany_collision .or. lany_closest) then call pl%rh2rb(system%cb) ! Update the central body barycenteric position vector to get us out of DH and into bary do k = 1, nenc + if (.not.lcollision(k) .and. .not. self%lclosest(k)) cycle i = self%index1(k) j = self%index2(k) - if (lcollision(k)) self%status(k) = COLLISION - self%tcollision(k) = t - self%x1(:,k) = pl%rh(:,i) + system%cb%rb(:) + self%r1(:,k) = pl%rh(:,i) + system%cb%rb(:) self%v1(:,k) = pl%vb(:,i) + if (lcollision(k)) then + self%status(k) = COLLISION + self%tcollision(k) = t + end if if (isplpl) then - self%x2(:,k) = pl%rh(:,j) + system%cb%rb(:) + self%r2(:,k) = pl%rh(:,j) + system%cb%rb(:) self%v2(:,k) = pl%vb(:,j) if (lcollision(k)) then ! Check to see if either of these bodies has been involved with a collision before, and if so, make this a collider pair @@ -364,7 +368,7 @@ module subroutine symba_collision_check_encounter(self, system, param, t, dt, ir call pl%info(j)%set_value(status="COLLISION", discard_time=t, discard_rh=pl%rh(:,j), discard_vh=pl%vh(:,j)) end if else - self%x2(:,k) = tp%rh(:,j) + system%cb%rb(:) + self%r2(:,k) = tp%rh(:,j) + system%cb%rb(:) self%v2(:,k) = tp%vb(:,j) if (lcollision(k)) then tp%status(j) = DISCARDED_PLR @@ -392,7 +396,7 @@ module subroutine symba_collision_check_encounter(self, system, param, t, dt, ir end if ! Take snapshots of pairs of bodies at close approach (but not collision) if requested - if (param%lenc_save_closest .and. any(lclosest(:))) call param%encounter_history%take_snapshot(param, system, t, "closest") + if (lany_closest) call param%encounter_history%take_snapshot(param, system, t, "closest") end select diff --git a/src/symba/symba_encounter_check.f90 b/src/symba/symba_encounter_check.f90 index dd60f2e00..f016af9d9 100644 --- a/src/symba/symba_encounter_check.f90 +++ b/src/symba/symba_encounter_check.f90 @@ -65,8 +65,8 @@ module function symba_encounter_check_pl(self, param, system, dt, irec) result(l plplenc_list%id2(k) = pl%id(j) plplenc_list%status(k) = ACTIVE plplenc_list%level(k) = irec - plplenc_list%x1(:,k) = pl%rh(:,i) - plplenc_list%x2(:,k) = pl%rh(:,j) + plplenc_list%r1(:,k) = pl%rh(:,i) + plplenc_list%r2(:,k) = pl%rh(:,j) plplenc_list%v1(:,k) = pl%vb(:,i) - cb%vb(:) plplenc_list%v2(:,k) = pl%vb(:,j) - cb%vb(:) pl%lencounter(i) = .true. diff --git a/src/symba/symba_util.f90 b/src/symba/symba_util.f90 index 74b555bce..06d75bac8 100644 --- a/src/symba/symba_util.f90 +++ b/src/symba/symba_util.f90 @@ -716,8 +716,8 @@ module subroutine symba_util_rearray_pl(self, system, param) system%plplenc_list%lvdotr(k) = plplenc_old%lvdotr(k) system%plplenc_list%lclosest(k) = plplenc_old%lclosest(k) system%plplenc_list%status(k) = plplenc_old%status(k) - system%plplenc_list%x1(:,k) = plplenc_old%x1(:,k) - system%plplenc_list%x2(:,k) = plplenc_old%x2(:,k) + system%plplenc_list%r1(:,k) = plplenc_old%r1(:,k) + system%plplenc_list%r2(:,k) = plplenc_old%r2(:,k) system%plplenc_list%v1(:,k) = plplenc_old%v1(:,k) system%plplenc_list%v2(:,k) = plplenc_old%v2(:,k) system%plplenc_list%tcollision(k) = plplenc_old%tcollision(k) @@ -727,8 +727,8 @@ module subroutine symba_util_rearray_pl(self, system, param) system%plplenc_list%lvdotr(k) = plplenc_old%lvdotr(k) system%plplenc_list%lclosest(k) = plplenc_old%lclosest(k) system%plplenc_list%status(k) = plplenc_old%status(k) - system%plplenc_list%x1(:,k) = plplenc_old%x2(:,k) - system%plplenc_list%x2(:,k) = plplenc_old%x1(:,k) + system%plplenc_list%r1(:,k) = plplenc_old%r2(:,k) + system%plplenc_list%r2(:,k) = plplenc_old%r1(:,k) system%plplenc_list%v1(:,k) = plplenc_old%v2(:,k) system%plplenc_list%v2(:,k) = plplenc_old%v1(:,k) system%plplenc_list%tcollision(k) = plplenc_old%tcollision(k) @@ -758,8 +758,8 @@ module subroutine symba_util_rearray_pl(self, system, param) system%plplenc_list%tcollision(1:nencmin) = pack(system%plplenc_list%tcollision(1:nenc_old), lmask(1:nenc_old)) system%plplenc_list%level(1:nencmin) = pack(system%plplenc_list%level(1:nenc_old), lmask(1:nenc_old)) do i = 1, NDIM - system%plplenc_list%x1(i, 1:nencmin) = pack(system%plplenc_list%x1(i, 1:nenc_old), lmask(1:nenc_old)) - system%plplenc_list%x2(i, 1:nencmin) = pack(system%plplenc_list%x2(i, 1:nenc_old), lmask(1:nenc_old)) + system%plplenc_list%r1(i, 1:nencmin) = pack(system%plplenc_list%r1(i, 1:nenc_old), lmask(1:nenc_old)) + system%plplenc_list%r2(i, 1:nencmin) = pack(system%plplenc_list%r2(i, 1:nenc_old), lmask(1:nenc_old)) system%plplenc_list%v1(i, 1:nencmin) = pack(system%plplenc_list%v1(i, 1:nenc_old), lmask(1:nenc_old)) system%plplenc_list%v2(i, 1:nencmin) = pack(system%plplenc_list%v2(i, 1:nenc_old), lmask(1:nenc_old)) end do