From 3a61a4909c4ca1d10aadc11cc561fdd0c4d2a527 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Tue, 17 Aug 2021 15:04:05 -0400 Subject: [PATCH] Fixed problem in which kinship properties were not reset during a rearray. Also made rearray able to handle the case wher a body is discarded, but no new ones are added --- src/symba/symba_discard.f90 | 6 +- src/symba/symba_util.f90 | 114 ++++++++++++++++++++---------------- 2 files changed, 67 insertions(+), 53 deletions(-) diff --git a/src/symba/symba_discard.f90 b/src/symba/symba_discard.f90 index f6f7c0e6b..bf70f15dd 100644 --- a/src/symba/symba_discard.f90 +++ b/src/symba/symba_discard.f90 @@ -176,10 +176,6 @@ subroutine symba_discard_nonplpl(pl, system, param) ! First check for collisions with the central body associate(npl => pl%nbody, cb => system%cb) if (npl == 0) return - if ((param%rmin >= 0.0_DP) .or. (param%rmax >= 0.0_DP) .or. & - (param%rmaxu >= 0.0_DP) .or. ((param%qmin >= 0.0_DP) .and. (param%qmin_coord == "BARY"))) then - call pl%h2b(cb) - end if if ((param%rmin >= 0.0_DP) .or. (param%rmax >= 0.0_DP) .or. (param%rmaxu >= 0.0_DP)) then call symba_discard_cb_pl(pl, system, param) end if @@ -290,6 +286,8 @@ module subroutine symba_discard_pl(self, system, param) select type(param) class is (symba_parameters) associate(pl => self, plplenc_list => system%plplenc_list, plplcollision_list => system%plplcollision_list) + call pl%vb2vh(system%cb) + call pl%xh2xb(system%cb) call plplenc_list%write(pl, pl, param) call symba_discard_nonplpl(self, system, param) diff --git a/src/symba/symba_util.f90 b/src/symba/symba_util.f90 index e6ac2b6be..a80c8e646 100644 --- a/src/symba/symba_util.f90 +++ b/src/symba/symba_util.f90 @@ -424,25 +424,25 @@ module subroutine symba_util_rearray_pl(self, system, param) associate(pl => self, pl_adds => system%pl_adds) - allocate(tmp, mold=self) + ! Deallocate any temporary variables + if (allocated(pl%xbeg)) deallocate(pl%xbeg) + if (allocated(pl%xend)) deallocate(pl%xend) + ! Remove the discards and destroy the list, as the system already tracks pl_discards elsewhere allocate(lmask, source=pl%ldiscard(:)) lmask(:) = lmask(:) .or. pl%status(:) == INACTIVE + allocate(tmp, mold=self) call pl%spill(tmp, lspill_list=lmask, ldestructive=.true.) - deallocate(lmask) call tmp%setup(0,param) deallocate(tmp) + deallocate(lmask) - ! Deallocate any temporary variables - if (allocated(pl%xbeg)) deallocate(pl%xbeg) - if (allocated(pl%xend)) deallocate(pl%xend) + ! Store the original plplenc list so we don't remove any of the original encounters + allocate(plplenc_old, source=system%plplenc_list) + call plplenc_old%copy(system%plplenc_list) ! Add in any new bodies if (pl_adds%nbody > 0) then - ! First store the original plplenc list so we don't remove any of the original encounters - allocate(plplenc_old, source=system%plplenc_list) - call plplenc_old%copy(system%plplenc_list) - ! Append the adds to the main pl object call pl%append(pl_adds, lsource_mask=[(.true., i=1, pl_adds%nbody)]) @@ -450,49 +450,65 @@ module subroutine symba_util_rearray_pl(self, system, param) lmask(:) = pl%status(1:pl%nbody) == NEW_PARTICLE call symba_io_dump_particle_info(system, param, plidx=pack([(i, i=1, pl%nbody)], lmask)) - where(pl%status(:) /= INACTIVE) pl%status(:) = ACTIVE - call pl%sort("mass", ascending=.false.) - pl%lmtiny(:) = pl%Gmass(:) > param%GMTINY - pl%nplm = count(pl%lmtiny(:)) - - ! Reindex the new list ofbodies - call pl%eucl_index() - - ! Re-build the encounter list - lencounter = pl%encounter_check(system, param%dt, 0) - select type(tp => system%tp) - class is (symba_tp) - lencounter = tp%encounter_check(system, param%dt, 0) - end select - - associate(idnew1 => system%plplenc_list%id1, idnew2 => system%plplenc_list%id2, idold1 => plplenc_old%id1, idold2 => plplenc_old%id2) - do k = 1, system%plplenc_list%nenc - if ((idnew1(k) == idold1(k)) .and. (idnew2(k) == idold2(k))) 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%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%v1(:,k) = plplenc_old%v1(:,k) - system%plplenc_list%v2(:,k) = plplenc_old%v2(:,k) - system%plplenc_list%t(k) = plplenc_old%t(k) - system%plplenc_list%level(k) = plplenc_old%level(k) - else if (((idnew1(k) == idold2(k)) .and. (idnew2(k) == idold1(k)))) 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%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%v1(:,k) = plplenc_old%v2(:,k) - system%plplenc_list%v2(:,k) = plplenc_old%v1(:,k) - system%plplenc_list%t(k) = plplenc_old%t(k) - system%plplenc_list%level(k) = plplenc_old%level(k) - end if - end do - end associate + deallocate(lmask) end if + ! Reset all of the status flags for this body + where(pl%status(:) /= INACTIVE) + pl%status(:) = ACTIVE + pl%ldiscard(:) = .false. + pl%lcollision(:) = .false. + pl%lmtiny(:) = pl%Gmass(:) > param%GMTINY + pl%lmask(:) = .true. + elsewhere + pl%ldiscard(:) = .true. + pl%lmask(:) = .false. + end where + + pl%nplm = count(pl%lmtiny(1:pl%nbody) .and. pl%lmask(1:pl%nbody)) + + ! Reindex the new list of bodies + call pl%sort("mass", ascending=.false.) + call pl%eucl_index() + + ! Reset the kinship trackers + pl%kin(1:pl%nbody)%nchild = 0 + pl%kin(1:pl%nbody)%parent = [(i, i=1, pl%nbody)] + + ! Re-build the encounter list + lencounter = pl%encounter_check(system, param%dt, 0) + select type(tp => system%tp) + class is (symba_tp) + lencounter = tp%encounter_check(system, param%dt, 0) + end select + + associate(idnew1 => system%plplenc_list%id1, idnew2 => system%plplenc_list%id2, idold1 => plplenc_old%id1, idold2 => plplenc_old%id2) + do k = 1, system%plplenc_list%nenc + if ((idnew1(k) == idold1(k)) .and. (idnew2(k) == idold2(k))) 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%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%v1(:,k) = plplenc_old%v1(:,k) + system%plplenc_list%v2(:,k) = plplenc_old%v2(:,k) + system%plplenc_list%t(k) = plplenc_old%t(k) + system%plplenc_list%level(k) = plplenc_old%level(k) + else if (((idnew1(k) == idold2(k)) .and. (idnew2(k) == idold1(k)))) 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%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%v1(:,k) = plplenc_old%v2(:,k) + system%plplenc_list%v2(:,k) = plplenc_old%v1(:,k) + system%plplenc_list%t(k) = plplenc_old%t(k) + system%plplenc_list%level(k) = plplenc_old%level(k) + end if + end do + end associate + end associate return