Fixed bugs that occur when number of massive bodies drops to 0 in SyMBA

src/helio/helio_kick.f90

@@ -74,10 +74,12 @@ module subroutine helio_kick_getacch_tp(self, nbody_system, param, t, lbeg)
 
       associate(tp => self, cb => nbody_system%cb, pl => nbody_system%pl, npl => nbody_system%pl%nbody)
          nbody_system%lbeg = lbeg
-         if (nbody_system%lbeg) then
-            call tp%accel_int(param, pl%Gmass(1:npl), pl%rbeg(:,1:npl), npl)
-         else
-            call tp%accel_int(param, pl%Gmass(1:npl), pl%rend(:,1:npl), npl)
+         if (npl > 0) then
+            if (nbody_system%lbeg) then
+               call tp%accel_int(param, pl%Gmass(1:npl), pl%rbeg(:,1:npl), npl)
+            else
+               call tp%accel_int(param, pl%Gmass(1:npl), pl%rend(:,1:npl), npl)
+            end if
          end if
          if (param%lnon_spherical_cb) call tp%accel_non_spherical_cb(nbody_system)
          if (param%lextra_force) call tp%accel_user(nbody_system, param, t, lbeg)

src/swiftest/swiftest_discard.f90

@@ -27,8 +27,8 @@ module subroutine swiftest_discard_system(self, param)
       class(swiftest_pl), allocatable :: plsub
       class(swiftest_tp), allocatable :: tpsub
 
