diff --git a/docs/src/io.f90 b/docs/src/io.f90 index 42cc8ddd9..de1226951 100644 --- a/docs/src/io.f90 +++ b/docs/src/io.f90 @@ -1211,16 +1211,20 @@ module subroutine io_write_discard(self, param) associate(tp_discards => self%tp_discards, nsp => self%tp_discards%nbody, pl => self%pl, npl => self%pl%nbody) if (nsp == 0) return - select case(param%out_stat) - case('APPEND') + if (lfirst) then + select case(param%out_stat) + case('APPEND') + open(unit = LUN, file = param%discard_out, status = 'OLD', position = 'APPEND', form = 'FORMATTED', iostat = ierr) + case('NEW', 'REPLACE', 'UNKNOWN') + open(unit = LUN, file = param%discard_out, status = param%out_stat, form = 'FORMATTED', iostat = ierr) + case default + write(*,*) 'Invalid status code for OUT_STAT: ',trim(adjustl(param%out_stat)) + call util_exit(FAILURE) + end select + lfirst = .false. + else open(unit = LUN, file = param%discard_out, status = 'OLD', position = 'APPEND', form = 'FORMATTED', iostat = ierr) - case('NEW', 'REPLACE', 'UNKNOWN') - open(unit = LUN, file = param%discard_out, status = param%out_stat, form = 'FORMATTED', iostat = ierr) - case default - write(*,*) 'Invalid status code for OUT_STAT: ',trim(adjustl(param%out_stat)) - call util_exit(FAILURE) - end select - lfirst = .false. + end if if (param%lgr) call tp_discards%pv2v(param) write(LUN, HDRFMT) param%t, nsp, param%lbig_discard diff --git a/examples/rmvs_swifter_comparison/mars_ejecta/param.swifter.in b/examples/rmvs_swifter_comparison/mars_ejecta/param.swifter.in index f4035c4c0..3ffb6a2ee 100644 --- a/examples/rmvs_swifter_comparison/mars_ejecta/param.swifter.in +++ b/examples/rmvs_swifter_comparison/mars_ejecta/param.swifter.in @@ -14,7 +14,7 @@ IN_TYPE ASCII BIN_OUT bin.swifter.dat OUT_TYPE REAL8 ! double precision real output OUT_FORM XV ! osculating element output -OUT_STAT NEW +OUT_STAT UNKNOWN J2 0.0 ! no J2 term J4 0.0 ! no J4 term CHK_CLOSE yes ! check for planetary close encounters diff --git a/examples/rmvs_swifter_comparison/mars_ejecta/param.swiftest.in b/examples/rmvs_swifter_comparison/mars_ejecta/param.swiftest.in index 7df10c4d0..943e45ca3 100644 --- a/examples/rmvs_swifter_comparison/mars_ejecta/param.swiftest.in +++ b/examples/rmvs_swifter_comparison/mars_ejecta/param.swiftest.in @@ -24,3 +24,6 @@ GR no MU2KG 1.988409870698051e+30 DU2M 149597870700.0 TU2S 86400.0000 +RHILL_PRESENT no +DISCARD_OUT discard.swiftest.out +BIG_DISCARD yes ! output all planets if anything discarded diff --git a/src/discard/discard.f90 b/src/discard/discard.f90 index c71790c2d..3d65a235d 100644 --- a/src/discard/discard.f90 +++ b/src/discard/discard.f90 @@ -12,15 +12,27 @@ module subroutine discard_system(self, param) class(swiftest_nbody_system), intent(inout) :: self !! Swiftest system object class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters ! Internals - logical :: lany_discards + logical :: lany_discards, lpl_discards, ltp_discards, lpl_check, ltp_check - associate(system => self, tp => self%tp, pl => self%pl) - lany_discards = .false. - call pl%discard(system, param) - call tp%discard(system, param) - if (tp%nbody > 0) lany_discards = lany_discards .or. any(tp%ldiscard(:)) - if (pl%nbody > 0) lany_discards = lany_discards .or. any(pl%ldiscard(:)) - if (lany_discards) call system%write_discard(param) + lpl_check = allocated(self%pl_discards) + ltp_check = allocated(self%tp_discards) + + associate(system => self, tp => self%tp, pl => self%pl, tp_discards => self%tp_discards, pl_discards => self%tp_discards) + lpl_discards = .false. + ltp_discards = .false. + if (lpl_check) then + call pl%discard(system, param) + lpl_discards = (pl_discards%nbody > 0) + end if + + if (ltp_check) then + call tp%discard(system, param) + ltp_discards = (tp_discards%nbody > 0) + end if + + if (lpl_discards .or. ltp_discards) call system%write_discard(param) + if (lpl_check) call pl_discards%setup(0,param) + if (ltp_check) call tp_discards%setup(0,param) end associate return diff --git a/src/modules/swiftest_classes.f90 b/src/modules/swiftest_classes.f90 index 2455e77f2..25c18295a 100644 --- a/src/modules/swiftest_classes.f90 +++ b/src/modules/swiftest_classes.f90 @@ -275,6 +275,7 @@ module swiftest_classes class(swiftest_pl), allocatable :: pl !! Massive body data structure class(swiftest_tp), allocatable :: tp !! Test particle data structure class(swiftest_tp), allocatable :: tp_discards !! Discarded test particle data structure + class(swiftest_pl), allocatable :: pl_discards !! Discarded massive body particle data structure real(DP) :: Gmtot = 0.0_DP !! Total system mass - used for barycentric coordinate conversion real(DP) :: ke_orbit = 0.0_DP !! System orbital kinetic energy real(DP) :: ke_spin = 0.0_DP !! System spin kinetic energy diff --git a/src/util/util_spill.f90 b/src/util/util_spill.f90 index 9acc6ae93..9c76ff5e9 100644 --- a/src/util/util_spill.f90 +++ b/src/util/util_spill.f90 @@ -183,7 +183,7 @@ module subroutine util_spill_body(self, discards, lspill_list, ldestructive) ! 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 = count(.not.lspill_list(:)) + keeps%nbody = keeps%nbody - discards%nbody if (keeps%nbody > size(keeps%status)) keeps%status(keeps%nbody+1:size(keeps%status)) = INACTIVE end associate diff --git a/src/whm/whm_setup.f90 b/src/whm/whm_setup.f90 index cbf36cc90..eaed16c14 100644 --- a/src/whm/whm_setup.f90 +++ b/src/whm/whm_setup.f90 @@ -82,7 +82,7 @@ module subroutine whm_setup_initialize_system(self, param) call self%pl%sort("ir3h", ascending=.false.) ! Make sure that the discard list gets allocated initially - call self%tp_discards%setup(self%tp%nbody, param) + call self%tp_discards%setup(0, param) call self%pl%set_mu(self%cb) call self%tp%set_mu(self%cb) if (param%lgr) then