diff --git a/src/modules/swiftest_classes.f90 b/src/modules/swiftest_classes.f90 index faffcfb62..746781117 100644 --- a/src/modules/swiftest_classes.f90 +++ b/src/modules/swiftest_classes.f90 @@ -355,7 +355,7 @@ module swiftest_classes real(DP), dimension(:), allocatable :: peri !! Perihelion distance real(DP), dimension(:), allocatable :: atp !! Semimajor axis following perihelion passage integer(I4B), dimension(:,:), allocatable :: k_pltp !! Index array used to convert flattened the body-body comparison upper triangular matrix - integer(I8B) :: npltp !! Number of pl-tp comparisons in the flattened upper triangular matrix + integer(I8B) :: npltp !! Number of pl-tp comparisons in the flattened upper triangular matrix !! Note to developers: If you add components to this class, be sure to update methods and subroutines that traverse the !! component list, such as setup_tp and util_spill_tp contains diff --git a/src/symba/symba_collision.f90 b/src/symba/symba_collision.f90 index e1bbe02da..384102ffc 100644 --- a/src/symba/symba_collision.f90 +++ b/src/symba/symba_collision.f90 @@ -539,7 +539,6 @@ module function symba_collision_check_encounter(self, system, param, t, dt, irec ! 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%ldiscard([i, j]) = .true. pl%status([i, j]) = COLLISION call pl%info(i)%set_value(status="COLLISION", discard_time=t, discard_xh=pl%xh(:,i), discard_vh=pl%vh(:,i)) call pl%info(j)%set_value(status="COLLISION", discard_time=t, discard_xh=pl%xh(:,j), discard_vh=pl%vh(:,j)) @@ -922,15 +921,9 @@ subroutine symba_collision_mergeaddsub(system, param, family, id_frag, Ip_frag, plnew%vb(:, 1:nfrag) = vb_frag(:, 1:nfrag) call pl%vb2vh(cb) call pl%xh2xb(cb) - write(54,*) "Fragment properties" - write(54,*) "xbcb : ", cb%xb(:) - write(54,*) "vbcb : ", cb%vb(:) do i = 1, nfrag plnew%xh(:,i) = xb_frag(:, i) - cb%xb(:) plnew%vh(:,i) = vb_frag(:, i) - cb%vb(:) - write(54,*) "index, id: ", i, plnew%id(i) - write(54,*) "xb : ", xb_frag(:,i) - write(54,*) "vb : ", vb_frag(:,i) end do plnew%mass(1:nfrag) = m_frag(1:nfrag) plnew%Gmass(1:nfrag) = param%GU * m_frag(1:nfrag) diff --git a/src/symba/symba_util.f90 b/src/symba/symba_util.f90 index 91f9549ab..4a40100c0 100644 --- a/src/symba/symba_util.f90 +++ b/src/symba/symba_util.f90 @@ -409,7 +409,7 @@ module subroutine symba_util_rearray_pl(self, system, param) class(symba_parameters), intent(inout) :: param !! Current run configuration parameters ! Internals class(symba_pl), allocatable :: tmp !! The discarded body list. - integer(I4B) :: i, j, k, npl, nadd, nencmin, idnew1, idnew2, idold1, idold2 + integer(I4B) :: i, j, k, npl, nadd, nencmin, nenc_old, idnew1, idnew2, idold1, idold2 logical, dimension(:), allocatable :: lmask, ldump_mask class(symba_plplenc), allocatable :: plplenc_old logical :: lencounter @@ -426,7 +426,7 @@ module subroutine symba_util_rearray_pl(self, system, param) ! Remove the discards and destroy the list, as the system already tracks pl_discards elsewhere allocate(lmask(npl)) - lmask(1:npl) = pl%ldiscard(1:npl) .or. pl%status(1:npl) == INACTIVE + lmask(1:npl) = pl%ldiscard(1:npl) allocate(tmp, mold=self) call pl%spill(tmp, lspill_list=lmask, ldestructive=.true.) npl = pl%nbody @@ -435,8 +435,11 @@ module subroutine symba_util_rearray_pl(self, system, param) deallocate(lmask) ! 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) + nenc_old = system%plplenc_list%nenc + if (nenc_old > 0) then + allocate(plplenc_old, source=system%plplenc_list) + call plplenc_old%copy(system%plplenc_list) + end if ! Add in any new bodies if (nadd > 0) then @@ -502,34 +505,36 @@ module subroutine symba_util_rearray_pl(self, system, param) call move_alloc(nplenc_orig_pl, pl%nplenc) ! Re-index the encounter list as the index values may have changed - nencmin = min(system%plplenc_list%nenc, plplenc_old%nenc) - do k = 1, nencmin - idnew1 = system%plplenc_list%id1(k) - idnew2 = system%plplenc_list%id2(k) - idold1 = plplenc_old%id1(k) - idold2 = plplenc_old%id2(k) - 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%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 == 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%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 + if (nenc_old > 0) then + nencmin = min(system%plplenc_list%nenc, plplenc_old%nenc) + do k = 1, nencmin + idnew1 = system%plplenc_list%id1(k) + idnew2 = system%plplenc_list%id2(k) + idold1 = plplenc_old%id1(k) + idold2 = plplenc_old%id2(k) + 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%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 == 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%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 if end associate return