diff --git a/src/helio/helio_kick.f90 b/src/helio/helio_kick.f90 index 6740e0eb0..7b7adee40 100644 --- a/src/helio/helio_kick.f90 +++ b/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) diff --git a/src/swiftest/swiftest_discard.f90 b/src/swiftest/swiftest_discard.f90 index e9a19e374..19f5cee27 100644 --- a/src/swiftest/swiftest_discard.f90 +++ b/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 diff --git a/src/swiftest/swiftest_obl.f90 b/src/swiftest/swiftest_obl.f90 index c2ecbb51f..68bd64dd0 100644 --- a/src/swiftest/swiftest_obl.f90 +++ b/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) diff --git a/src/swiftest/swiftest_util.f90 b/src/swiftest/swiftest_util.f90 index 469163b13..e1bfbb17b 100644 --- a/src/swiftest/swiftest_util.f90 +++ b/src/swiftest/swiftest_util.f90 @@ -1213,32 +1213,36 @@ 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(:)) @@ -1246,38 +1250,48 @@ module subroutine swiftest_util_get_energy_and_momentum_system(self, param) ! 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 diff --git a/src/symba/symba_discard.f90 b/src/symba/symba_discard.f90 index 7c95da30a..2d22ef299 100644 --- a/src/symba/symba_discard.f90 +++ b/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 diff --git a/src/symba/symba_encounter_check.f90 b/src/symba/symba_encounter_check.f90 index 9dec2072c..768349e7e 100644 --- a/src/symba/symba_encounter_check.f90 +++ b/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) diff --git a/src/symba/symba_util.f90 b/src/symba/symba_util.f90 index 38b13978a..e54a96f5c 100644 --- a/src/symba/symba_util.f90 +++ b/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