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

Commit

Permalink
Moved interaction acceleration to the main kick submodule, as the alg…
Browse files Browse the repository at this point in the history 
…orithm is shared by multiple integrators
  • Loading branch information
@daminton
daminton committed Jul 23, 2021
1 parent 8070dc5 commit 4efd05b
Show file tree
Hide file tree
Showing 6 changed files with 112 additions and 207 deletions.
4 changes: 2 additions & 2 deletions examples/whm_swifter_comparison/swiftest_vs_swifter.ipynb
View file

Large diffs are not rendered by default.

27 changes: 0 additions & 27 deletions src/eucl/eucl.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -20,9 +20,7 @@ module subroutine eucl_dist_index_plpl(self)
associate(nplpl => self%nplpl)
nplpl = (npl * (npl - 1) / 2) ! number of entries in a strict lower triangle, nplm x npl, minus first column
if (allocated(self%k_eucl)) deallocate(self%k_eucl) ! Reset the index array if it's been set previously
if (allocated(self%irij3)) deallocate(self%irij3)
allocate(self%k_eucl(2, nplpl))
allocate(self%irij3(nplpl))
do i = 1, npl
counter = (i - 1_I8B) * npl - i * (i - 1_I8B) / 2_I8B + 1_I8B
do j = i + 1_I8B, npl
Expand All @@ -49,29 +47,4 @@ module subroutine eucl_dist_index_pltp(self, pl)
class(swiftest_pl), intent(inout) :: pl !! Swiftest massive body object
end subroutine eucl_dist_index_pltp

module subroutine eucl_irij3_plpl(self)
!! author: Jacob R. Elliott and David A. Minton
!!
!! Efficient parallel loop-blocking algrorithm for evaluating the Euclidean distance matrix for planet-planet
implicit none
! Arguments
class(swiftest_pl), intent(inout) :: self !! Swiftest massive body object
! Internals
integer(I4B) :: k, i, j
real(DP), dimension(NDIM) :: dx
real(DP) :: rji2

associate(k_eucl => self%k_eucl, xh => self%xh, irij3 => self%irij3, nk => self%nplpl)
do k = 1, nk
i = k_eucl(1, k)
j = k_eucl(2, k)
dx(:) = xh(:, j) - xh(:, i)
rji2 = dot_product(dx(:), dx(:))
irij3(k) = 1.0_DP / (rji2 * sqrt(rji2))
end do
end associate

return
end subroutine eucl_irij3_plpl

end submodule s_eucl
91 changes: 8 additions & 83 deletions src/helio/helio_getacch.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -17,7 +17,8 @@ module subroutine helio_getacch_pl(self, system, param, t, lbeg)
logical, optional, intent(in) :: lbeg !! Optional argument that determines whether or not this is the beginning or end of the step

associate(cb => system%cb, pl => self, npl => self%nbody)
call helio_getacch_int_pl(pl, t)
pl%ah(:,:) = 0.0_DP
call pl%accel_int()
if (param%loblatecb) then
cb%aoblbeg = cb%aobl
call pl%accel_obl(system)
Expand Down Expand Up @@ -49,95 +50,19 @@ module subroutine helio_getacch_tp(self, system, param, t, lbeg)
class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters
real(DP), intent(in) :: t !! Current time
logical, optional, intent(in) :: lbeg !! Optional argument that determines whether or not this is the beginning or end of the step
! Internals
logical, save :: lmalloc = .true.
integer(I4B) :: i
real(DP) :: r2, mu
real(DP), dimension(:), allocatable, save :: irh, irht

associate(tp => self, ntp => self%nbody, cb => system%cb, npl => system%pl%nbody)
associate(tp => self, ntp => self%nbody, cb => system%cb, pl => system%pl, npl => system%pl%nbody)
if (present(lbeg)) system%lbeg = lbeg
call helio_getacch_int_tp(tp, system, param, t)
if (lbeg) then
call tp%accel_int(pl%Gmass(:), pl%xbeg(:,:), npl)
else
call tp%accel_int(pl%Gmass(:), pl%xend(:,:), npl)
end if
if (param%loblatecb) call tp%accel_obl(system)
if (param%lextra_force) call tp%accel_user(system, param, t)
!if (param%lgr) call tp%gr_accel(param)
end associate
return
end subroutine helio_getacch_tp

