From 9196c53e45dddeacda93a8d100ddeac178745ffe Mon Sep 17 00:00:00 2001 From: David A Minton Date: Tue, 31 Aug 2021 14:07:07 -0400 Subject: [PATCH] added more explicit setup/destruct calls for encounter lists and moved the kinship reset higher up in the order of operations, as part of the family consolidation step. --- src/orbel/orbel.f90 | 6 +++--- src/symba/symba_collision.f90 | 25 ++++++++++++++++--------- src/symba/symba_step.f90 | 12 +++++++++--- 3 files changed, 28 insertions(+), 15 deletions(-) diff --git a/src/orbel/orbel.f90 b/src/orbel/orbel.f90 index f01c97529..620e0f839 100644 --- a/src/orbel/orbel.f90 +++ b/src/orbel/orbel.f90 @@ -931,12 +931,12 @@ module pure subroutine orbel_xv2el(mu, x, v, a, e, inc, capom, omega, capm) capom = 0.0_DP omega = 0.0_DP capm = 0.0_DP - r = sqrt(dot_product(x(:), x(:))) + r = .mag. x(:) v2 = dot_product(v(:), v(:)) hvec = x(:) .cross. v(:) h2 = dot_product(hvec(:), hvec(:)) - h = sqrt(h2) - if (h2 == 0.0_DP) return + h = .mag. hvec(:) + if (h2 <= 10 * tiny(0.0_DP)) return rdotv = dot_product(x(:), v(:)) energy = 0.5_DP * v2 - mu / r fac = hvec(3) / h diff --git a/src/symba/symba_collision.f90 b/src/symba/symba_collision.f90 index ea93948f3..39d09970a 100644 --- a/src/symba/symba_collision.f90 +++ b/src/symba/symba_collision.f90 @@ -90,7 +90,6 @@ module function symba_collision_casedisruption(system, param, family, x, v, mass pl%status(family(:)) = status pl%ldiscard(family(:)) = .false. pl%lcollision(family(:)) = .false. - call pl%reset_kinship(family(:)) end select else ! Populate the list of new bodies @@ -202,7 +201,6 @@ module function symba_collision_casehitandrun(system, param, family, x, v, mass, pl%status(family(:)) = ACTIVE pl%ldiscard(family(:)) = .false. pl%lcollision(family(:)) = .false. - call pl%reset_kinship(family(:)) end select else status = HIT_AND_RUN_DISRUPT @@ -410,7 +408,6 @@ module function symba_collision_casesupercatastrophic(system, param, family, x, pl%status(family(:)) = status pl%ldiscard(family(:)) = .false. pl%lcollision(family(:)) = .false. - call pl%reset_kinship(family(:)) end select else ! Populate the list of new bodies @@ -657,6 +654,7 @@ function symba_collision_consolidate_familes(pl, cb, param, idx_parent, family, if (idx_parent(1) == idx_parent(2)) then if (nchild(1) == 0) then ! There is only one valid body recorded in this pair (this could happen due to restructuring of the kinship relationships, though it should be rare) lflag = .false. + call pl%reset_kinship([idx_parent(1)]) return end if idx_parent(2) = pl%kin(idx_parent(1))%child(nchild(1)) @@ -744,6 +742,9 @@ function symba_collision_consolidate_familes(pl, cb, param, idx_parent, family, end do lflag = .true. + ! Destroy the kinship relationships for all members of this family + call pl%reset_kinship(family(:)) + return end function symba_collision_consolidate_familes @@ -921,14 +922,21 @@ subroutine symba_collision_mergeaddsub(system, param, family, id_frag, Ip_frag, plnew%vb(:, 1:nfrag) = vb_frag(:, 1:nfrag) call pl%vb2vh(cb) call pl%xh2xb(cb) + ! write(54,*) "Fragment properties" + ! write(54,*) "xbcb : ", cb%xb(:) + ! write(54,*) "vbcb : ", cb%vb(:) do i = 1, nfrag plnew%xh(:,i) = xb_frag(:, i) - cb%xb(:) plnew%vh(:,i) = vb_frag(:, i) - cb%vb(:) + ! write(54,*) "index, id: ", i, plnew%id(i) + ! write(54,*) "xb : ", xb_frag(:,i) + ! write(54,*) "vb : ", vb_frag(:,i) end do plnew%mass(1:nfrag) = m_frag(1:nfrag) plnew%Gmass(1:nfrag) = param%GU * m_frag(1:nfrag) plnew%radius(1:nfrag) = rad_frag(1:nfrag) plnew%density(1:nfrag) = m_frag(1:nfrag) / rad_frag(1:nfrag) + call plnew%set_rhill(cb) select case(status) case(DISRUPTION) @@ -993,8 +1001,6 @@ subroutine symba_collision_mergeaddsub(system, param, family, id_frag, Ip_frag, plnew%tlag = pl%tlag(ibiggest) end if - call plnew%set_mu(cb) - call plnew%set_rhill(cb) !Copy over or set integration parameters for new bodies plnew%lcollision(1:nfrag) = .false. plnew%ldiscard(1:nfrag) = .false. @@ -1003,8 +1009,8 @@ subroutine symba_collision_mergeaddsub(system, param, family, id_frag, Ip_frag, ! Append the new merged body to the list nstart = pl_adds%nbody + 1 - nend = pl_adds%nbody + plnew%nbody - call pl_adds%append(plnew, lsource_mask=[(.true., i=1, plnew%nbody)]) + nend = pl_adds%nbody + nfrag + call pl_adds%append(plnew, lsource_mask=[(.true., i=1, nfrag)]) ! Record how many bodies were added in this event pl_adds%ncomp(nstart:nend) = plnew%nbody @@ -1012,11 +1018,11 @@ subroutine symba_collision_mergeaddsub(system, param, family, id_frag, Ip_frag, pl%status(family(:)) = MERGED pl%ldiscard(family(:)) = .true. pl%lcollision(family(:)) = .true. - call pl%reset_kinship(family(:)) lmask(:) = .false. lmask(family(:)) = .true. call plnew%setup(0, param) + deallocate(plnew) allocate(plsub, mold=pl) call pl%spill(plsub, lmask, ldestructive=.false.) @@ -1028,7 +1034,8 @@ subroutine symba_collision_mergeaddsub(system, param, family, id_frag, Ip_frag, ! Record how many bodies were subtracted in this event pl_discards%ncomp(nstart:nend) = nfamily - deallocate(plnew) + call plsub%setup(0, param) + deallocate(plsub) end associate end select end select diff --git a/src/symba/symba_step.f90 b/src/symba/symba_step.f90 index 7ab8113c0..9f4508550 100644 --- a/src/symba/symba_step.f90 +++ b/src/symba/symba_step.f90 @@ -226,7 +226,7 @@ module subroutine symba_step_reset_system(self, param) class(symba_nbody_system), intent(inout) :: self !! SyMBA nbody system object class(symba_parameters), intent(in) :: param !! Current run configuration parameters with SyMBA additions ! Internals - integer(I4B) :: i + integer(I4B) :: i, nenc_old associate(system => self) select type(pl => system%pl) @@ -245,8 +245,11 @@ module subroutine symba_step_reset_system(self, param) pl%lcollision(1:npl) = .false. pl%ldiscard(1:npl) = .false. pl%lmask(1:npl) = .true. + nenc_old = system%plplenc_list%nenc + call system%plplenc_list%setup(0) + call system%plplenc_list%setup(nenc_old) system%plplenc_list%nenc = 0 - system%plplcollision_list%nenc = 0 + call system%plplcollision_list%setup(0) end if if (ntp > 0) then @@ -254,7 +257,10 @@ module subroutine symba_step_reset_system(self, param) tp%levelg(1:ntp) = -1 tp%levelm(1:ntp) = -1 tp%lmask(1:ntp) = .true. - pl%ldiscard(1:npl) = .false. + tp%ldiscard(1:npl) = .false. + nenc_old = system%pltpenc_list%nenc + call system%pltpenc_list%setup(0) + call system%pltpenc_list%setup(nenc_old) system%pltpenc_list%nenc = 0 end if