-      lpl_check = allocated(self%pl_discards)
-      ltp_check = allocated(self%tp_discards)
+      lpl_check = allocated(self%pl_discards) .and. self%pl%nbody > 0
+      ltp_check = allocated(self%tp_discards) .and. self%tp%nbody > 0
 
       associate(nbody_system => self,tp => self%tp,pl => self%pl,tp_discards => self%tp_discards,pl_discards => self%pl_discards, &
                npl => self%pl%nbody, ntp => self%tp%nbody, t => self%t, collision_history => self%collision_history, &
@@ -38,14 +38,14 @@ module subroutine swiftest_discard_system(self, param)
          if (lpl_check .and. pl%nbody > 0) then
             pl%ldiscard = pl%status(:) /= ACTIVE
             call pl%discard(nbody_system, param)
-            lpl_discards = any(pl%ldiscard(1:npl))
+            if (npl > 0) lpl_discards = any(pl%ldiscard(1:npl))
          end if
 
          if (ltp_check .and. tp%nbody > 0) then
             tp%ldiscard = tp%status(:) /= ACTIVE
             call tp%discard(nbody_system, param)
-            ltp_discards = any(tp%ldiscard(1:ntp))
-            lpl_discards = any(pl%ldiscard(1:npl))
+            if (ntp > 0) ltp_discards = any(tp%ldiscard(1:ntp))
+            if (npl > 0) lpl_discards = any(pl%ldiscard(1:npl))
          end if
 
          if (ltp_discards.or.lpl_discards) then

src/swiftest/swiftest_obl.f90

@@ -313,6 +313,7 @@ module subroutine swiftest_obl_pot_system(self)
 
       associate(nbody_system => self, pl => self%pl, cb => self%cb)
          npl = self%pl%nbody
+         if (npl == 0) return
          if (.not. any(pl%lmask(1:npl))) return
 #ifdef DOCONLOC
          do concurrent (i = 1:npl, pl%lmask(i)) shared(cb,pl,oblpot_arr)

src/swiftest/swiftest_util.f90

@@ -1213,71 +1213,85 @@ module subroutine swiftest_util_get_energy_and_momentum_system(self, param)
          nbody_system%ke_orbit = 0.0_DP
          nbody_system%ke_spin = 0.0_DP
 
-         kepl(:) = 0.0_DP
-         Lplorbit(:,:) = 0.0_DP
-         Lplspin(:,:) = 0.0_DP
-
-         pl%lmask(1:npl) = pl%status(1:npl) /= INACTIVE
-
-         nbody_system%GMtot = cb%Gmass + sum(pl%Gmass(1:npl), pl%lmask(1:npl)) 
+         nbody_system%GMtot = cb%Gmass
+         if (npl > 0) then
+            kepl(:) = 0.0_DP
+            Lplorbit(:,:) = 0.0_DP
+            Lplspin(:,:) = 0.0_DP
+            pl%lmask(1:npl) = pl%status(1:npl) /= INACTIVE
+            nbody_system%GMtot = nbody_system%GMtot + sum(pl%Gmass(1:npl), pl%lmask(1:npl)) 
+         end if
+
          kecb = cb%mass * dot_product(cb%vb(:), cb%vb(:))
          nbody_system%be_cb = -3*cb%Gmass * cb%mass / (5 * cb%radius) 
          Lcborbit(:) = cb%mass * (cb%rb(:) .cross. cb%vb(:))
+         if (npl > 0) then
 
 #ifdef DOCONLOC
-         do concurrent (i = 1:npl, pl%lmask(i)) shared(pl,Lplorbit,kepl,npl) local(h) 
+            do concurrent (i = 1:npl, pl%lmask(i)) shared(pl,Lplorbit,kepl,npl) local(h) 
 #else
-         do concurrent (i = 1:npl, pl%lmask(i))
+            do concurrent (i = 1:npl, pl%lmask(i))
 #endif
-            h(1) = pl%rb(2,i) * pl%vb(3,i) - pl%rb(3,i) * pl%vb(2,i)
-            h(2) = pl%rb(3,i) * pl%vb(1,i) - pl%rb(1,i) * pl%vb(3,i)
-            h(3) = pl%rb(1,i) * pl%vb(2,i) - pl%rb(2,i) * pl%vb(1,i)
-
+               h(1) = pl%rb(2,i) * pl%vb(3,i) - pl%rb(3,i) * pl%vb(2,i)
+               h(2) = pl%rb(3,i) * pl%vb(1,i) - pl%rb(1,i) * pl%vb(3,i)
+               h(3) = pl%rb(1,i) * pl%vb(2,i) - pl%rb(2,i) * pl%vb(1,i)
+   
             ! Angular momentum from orbit 
-            Lplorbit(:,i) = pl%mass(i) * h(:)
+               Lplorbit(:,i) = pl%mass(i) * h(:)
 
-            ! Kinetic energy from orbit
-            kepl(i) = pl%mass(i) * dot_product(pl%vb(:,i), pl%vb(:,i)) 
-         end do
+               ! Kinetic energy from orbit
+               kepl(i) = pl%mass(i) * dot_product(pl%vb(:,i), pl%vb(:,i)) 
+            end do
+         end if
 
          if (param%lrotation) then
             kespincb = cb%mass * cb%Ip(3) * cb%radius**2 * dot_product(cb%rot(:), cb%rot(:))
 
             ! For simplicity, we always assume that the rotation pole is the 3rd principal axis
             Lcbspin(:) = cb%Ip(3) * cb%mass * cb%radius**2 * cb%rot(:)
 
+            if (npl > 0) then
 #ifdef DOCONLOC
-            do concurrent (i = 1:npl, pl%lmask(i)) shared(pl,Lplspin,kespinpl)
+               do concurrent (i = 1:npl, pl%lmask(i)) shared(pl,Lplspin,kespinpl)
 #else
-            do concurrent (i = 1:npl, pl%lmask(i))
+               do concurrent (i = 1:npl, pl%lmask(i))
 #endif
-               ! Currently we assume that the rotation pole is the 3rd principal axis
-               ! Angular momentum from spin
-               Lplspin(:,i) = pl%mass(i) * pl%Ip(3,i) * pl%radius(i)**2 * pl%rot(:,i)
+                  ! Currently we assume that the rotation pole is the 3rd principal axis
+                  ! Angular momentum from spin
+                  Lplspin(:,i) = pl%mass(i) * pl%Ip(3,i) * pl%radius(i)**2 * pl%rot(:,i)
 
-               ! Kinetic energy from spin
-               kespinpl(i) = pl%mass(i) * pl%Ip(3,i) * pl%radius(i)**2 * dot_product(pl%rot(:,i), pl%rot(:,i))
-            end do
+                  ! Kinetic energy from spin
+                  kespinpl(i) = pl%mass(i) * pl%Ip(3,i) * pl%radius(i)**2 * dot_product(pl%rot(:,i), pl%rot(:,i))
+               end do
 
-            nbody_system%ke_spin = 0.5_DP * (kespincb + sum(kespinpl(1:npl), pl%lmask(1:npl)))
+               nbody_system%ke_spin = 0.5_DP * (kespincb + sum(kespinpl(1:npl), pl%lmask(1:npl)))
+            else
+               nbody_system%ke_spin = 0.5_DP * kespincb
+            end if
 
+            if (npl > 0) then
 #ifdef DOCONLOC
-            do concurrent (j = 1:NDIM) shared(nbody_system,pl,Lplspin,Lcbspin)
+               do concurrent (j = 1:NDIM) shared(nbody_system,pl,Lplspin,Lcbspin)
 #else
-            do concurrent (j = 1:NDIM)
+               do concurrent (j = 1:NDIM)
 #endif
-               nbody_system%L_spin(j) = Lcbspin(j) + sum(Lplspin(j,1:npl), pl%lmask(1:npl))
-            end do
+                  nbody_system%L_spin(j) = Lcbspin(j) + sum(Lplspin(j,1:npl), pl%lmask(1:npl))
+               end do
+            else
+               nbody_system%L_spin(:) = Lcbspin(:)
+            end if
          else
             nbody_system%ke_spin = 0.0_DP
             nbody_system%L_spin(:) = 0.0_DP
          end if
-
-         if (param%lflatten_interactions) then
-            call swiftest_util_get_potential_energy(npl, pl%nplpl, pl%k_plpl, pl%lmask, cb%Gmass, pl%Gmass, pl%mass, pl%rb, &
-                                                    nbody_system%pe)
-         else
-            call swiftest_util_get_potential_energy(npl, pl%lmask, cb%Gmass, pl%Gmass, pl%mass, pl%rb, nbody_system%pe)
+
+         if (npl > 0) then
+            if (param%lflatten_interactions) then
+               call swiftest_util_get_potential_energy(npl, pl%nplpl, pl%k_plpl, pl%lmask, cb%Gmass, pl%Gmass, pl%mass, pl%rb, &
+                                                      nbody_system%pe)
+            else
+               call swiftest_util_get_potential_energy(npl, pl%lmask, cb%Gmass, pl%Gmass, pl%mass, pl%rb, nbody_system%pe)
+            end if
          end if
 
          ! Potential energy from the oblateness term
@@ -1286,16 +1300,21 @@ module subroutine swiftest_util_get_energy_and_momentum_system(self, param)
             nbody_system%pe = nbody_system%pe + nbody_system%oblpot
          end if
 
-         nbody_system%ke_orbit = 0.5_DP * (kecb + sum(kepl(1:npl), pl%lmask(1:npl)))
+         if (npl > 0) then
+            nbody_system%ke_orbit = 0.5_DP * (kecb + sum(kepl(1:npl), pl%lmask(1:npl)))
 #ifdef DOCONLOC
-         do concurrent (j = 1:NDIM) shared(nbody_system,pl,Lcborbit,Lplorbit,npl)
+            do concurrent (j = 1:NDIM) shared(nbody_system,pl,Lcborbit,Lplorbit,npl)
 #else  
-         do concurrent (j = 1:NDIM)
+            do concurrent (j = 1:NDIM)
 #endif
-            nbody_system%L_orbit(j) = Lcborbit(j) + sum(Lplorbit(j,1:npl), pl%lmask(1:npl)) 
-         end do
+               nbody_system%L_orbit(j) = Lcborbit(j) + sum(Lplorbit(j,1:npl), pl%lmask(1:npl)) 
+            end do
+         else
+            nbody_system%ke_orbit = 0.5_DP * kecb
+            nbody_system%L_orbit(:) = Lcborbit(:)
+         end if
 
-         if ((param%lclose)) then
+         if ((param%lclose .and. (npl > 0))) then
             nbody_system%be = sum(-3*pl%Gmass(1:npl)*pl%mass(1:npl)/(5*pl%radius(1:npl)), pl%lmask(1:npl)) 
          else
             nbody_system%be = 0.0_DP

src/symba/symba_discard.f90

@@ -48,11 +48,13 @@ subroutine symba_discard_cb_pl(pl, nbody_system, param)
                   write(message, *) trim(adjustl(pl%info(i)%name)) // " (" // trim(adjustl(idstr)) // ")" // &
                                     " too far from the central body at t = " // trim(adjustl(timestr))
                   call swiftest_io_log_one_message(COLLISION_LOG_OUT, "")
-                  call swiftest_io_log_one_message(COLLISION_LOG_OUT, "***********************************************************" // &
-                                                           "***********************************************************")
+                  call swiftest_io_log_one_message(COLLISION_LOG_OUT, &
+                                                   "***********************************************************" // &
+                                                   "***********************************************************")
                   call swiftest_io_log_one_message(COLLISION_LOG_OUT, message)
-                  call swiftest_io_log_one_message(COLLISION_LOG_OUT, "***********************************************************" // &
-                                                           "***********************************************************")
+                  call swiftest_io_log_one_message(COLLISION_LOG_OUT, & 
+                                                   "***********************************************************" // &
+                                                   "***********************************************************")
                   call swiftest_io_log_one_message(COLLISION_LOG_OUT, "")
                   call pl%info(i)%set_value(status="DISCARDED_RMAX", discard_time=nbody_system%t, discard_rh=pl%rh(:,i), &
                                             discard_vh=pl%vh(:,i))
@@ -65,11 +67,13 @@ subroutine symba_discard_cb_pl(pl, nbody_system, param)
                   write(message, *) trim(adjustl(pl%info(i)%name)) // " ("  // trim(adjustl(idstr)) // ")" // &
                                     " too close to the central body at t = " // trim(adjustl(timestr))
                   call swiftest_io_log_one_message(COLLISION_LOG_OUT, "")
-                  call swiftest_io_log_one_message(COLLISION_LOG_OUT, "************************************************************" // &
-                                                           "************************************************************")
+                  call swiftest_io_log_one_message(COLLISION_LOG_OUT, &
+                                                    "************************************************************" // &
+                                                    "************************************************************")
                   call swiftest_io_log_one_message(COLLISION_LOG_OUT, message)
-                  call swiftest_io_log_one_message(COLLISION_LOG_OUT, "************************************************************" // &
-                                                           "************************************************************")
+                  call swiftest_io_log_one_message(COLLISION_LOG_OUT, &
+                                                   "************************************************************" // &
+                                                   "************************************************************")
                   call swiftest_io_log_one_message(COLLISION_LOG_OUT, "")
                   call pl%info(i)%set_value(status="DISCARDED_RMIN", discard_time=nbody_system%t, discard_rh=pl%rh(:,i), &
                                             discard_vh=pl%vh(:,i), discard_body_id=cb%id)
@@ -86,11 +90,13 @@ subroutine symba_discard_cb_pl(pl, nbody_system, param)
                      write(message, *) trim(adjustl(pl%info(i)%name)) // " (" // trim(adjustl(idstr)) // ")" // &
                                        " is unbound and too far from barycenter at t = " // trim(adjustl(timestr))
                      call swiftest_io_log_one_message(COLLISION_LOG_OUT, "")
-                     call swiftest_io_log_one_message(COLLISION_LOG_OUT, "************************************************************" // &
-                                                              "************************************************************")
+                     call swiftest_io_log_one_message(COLLISION_LOG_OUT, &
+                                                      "************************************************************" // &
+                                                      "************************************************************")
                      call swiftest_io_log_one_message(COLLISION_LOG_OUT, message)
-                     call swiftest_io_log_one_message(COLLISION_LOG_OUT, "************************************************************" // &
-                                                              "************************************************************")
+                     call swiftest_io_log_one_message(COLLISION_LOG_OUT, &
+                                                      "************************************************************" // &
+                                                      "************************************************************")
                      call swiftest_io_log_one_message(COLLISION_LOG_OUT, "")
                      call pl%info(i)%set_value(status="DISCARDED_RMAXU", discard_time=nbody_system%t, discard_rh=pl%rh(:,i), &
                                                discard_vh=pl%vh(:,i))
@@ -233,16 +239,6 @@ subroutine symba_discard_nonplpl(pl, nbody_system, param)
             call symba_discard_cb_pl(pl, nbody_system, param)
          end if
          if (param%qmin >= 0.0_DP) call symba_discard_peri_pl(pl, nbody_system, param)
-         ! if (any(pl%ldiscard(1:npl))) then
-         !    ldiscard(1:npl) = pl%ldiscard(1:npl)
-
-         !    allocate(plsub, mold=pl)
-         !    call pl%spill(plsub, ldiscard, ldestructive=.false.)
-         !    nsub = plsub%nbody
-         !    nstart = pl_discards%nbody + 1
-         !    nend = pl_discards%nbody + nsub
-         !    call pl_discards%append(plsub, lsource_mask=[(.true., i = 1, nsub)])
-         ! end if
       end associate
 
       return
@@ -266,7 +262,8 @@ subroutine symba_discard_nonplpl_conservation(pl, nbody_system, param)
       integer(I4B), dimension(:), allocatable :: discard_index_list
 
       associate(npl => pl%nbody)
-         discard_l_pl(1:npl) = pl%ldiscard(1:npl) .and. .not. pl%lcollision(1:npl) ! These are bodies that are discarded but not flagged as pl-pl collision
+         discard_l_pl(1:npl) = pl%ldiscard(1:npl) .and. .not. pl%lcollision(1:npl) ! These are bodies that are discarded but not 
+                                                                                   ! flagged as pl-pl collision
          ndiscard = count(discard_l_pl(:)) 
          allocate(discard_index_list(ndiscard))
          discard_index_list(:) = pack([(i, i = 1, npl)], discard_l_pl(1:npl))
@@ -341,7 +338,8 @@ end subroutine symba_discard_peri_pl
    module subroutine symba_discard_pl(self, nbody_system, param)
       !! author: David A. Minton
       !!
-      !! Call the various flavors of discards for massive bodies in SyMBA runs, including discards due to colliding with the central body or escaping the nbody_system
+      !! Call the various flavors of discards for massive bodies in SyMBA runs, including discards due to colliding with the central
+      !! body or escaping the nbody_system
       implicit none
       ! Arguments
       class(symba_pl),              intent(inout) :: self   !! SyMBA test particle object

src/symba/symba_encounter_check.f90

@@ -256,7 +256,7 @@ module function symba_encounter_check_tp(self, param, nbody_system, dt, irec) re
       integer(I4B), dimension(:),   allocatable :: index1, index2
 
       lany_encounter = .false.
-      if (self%nbody == 0) return
+      if (self%nbody == 0 .or. nbody_system%pl%nbody == 0) return
 
       associate(tp => self, ntp => self%nbody, pl => nbody_system%pl, npl => nbody_system%pl%nbody, cb => nbody_system%cb)
          call pl%set_renc(irec)

src/symba/symba_util.f90

@@ -254,6 +254,7 @@ module subroutine symba_util_set_renc(self, scale)
       integer(I4B) :: i
       real(DP)     :: rshell_irec
 
+      if (self%nbody == 0) return
       associate(pl => self, npl => self%nbody)
          rshell_irec = 1._DP
          do i = 1, scale