From 43f6b4cd1e81cc1a45535ed1ce34296095dac61e Mon Sep 17 00:00:00 2001 From: David A Minton Date: Sat, 24 Jul 2021 16:00:00 -0400 Subject: [PATCH] Added encounter checks for the plpl and pltp encounter lists --- src/modules/symba_classes.f90 | 20 ++++++ src/symba/symba_encounter_check.f90 | 96 +++++++++++++++++++++++++++++ src/symba/symba_step.f90 | 14 +---- 3 files changed, 119 insertions(+), 11 deletions(-) diff --git a/src/modules/symba_classes.f90 b/src/modules/symba_classes.f90 index 837394571..877f5b585 100644 --- a/src/modules/symba_classes.f90 +++ b/src/modules/symba_classes.f90 @@ -121,6 +121,7 @@ module symba_classes integer(I4B), dimension(:), allocatable :: index1 !! position of the planet in encounter integer(I4B), dimension(:), allocatable :: index2 !! position of the test particle in encounter contains + procedure, public :: encounter_check => symba_encounter_check_pltpenc !! Checks if massive bodies are going through close encounters with each other procedure, public :: setup => symba_setup_pltpenc !! A constructor that sets the number of encounters and allocates and initializes all arrays procedure, public :: copy => symba_util_copy_pltpenc !! Copies all elements of one pltpenc list to another procedure, public :: resize => symba_util_resize_pltpenc !! Checks the current size of the pltpenc_list against the required size and extends it by a factor of 2 more than requested if it is too small @@ -136,6 +137,7 @@ module symba_classes real(DP), dimension(:,:), allocatable :: vb1 !! the barycentric velocity of parent 1 in encounter real(DP), dimension(:,:), allocatable :: vb2 !! the barycentric velocity of parent 2 in encounter contains + procedure, public :: encounter_check => symba_encounter_check_plplenc !! Checks if massive bodies are going through close encounters with each other procedure, public :: setup => symba_setup_plplenc !! A constructor that sets the number of encounters and allocates and initializes all arrays procedure, public :: copy => symba_util_copy_plplenc !! Copies all elements of one plplenc list to another end type symba_plplenc @@ -192,6 +194,24 @@ module function symba_encounter_check_pl(self, system, dt, irec) result(lany_enc logical :: lany_encounter !! Returns true if there is at least one close encounter end function symba_encounter_check_pl + module function symba_encounter_check_plplenc(self, system, dt, irec) result(lany_encounter) + implicit none + class(symba_plplenc), intent(inout) :: self !! SyMBA pl-pl encounter list object + class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object + real(DP), intent(in) :: dt !! step size + integer(I4B), intent(in) :: irec !! Current recursion level + logical :: lany_encounter !! Returns true if there is at least one close encounter + end function symba_encounter_check_plplenc + + module function symba_encounter_check_pltpenc(self, system, dt, irec) result(lany_encounter) + implicit none + class(symba_pltpenc), intent(inout) :: self !! SyMBA pl-pl encounter list object + class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object + real(DP), intent(in) :: dt !! step size + integer(I4B), intent(in) :: irec !! Current recursion level + logical :: lany_encounter !! Returns true if there is at least one close encounter + end function symba_encounter_check_pltpenc + module function symba_encounter_check_tp(self, system, dt, irec) result(lany_encounter) implicit none class(symba_tp), intent(inout) :: self !! SyMBA test particle object diff --git a/src/symba/symba_encounter_check.f90 b/src/symba/symba_encounter_check.f90 index d76a07af5..91feaafd1 100644 --- a/src/symba/symba_encounter_check.f90 +++ b/src/symba/symba_encounter_check.f90 @@ -50,6 +50,102 @@ module function symba_encounter_check_pl(self, system, dt, irec) result(lany_enc return end function symba_encounter_check_pl + module function symba_encounter_check_plplenc(self, system, dt, irec) result(lany_encounter) + !! author: David A. Minton + !! + !! Check for an encounter between massive bodies in the plplenc list. + !! + !! Adapted from portions of David E. Kaufmann's Swifter routine: symba_step_recur.f90 + implicit none + ! Arguments + class(symba_plplenc), intent(inout) :: self !! SyMBA pl-pl encounter list object + class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object + real(DP), intent(in) :: dt !! step size + integer(I4B), intent(in) :: irec !! Current recursion level + logical :: lany_encounter !! Returns true if there is at least one close encounter + ! Internals + integer(I4B) :: i + real(DP), dimension(NDIM) :: xr, vr + logical :: lencounter + real(DP) :: rlim2, rji2 + + lany_encounter = .false. + associate(plplenc_list => self, nplplenc => self%nenc) + select type(pl => system%pl) + class is (symba_pl) + do i = 1, nplplenc + associate(index_i => plplenc_list%index1(i), index_j => plplenc_list%index2(i)) + xr(:) = pl%xh(:,index_j) - pl%xh(:,index_i) + vr(:) = pl%vb(:,index_j) - pl%vb(:,index_i) + call symba_encounter_check_one(xr(1), xr(2), xr(3), vr(1), vr(2), vr(3), pl%rhill(index_i), pl%rhill(index_j), dt, irec, lencounter, plplenc_list%lvdotr(i)) + if (lencounter) then + rlim2 = (pl%radius(index_i) + pl%radius(index_j))**2 + rji2 = dot_product(xr(:), xr(:))! Check to see if these are physically overlapping bodies first, which we should ignore + if (rji2 > rlim2) then + lany_encounter = .true. + pl%levelg(index_i) = irec + pl%levelm(index_i) = MAX(irec, pl%levelm(index_i)) + pl%levelg(index_j) = irec + pl%levelm(index_j) = MAX(irec, pl%levelm(index_j)) + plplenc_list%level(i) = irec + end if + end if + end associate + end do + end select + end associate + return + end function symba_encounter_check_plplenc + + module function symba_encounter_check_pltpenc(self, system, dt, irec) result(lany_encounter) + !! author: David A. Minton + !! + !! Check for an encounter between test particles and massive bodies in the pltpenc list. + !! + !! Adapted from portions of David E. Kaufmann's Swifter routine: symba_step_recur.f90 + implicit none + ! Arguments + class(symba_pltpenc), intent(inout) :: self !! SyMBA pl-pl encounter list object + class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object + real(DP), intent(in) :: dt !! step size + integer(I4B), intent(in) :: irec !! Current recursion level + logical :: lany_encounter !! Returns true if there is at least one close encounter + ! Internals + integer(I4B) :: i + real(DP), dimension(NDIM) :: xr, vr + logical :: lencounter + real(DP) :: rlim2, rji2 + + lany_encounter = .false. + associate(pltpenc_list => self, npltpenc => self%nenc) + select type(pl => system%pl) + class is (symba_pl) + select type(tp => system%tp) + class is (symba_tp) + do i = 1, npltpenc + associate(index_i => pltpenc_list%index1(i), index_j => pltpenc_list%index2(i)) + xr(:) = tp%xh(:,index_j) - pl%xh(:,index_i) + vr(:) = tp%vb(:,index_j) - pl%vb(:,index_i) + call symba_encounter_check_one(xr(1), xr(2), xr(3), vr(1), vr(2), vr(3), pl%rhill(index_i), 0.0_DP, dt, irec, lencounter, pltpenc_list%lvdotr(i)) + if (lencounter) then + rlim2 = (pl%radius(index_i))**2 + rji2 = dot_product(xr(:), xr(:))! Check to see if these are physically overlapping bodies first, which we should ignore + if (rji2 > rlim2) then + lany_encounter = .true. + pl%levelg(index_i) = irec + pl%levelm(index_i) = MAX(irec, pl%levelm(index_i)) + tp%levelg(index_j) = irec + tp%levelm(index_j) = MAX(irec, tp%levelm(index_j)) + pltpenc_list%level(i) = irec + end if + end if + end associate + end do + end select + end select + end associate + end function symba_encounter_check_pltpenc + module function symba_encounter_check_tp(self, system, dt, irec) result(lany_encounter) !! author: David A. Minton !! diff --git a/src/symba/symba_step.f90 b/src/symba/symba_step.f90 index 87cdab444..cc7090065 100644 --- a/src/symba/symba_step.f90 +++ b/src/symba/symba_step.f90 @@ -112,12 +112,12 @@ module recursive subroutine symba_step_recur_system(self, param, t, dt, ireci) real(DP), intent(in) :: dt !! Current stepsize integer(I4B), value, intent(in) :: ireci !! input recursion level ! Internals - integer(I4B) :: i, j, irecp, icflg, nloops + integer(I4B) :: i, j, irecp, nloops real(DP) :: dtl, dth, sgn real(DP), dimension(NDIM) :: xr, vr logical :: lencounter - associate(pl => self%pl, tp => self%tp, plplenc_list => self%plplenc_list, nplplenc => self%plplenc_list%nenc, pltpenc_list => self%pltpenc_list, npltpenc => self%pltpenc_list%nenc) + associate(system => self, pl => self%pl, tp => self%tp, plplenc_list => self%plplenc_list, pltpenc_list => self%pltpenc_list) dtl = param%dt / (NTENC**ireci) dth = 0.5_DP * dtl IF (dtl / param%dt < VSMALL) THEN @@ -133,15 +133,7 @@ module recursive subroutine symba_step_recur_system(self, param, t, dt, ireci) nloops = NTENC end if do j = 1, nloops - icflg = 0 - do i = 1, nplplenc - associate(index_i => plplenc_list%index1(i), index_j => plplenc_list%index2(i)) - xr(:) = pl%xh(:,index_j) - pl%xh(:,index_i) - vr(:) = pl%vb(:,index_j) - pl%vb(:,index_i) - call symba_encounter_check_one(xr(1), xr(2), xr(3), vr(1), vr(2), vr(3), pl%rhill(index_i), pl%rhill(index_j), dtl, irecp, lencounter, plplenc_list%lvdotr(i)) - end associate - end do - + lencounter = plplenc_list%encounter_check(system, dtl, irecp) .or. pltpenc_list%encounter_check(system, dtl, irecp) end do end associate