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

Commit

Permalink
Converted symba_kick to OOF
Browse files Browse the repository at this point in the history
  • Loading branch information
@daminton
daminton committed Jul 25, 2021
1 parent 9fd453f commit 2dacd1b
Show file tree
Hide file tree
Showing 2 changed files with 95 additions and 39 deletions.
19 changes: 0 additions & 19 deletions src/modules/symba_classes.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -138,7 +138,6 @@ module symba_classes
real(DP), dimension(:,:), allocatable :: vb1 !! the barycentric velocity of parent 1 in encounter
real(DP), dimension(:,:), allocatable :: vb2 !! the barycentric velocity of parent 2 in encounter
contains
procedure, public :: kick => symba_kick_plplenc !! Kick barycentric velocities of massive bodies within SyMBA recursion
procedure, public :: setup => symba_setup_plplenc !! A constructor that sets the number of encounters and allocates and initializes all arrays
procedure, public :: copy => symba_util_copy_plplenc !! Copies all elements of one plplenc list to another
end type symba_plplenc
Expand Down Expand Up @@ -195,15 +194,6 @@ module function symba_encounter_check_pl(self, system, dt, irec) result(lany_enc
logical :: lany_encounter !! Returns true if there is at least one close encounter
end function symba_encounter_check_pl

module function symba_encounter_check_plplenc(self, system, dt, irec) result(lany_encounter)
implicit none
class(symba_plplenc), intent(inout) :: self !! SyMBA pl-pl encounter list object
class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object
real(DP), intent(in) :: dt !! step size
integer(I4B), intent(in) :: irec !! Current recursion level
logical :: lany_encounter !! Returns true if there is at least one close encounter
end function symba_encounter_check_plplenc

module function symba_encounter_check_pltpenc(self, system, dt, irec) result(lany_encounter)
implicit none
class(symba_pltpenc), intent(inout) :: self !! SyMBA pl-pl encounter list object
Expand All @@ -222,15 +212,6 @@ module function symba_encounter_check_tp(self, system, dt, irec) result(lany_enc
logical :: lany_encounter !! Returns true if there is at least one close encounter
end function symba_encounter_check_tp

module subroutine symba_kick_plplenc(self, system, dt, irec, sgn)
implicit none
class(symba_plplenc), intent(in) :: self !! SyMBA pl-pl encounter list object
class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object
real(DP), intent(in) :: dt !! step size
integer(I4B), intent(in) :: irec !! Current recursion level
integer(I4B), intent(in) :: sgn !! sign to be applied to acceleration
end subroutine symba_kick_plplenc

module subroutine symba_kick_pltpenc(self, system, dt, irec, sgn)
implicit none
class(symba_pltpenc), intent(in) :: self !! SyMBA pl-tp encounter list object
Expand Down
115 changes: 95 additions & 20 deletions src/symba/symba_kick.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,34 +2,109 @@
use swiftest
contains

module subroutine symba_kick_plplenc(self, system, dt, irec, sgn)
!! author: David A. Minton
!!
!! !! Kick barycentric velocities of massive bodies within SyMBA recursion.
!!
!! Adapted from David E. Kaufmann's Swifter routine: symba_kick.f90
!! Adapted from Hal Levison's Swift routine symba5_kick.f
implicit none
class(symba_plplenc), intent(in) :: self !! SyMBA pl-pl encounter list object
class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object
real(DP), intent(in) :: dt !! step size
integer(I4B), intent(in) :: irec !! Current recursion level
integer(I4B), intent(in) :: sgn !! sign to be applied to acceleration
end subroutine symba_kick_plplenc

module subroutine symba_kick_pltpenc(self, system, dt, irec, sgn)
!! author: David A. Minton
!!
!! !! Kick barycentric velocities of active test particles within SyMBA recursion.
!! Kick barycentric velocities of massive bodies and ACTIVE test particles within SyMBA recursion.
!! Note: This method works for the polymorphic symba_pltpenc and symba_plplenc types
!!
!! Adapted from David E. Kaufmann's Swifter routine: symba_kick.f90
!! Adapted from Hal Levison's Swift routine symba5_kick.f
implicit none
class(symba_pltpenc), intent(in) :: self !! SyMBA pl-tp encounter list object
! Arguments
class(symba_pltpenc), intent(in) :: self !! SyMBA pl-tp encounter list object
class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object
real(DP), intent(in) :: dt !! step size
integer(I4B), intent(in) :: irec !! Current recursion level
integer(I4B), intent(in) :: sgn !! sign to be applied to acceleration
real(DP), intent(in) :: dt !! step size
integer(I4B), intent(in) :: irec !! Current recursion level
integer(I4B), intent(in) :: sgn !! sign to be applied to acceleration
! Internals
integer(I4B) :: i, irm1, irecl
real(DP) :: r, rr, ri, ris, rim1, r2, ir3, fac, faci, facj
real(DP), dimension(NDIM) :: dx
logical :: isplpl, lgoodlevel

select type(self)
class is (symba_plplenc)
isplpl = .true.
class is (symba_pltpenc)
isplpl = .false.
end select
select type(pl => system%pl)
class is (symba_pl)
select type(tp => system%tp)
class is (symba_tp)
irm1 = irec - 1
if (sgn < 0) then
irecl = irec - 1
else
irecl = irec
end if
do i = 1, self%nenc
associate(index_i => self%index1(i), index_j => self%index2(i))
if (isplpl) then
pl%ah(:,index_i) = 0.0_DP
pl%ah(:,index_j) = 0.0_DP
else
tp%ah(:,index_j) = 0.0_DP
end if
if (isplpl) then
lgoodlevel = (pl%levelg(index_i) >= irm1) .and. (pl%levelg(index_j) >= irm1)
else
lgoodlevel = (pl%levelg(index_i) >= irm1) .and. (tp%levelg(index_j) >= irm1)
end if
if ((self%status(i) == ACTIVE) .and. lgoodlevel) then
if (isplpl) then
ri = ((pl%rhill(index_i) + pl%rhill(index_j))**2) * (RHSCALE**2) * (RSHELL**(2*irecl))
rim1 = ri * (RSHELL**2)
dx(:) = pl%xh(:,index_j) - pl%xh(:,index_i)
else
ri = ((pl%rhill(index_i))**2) * (RHSCALE**2) * (RSHELL**(2*irecl))
rim1 = ri * (RSHELL**2)
dx(:) = tp%xh(:,index_j) - pl%xh(:,index_i)
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%mass(index_i)
if (isplpl) then
facj = fac * pl%mass(index_j)
pl%ah(:,index_i) = pl%ah(:,index_i) + facj*dx(:)
pl%ah(:,index_j) = pl%ah(:,index_j) - faci*dx(:)
else
tp%ah(:,index_j) = tp%ah(:,index_j) - faci*dx(:)
end if
end if
end associate
end do
if (isplpl) then
do i = 1, self%nenc
associate(index_i => self%index1(i), index_j => self%index2(i))
pl%vb(:,index_i) = pl%vb(:,index_i) + sgn * dt * pl%ah(:,index_i)
pl%vb(:,index_j) = pl%vb(:,index_j) + sgn * dt * pl%ah(:,index_j)
pl%ah(:,index_i) = 0.0_DP
pl%ah(:,index_j) = 0.0_DP
end associate
end do
else
where(tp%status(self%index2(1:self%nenc)) == ACTIVE)
tp%vb(1,self%index2(:)) = tp%vb(1,self%index2(:)) + sgn * dt * tp%ah(1,self%index2(:))
tp%vb(2,self%index2(:)) = tp%vb(2,self%index2(:)) + sgn * dt * tp%ah(2,self%index2(:))
tp%vb(3,self%index2(:)) = tp%vb(3,self%index2(:)) + sgn * dt * tp%ah(3,self%index2(:))
end where
tp%ah(:,self%index2(1:self%nenc)) = 0.0_DP
end if
end select
end select
return
end subroutine symba_kick_pltpenc

end submodule s_symba_kick

0 comments on commit 2dacd1b

Please sign in to comment.