From d18ffb5fda8b22dbb4567876e23bc7f7e89ff99d Mon Sep 17 00:00:00 2001 From: David A Minton Date: Sat, 24 Jul 2021 15:24:14 -0400 Subject: [PATCH] Fixed encounter check and fleshing out recursive step subroutine --- src/modules/symba_classes.f90 | 8 +++ src/symba/symba_encounter_check.f90 | 93 +++++++++++++---------------- src/symba/symba_step.f90 | 28 ++++++--- src/symba/symba_util.f90 | 1 + 4 files changed, 71 insertions(+), 59 deletions(-) diff --git a/src/modules/symba_classes.f90 b/src/modules/symba_classes.f90 index 4edad0767..837394571 100644 --- a/src/modules/symba_classes.f90 +++ b/src/modules/symba_classes.f90 @@ -175,6 +175,14 @@ module subroutine symba_discard_tp(self, system, param) class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters end subroutine symba_discard_tp + module pure elemental subroutine symba_encounter_check_one(xr, yr, zr, vxr, vyr, vzr, rhill1, rhill2, dt, irec, lencounter, lvdotr) + implicit none + real(DP), intent(in) :: xr, yr, zr, vxr, vyr, vzr + real(DP), intent(in) :: rhill1, rhill2, dt + integer(I4B), intent(in) :: irec + logical, intent(out) :: lencounter, lvdotr + end subroutine symba_encounter_check_one + module function symba_encounter_check_pl(self, system, dt, irec) result(lany_encounter) implicit none class(symba_pl), intent(inout) :: self !! SyMBA test particle object diff --git a/src/symba/symba_encounter_check.f90 b/src/symba/symba_encounter_check.f90 index e188e57ac..d76a07af5 100644 --- a/src/symba/symba_encounter_check.f90 +++ b/src/symba/symba_encounter_check.f90 @@ -2,6 +2,10 @@ use swiftest contains module function symba_encounter_check_pl(self, system, dt, irec) result(lany_encounter) + !! author: David A. Minton + !! + !! Check for an encounter between massive bodies. + !! implicit none ! Arguments class(symba_pl), intent(inout) :: self !! SyMBA test particle object @@ -11,7 +15,6 @@ module function symba_encounter_check_pl(self, system, dt, irec) result(lany_enc ! Result logical :: lany_encounter !! Returns true if there is at least one close encounter ! Internals - real(DP) :: r2crit, vdotr, r2, v2, tmin, r2min, term2 integer(I4B) :: nenc_old integer(I8B) :: k real(DP), dimension(NDIM) :: xr, vr @@ -21,34 +24,11 @@ module function symba_encounter_check_pl(self, system, dt, irec) result(lany_enc allocate(lencounter(nplpl), loc_lvdotr(nplpl)) lencounter(:) = .false. - term2 = RHSCALE * (RSHELL**irec) - do k = 1, nplpl associate(i => pl%k_plpl(1, k), j => pl%k_plpl(2, k)) xr(:) = pl%xh(:, j) - pl%xh(:, i) - r2 = dot_product(xr(:), xr(:)) - r2crit = ((pl%rhill(i) + pl%rhill(i)) * term2)**2 vr(:) = pl%vh(:, j) - pl%vh(:, i) - vdotr = dot_product(vr(:), xr(:)) - if (r2 < r2crit) then - lencounter(k) = .true. - loc_lvdotr(k) = (vdotr < 0.0_DP) - else - if (vdotr < 0.0_DP) then - v2 = dot_product(vr(:), vr(:)) - tmin = -vdotr / v2 - if (tmin < dt) then - r2min = r2 - vdotr * vdotr / v2 - else - r2min = r2 + 2 * vdotr * dt + v2 * dt * dt - end if - r2min = min(r2min, r2) - if (r2min <= r2crit) then - lencounter(k) = .true. - loc_lvdotr(k) = (vdotr < 0.0_DP) - end if - end if - end if + call symba_encounter_check_one(xr(1), xr(2), xr(3), vr(1), vr(2), vr(3), pl%rhill(i), pl%rhill(j), dt, irec, lencounter(k), loc_lvdotr(k)) end associate end do @@ -71,6 +51,10 @@ module function symba_encounter_check_pl(self, system, dt, irec) result(lany_enc end function symba_encounter_check_pl module function symba_encounter_check_tp(self, system, dt, irec) result(lany_encounter) + !! author: David A. Minton + !! + !! Check for an encounter between test particles and massive bodies. + !! implicit none ! Arguments class(symba_tp), intent(inout) :: self !! SyMBA test particle object @@ -89,34 +73,11 @@ module function symba_encounter_check_tp(self, system, dt, irec) result(lany_enc allocate(lencounter(npl, ntp), loc_lvdotr(npl, ntp)) lencounter(:,:) = .false. - term2 = RHSCALE * (RSHELL**irec) - do j = 1, ntp do i = 1, npl xr(:) = tp%xh(:, j) - pl%xh(:, i) - r2 = dot_product(xr(:), xr(:)) - r2crit = (pl%rhill(i) * term2)**2 vr(:) = tp%vh(:, j) - pl%vh(:, i) - vdotr = dot_product(vr(:), xr(:)) - if (r2 < r2crit) then - lencounter(i,j) = .true. - loc_lvdotr(i,j) = (vdotr < 0.0_DP) - else - if (vdotr < 0.0_DP) then - v2 = dot_product(vr(:), vr(:)) - tmin = -vdotr / v2 - if (tmin < dt) then - r2min = r2 - vdotr * vdotr / v2 - else - r2min = r2 + 2 * vdotr * dt + v2 * dt * dt - end if - r2min = min(r2min, r2) - if (r2min <= r2crit) then - lencounter(i,j) = .true. - loc_lvdotr(i,j) = (vdotr < 0.0_DP) - end if - end if - end if + call symba_encounter_check_one(xr(1), xr(2), xr(3), vr(1), vr(2), vr(3), pl%rhill(i), 0.0_DP, dt, irec, lencounter(i,j), loc_lvdotr(i,j)) end do end do @@ -128,11 +89,8 @@ module function symba_encounter_check_tp(self, system, dt, irec) result(lany_enc pltpenc_list%status(nenc_old+1:nenc) = ACTIVE pltpenc_list%level(nenc_old+1:nenc) = irec pltpenc_list%lvdotr(nenc_old+1:nenc) = pack(loc_lvdotr(:,:), lencounter(:,:)) - !********************************************************************************************************* - ! This needs to be tested pltpenc_list%index1(nenc_old+1:nenc) = pack(spread([(i, i = 1, npl)], dim=2, ncopies=ntp), lencounter(:,:)) pltpenc_list%index2(nenc_old+1:nenc) = pack(spread([(j, j = 1, ntp)], dim=1, ncopies=npl), lencounter(:,:)) - !********************************************************************************************************* select type(pl) class is (symba_pl) pl%lencounter(pltpenc_list%index1(nenc_old+1:nenc)) = .true. @@ -143,4 +101,35 @@ module function symba_encounter_check_tp(self, system, dt, irec) result(lany_enc return end function symba_encounter_check_tp + module pure elemental subroutine symba_encounter_check_one(xr, yr, zr, vxr, vyr, vzr, rhill1, rhill2, dt, irec, lencounter, lvdotr) + !! author: David A. Minton + !! + !! Check for an encounter. + !! + !! Adapted from David E. Kaufmann's Swifter routine: symba_chk.f90 + !! Adapted from Hal Levison's Swift routine symba5_chk.f + implicit none + ! Arguments + real(DP), intent(in) :: xr, yr, zr, vxr, vyr, vzr + real(DP), intent(in) :: rhill1, rhill2, dt + integer(I4B), intent(in) :: irec + logical, intent(out) :: lencounter, lvdotr + ! Internals + integer(I4B) :: iflag + real(DP) :: r2, v2, rcrit, r2crit, vdotr + + lencounter = .false. + rcrit = (rhill1 + rhill2)*RHSCALE*(RSHELL**(irec)) + r2crit = rcrit**2 + r2 = xr**2 + yr**2 + zr**2 + v2 = vxr**2 + vyr**2 + vzr**2 + vdotr = xr * vxr + yr * vyr + zr * vzr + iflag = rmvs_chk_ind(r2, v2, vdotr, dt, r2crit) + if (iflag /= 0) lencounter = .true. + lvdotr = (vdotr < 0.0_DP) + + return + end subroutine symba_encounter_check_one + + end submodule s_symba_encounter_check \ No newline at end of file diff --git a/src/symba/symba_step.f90 b/src/symba/symba_step.f90 index 3f042fbd6..87cdab444 100644 --- a/src/symba/symba_step.f90 +++ b/src/symba/symba_step.f90 @@ -52,7 +52,6 @@ module subroutine symba_step_interp_system(self, param, t, dt) real(DP), intent(in) :: dt !! Current stepsize ! Internals real(DP) :: dth !! Half step size - integer(I4B) :: irec !! Recursion level dth = 0.5_DP * dt associate(system => self) @@ -76,8 +75,8 @@ module subroutine symba_step_interp_system(self, param, t, dt) call pl%drift(system, param, dt, pl%status(:) == ACTIVE) call tp%drift(system, param, dt, tp%status(:) == ACTIVE) - irec = 0 - call system%recursive_step(param, t, dt, irec) + + call system%recursive_step(param, t, dt, 0) call pl%set_beg_end(xend = pl%xh) call pl%accel(system, param, t + dt) @@ -113,10 +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, index_i, index_j, index_pl, index_tp - real(DP) :: dtl, dth,sgn + integer(I4B) :: i, j, irecp, icflg, nloops + real(DP) :: dtl, dth, sgn + real(DP), dimension(NDIM) :: xr, vr + logical :: lencounter - associate(plplenc_list => self%plplenc_list, pltpenc_list => self%pltpenc_list) + 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) dtl = param%dt / (NTENC**ireci) dth = 0.5_DP * dtl IF (dtl / param%dt < VSMALL) THEN @@ -127,8 +128,21 @@ module recursive subroutine symba_step_recur_system(self, param, t, dt, ireci) END IF irecp = ireci + 1 if (ireci == 0) then - icflg = 0 + nloops = 1 + else + 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 + + end do end associate end subroutine symba_step_recur_system diff --git a/src/symba/symba_util.f90 b/src/symba/symba_util.f90 index 81b351e65..c356d686b 100644 --- a/src/symba/symba_util.f90 +++ b/src/symba/symba_util.f90 @@ -59,6 +59,7 @@ module subroutine symba_util_resize_pltpenc(self, nrequested) allocate(enc_temp, source=self) call self%setup(2 * nrequested) call self%copy(enc_temp) + self%nenc = nrequested deallocate(enc_temp) return end subroutine symba_util_resize_pltpenc