From 72904e2aad485f8633843fa528d199689ac4fc6d Mon Sep 17 00:00:00 2001 From: David A Minton Date: Tue, 14 Dec 2021 12:36:21 -0500 Subject: [PATCH 1/2] Made the planetesimal mass in the symba_gr_test initial conditions generator a fraction of Mercury's mass --- examples/symba_gr_test/init_cond.py | 2 +- examples/symba_gr_test/pl.swiftest.in | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/examples/symba_gr_test/init_cond.py b/examples/symba_gr_test/init_cond.py index e59294190..a88ed26fb 100755 --- a/examples/symba_gr_test/init_cond.py +++ b/examples/symba_gr_test/init_cond.py @@ -61,7 +61,7 @@ GMcb = sim.ds.isel(id=0)['Gmass'].values GU = swiftest.GC / (sim.param['DU2M']**3 / (sim.param['MU2KG'] * sim.param['TU2S']**2)) dens = 3000.0 / (sim.param['MU2KG'] / sim.param['DU2M']**3) # Assume a bulk density of 3 g/cm^3 -GM_pl = 1e-11 +GM_pl = sim.ds.sel(id=1)['Gmass'].values * 0.01 M_pl = GM_pl / GU R_pl = (3 * M_pl / (4 * np.pi * dens))**(1.0 / 3.0) Rh_pl = Me_a * (GM_pl / (3 * GMcb))**(1.0/3.0) diff --git a/examples/symba_gr_test/pl.swiftest.in b/examples/symba_gr_test/pl.swiftest.in index cd36bda61..1e4f398e9 100644 --- a/examples/symba_gr_test/pl.swiftest.in +++ b/examples/symba_gr_test/pl.swiftest.in @@ -31,7 +31,7 @@ Neptune 0.0020336100526728302319 0.7749718408665498732 0.000164587904124493665 30.038959912561129073 0.008955570159296157365 1.7711193542961420899 131.82211597488270627 284.47484279411258967 308.45137222693909962 -Planetesimal 9.999999999999999395e-12 1.698243864025464137e-05 -2.2876635862715332271e-07 +Planetesimal 6.553709809565314101e-08 0.00031780626665496970825 +4.281092009908029071e-06 0.38709858430953941744 0.20562340108970009189 7.0033025080013837638 48.296118373786072198 29.20442403952453958 339.33948746828792764 From 2fbd36738fce58178f8fbe38f64bff0bc7e875a3 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Tue, 14 Dec 2021 14:03:01 -0500 Subject: [PATCH 2/2] Fixed bugs related to not properly removing pl-tp encounters from the list after a pl-tp collision takes place. --- src/symba/symba_collision.f90 | 90 ++++++++++++++++++----------------- src/symba/symba_kick.f90 | 9 ++-- 2 files changed, 53 insertions(+), 46 deletions(-) diff --git a/src/symba/symba_collision.f90 b/src/symba/symba_collision.f90 index 8b429a83d..fd0d569d1 100644 --- a/src/symba/symba_collision.f90 +++ b/src/symba/symba_collision.f90 @@ -275,7 +275,7 @@ module function symba_collision_check_encounter(self, system, param, t, dt, irec real(DP) :: rlim, Gmtot logical :: isplpl character(len=STRMAX) :: timestr, idstri, idstrj, message - + class(symba_encounter), allocatable :: tmp lany_collision = .false. if (self%nenc == 0) return @@ -324,54 +324,58 @@ module function symba_collision_check_encounter(self, system, param, t, dt, irec end do end if - if (any(lcollision(1:nenc))) call pl%xh2xb(system%cb) ! Update the central body barycenteric position vector to get us out of DH and into bary - do k = 1, nenc - i = self%index1(k) - j = self%index2(k) - if (lcollision(k)) self%status(k) = COLLISION - self%t(k) = t - self%x1(:,k) = pl%xh(:,i) + system%cb%xb(:) - self%v1(:,k) = pl%vb(:,i) - if (isplpl) then - self%x2(:,k) = pl%xh(:,j) + system%cb%xb(:) - self%v2(:,k) = pl%vb(:,j) - if (lcollision(k)) then - ! Check to see if either of these bodies has been involved with a collision before, and if so, make this a collisional colliders%idx - if (pl%lcollision(i) .or. pl%lcollision(j)) call pl%make_colliders([i,j]) - - ! Set the collision flag for these to bodies to true in case they become involved in another collision later in the step - pl%lcollision([i, j]) = .true. - pl%status([i, j]) = COLLISION - call pl%info(i)%set_value(status="COLLISION", discard_time=t, discard_xh=pl%xh(:,i), discard_vh=pl%vh(:,i)) - call pl%info(j)%set_value(status="COLLISION", discard_time=t, discard_xh=pl%xh(:,j), discard_vh=pl%vh(:,j)) - end if - else - self%x2(:,k) = tp%xh(:,j) + system%cb%xb(:) - self%v2(:,k) = tp%vb(:,j) - if (lcollision(k)) then - tp%status(j) = DISCARDED_PLR - tp%ldiscard(j) = .true. - write(idstri, *) pl%id(i) - write(idstrj, *) tp%id(j) - write(timestr, *) t - call tp%info(j)%set_value(status="DISCARDED_PLR", discard_time=t, discard_xh=tp%xh(:,j), discard_vh=tp%vh(:,j)) - write(message, *) "Particle " // trim(adjustl(tp%info(j)%name)) // " (" // trim(adjustl(idstrj)) // ")" & - // " collided with massive body " // trim(adjustl(pl%info(i)%name)) // " (" // trim(adjustl(idstri)) // ")" & - // " at t = " // trim(adjustl(timestr)) - call io_log_one_message(FRAGGLE_LOG_OUT, message) + lany_collision = any(lcollision(:)) + if (lany_collision) then + call pl%xh2xb(system%cb) ! Update the central body barycenteric position vector to get us out of DH and into bary + do k = 1, nenc + i = self%index1(k) + j = self%index2(k) + if (lcollision(k)) self%status(k) = COLLISION + self%t(k) = t + self%x1(:,k) = pl%xh(:,i) + system%cb%xb(:) + self%v1(:,k) = pl%vb(:,i) + if (isplpl) then + self%x2(:,k) = pl%xh(:,j) + system%cb%xb(:) + self%v2(:,k) = pl%vb(:,j) + if (lcollision(k)) then + ! Check to see if either of these bodies has been involved with a collision before, and if so, make this a collisional colliders%idx + if (pl%lcollision(i) .or. pl%lcollision(j)) call pl%make_colliders([i,j]) + + ! Set the collision flag for these to bodies to true in case they become involved in another collision later in the step + pl%lcollision([i, j]) = .true. + pl%status([i, j]) = COLLISION + call pl%info(i)%set_value(status="COLLISION", discard_time=t, discard_xh=pl%xh(:,i), discard_vh=pl%vh(:,i)) + call pl%info(j)%set_value(status="COLLISION", discard_time=t, discard_xh=pl%xh(:,j), discard_vh=pl%vh(:,j)) + end if + else + self%x2(:,k) = tp%xh(:,j) + system%cb%xb(:) + self%v2(:,k) = tp%vb(:,j) + if (lcollision(k)) then + tp%status(j) = DISCARDED_PLR + tp%ldiscard(j) = .true. + write(idstri, *) pl%id(i) + write(idstrj, *) tp%id(j) + write(timestr, *) t + call tp%info(j)%set_value(status="DISCARDED_PLR", discard_time=t, discard_xh=tp%xh(:,j), discard_vh=tp%vh(:,j)) + write(message, *) "Particle " // trim(adjustl(tp%info(j)%name)) // " (" // trim(adjustl(idstrj)) // ")" & + // " collided with massive body " // trim(adjustl(pl%info(i)%name)) // " (" // trim(adjustl(idstri)) // ")" & + // " at t = " // trim(adjustl(timestr)) + call io_log_one_message(FRAGGLE_LOG_OUT, message) + end if end if - end if - end do + end do + end if end select end select - lany_collision = any(lcollision(:)) - - ! Extract the pl-pl encounter list and return the plplcollision_list + ! Extract the pl-pl or pl-tpencounter list and return the plplcollision_list if (lany_collision) then - select type(plplenc_list => self) + select type(self) class is (symba_plplenc) - call plplenc_list%extract_collisions(system, param) + call self%extract_collisions(system, param) + class default + allocate(tmp, mold=self) + call self%spill(tmp, lcollision, ldestructive=.true.) ! Remove this encounter pair from the encounter list end select end if diff --git a/src/symba/symba_kick.f90 b/src/symba/symba_kick.f90 index fefad67e7..33dc6e804 100644 --- a/src/symba/symba_kick.f90 +++ b/src/symba/symba_kick.f90 @@ -163,7 +163,6 @@ module subroutine symba_kick_encounter(self, system, dt, irec, sgn) if (self%nenc == 0) return - select type(self) class is (symba_plplenc) isplpl = .true. @@ -175,8 +174,12 @@ module subroutine symba_kick_encounter(self, system, dt, irec, sgn) select type(tp => system%tp) class is (symba_tp) associate(npl => pl%nbody, ntp => tp%nbody, nenc => self%nenc) - if (npl > 0) pl%lmask(1:npl) = pl%status(1:npl) /= INACTIVE - if (ntp > 0) tp%lmask(1:ntp) = tp%status(1:ntp) /= INACTIVE + if (npl == 0) return + pl%lmask(1:npl) = pl%status(1:npl) /= INACTIVE + if (.not. isplpl) then + if (ntp == 0) return + tp%lmask(1:ntp) = tp%status(1:ntp) /= INACTIVE + end if allocate(lgoodlevel(nenc)) irm1 = irec - 1