subroutine helio_getacch_int_pl(pl, t)
!! author: David A. Minton
!!
!! Compute direct cross term heliocentric accelerations of massive bodiese
!!
!! Adapted from David E. Kaufmann's Swifter routine helio_getacch_int.f90
!! Adapted from Hal Levison's Swift routine getacch_ah3.f
implicit none
! Arguments
class(helio_pl), intent(inout) :: pl !! Helio massive body particle data structure
real(DP), intent(in) :: t !! Current time
! Internals
integer(I4B) :: i, j
real(DP) :: rji2, irij3, faci, facj
real(DP), dimension(NDIM) :: dx

associate(npl => pl%nbody)
pl%ah(:,:) = 0.0_DP
do i = 1, npl - 1
do j = i + 1, npl
dx(:) = pl%xh(:,j) - pl%xh(:,i)
rji2 = dot_product(dx(:), dx(:))
irij3 = 1.0_DP / (rji2 * sqrt(rji2))
faci = pl%Gmass(i) * irij3
facj = pl%Gmass(j) * irij3
pl%ah(:,i) = pl%ah(:,i) + facj * dx(:)
pl%ah(:,j) = pl%ah(:,j) - faci * dx(:)
end do
end do
end associate

return
end subroutine helio_getacch_int_pl

subroutine helio_getacch_int_tp(tp, system, param, t)
!! author: David A. Minton
!!
!! Compute direct cross term heliocentric accelerations of test particles
!!
!! Adapted from David E. Kaufmann's Swifter routine helio_getacch_int_tp.f90
!! Adapted from Hal Levison's Swift routine getacch_ah3_tp.f
implicit none
! Arguments
class(helio_tp), intent(inout) :: tp !! Helio test particle object
class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object
class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters
real(DP), intent(in) :: t !! Current times
! Internals
integer(I4B) :: i, j
real(DP) :: r2, fac
real(DP), dimension(NDIM) :: dx
real(DP), dimension(:, :), allocatable :: xhp

associate(ntp => tp%nbody, pl => system%pl, npl => system%pl%nbody, lbeg => system%lbeg)
if (lbeg) then
allocate(xhp, source=pl%xbeg)
else
allocate(xhp, source=pl%xend)
end if

tp%ah(:,:) = 0.0_DP
do i = 1, ntp
if (tp%status(i) == ACTIVE) then
do j = 1, npl
dx(:) = tp%xh(:,i) - xhp(:,j)
r2 = dot_product(dx(:), dx(:))
fac = pl%Gmass(j) / (r2 * sqrt(r2))
tp%ah(:,i) = tp%ah(:,i) - fac * dx(:)
end do
end if
end do
end associate
return
end subroutine helio_getacch_int_tp

end submodule s_helio_getacch
72 changes: 70 additions & 2 deletions src/kick/kick.f90
View file
Original file line number Diff line number Diff line change
@@ -1,7 +1,73 @@
submodule(swiftest_classes) s_kick
use swiftest
contains
module subroutine kickvh_body(self, dt)
module pure subroutine kick_getacch_int_pl(self)
!! author: David A. Minton
!!
!! Compute direct cross (third) term heliocentric accelerations of massive bodies
!!
!! Adapted from Hal Levison's Swift routine getacch_ah3.f
!! Adapted from David E. Kaufmann's Swifter routine whm_getacch_ah3.f90 and helio_getacch_int.f90
implicit none
! Arguments
class(swiftest_pl), intent(inout) :: self
! Internals
integer(I4B) :: k
real(DP) :: rji2, irij3, faci, facj
real(DP), dimension(NDIM) :: dx

associate(pl => self, npl => self%nbody, nplpl => self%nplpl)
do k = 1, nplpl
associate(i => pl%k_eucl(1, k), j => pl%k_eucl(2, k))
dx(:) = pl%xh(:, j) - pl%xh(:, i)
rji2 = dot_product(dx(:), dx(:))
irij3 = 1.0_DP / (rji2 * sqrt(rji2))
faci = pl%Gmass(i) * irij3
facj = pl%Gmass(j) * irij3
pl%ah(:, i) = pl%ah(:, i) + facj * dx(:)
pl%ah(:, j) = pl%ah(:, j) - faci * dx(:)
end associate
end do
end associate

return
end subroutine kick_getacch_int_pl

module pure subroutine kick_getacch_int_tp(self, GMpl, xhp, npl)
!! author: David A. Minton
!!
!! Compute direct cross (third) term heliocentric accelerations of test particles by massive bodies
!!
!! Adapted from Hal Levison's Swift routine getacch_ah3_tp.f
!! Adapted from David E. Kaufmann's Swifter routine whm_getacch_ah3.f90 and helio_getacch_int_tp.f90
implicit none
! Arguments
class(swiftest_tp), intent(inout) :: self !! Swiftest test particle
real(DP), dimension(:), intent(in) :: GMpl !! Massive body masses
real(DP), dimension(:,:), intent(in) :: xhp !! Massive body position vectors
integer(I4B), intent(in) :: npl !! Number of active massive bodies
! Internals
integer(I4B) :: i, j
real(DP) :: rji2, irij3, fac
real(DP), dimension(NDIM) :: dx, acc

