From 70bfb25a8639475970bceac6787518850721c625 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Wed, 18 Aug 2021 16:48:59 -0400 Subject: [PATCH] Fixed some indexing issues with SyMBA recursion levels and tried to make the kick routine cleaner --- src/symba/symba_kick.f90 | 138 ++++++++++++++++++++++++--------------- src/symba/symba_step.f90 | 4 +- 2 files changed, 86 insertions(+), 56 deletions(-) diff --git a/src/symba/symba_kick.f90 b/src/symba/symba_kick.f90 index 1eabf9b6d..b531e447a 100644 --- a/src/symba/symba_kick.f90 +++ b/src/symba/symba_kick.f90 @@ -109,10 +109,12 @@ module subroutine symba_kick_encounter(self, system, dt, irec, sgn) integer(I4B), intent(in) :: irec !! Current recursion level integer(I4B), intent(in) :: sgn !! sign to be applied to acceleration ! Internals - integer(I4B) :: k, irm1, irecl + integer(I4B) :: i, irm1, irecl, ngood real(DP) :: r, rr, ri, ris, rim1, r2, ir3, fac, faci, facj real(DP), dimension(NDIM) :: dx - logical :: isplpl, lgoodlevel + logical :: isplpl + logical, dimension(self%nenc) :: lgoodlevel + integer(I4B), dimension(:), allocatable :: good_idx if (self%nenc == 0) return @@ -138,63 +140,91 @@ module subroutine symba_kick_encounter(self, system, dt, irec, sgn) irecl = irec end if - if (isplpl) then - pl%ah(:,ind1(1:nenc)) = 0.0_DP - pl%ah(:,ind2(1:nenc)) = 0.0_DP - else - tp%ah(:,ind2(1:nenc)) = 0.0_DP - end if - - do k = 1, nenc + do i = 1, nenc if (isplpl) then - lgoodlevel = (pl%levelg(ind1(k)) >= irm1) .and. (pl%levelg(ind2(k)) >= irm1) + lgoodlevel(i) = (pl%levelg(ind1(i)) >= irm1) .and. (pl%levelg(ind2(i)) >= irm1) else - lgoodlevel = (pl%levelg(ind1(k)) >= irm1) .and. (tp%levelg(ind2(k)) >= irm1) - end if - if ((self%status(k) /= INACTIVE) .and. lgoodlevel) then - if (isplpl) then - ri = ((pl%rhill(ind1(k)) + pl%rhill(ind2(k)))**2) * (RHSCALE**2) * (RSHELL**(2*irecl)) - rim1 = ri * (RSHELL**2) - dx(:) = pl%xh(:,ind2(k)) - pl%xh(:,ind1(k)) - else - ri = ((pl%rhill(ind1(k)))**2) * (RHSCALE**2) * (RSHELL**(2*irecl)) - rim1 = ri * (RSHELL**2) - dx(:) = tp%xh(:,ind2(k)) - pl%xh(:,ind1(k)) - end if - r2 = dot_product(dx(:), dx(:)) - if (r2 < rim1) then - fac = 0.0_DP - else if (r2 < ri) then - ris = sqrt(ri) - r = sqrt(r2) - rr = (ris - r) / (ris * (1.0_DP - RSHELL)) - fac = (r2**(-1.5_DP)) * (1.0_DP - 3 * (rr**2) + 2 * (rr**3)) - else - ir3 = 1.0_DP / (r2 * sqrt(r2)) - fac = ir3 - end if - faci = fac * pl%Gmass(ind1(k)) - if (isplpl) then - facj = fac * pl%Gmass(ind2(k)) - pl%ah(:, ind1(k)) = pl%ah(:, ind1(k)) + facj * dx(:) - pl%ah(:, ind2(k)) = pl%ah(:, ind2(k)) - faci * dx(:) - else - tp%ah(:, ind2(k)) = tp%ah(:, ind2(k)) - faci * dx(:) - end if + lgoodlevel(i) = (pl%levelg(ind1(i)) >= irm1) .and. (tp%levelg(ind2(i)) >= irm1) end if + lgoodlevel(i) = (self%status(i) == ACTIVE) .and. lgoodlevel(i) end do - if (isplpl) then - do k = 1, nenc - pl%vb(:,ind1(k)) = pl%vb(:,ind1(k)) + sgn * dt * pl%ah(:,ind1(k)) - pl%vb(:,ind2(k)) = pl%vb(:,ind2(k)) + sgn * dt * pl%ah(:,ind2(k)) - pl%ah(:,ind1(k)) = 0.0_DP - pl%ah(:,ind2(k)) = 0.0_DP - end do - else - do k = 1, self%nenc - tp%vb(:,ind2(k)) = tp%vb(:,ind2(k)) + sgn * dt * tp%ah(:,ind2(k)) - tp%ah(:,ind2(k)) = 0.0_DP + ngood = count(lgoodlevel(:)) + if (ngood > 0) then + allocate(good_idx(ngood)) + good_idx(:) = pack([(i, i = 1, nenc)], lgoodlevel(:)) + + if (isplpl) then + do i = 1, ngood + associate(i1 => ind1(good_idx(i)), i2 => ind2(good_idx(i))) + pl%ah(:,i1) = 0.0_DP + pl%ah(:,i2) = 0.0_DP + end associate + end do + else + do i = 1, ngood + associate(i2 => ind2(good_idx(i))) + tp%ah(:,i2) = 0.0_DP + end associate + end do + end if + + do i = 1, ngood + associate(k => good_idx(i), i1 => ind1(good_idx(i)), i2 => ind2(good_idx(i))) + if (isplpl) then + ri = ((pl%rhill(i1) + pl%rhill(i2))**2) * (RHSCALE**2) * (RSHELL**(2*irecl)) + rim1 = ri * (RSHELL**2) + dx(:) = pl%xh(:,i2) - pl%xh(:,i1) + else + ri = ((pl%rhill(i1))**2) * (RHSCALE**2) * (RSHELL**(2*irecl)) + rim1 = ri * (RSHELL**2) + dx(:) = tp%xh(:,i2) - pl%xh(:,i1) + end if + r2 = dot_product(dx(:), dx(:)) + if (r2 < rim1) then + fac = 0.0_DP + lgoodlevel(k) = .false. + cycle + end if + if (r2 < ri) then + ris = sqrt(ri) + r = sqrt(r2) + rr = (ris - r) / (ris * (1.0_DP - RSHELL)) + fac = (r2**(-1.5_DP)) * (1.0_DP - 3 * (rr**2) + 2 * (rr**3)) + else + ir3 = 1.0_DP / (r2 * sqrt(r2)) + fac = ir3 + end if + faci = fac * pl%Gmass(i1) + if (isplpl) then + facj = fac * pl%Gmass(i2) + pl%ah(:, i1) = pl%ah(:, i1) + facj * dx(:) + pl%ah(:, i2) = pl%ah(:, i2) - faci * dx(:) + else + tp%ah(:, i2) = tp%ah(:, i2) - faci * dx(:) + end if + end associate end do + ngood = count(lgoodlevel(:)) + if (ngood == 0) return + good_idx(1:ngood) = pack([(i, i = 1, nenc)], lgoodlevel(:)) + + if (isplpl) then + do i = 1, ngood + associate(i1 => ind1(good_idx(i)), i2 => ind2(good_idx(i))) + pl%vb(:,i1) = pl%vb(:,i1) + sgn * dt * pl%ah(:,i1) + pl%vb(:,i2) = pl%vb(:,i2) + sgn * dt * pl%ah(:,i2) + pl%ah(:,ind1(1:nenc)) = 0.0_DP + pl%ah(:,ind2(1:nenc)) = 0.0_DP + end associate + end do + else + do i = 1, ngood + associate(i2 => ind2(good_idx(i))) + tp%vb(:,i2) = tp%vb(:,i2) + sgn * dt * tp%ah(:,i2) + tp%ah(:,i2) = 0.0_DP + end associate + end do + end if end if end associate end select diff --git a/src/symba/symba_step.f90 b/src/symba/symba_step.f90 index 249c73a82..1988b4fbc 100644 --- a/src/symba/symba_step.f90 +++ b/src/symba/symba_step.f90 @@ -243,8 +243,8 @@ module subroutine symba_step_reset_system(self, param) end do pl%nplenc(:) = 0 pl%ntpenc(:) = 0 - pl%levelg(:) = 0 - pl%levelm(:) = 0 + pl%levelg(:) = -1 + pl%levelm(:) = -1 pl%lencounter = .false. pl%lcollision = .false. system%plplenc_list%nenc = 0