From 64909af626beeb85141c87f367ecf95b0ba3b2cf Mon Sep 17 00:00:00 2001 From: David A Minton Date: Tue, 13 Dec 2022 11:12:23 -0500 Subject: [PATCH] More restructruing anc getting close approach saving infrastructure in place --- python/swiftest/swiftest/io.py | 11 +- src/encounter/encounter_setup.f90 | 14 +- src/encounter/encounter_util.f90 | 219 ++++++++++++++++-------------- src/modules/encounter_classes.f90 | 3 +- src/modules/swiftest_classes.f90 | 2 +- src/modules/symba_classes.f90 | 8 +- src/setup/setup.f90 | 2 +- src/symba/symba_collision.f90 | 180 ++++++++++++------------ src/symba/symba_io.f90 | 4 +- src/symba/symba_step.f90 | 6 +- src/symba/symba_util.f90 | 11 +- 11 files changed, 237 insertions(+), 223 deletions(-) diff --git a/python/swiftest/swiftest/io.py b/python/swiftest/swiftest/io.py index 029672ec7..c002978b9 100644 --- a/python/swiftest/swiftest/io.py +++ b/python/swiftest/swiftest/io.py @@ -900,15 +900,14 @@ def string_converter(da): ------- da : xarray dataset with the strings cleaned up """ - if da.dtype == np.dtype(object): - da = da.astype(' self) call util_spill(keeps%lvdotr, discards%lvdotr, lspill_list, ldestructive) + call util_spill(keeps%lclosest, discards%lclosest, 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) @@ -398,7 +403,7 @@ end subroutine encounter_util_spill_list - subroutine encounter_util_save_collision(param, snapshot) + subroutine encounter_util_save_collision(collision_history, snapshot) !! author: David A. Minton !! !! Checks the current size of the encounter storage against the required size and extends it by a factor of 2 more than requested if it is too small. @@ -407,18 +412,18 @@ subroutine encounter_util_save_collision(param, snapshot) !! Memory usage grows by a factor of 2 each time it fills up, but no more. implicit none ! Arguments - class(symba_parameters), intent(inout) :: param !! SyMBA parameter object - class(encounter_snapshot), intent(in) :: snapshot !! Encounter snapshot object + type(collision_storage(*)), allocatable, intent(inout) :: collision_history !! Collision history object + class(encounter_snapshot), intent(in) :: snapshot !! Encounter snapshot object ! Internals type(collision_storage(nframes=:)), allocatable :: tmp integer(I4B) :: i, nnew, nold, nbig ! Advance the snapshot frame counter - param%collision_history%iframe = param%collision_history%iframe + 1 + collision_history%iframe = collision_history%iframe + 1 ! Check to make sure the current encounter_history object is big enough. If not, grow it by a factor of 2 - nnew = param%collision_history%iframe - nold = param%collision_history%nframes + nnew = collision_history%iframe + nold = collision_history%nframes if (nnew > nold) then nbig = nold @@ -426,24 +431,24 @@ subroutine encounter_util_save_collision(param, snapshot) nbig = nbig * 2 end do allocate(collision_storage(nbig) :: tmp) - tmp%iframe = param%collision_history%iframe - call move_alloc(param%collision_history%nc, tmp%nc) + tmp%iframe = collision_history%iframe + call move_alloc(collision_history%nc, tmp%nc) do i = 1, nold - if (allocated(param%collision_history%frame(i)%item)) call move_alloc(param%collision_history%frame(i)%item, tmp%frame(i)%item) + if (allocated(collision_history%frame(i)%item)) call move_alloc(collision_history%frame(i)%item, tmp%frame(i)%item) end do - deallocate(param%collision_history) - call move_alloc(tmp,param%collision_history) + deallocate(collision_history) + call move_alloc(tmp,collision_history) nnew = nbig end if - param%collision_history%frame(nnew) = snapshot + collision_history%frame(nnew) = snapshot return end subroutine encounter_util_save_collision - subroutine encounter_util_save_encounter(param, snapshot, t) + subroutine encounter_util_save_encounter(encounter_history, snapshot, t) !! author: David A. Minton !! !! Checks the current size of the encounter storage against the required size and extends it by a factor of 2 more than requested if it is too small. @@ -452,19 +457,19 @@ subroutine encounter_util_save_encounter(param, snapshot, t) !! Memory usage grows by a factor of 2 each time it fills up, but no more. implicit none ! Arguments - class(symba_parameters), intent(inout) :: param !! SyMBA parameter object - class(encounter_snapshot), intent(in) :: snapshot !! Encounter snapshot object - real(DP), intent(in) :: t !! The time of the snapshot + type(encounter_storage(*)), allocatable, intent(inout) :: encounter_history !! SyMBA encounter storage object + class(encounter_snapshot), intent(in) :: snapshot !! Encounter snapshot object + real(DP), intent(in) :: t !! The time of the snapshot ! Internals type(encounter_storage(nframes=:)), allocatable :: tmp integer(I4B) :: i, nnew, nold, nbig ! Advance the snapshot frame counter - param%encounter_history%iframe = param%encounter_history%iframe + 1 + encounter_history%iframe = encounter_history%iframe + 1 ! Check to make sure the current encounter_history object is big enough. If not, grow it by a factor of 2 - nnew = param%encounter_history%iframe - nold = param%encounter_history%nframes + nnew = encounter_history%iframe + nold = encounter_history%nframes if (nnew > nold) then nbig = nold @@ -472,20 +477,20 @@ subroutine encounter_util_save_encounter(param, snapshot, t) nbig = nbig * 2 end do allocate(encounter_storage(nbig) :: tmp) - tmp%iframe = param%encounter_history%iframe - call move_alloc(param%encounter_history%nc, tmp%nc) + tmp%iframe = encounter_history%iframe + call move_alloc(encounter_history%nc, tmp%nc) do i = 1, nold - if (allocated(param%encounter_history%frame(i)%item)) call move_alloc(param%encounter_history%frame(i)%item, tmp%frame(i)%item) + if (allocated(encounter_history%frame(i)%item)) call move_alloc(encounter_history%frame(i)%item, tmp%frame(i)%item) end do - deallocate(param%encounter_history) - call move_alloc(tmp,param%encounter_history) + deallocate(encounter_history) + call move_alloc(tmp,encounter_history) nnew = nbig end if ! Find out which time slot this belongs in by searching for an existing slot ! with the same value of time or the first available one - param%encounter_history%frame(nnew) = snapshot + encounter_history%frame(nnew) = snapshot return end subroutine encounter_util_save_encounter @@ -502,13 +507,11 @@ module subroutine encounter_util_snapshot_collision(self, param, system, t, arg) class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object to store real(DP), intent(in), optional :: t !! Time of snapshot if different from system time - character(*), intent(in), optional :: arg !! Optional argument (needed for extended storage type used in collision snapshots) + character(*), intent(in), optional :: arg !! "before": takes a snapshot just before the collision. "after" takes the snapshot just after the collision. ! Arguments class(fraggle_snapshot), allocatable :: snapshot type(symba_pl) :: pl character(len=:), allocatable :: stage - integer(I4B) :: i,j - if (present(arg)) then stage = arg @@ -539,10 +542,10 @@ module subroutine encounter_util_snapshot_collision(self, param, system, t, arg) allocate(fraggle_snapshot :: snapshot) allocate(snapshot%colliders, source=system%colliders) allocate(snapshot%fragments, source=system%fragments) - select type (param) + select type(param) class is (symba_parameters) - call encounter_util_save_collision(param,snapshot) - end select + call encounter_util_save_collision(param%collision_history,snapshot) + end select case default write(*,*) "encounter_util_snapshot_collision requies either 'before' or 'after' passed to 'arg'" end select @@ -571,87 +574,99 @@ module subroutine encounter_util_snapshot_encounter(self, param, system, t, arg) if (.not.present(t)) then write(*,*) "encounter_util_snapshot_encounter requires `t` to be passed" + return end if + + if (.not.present(arg)) then + write(*,*) "encounter_util_snapshot_encounter requires `arg` to be passed" + return + end if + 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(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 - - 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)) + select case(arg) + case("trajectory") + 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(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 (param%lrotation) then + 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 + + 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%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)) + 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 - 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 + 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 + 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 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,snapshot,t) 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 8e74913ed..f566a8070 100644 --- a/src/modules/encounter_classes.f90 +++ b/src/modules/encounter_classes.f90 @@ -22,6 +22,7 @@ module encounter_classes integer(I8B) :: nenc = 0 !! Total number of encounters logical :: lcollision !! Indicates if the encounter resulted in at least one collision real(DP) :: t !! Time of encounter + logical, dimension(:), allocatable :: lclosest !! indicates that thie pair of bodies is in currently at its closest approach point 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 @@ -329,7 +330,7 @@ module subroutine encounter_util_snapshot_collision(self, param, system, t, arg) class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object to store real(DP), intent(in), optional :: t !! Time of snapshot if different from system time - character(*), intent(in), optional :: arg !! Optional argument (needed for extended storage type used in collision snapshots) + character(*), intent(in), optional :: arg !! "before": takes a snapshot just before the collision. "after" takes the snapshot just after the collision. end subroutine encounter_util_snapshot_collision module subroutine encounter_util_snapshot_encounter(self, param, system, t, arg) diff --git a/src/modules/swiftest_classes.f90 b/src/modules/swiftest_classes.f90 index d86a742b2..2e32f8c1d 100644 --- a/src/modules/swiftest_classes.f90 +++ b/src/modules/swiftest_classes.f90 @@ -1705,7 +1705,7 @@ module subroutine util_snapshot_system(self, param, system, t, arg) class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object to store real(DP), intent(in), optional :: t !! Time of snapshot if different from system time - character(*), intent(in), optional :: arg !! Optional argument (needed for extended storage type used in collision snapshots) + character(*), intent(in), optional :: arg !! Optional argument (needed for extended storage type used in encounter snapshots) end subroutine util_snapshot_system end interface diff --git a/src/modules/symba_classes.f90 b/src/modules/symba_classes.f90 index ce95269e2..1d5a708b3 100644 --- a/src/modules/symba_classes.f90 +++ b/src/modules/symba_classes.f90 @@ -16,7 +16,7 @@ module symba_classes use swiftest_classes, only : swiftest_parameters, swiftest_base, swiftest_particle_info, swiftest_storage, 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, encounter_snapshot, encounter_storage, collision_storage + use encounter_classes, only : encounter_list, encounter_storage, collision_storage implicit none public @@ -31,8 +31,8 @@ module symba_classes integer(I4B), dimension(:), allocatable :: seed !! Random seeds for fragmentation modeling logical :: lfragmentation = .false. !! Do fragmentation modeling instead of simple merger. character(STRMAX) :: encounter_save = "NONE" !! Indicate if and how encounter data should be saved - logical :: lenc_trajectory_save = .false. !! Indicates that when encounters are saved, the full trajectory through recursion steps are saved - logical :: lenc_closest_save = .false. !! Indicates that when encounters are saved, the closest approach distance between pairs of bodies is saved + logical :: lenc_save_trajectory = .false. !! Indicates that when encounters are saved, the full trajectory through recursion steps are saved + logical :: lenc_save_closest = .false. !! Indicates that when encounters are saved, the closest approach distance between pairs of bodies is saved type(encounter_storage(nframes=:)), allocatable :: encounter_history !! Stores encounter history for later retrieval and saving to file type(collision_storage(nframes=:)), allocatable :: collision_history !! Stores encounter history for later retrieval and saving to file contains @@ -212,7 +212,7 @@ module subroutine symba_collision_check_encounter(self, system, param, t, dt, ir implicit none 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 + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters real(DP), intent(in) :: t !! current time real(DP), intent(in) :: dt !! step size integer(I4B), intent(in) :: irec !! Current recursion level diff --git a/src/setup/setup.f90 b/src/setup/setup.f90 index b840528a0..ef8558aef 100644 --- a/src/setup/setup.f90 +++ b/src/setup/setup.f90 @@ -75,7 +75,7 @@ module subroutine setup_construct_system(system, param) select type(param) class is (symba_parameters) - if (param%lenc_trajectory_save .or. param%lenc_closest_save) then + if (param%lenc_save_trajectory .or. param%lenc_save_closest) then allocate(encounter_storage :: param%encounter_history) associate (encounter_history => param%encounter_history) allocate(encounter_io_parameters :: encounter_history%nc) diff --git a/src/symba/symba_collision.f90 b/src/symba/symba_collision.f90 index 637efa412..9015d3403 100644 --- a/src/symba/symba_collision.f90 +++ b/src/symba/symba_collision.f90 @@ -274,7 +274,7 @@ module subroutine symba_collision_check_encounter(self, system, param, t, dt, ir ! Arguments 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 + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters real(DP), intent(in) :: t !! current time real(DP), intent(in) :: dt !! step size integer(I4B), intent(in) :: irec !! Current recursion level @@ -302,97 +302,99 @@ module subroutine symba_collision_check_encounter(self, system, param, t, dt, ir class is (symba_pl) select type(tp => system%tp) class is (symba_tp) - nenc = self%nenc - allocate(lmask(nenc)) - lmask(:) = ((self%status(1:nenc) == ACTIVE) .and. (pl%levelg(self%index1(1:nenc)) >= irec)) - if (isplpl) then - lmask(:) = lmask(:) .and. (pl%levelg(self%index2(1:nenc)) >= irec) - else - lmask(:) = lmask(:) .and. (tp%levelg(self%index2(1:nenc)) >= irec) - end if - if (.not.any(lmask(:))) return - - allocate(lcollision(nenc)) - lcollision(:) = .false. - allocate(lclosest(nenc)) - lclosest(:) = .false. - - if (isplpl) then - do concurrent(k = 1:nenc, lmask(k)) - i = self%index1(k) - j = self%index2(k) - xr(:) = pl%rh(:, i) - pl%rh(:, j) - 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)) - end do - else - do concurrent(k = 1:nenc, lmask(k)) - i = self%index1(k) - 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)) - end do - end if + select type (param) + class is (symba_parameters) + nenc = self%nenc + allocate(lmask(nenc)) + lmask(:) = ((self%status(1:nenc) == ACTIVE) .and. (pl%levelg(self%index1(1:nenc)) >= irec)) + if (isplpl) then + lmask(:) = lmask(:) .and. (pl%levelg(self%index2(1:nenc)) >= irec) + else + lmask(:) = lmask(:) .and. (tp%levelg(self%index2(1:nenc)) >= irec) + end if + if (.not.any(lmask(:))) return + + allocate(lcollision(nenc)) + lcollision(:) = .false. + allocate(lclosest(nenc)) + lclosest(:) = .false. + + if (isplpl) then + do concurrent(k = 1:nenc, lmask(k)) + i = self%index1(k) + j = self%index2(k) + xr(:) = pl%rh(:, i) - pl%rh(:, j) + 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)) + end do + else + do concurrent(k = 1:nenc, lmask(k)) + i = self%index1(k) + 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)) + end do + end if - lany_collision = any(lcollision(:)) - - if (lany_collision) 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 - 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%v1(:,k) = pl%vb(:,i) - if (isplpl) then - self%x2(:,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 collisional colliders%idx - if (pl%lcollision(i) .or. pl%lcollision(j)) call pl%make_colliders([i,j]) - - ! Set the collision flag for these to bodies to true in case they become involved in another collision later in the step - pl%lcollision([i, j]) = .true. - pl%status([i, j]) = COLLISION - call pl%info(i)%set_value(status="COLLISION", discard_time=t, discard_rh=pl%rh(:,i), discard_vh=pl%vh(:,i)) - 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%v2(:,k) = tp%vb(:,j) - if (lcollision(k)) then - tp%status(j) = DISCARDED_PLR - tp%ldiscard(j) = .true. - write(idstri, *) pl%id(i) - write(idstrj, *) tp%id(j) - write(timestr, *) t - call tp%info(j)%set_value(status="DISCARDED_PLR", discard_time=t, discard_rh=tp%rh(:,j), discard_vh=tp%vh(:,j)) - write(message, *) "Particle " // trim(adjustl(tp%info(j)%name)) // " (" // trim(adjustl(idstrj)) // ")" & - // " collided with massive body " // trim(adjustl(pl%info(i)%name)) // " (" // trim(adjustl(idstri)) // ")" & - // " at t = " // trim(adjustl(timestr)) - call io_log_one_message(FRAGGLE_LOG_OUT, message) + lany_collision = any(lcollision(:)) + + if (lany_collision) 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 + 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%v1(:,k) = pl%vb(:,i) + if (isplpl) then + self%x2(:,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 + if (pl%lcollision(i) .or. pl%lcollision(j)) call pl%make_colliders([i,j]) + + ! Set the collision flag for these to bodies to true in case they become involved in another collision later in the step + pl%lcollision([i, j]) = .true. + pl%status([i, j]) = COLLISION + call pl%info(i)%set_value(status="COLLISION", discard_time=t, discard_rh=pl%rh(:,i), discard_vh=pl%vh(:,i)) + 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%v2(:,k) = tp%vb(:,j) + if (lcollision(k)) then + tp%status(j) = DISCARDED_PLR + tp%ldiscard(j) = .true. + write(idstri, *) pl%id(i) + write(idstrj, *) tp%id(j) + write(timestr, *) t + call tp%info(j)%set_value(status="DISCARDED_PLR", discard_time=t, discard_rh=tp%rh(:,j), discard_vh=tp%vh(:,j)) + write(message, *) "Particle " // trim(adjustl(tp%info(j)%name)) // " (" // trim(adjustl(idstrj)) // ")" & + // " collided with massive body " // trim(adjustl(pl%info(i)%name)) // " (" // trim(adjustl(idstri)) // ")" & + // " at t = " // trim(adjustl(timestr)) + call io_log_one_message(FRAGGLE_LOG_OUT, message) + end if end if - end if - end do + end do - ! Extract the pl-pl or pl-tp encounter list and return the pl-pl or pl-tp collision_list - select type(self) - class is (symba_plplenc) - call self%extract_collisions(system, param) - class is (symba_pltpenc) - allocate(tmp, mold=self) - call self%spill(tmp, lcollision, ldestructive=.true.) ! Remove this encounter pair from the encounter list - end select - end if + ! Extract the pl-pl or pl-tp encounter list and return the pl-pl or pl-tp collision_list + select type(self) + class is (symba_plplenc) + call self%extract_collisions(system, param) + class is (symba_pltpenc) + allocate(tmp, mold=self) + call self%spill(tmp, lcollision, ldestructive=.true.) ! Remove this encounter pair from the encounter list + end select + end if - ! Take snapshots of pairs of bodies at close approach (but not collision) if requested - if (any(lclosest(:))) then + ! 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") - end if + end select end select end select @@ -905,7 +907,7 @@ module subroutine symba_resolve_collision_fragmentations(self, system, param) call system%colliders%regime(system%fragments, system, param) - if (param%lenc_trajectory_save) call collision_history%take_snapshot(param,system, t, "before") + if (param%lenc_save_trajectory) call collision_history%take_snapshot(param,system, t, "before") select case (system%fragments%regime) case (COLLRESOLVE_REGIME_DISRUPTION, COLLRESOLVE_REGIME_SUPERCATASTROPHIC) plplcollision_list%status(i) = symba_collision_casedisruption(system, param) @@ -917,7 +919,7 @@ module subroutine symba_resolve_collision_fragmentations(self, system, param) write(*,*) "Error in symba_collision, unrecognized collision regime" call util_exit(FAILURE) end select - if (param%lenc_trajectory_save) call collision_history%take_snapshot(param,system, t, "after") + if (param%lenc_save_trajectory) call collision_history%take_snapshot(param,system, t, "after") deallocate(system%colliders,system%fragments) end do end select diff --git a/src/symba/symba_io.f90 b/src/symba/symba_io.f90 index 46fa3dfd5..916fb1a4c 100644 --- a/src/symba/symba_io.f90 +++ b/src/symba/symba_io.f90 @@ -125,8 +125,8 @@ module subroutine symba_io_param_reader(self, unit, iotype, v_list, iostat, ioms return end if - param%lenc_trajectory_save = (param%encounter_save == "TRAJECTORY") - param%lenc_closest_save = (param%encounter_save == "CLOSEST") .or. param%lenc_trajectory_save ! Closest approaches are always saved when trajectories are saved + param%lenc_save_trajectory = (param%encounter_save == "TRAJECTORY") + param%lenc_save_closest = (param%encounter_save == "CLOSEST") .or. param%lenc_save_trajectory ! Closest approaches are always saved when trajectories are saved ! Call the base method (which also prints the contents to screen) call io_param_reader(param, unit, iotype, v_list, iostat, iomsg) diff --git a/src/symba/symba_step.f90 b/src/symba/symba_step.f90 index 68548bdbe..3b217305f 100644 --- a/src/symba/symba_step.f90 +++ b/src/symba/symba_step.f90 @@ -38,9 +38,9 @@ module subroutine symba_step_system(self, param, t, dt) call self%reset(param) lencounter = pl%encounter_check(param, self, dt, 0) .or. tp%encounter_check(param, self, dt, 0) if (lencounter) then - if (param%lenc_trajectory_save) call encounter_history%take_snapshot(param, self, t) + if (param%lenc_save_trajectory) call encounter_history%take_snapshot(param, self, t, "trajectory") call self%interp(param, t, dt) - if (param%lenc_trajectory_save) call encounter_history%take_snapshot(param, self, t+dt) + if (param%lenc_save_trajectory) call encounter_history%take_snapshot(param, self, t+dt, "trajectory") else self%irec = -1 call helio_step_system(self, param, t, dt) @@ -247,7 +247,7 @@ recursive module subroutine symba_step_recur_system(self, param, t, ireci) if (lplpl_collision) call plplenc_list%resolve_collision(system, param, t+dtl, dtl, ireci) if (lpltp_collision) call pltpenc_list%resolve_collision(system, param, t+dtl, dtl, ireci) end if - if (param%lenc_trajectory_save) call encounter_history%take_snapshot(param, self, t+dtl) + if (param%lenc_save_trajectory) call encounter_history%take_snapshot(param, self, t+dtl, "trajectory") call self%set_recur_levels(ireci) diff --git a/src/symba/symba_util.f90 b/src/symba/symba_util.f90 index 8874be49b..74b555bce 100644 --- a/src/symba/symba_util.f90 +++ b/src/symba/symba_util.f90 @@ -204,6 +204,8 @@ module subroutine symba_util_dealloc_encounter_list(self) if (allocated(self%level)) deallocate(self%level) if (allocated(self%tcollision)) deallocate(self%tcollision) + call self%encounter_list%dealloc() + return end subroutine symba_util_dealloc_encounter_list @@ -232,7 +234,7 @@ module subroutine symba_util_dealloc_merger(self) if (allocated(self%ncomp)) deallocate(self%ncomp) - call symba_util_dealloc_pl(self) + call self%symba_pl%dealloc() return end subroutine symba_util_dealloc_merger @@ -266,7 +268,7 @@ module subroutine symba_util_dealloc_pl(self) deallocate(self%kin) end if - call util_dealloc_pl(self) + call self%helio_pl%dealloc() return end subroutine symba_util_dealloc_pl @@ -284,7 +286,7 @@ module subroutine symba_util_dealloc_tp(self) if (allocated(self%levelg)) deallocate(self%levelg) if (allocated(self%levelm)) deallocate(self%levelm) - call util_dealloc_tp(self) + call self%helio_tp%dealloc() return end subroutine symba_util_dealloc_tp @@ -712,6 +714,7 @@ module subroutine symba_util_rearray_pl(self, system, param) if ((idnew1 == idold1) .and. (idnew2 == idold2)) then ! This is an encounter we already know about, so save the old information 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) @@ -722,6 +725,7 @@ module subroutine symba_util_rearray_pl(self, system, param) else if (((idnew1 == idold2) .and. (idnew2 == idold1))) then ! This is an encounter we already know about, but with the order reversed, so save the old information 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) @@ -749,6 +753,7 @@ module subroutine symba_util_rearray_pl(self, system, param) system%plplenc_list%id1(1:nencmin) = pack(system%plplenc_list%id1(1:nenc_old), lmask(1:nenc_old)) system%plplenc_list%id2(1:nencmin) = pack(system%plplenc_list%id2(1:nenc_old), lmask(1:nenc_old)) system%plplenc_list%lvdotr(1:nencmin) = pack(system%plplenc_list%lvdotr(1:nenc_old), lmask(1:nenc_old)) + system%plplenc_list%lclosest(1:nencmin) = pack(system%plplenc_list%lclosest(1:nenc_old), lmask(1:nenc_old)) system%plplenc_list%status(1:nencmin) = pack(system%plplenc_list%status(1:nenc_old), lmask(1:nenc_old)) 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))