From c1ca5fc2c8bc91ff5215a199f2882f16265978d3 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Mon, 2 Aug 2021 18:35:34 -0400 Subject: [PATCH] Added the central body discard code from the Fragmentation branch --- src/discard/discard.f90 | 6 +-- src/main/swiftest_driver.f90 | 2 +- src/symba/symba_discard.f90 | 89 +++++++++++++++++++++++++++++++++++- src/util/util_set.f90 | 2 +- 4 files changed, 93 insertions(+), 6 deletions(-) diff --git a/src/discard/discard.f90 b/src/discard/discard.f90 index 0c84c9e88..5f87ee752 100644 --- a/src/discard/discard.f90 +++ b/src/discard/discard.f90 @@ -65,7 +65,7 @@ module subroutine discard_tp(self, system, param) if (ntp > 0) call tp%h2b(cb) end if if ((param%rmin >= 0.0_DP) .or. (param%rmax >= 0.0_DP) .or. (param%rmaxu >= 0.0_DP)) then - if (ntp > 0) call discard_sun_tp(tp, system, param) + if (ntp > 0) call discard_cb_tp(tp, system, param) end if if (param%qmin >= 0.0_DP .and. ntp > 0) call discard_peri_tp(tp, system, param) if (param%lclose .and. ntp > 0) call discard_pl_tp(tp, system, param) @@ -76,7 +76,7 @@ module subroutine discard_tp(self, system, param) end subroutine discard_tp - subroutine discard_sun_tp(tp, system, param) + subroutine discard_cb_tp(tp, system, param) !! author: David A. Minton !! !! Check to see if test particles should be discarded based on their positions relative to the Sun @@ -126,7 +126,7 @@ subroutine discard_sun_tp(tp, system, param) end associate return - end subroutine discard_sun_tp + end subroutine discard_cb_tp subroutine discard_peri_tp(tp, system, param) diff --git a/src/main/swiftest_driver.f90 b/src/main/swiftest_driver.f90 index 805264c2c..55eb1bc89 100644 --- a/src/main/swiftest_driver.f90 +++ b/src/main/swiftest_driver.f90 @@ -67,7 +67,7 @@ program swiftest_driver t = t0 + iloop * dt - !> Evaluate any discards or mergers + !> Evaluate any discards or collisional outcomes call nbody_system%discard(param) !> If the loop counter is at the output cadence value, append the data file with a single frame diff --git a/src/symba/symba_discard.f90 b/src/symba/symba_discard.f90 index 3f8ada6fe..9a90d7ea1 100644 --- a/src/symba/symba_discard.f90 +++ b/src/symba/symba_discard.f90 @@ -8,8 +8,95 @@ module subroutine symba_discard_pl(self, system, param) class(symba_pl), intent(inout) :: self !! SyMBA test particle object class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters - + + ! First check for collisions with the central body + associate(pl => self, npl => self%nbody, cb => system%cb) + if (npl == 0) return + if ((param%rmin >= 0.0_DP) .or. (param%rmax >= 0.0_DP) .or. & + (param%rmaxu >= 0.0_DP) .or. ((param%qmin >= 0.0_DP) .and. (param%qmin_coord == "BARY"))) then + call pl%h2b(cb) + end if + if ((param%rmin >= 0.0_DP) .or. (param%rmax >= 0.0_DP) .or. (param%rmaxu >= 0.0_DP)) then + call symba_discard_cb_pl(pl, system, param) + end if + if (param%qmin >= 0.0_DP .and. npl > 0) call symba_discard_peri_pl(pl, system, param) + end associate + + select type(param) + class is (symba_parameters) + if (param%lfragmentation) then + + end if + + end select + return end subroutine symba_discard_pl + + subroutine symba_discard_cb_pl(pl, system, param) + !! author: David A. Minton + !! + !! Check to see if planets should be discarded based on their positions relative to the central body + !! + !! Adapted from David E. Kaufmann's Swifter routine: symba_discard_cb.f90 + !! Adapted from Hal Levison's Swift routine discard_massive5.f + implicit none + ! Arguments + class(symba_pl), intent(inout) :: pl !! SyMBA massive body object + class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object + class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters + ! Internals + integer(I4B) :: i, j + real(DP) :: energy, vb2, rb2, rh2, rmin2, rmax2, rmaxu2 + + associate(npl => pl%nbody, cb => system%cb) + call system%set_msys() + rmin2 = param%rmin**2 + rmax2 = param%rmax*2 + rmaxu2 = param%rmaxu**2 + do i = 1, npl + if (pl%status(i) == ACTIVE) then + rh2 = dot_product(pl%xh(:,i), pl%xh(:,i)) + if ((param%rmax >= 0.0_DP) .and. (rh2 > rmax2)) then + pl%ldiscard(i) = .true. + pl%status(i) = DISCARDED_RMAX + write(*, *) "Massive body ", pl%id(i), " too far from the central body at t = ", param%t + else if ((param%rmin >= 0.0_DP) .and. (rh2 < rmin2)) then + pl%ldiscard(i) = .true. + pl%status(i) = DISCARDED_RMIN + write(*, *) "Massive body ", pl%id(i), " too close to the central body at t = ", param%t + else if (param%rmaxu >= 0.0_DP) then + rb2 = dot_product(pl%xb(:,i), pl%xb(:,i)) + vb2 = dot_product(pl%vb(:,i), pl%vb(:,i)) + energy = 0.5_DP * vb2 - system%msys / sqrt(rb2) + if ((energy > 0.0_DP) .and. (rb2 > rmaxu2)) then + pl%ldiscard(i) = .true. + pl%status(i) = DISCARDED_RMAXU + write(*, *) "Massive body ", pl%id(i), " is unbound and too far from barycenter at t = ", param%t + end if + end if + end if + end do + end associate + + return + end subroutine symba_discard_cb_pl + + + subroutine symba_discard_peri_pl(pl, system, param) + !! author: David A. Minton + !! + !! Check to see if a test particle should be discarded because its perihelion distance becomes too small + !! + !! Adapted from David E. Kaufmann's Swifter routine: discard_peri.f90 + !! Adapted from Hal Levison's Swift routine discard_peri.f + implicit none + ! Arguments + class(symba_pl), intent(inout) :: pl !! SyMBA massive body object + class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object + class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters + return + end subroutine symba_discard_peri_pl + end submodule s_symba_discard \ No newline at end of file diff --git a/src/util/util_set.f90 b/src/util/util_set.f90 index c401cb0ce..a1c4075b6 100644 --- a/src/util/util_set.f90 +++ b/src/util/util_set.f90 @@ -62,7 +62,7 @@ module subroutine util_set_msys(self) ! Arguments class(swiftest_nbody_system), intent(inout) :: self !! Swiftest nobdy system object - self%msys = self%cb%mass + sum(self%pl%mass(1:self%pl%nbody)) + self%msys = self%cb%Gmass + sum(self%pl%Gmass(1:self%pl%nbody), self%pl%status(1:self%pl%nbody) /= INACTIVE) return end subroutine util_set_msys