From 8a3c43cf198daa7f0b3fccbdd19a412eab7e6a9c Mon Sep 17 00:00:00 2001 From: David A Minton Date: Wed, 4 Aug 2021 10:03:08 -0400 Subject: [PATCH] Refactored msys to system%Gmtot. Added back the massive body pericenter check (originally symba_peri.f90 and symba_discard_peri_pl.f90 from Swifter. Now symba_util_peri_pl bound to the get_peri method call, and symba_discard_peri_pl --- src/discard/discard.f90 | 4 +- src/helio/helio_coord.f90 | 6 +- src/modules/swiftest_classes.f90 | 2 +- src/modules/symba_classes.f90 | 9 ++ src/symba/symba_discard.f90 | 147 ++++++++++++++++++------------- src/symba/symba_util.f90 | 87 ++++++++++++++++++ src/util/util_coord.f90 | 10 +-- src/util/util_peri.f90 | 2 +- src/util/util_set.f90 | 2 +- 9 files changed, 197 insertions(+), 72 deletions(-) diff --git a/src/discard/discard.f90 b/src/discard/discard.f90 index f4244145a..7dbbd95c2 100644 --- a/src/discard/discard.f90 +++ b/src/discard/discard.f90 @@ -94,7 +94,7 @@ subroutine discard_cb_tp(tp, system, param) integer(I4B) :: i real(DP) :: energy, vb2, rb2, rh2, rmin2, rmax2, rmaxu2 - associate(ntp => tp%nbody, cb => system%cb, t => param%t, msys => system%msys) + associate(ntp => tp%nbody, cb => system%cb, t => param%t, Gmtot => system%Gmtot) rmin2 = max(param%rmin * param%rmin, cb%radius * cb%radius) rmax2 = param%rmax**2 rmaxu2 = param%rmaxu**2 @@ -114,7 +114,7 @@ subroutine discard_cb_tp(tp, system, param) else if (param%rmaxu >= 0.0_DP) then rb2 = dot_product(tp%xb(:, i), tp%xb(:, i)) vb2 = dot_product(tp%vb(:, i), tp%vb(:, i)) - energy = 0.5_DP * vb2 - msys / sqrt(rb2) + energy = 0.5_DP * vb2 - Gmtot / sqrt(rb2) if ((energy > 0.0_DP) .and. (rb2 > rmaxu2)) then tp%status(i) = DISCARDED_RMAXU write(*, *) "Particle ", tp%id(i), " is unbound and too far from barycenter at t = ", t diff --git a/src/helio/helio_coord.f90 b/src/helio/helio_coord.f90 index 0e58a3ab6..f40781810 100644 --- a/src/helio/helio_coord.f90 +++ b/src/helio/helio_coord.f90 @@ -68,14 +68,14 @@ module subroutine helio_coord_vh2vb_pl(self, cb) class(swiftest_cb), intent(inout) :: cb !! Swiftest central body object ! Internals integer(I4B) :: i - real(DP) :: msys + real(DP) :: Gmtot if (self%nbody == 0) return associate(pl => self, npl => self%nbody) - msys = cb%Gmass + sum(pl%Gmass(1:npl)) + Gmtot = cb%Gmass + sum(pl%Gmass(1:npl)) do i = 1, NDIM - cb%vb(i) = -sum(pl%Gmass(1:npl) * pl%vh(i, 1:npl)) / msys + cb%vb(i) = -sum(pl%Gmass(1:npl) * pl%vh(i, 1:npl)) / Gmtot pl%vb(i, 1:npl) = pl%vh(i, 1:npl) + cb%vb(i) end do end associate diff --git a/src/modules/swiftest_classes.f90 b/src/modules/swiftest_classes.f90 index fdb2360ee..c30947d68 100644 --- a/src/modules/swiftest_classes.f90 +++ b/src/modules/swiftest_classes.f90 @@ -259,7 +259,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 - real(DP) :: msys = 0.0_DP !! Total system mass - used for barycentric coordinate conversion + real(DP) :: Gmtot = 0.0_DP !! Total system mass - used for barycentric coordinate conversion real(DP) :: ke = 0.0_DP !! System kinetic energy real(DP) :: pe = 0.0_DP !! System potential energy real(DP) :: te = 0.0_DP !! System total energy diff --git a/src/modules/symba_classes.f90 b/src/modules/symba_classes.f90 index 92caa0bb3..c60b2a838 100644 --- a/src/modules/symba_classes.f90 +++ b/src/modules/symba_classes.f90 @@ -94,6 +94,7 @@ module symba_classes procedure :: setup => symba_setup_pl !! Constructor method - Allocates space for number of particle procedure :: append => symba_util_append_pl !! Appends elements from one structure to another procedure :: fill => symba_util_fill_pl !! "Fills" bodies from one object into another depending on the results of a mask (uses the UNPACK intrinsic) + procedure :: get_peri => symba_util_peri_pl !! Determine system pericenter passages for massive bodies procedure :: resize => symba_util_resize_pl !! Checks the current size of a Swiftest body against the requested size and resizes it if it is too small. procedure :: sort => symba_util_sort_pl !! Sorts body arrays by a sortable componen procedure :: rearrange => symba_util_sort_rearrange_pl !! Rearranges the order of array elements of body based on an input index array. Used in sorting methods @@ -479,6 +480,14 @@ module subroutine symba_util_fill_tp(self, inserts, lfill_list) class(swiftest_body), intent(in) :: inserts !! Inserted object logical, dimension(:), intent(in) :: lfill_list !! Logical array of bodies to merge into the keeps end subroutine symba_util_fill_tp + + module subroutine symba_util_peri_pl(self, system, param) + use swiftest_classes, only : swiftest_nbody_system, swiftest_parameters + implicit none + class(symba_pl), intent(inout) :: self !! 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 + end subroutine symba_util_peri_pl end interface interface util_resize diff --git a/src/symba/symba_discard.f90 b/src/symba/symba_discard.f90 index 7d15d996b..a84e158d3 100644 --- a/src/symba/symba_discard.f90 +++ b/src/symba/symba_discard.f90 @@ -2,62 +2,6 @@ use swiftest contains - module subroutine symba_discard_pl(self, system, param) - !! author: David A. Minton - !! - !! Call the various flavors of discards for massive bodies in SyMBA runs, including discards due to colling with the central body, - !! escaping the system, or colliding with each other. - implicit none - ! Arguments - 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 - - call symba_discard_nonplpl(self, system, param) - - select type(param) - class is (symba_parameters) - if (param%lfragmentation) then - - end if - - end select - - return - end subroutine symba_discard_pl - - subroutine symba_discard_nonplpl(pl, system, param) - !! author: David A. Minton - !! - !! Check to see if planets should be discarded based on their positions or because they are unbound - !s - !! - !! Adapted from David E. Kaufmann's Swifter routine: symba_discard_pl.f90 - !! Adapted from Hal Levison's Swift routine discard_massive5.f - implicit none - ! Arguments - class(symba_pl), intent(inout) :: pl !! 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(npl => pl%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 - - return - end subroutine symba_discard_nonplpl - - - subroutine symba_discard_cb_pl(pl, system, param) !! author: David A. Minton !! @@ -98,7 +42,7 @@ subroutine symba_discard_cb_pl(pl, system, param) 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) + energy = 0.5_DP * vb2 - system%Gmtot / sqrt(rb2) if ((energy > 0.0_DP) .and. (rb2 > rmaxu2)) then pl%ldiscard(i) = .true. pl%lcollision(i) = .false. @@ -114,19 +58,104 @@ subroutine symba_discard_cb_pl(pl, system, param) end subroutine symba_discard_cb_pl + subroutine symba_discard_nonplpl(pl, system, param) + !! author: David A. Minton + !! + !! Check to see if planets should be discarded based on their positions or because they are unbound + !s + !! + !! Adapted from David E. Kaufmann's Swifter routine: symba_discard_pl.f90 + !! Adapted from Hal Levison's Swift routine discard_massive5.f + implicit none + ! Arguments + class(symba_pl), intent(inout) :: pl !! 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(npl => pl%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 + + return + end subroutine symba_discard_nonplpl + + 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 + !! Adapted from David E. Kaufmann's Swifter routine: symba_discard_peri_pl.f90 + !! Adapted from Hal Levison's Swift routine discard_mass_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 + ! Internals + logical, save :: lfirst = .true. + logical :: lfirst_orig + integer(I4B) :: i + + + lfirst_orig = pl%lfirst + pl%lfirst = lfirst + if (lfirst) then + call pl%get_peri(system, param) + lfirst = .false. + else + call pl%get_peri(system, param) + do i = 1, pl%nbody + if (pl%status(i) == ACTIVE) then + if ((pl%isperi(i) == 0) .and. (pl%nplenc(i)== 0)) then + if ((pl%atp(i) >= param%qmin_alo) .and. (pl%atp(i) <= param%qmin_ahi) .and. (pl%peri(i) <= param%qmin)) then + pl%ldiscard(i) = .true. + pl%lcollision(i) = .false. + pl%status(i) = DISCARDED_PERI + write(*, *) "Particle ", pl%id(i), " perihelion distance too small at t = ", param%t + end if + end if + end if + end do + end if + pl%lfirst = lfirst_orig + return + end subroutine symba_discard_peri_pl + + module subroutine symba_discard_pl(self, system, param) + !! author: David A. Minton + !! + !! Call the various flavors of discards for massive bodies in SyMBA runs, including discards due to colling with the central body, + !! escaping the system, or colliding with each other. + implicit none + ! Arguments + 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 + + call symba_discard_nonplpl(self, system, param) + + select type(param) + class is (symba_parameters) + if (param%lfragmentation) then + + end if + + end select + + return + end subroutine symba_discard_pl + end submodule s_symba_discard \ No newline at end of file diff --git a/src/symba/symba_util.f90 b/src/symba/symba_util.f90 index feda27a07..7050bf050 100644 --- a/src/symba/symba_util.f90 +++ b/src/symba/symba_util.f90 @@ -271,6 +271,92 @@ module subroutine symba_util_fill_tp(self, inserts, lfill_list) end subroutine symba_util_fill_tp + module subroutine symba_util_peri_pl(self, system, param) + !! author: David A. Minton + !! + !! Determine system pericenter passages for planets in SyMBA + !! + !! Adapted from David E. Kaufmann's Swifter routine: symba_peri.f90 + !! Adapted from Hal Levison's Swift routine util_mass_peri.f + implicit none + ! Arguments + class(symba_pl), intent(inout) :: self !! 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 + real(DP) :: vdotr, e + + associate(pl => self, npl => self%nbody) + if (pl%lfirst) then + if (param%qmin_coord == "HELIO") then + do i = 1, npl + if (pl%status(i) == ACTIVE) then + vdotr = dot_product(pl%xh(:,i), pl%vh(:,i)) + if (vdotr > 0.0_DP) then + pl%isperi(i) = 1 + else + pl%isperi(i) = -1 + end if + end if + end do + else + do i = 1, npl + if (pl%status(i) == ACTIVE) then + vdotr = dot_product(pl%xb(:,i), pl%vb(:,i)) + if (vdotr > 0.0_DP) then + pl%isperi(i) = 1 + else + pl%isperi(i) = -1 + end if + end if + end do + end if + else + if (param%qmin_coord == "HELIO") then + do i = 1, npl + if (pl%status(i) == ACTIVE) then + vdotr = dot_product(pl%xh(:,i), pl%vh(:,i)) + if (pl%isperi(i) == -1) then + if (vdotr >= 0.0_DP) then + pl%isperi(i) = 0 + CALL orbel_xv2aeq(pl%mu(i), pl%xh(:,i), pl%vh(:,i), pl%atp(i), e, pl%peri(i)) + end if + else + if (vdotr > 0.0_DP) then + pl%isperi(i) = 1 + else + pl%isperi(i) = -1 + end if + end if + end if + end do + else + do i = 1, npl + if (pl%status(i) == ACTIVE) then + vdotr = dot_product(pl%xb(:,i), pl%vb(:,i)) + if (pl%isperi(i) == -1) then + if (vdotr >= 0.0_DP) then + pl%isperi(i) = 0 + CALL orbel_xv2aeq(system%Gmtot, pl%xb(:,i), pl%vb(:,i), pl%atp(i), e, pl%peri(i)) + end if + else + if (vdotr > 0.0_DP) then + pl%isperi(i) = 1 + else + pl%isperi(i) = -1 + end if + end if + end if + end do + end if + end if + end associate + + return + end subroutine symba_util_peri_pl + + module subroutine symba_util_resize_arr_info(arr, nnew) !! author: David A. Minton !! @@ -385,6 +471,7 @@ module subroutine symba_util_resize_tp(self, nnew) return end subroutine symba_util_resize_tp + module subroutine symba_util_sort_pl(self, sortby, ascending) !! author: David A. Minton !! diff --git a/src/util/util_coord.f90 b/src/util/util_coord.f90 index ab5fe85df..c10dbace7 100644 --- a/src/util/util_coord.f90 +++ b/src/util/util_coord.f90 @@ -15,21 +15,21 @@ module subroutine util_coord_h2b_pl(self, cb) class(swiftest_cb), intent(inout) :: cb !! Swiftest central body object ! Internals integer(I4B) :: i - real(DP) :: msys + real(DP) :: Gmtot real(DP), dimension(NDIM) :: xtmp, vtmp if (self%nbody == 0) return associate(pl => self, npl => self%nbody) - msys = cb%Gmass + Gmtot = cb%Gmass xtmp(:) = 0.0_DP vtmp(:) = 0.0_DP do i = 1, npl - msys = msys + pl%Gmass(i) + Gmtot = Gmtot + pl%Gmass(i) xtmp(:) = xtmp(:) + pl%Gmass(i) * pl%xh(:,i) vtmp(:) = vtmp(:) + pl%Gmass(i) * pl%vh(:,i) end do - cb%xb(:) = -xtmp(:) / msys - cb%vb(:) = -vtmp(:) / msys + cb%xb(:) = -xtmp(:) / Gmtot + cb%vb(:) = -vtmp(:) / Gmtot do i = 1, npl pl%xb(:,i) = pl%xh(:,i) + cb%xb(:) pl%vb(:,i) = pl%vh(:,i) + cb%vb(:) diff --git a/src/util/util_peri.f90 b/src/util/util_peri.f90 index 407ee5097..66f2254e1 100644 --- a/src/util/util_peri.f90 +++ b/src/util/util_peri.f90 @@ -45,7 +45,7 @@ module subroutine util_peri_tp(self, system, param) if (tp%isperi(i) == -1) then if (vdotr(i) >= 0.0_DP) then tp%isperi(i) = 0 - call orbel_xv2aeq(system%msys, tp%xb(:, i), tp%vb(:, i), tp%atp(i), e, tp%peri(i)) + call orbel_xv2aeq(system%Gmtot, tp%xb(:, i), tp%vb(:, i), tp%atp(i), e, tp%peri(i)) end if else if (vdotr(i) > 0.0_DP) then diff --git a/src/util/util_set.f90 b/src/util/util_set.f90 index a1c4075b6..86e021ab6 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%Gmass + sum(self%pl%Gmass(1:self%pl%nbody), self%pl%status(1:self%pl%nbody) /= INACTIVE) + self%Gmtot = 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