diff --git a/src/io/io.f90 b/src/io/io.f90 index 637d1705f..e01efe4ec 100644 --- a/src/io/io.f90 +++ b/src/io/io.f90 @@ -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 @@ -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 @@ -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 @@ -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) diff --git a/src/symba/symba_collision.f90 b/src/symba/symba_collision.f90 index 922cadc09..ae63edcba 100644 --- a/src/symba/symba_collision.f90 +++ b/src/symba/symba_collision.f90 @@ -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 diff --git a/src/util/util_append.f90 b/src/util/util_append.f90 index ecefb1572..ffe576c85 100644 --- a/src/util/util_append.f90 +++ b/src/util/util_append.f90 @@ -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 @@ -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 diff --git a/src/util/util_spill.f90 b/src/util/util_spill.f90 index 8f9567ef9..16039c4da 100644 --- a/src/util/util_spill.f90 +++ b/src/util/util_spill.f90 @@ -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) @@ -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 @@ -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) @@ -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