associate(tp => self, ntp => self%nbody)
do concurrent(i = 1:ntp, tp%status(i) == ACTIVE)
acc(:) = 0.0_DP
do j = 1, npl
dx(:) = tp%xh(:,i) - xhp(:, j)
rji2 = dot_product(dx(:), dx(:))
irij3 = 1.0_DP / (rji2 * sqrt(rji2))
fac = GMpl(j) * irij3
acc(:) = acc(:) - fac * dx(:)
end do
tp%ah(:, i) = tp%ah(:, i) + acc(:)
end do
end associate
return
end subroutine kick_getacch_int_tp

module subroutine kick_vh_body(self, dt)
!! author: David A. Minton
!!
!! Kick heliocentric velocities of bodies
Expand All @@ -23,6 +89,8 @@ module subroutine kickvh_body(self, dt)
end associate

return
end subroutine kickvh_body
end subroutine kick_vh_body



end submodule s_kick
35 changes: 21 additions & 14 deletions src/modules/swiftest_classes.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -8,12 +8,12 @@ module swiftest_classes
private
public :: discard_pl, discard_system, discard_tp
public :: drift_body, drift_one
public :: eucl_dist_index_plpl, eucl_dist_index_pltp, eucl_irij3_plpl
public :: eucl_dist_index_plpl, eucl_dist_index_pltp
public :: gr_getaccb_ns_body, gr_p4_pos_kick, gr_pseudovel2vel, gr_vel2pseudovel
public :: io_dump_param, io_dump_swiftest, io_dump_system, io_get_args, io_get_token, io_param_reader, io_param_writer, io_read_body_in, &
io_read_cb_in, io_read_param_in, io_read_frame_body, io_read_frame_cb, io_read_frame_system, &
io_toupper, io_write_discard, io_write_encounter, io_write_frame_body, io_write_frame_cb, io_write_frame_system
public :: kickvh_body
public :: kick_getacch_int_pl, kick_vh_body
public :: obl_acc_body, obl_acc_pl, obl_acc_tp
public :: orbel_el2xv_vec, orbel_xv2el_vec, orbel_scget, orbel_xv2aeq, orbel_xv2aqt
public :: setup_body, setup_construct_system, setup_initialize_system, setup_pl, setup_tp
Expand Down Expand Up @@ -179,7 +179,7 @@ module swiftest_classes
procedure, public :: initialize => io_read_body_in !! Read in body initial conditions from a file
procedure, public :: read_frame => io_read_frame_body !! I/O routine for writing out a single frame of time-series data for the central body
procedure, public :: write_frame => io_write_frame_body !! I/O routine for writing out a single frame of time-series data for the central body
procedure, public :: kick => kickvh_body !! Kicks the heliocentric velocities
procedure, public :: kick => kick_vh_body !! Kicks the heliocentric velocities
procedure, public :: accel_obl => obl_acc_body !! Compute the barycentric accelerations of bodies due to the oblateness of the central body
procedure, public :: el2xv => orbel_el2xv_vec !! Convert orbital elements to position and velocity vectors
procedure, public :: xv2el => orbel_xv2el_vec !! Convert position and velocity vectors to orbital elements
Expand All @@ -201,7 +201,6 @@ module swiftest_classes
real(DP), dimension(:), allocatable :: Gmass !! Mass gravitational term G * mass (units GU * MU)
real(DP), dimension(:), allocatable :: rhill !! Body mass (units MU)
real(DP), dimension(:), allocatable :: radius !! Body radius (units DU)
real(DP), dimension(:), allocatable :: irij3 !! 1.0_DP / (rji2 * sqrt(rji2)) where rji2 is the square of the Euclidean distance
real(DP), dimension(:,:), allocatable :: xbeg !! Position at beginning of step
real(DP), dimension(:,:), allocatable :: xend !! Position at end of step
real(DP), dimension(:,:), allocatable :: vbeg !! Velocity at beginning of step
Expand All @@ -219,10 +218,10 @@ module swiftest_classes
! These are concrete because they are the same implemenation for all integrators
procedure, public :: discard => discard_pl !! Placeholder method for discarding massive bodies
procedure, public :: eucl_index => eucl_dist_index_plpl !! Sets up the (i, j) -> k indexing used for the single-loop blocking Euclidean distance matrix
procedure, public :: eucl_irij3 => eucl_irij3_plpl !! Parallelized single loop blocking for Euclidean distance matrix calcualtion
procedure, public :: accel_int => kick_getacch_int_pl !! Compute direct cross (third) term heliocentric accelerations of massive bodies
procedure, public :: accel_obl => obl_acc_pl !! Compute the barycentric accelerations of bodies due to the oblateness of the central body
procedure, public :: accel_tides => tides_getacch_pl !! Compute the accelerations of bodies due to tidal interactions with the central body
procedure, public :: setup => setup_pl !! A base constructor that sets the number of bodies and allocates and initializes all arrays
procedure, public :: accel_tides => tides_getacch_pl !! Compute the accelerations of bodies due to tidal interactions with the central body
procedure, public :: set_mu => util_set_mu_pl !! Method used to construct the vectorized form of the central body mass
procedure, public :: set_rhill => util_set_rhill !! Calculates the Hill's radii for each body
procedure, public :: h2b => util_coord_h2b_pl !! Convert massive bodies from heliocentric to barycentric coordinates (position and velocity)
Expand All @@ -241,7 +240,6 @@ module swiftest_classes
integer(I4B), dimension(:), allocatable :: isperi !! Perihelion passage flag
real(DP), dimension(:), allocatable :: peri !! Perihelion distance
real(DP), dimension(:), allocatable :: atp !! Semimajor axis following perihelion passage
real(DP), dimension(:, :), allocatable :: irij3 !! 1.0_DP / (rji2 * sqrt(rji2)) where rji2 is the square of the Euclidean distance betwen each pl-tp
!! Note to developers: If you add components to this class, be sure to update methods and subroutines that traverse the
!! component list, such as setup_tp and util_spill_tp
contains
Expand All @@ -250,6 +248,7 @@ module swiftest_classes
! These are concrete because they are the same implemenation for all integrators
procedure, public :: discard => discard_tp !! Check to see if test particles should be discarded based on their positions relative to the massive bodies
procedure, public :: eucl_index => eucl_dist_index_pltp !! Sets up the (i, j) -> k indexing used for the single-loop blocking Euclidean distance matrix
procedure, public :: accel_int => kick_getacch_int_tp !! Compute direct cross (third) term heliocentric accelerations of test particles by massive bodies
procedure, public :: accel_obl => obl_acc_tp !! Compute the barycentric accelerations of bodies due to the oblateness of the central body
procedure, public :: setup => setup_tp !! A base constructor that sets the number of bodies and
procedure, public :: set_mu => util_set_mu_tp !! Method used to construct the vectorized form of the central body mass
Expand Down Expand Up @@ -412,11 +411,6 @@ module subroutine eucl_dist_index_pltp(self, pl)
class(swiftest_pl), intent(inout) :: pl !! Swiftest massive body object
end subroutine

