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

Commit

Permalink
Changed the way discards are handeled and ensured that non-distructiv…
Browse files Browse the repository at this point in the history 
…e spill method doesn't change the number of bodies/encounters.
  • Loading branch information
@daminton
daminton committed Aug 30, 2021
1 parent d4e1776 commit 5e13952
Show file tree
Hide file tree
Showing 4 changed files with 68 additions and 50 deletions.
7 changes: 4 additions & 3 deletions src/io/io.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -67,14 +67,15 @@ module subroutine io_conservation_report(self, param, lterminal)
Ecoll_error = param%Ecollisions / abs(param%Eorbit_orig)
Etotal_error = (Eorbit_now - param%Ecollisions - param%Eorbit_orig - param%Euntracked) / abs(param%Eorbit_orig)
Merror = (GMtot_now - param%GMtot_orig) / param%GMtot_orig
write(*, EGYTERMFMT) Lerror, Ecoll_error, Etotal_error, Merror
if (Merror < -10 * epsilon(Merror)) then
write(*,*) 'Mass loss! Halting!'
call pl%xv2el(cb)
call param%nciu%open(param)
call pl%write_frame(param%nciu, param)
call param%nciu%close(param)
call util_exit(FAILURE)
end if
write(*, EGYTERMFMT) Lerror, Ecoll_error, Etotal_error, Merror
end if
end associate

Expand Down Expand Up @@ -1868,7 +1869,6 @@ module subroutine io_write_frame_system(self, param)
else
open(unit = iu, file = param%outfile, status = 'OLD', position = 'APPEND', form = 'UNFORMATTED', err = 667, iomsg = errmsg)
end if
call self%write_hdr(iu, param)
else if ((param%out_type == NETCDF_FLOAT_TYPE) .or. (param%out_type == NETCDF_DOUBLE_TYPE)) then

if (lfirst) then
Expand Down Expand Up @@ -1904,7 +1904,6 @@ module subroutine io_write_frame_system(self, param)
else
call param%nciu%open(param)
end if
call self%write_hdr(param%nciu, param)
end if

if (param%lgr) then
Expand All @@ -1919,11 +1918,13 @@ module subroutine io_write_frame_system(self, param)

