Skip to content
This repository was archived by the owner on Aug 28, 2024. It is now read-only.

Commit

Permalink
Fixed encounter check and fleshing out recursive step subroutine
Browse files Browse the repository at this point in the history
  • Loading branch information
@daminton
daminton committed Jul 24, 2021
1 parent eb8fbe0 commit d18ffb5
Show file tree
Hide file tree
Showing 4 changed files with 71 additions and 59 deletions.
8 changes: 8 additions & 0 deletions src/modules/symba_classes.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
93 changes: 41 additions & 52 deletions src/symba/symba_encounter_check.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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
Expand All @@ -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

Expand All @@ -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
Expand All @@ -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

Expand All @@ -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.
Expand All @@ -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
28 changes: 21 additions & 7 deletions src/symba/symba_step.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand All @@ -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)
Expand Down Expand Up @@ -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
Expand All @@ -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
Expand Down
1 change: 1 addition & 0 deletions src/symba/symba_util.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down

0 comments on commit d18ffb5

Please sign in to comment.