module subroutine eucl_irij3_plpl(self)
implicit none
class(swiftest_pl), intent(inout) :: self !! Swiftest massive body object
end subroutine eucl_irij3_plpl

module pure subroutine gr_getaccb_ns_body(self, system, param)
implicit none
class(swiftest_body), intent(inout) :: self !! Swiftest generic body object
Expand Down Expand Up @@ -606,11 +600,24 @@ module subroutine io_write_frame_system(self, iu, param)
class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters
end subroutine io_write_frame_system

module subroutine kickvh_body(self, dt)
module pure subroutine kick_getacch_int_pl(self)
implicit none
class(swiftest_pl), intent(inout) :: self
end subroutine kick_getacch_int_pl

module pure subroutine kick_getacch_int_tp(self, GMpl, xhp, npl)
implicit none
class(swiftest_tp), intent(inout) :: self !! Swiftest test particle
real(DP), dimension(:), intent(in) :: GMpl !! Massive body masses
real(DP), dimension(:,:), intent(in) :: xhp !! Massive body position vectors
integer(I4B), intent(in) :: npl !! Number of active massive bodies
end subroutine kick_getacch_int_tp

module subroutine kick_vh_body(self, dt)
implicit none
class(swiftest_body), intent(inout) :: self !! Swiftest body object
real(DP), intent(in) :: dt !! Stepsize
end subroutine kickvh_body
end subroutine kick_vh_body

module subroutine obl_acc_body(self, system)
implicit none
Expand Down
Loading

0 comments on commit 4efd05b

Please sign in to comment.