Fixed discard particle number issue and added new checks for discards

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

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

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         

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

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

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

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