! Write out each data type frame
if ((param%out_type == REAL4_TYPE) .or. (param%out_type == REAL8_TYPE)) then
call self%write_hdr(iu, param)
call cb%write_frame(iu, param)
call pl%write_frame(iu, param)
call tp%write_frame(iu, param)
close(iu, err = 667, iomsg = errmsg)
else if ((param%out_type == NETCDF_FLOAT_TYPE) .or. (param%out_type == NETCDF_DOUBLE_TYPE)) then
call self%write_hdr(param%nciu, param)
call cb%write_frame(param%nciu, param)
call pl%write_frame(param%nciu, param)
call tp%write_frame(param%nciu, param)
Expand Down
10 changes: 6 additions & 4 deletions src/symba/symba_collision.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -989,18 +989,20 @@ subroutine symba_collision_mergeaddsub(system, param, family, id_frag, Ip_frag,
call pl_adds%append(plnew, lsource_mask=[(.true., i=1, plnew%nbody)])
! Record how many bodies were added in this event
pl_adds%ncomp(nstart:nend) = plnew%nbody
call plnew%setup(0, param)

! Add the discarded bodies to the discard list
pl%status(family(:)) = MERGED
lmask(:) = .false.
lmask(family(:)) = .true.
pl%ldiscard(family(:)) = .true.
pl%lcollision(family(:)) = .true.
lmask(:) = .false.
lmask(family(:)) = .true.

call plnew%setup(0, param)
call pl%spill(plnew, lmask, ldestructive=.false.)

nstart = pl_discards%nbody + 1
nend = pl_discards%nbody + nfamily
call pl_discards%append(pl, lmask)
call pl_discards%append(plnew, lsource_mask=[(.true., i = 1, nfamily)])

! Record how many bodies were subtracted in this event
pl_discards%ncomp(nstart:nend) = nfamily
Expand Down
85 changes: 46 additions & 39 deletions src/util/util_append.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -187,31 +187,36 @@ module subroutine util_append_body(self, source, lsource_mask)
class(swiftest_body), intent(inout) :: self !! Swiftest body object
class(swiftest_body), intent(in) :: source !! Source object to append
logical, dimension(:), intent(in) :: lsource_mask !! Logical mask indicating which elements to append to
! Internals
integer(I4B) :: nold, nsrc, nnew

nold = self%nbody
nsrc = source%nbody
nnew = count(lsource_mask(1:nsrc))

associate(nold => self%nbody, nsrc => source%nbody)
call util_append(self%info, source%info, nold, nsrc, lsource_mask)
call util_append(self%id, source%id, nold, nsrc, lsource_mask)
call util_append(self%status, source%status, nold, nsrc, lsource_mask)
call util_append(self%ldiscard, source%ldiscard, nold, nsrc, lsource_mask)
call util_append(self%lmask, source%lmask, nold, nsrc, lsource_mask)
call util_append(self%mu, source%mu, nold, nsrc, lsource_mask)
call util_append(self%xh, source%xh, nold, nsrc, lsource_mask)
call util_append(self%vh, source%vh, nold, nsrc, lsource_mask)
call util_append(self%xb, source%xb, nold, nsrc, lsource_mask)
call util_append(self%vb, source%vb, nold, nsrc, lsource_mask)
call util_append(self%ah, source%ah, nold, nsrc, lsource_mask)
call util_append(self%aobl, source%aobl, nold, nsrc, lsource_mask)
call util_append(self%atide, source%atide, nold, nsrc, lsource_mask)
call util_append(self%agr, source%agr, nold, nsrc, lsource_mask)
call util_append(self%ir3h, source%ir3h, nold, nsrc, lsource_mask)
call util_append(self%a, source%a, nold, nsrc, lsource_mask)
call util_append(self%e, source%e, nold, nsrc, lsource_mask)
call util_append(self%inc, source%inc, nold, nsrc, lsource_mask)
call util_append(self%capom, source%capom, nold, nsrc, lsource_mask)
call util_append(self%omega, source%omega, nold, nsrc, lsource_mask)
call util_append(self%capm, source%capm, nold, nsrc, lsource_mask)
self%nbody = nold + count(lsource_mask(:))
end associate
call util_append(self%info, source%info, nold, nsrc, lsource_mask)
call util_append(self%id, source%id, nold, nsrc, lsource_mask)
call util_append(self%status, source%status, nold, nsrc, lsource_mask)
call util_append(self%ldiscard, source%ldiscard, nold, nsrc, lsource_mask)
call util_append(self%lmask, source%lmask, nold, nsrc, lsource_mask)
call util_append(self%mu, source%mu, nold, nsrc, lsource_mask)
call util_append(self%xh, source%xh, nold, nsrc, lsource_mask)
call util_append(self%vh, source%vh, nold, nsrc, lsource_mask)
call util_append(self%xb, source%xb, nold, nsrc, lsource_mask)
call util_append(self%vb, source%vb, nold, nsrc, lsource_mask)
call util_append(self%ah, source%ah, nold, nsrc, lsource_mask)
call util_append(self%aobl, source%aobl, nold, nsrc, lsource_mask)
call util_append(self%atide, source%atide, nold, nsrc, lsource_mask)
call util_append(self%agr, source%agr, nold, nsrc, lsource_mask)
call util_append(self%ir3h, source%ir3h, nold, nsrc, lsource_mask)
call util_append(self%a, source%a, nold, nsrc, lsource_mask)
call util_append(self%e, source%e, nold, nsrc, lsource_mask)
call util_append(self%inc, source%inc, nold, nsrc, lsource_mask)
call util_append(self%capom, source%capom, nold, nsrc, lsource_mask)
call util_append(self%omega, source%omega, nold, nsrc, lsource_mask)
call util_append(self%capm, source%capm, nold, nsrc, lsource_mask)

self%nbody = nold + nnew

return
end subroutine util_append_body
Expand All @@ -227,21 +232,23 @@ module subroutine util_append_encounter(self, source, lsource_mask)
class(swiftest_encounter), intent(inout) :: self !! Swiftest encounter list object
class(swiftest_encounter), intent(in) :: source !! Source object to append
logical, dimension(:), intent(in) :: lsource_mask !! Logical mask indicating which elements to append to

associate(nold => self%nenc, nsrc => source%nenc)
call util_append(self%lvdotr, source%lvdotr, nold, nsrc, lsource_mask)
call util_append(self%status, source%status, nold, nsrc, lsource_mask)
call util_append(self%index1, source%index1, nold, nsrc, 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%v1, source%v1, nold, nsrc, lsource_mask)
call util_append(self%v2, source%v2, nold, nsrc, lsource_mask)
call util_append(self%t, source%t, nold, nsrc, lsource_mask)
self%nenc = nold + count(lsource_mask(:))
end associate
! Internals
integer(I4B) :: nold, nsrc

nold = self%nenc
nsrc = source%nenc
call util_append(self%lvdotr, source%lvdotr, nold, nsrc, lsource_mask)
call util_append(self%status, source%status, nold, nsrc, lsource_mask)
call util_append(self%index1, source%index1, nold, nsrc, 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%v1, source%v1, nold, nsrc, lsource_mask)
call util_append(self%v2, source%v2, nold, nsrc, lsource_mask)
call util_append(self%t, source%t, nold, nsrc, lsource_mask)
self%nenc = nold + count(lsource_mask(1:nsrc))

return
end subroutine util_append_encounter
Expand Down
16 changes: 12 additions & 4 deletions src/util/util_spill.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -318,10 +318,12 @@ module subroutine util_spill_body(self, discards, lspill_list, ldestructive)
logical, dimension(:), intent(in) :: lspill_list !! Logical array of bodies to spill into the discards
logical, intent(in) :: ldestructive !! Logical flag indicating whether or not this operation should alter body by removing the discard list
! Internals
integer(I4B) :: nbody_old

! For each component, pack the discarded bodies into the discard object and do the inverse with the keeps
!> Spill all the common components
associate(keeps => self)

call util_spill(keeps%id, discards%id, lspill_list, ldestructive)
call util_spill(keeps%info, discards%info, lspill_list, ldestructive)
call util_spill(keeps%status, discards%status, lspill_list, ldestructive)
Expand All @@ -343,10 +345,12 @@ module subroutine util_spill_body(self, discards, lspill_list, ldestructive)
call util_spill(keeps%omega, discards%omega, lspill_list, ldestructive)
call util_spill(keeps%capm, discards%capm, lspill_list, ldestructive)

nbody_old = keeps%nbody

! This is the base class, so will be the last to be called in the cascade.
! Therefore we need to set the nbody values for both the keeps and discareds
discards%nbody = count(lspill_list(:))
keeps%nbody = keeps%nbody - discards%nbody
discards%nbody = count(lspill_list(1:nbody_old))
if (ldestructive) keeps%nbody = nbody_old- discards%nbody
end associate

return
Expand All @@ -363,6 +367,8 @@ module subroutine util_spill_encounter(self, discards, lspill_list, ldestructive
class(swiftest_encounter), intent(inout) :: discards !! Discarded object
logical, dimension(:), intent(in) :: lspill_list !! Logical array of bodies to spill into the discards
logical, intent(in) :: ldestructive !! Logical flag indicating whether or not this operation should alter body by removing the discard list
! Internals
integer(I4B) :: nenc_old

associate(keeps => self)
call util_spill(keeps%lvdotr, discards%lvdotr, lspill_list, ldestructive)
Expand All @@ -378,10 +384,12 @@ module subroutine util_spill_encounter(self, discards, lspill_list, ldestructive
call util_spill(keeps%v2, discards%v2, lspill_list, ldestructive)
call util_spill(keeps%t, discards%t, lspill_list, ldestructive)

nenc_old = keeps%nenc

! This is the base class, so will be the last to be called in the cascade.
! Therefore we need to set the nenc values for both the keeps and discareds
discards%nenc = count(lspill_list(:))
keeps%nenc = count(.not.lspill_list(:))
discards%nenc = count(lspill_list(1:nenc_old))
if (ldestructive) keeps%nenc = nenc_old - discards%nenc
end associate

return
Expand Down

0 comments on commit 5e13952

Please sign in to comment.