Skip to content
This repository was archived by the owner on Aug 28, 2024. It is now read-only.

Commit

Permalink
Fixed issues involving discard flags and bookkeeping of discards/adds
Browse files Browse the repository at this point in the history
  • Loading branch information
@daminton
daminton committed Sep 1, 2021
1 parent 3f292fa commit 9b5b3c3
Show file tree
Hide file tree
Showing 3 changed files with 38 additions and 40 deletions.
2 changes: 1 addition & 1 deletion src/modules/swiftest_classes.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
7 changes: 0 additions & 7 deletions src/symba/symba_collision.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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))
Expand Down Expand Up @@ -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)
Expand Down
69 changes: 37 additions & 32 deletions src/symba/symba_util.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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
Expand All @@ -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
Expand Down Expand Up @@ -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
Expand Down

0 comments on commit 9b5b3c3

Please sign in to comment.