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

Commit

Permalink
Completed a major restructuring to streamline the beginning/ending po…
Browse files Browse the repository at this point in the history 
…sition and acceleration saving and retreival system for test particle accelerations
  • Loading branch information
@daminton
daminton committed Jul 8, 2021
1 parent ae70ca1 commit 3b0041f
Show file tree
Hide file tree
Showing 22 changed files with 867 additions and 763 deletions.
10 changes: 5 additions & 5 deletions examples/rmvs_swifter_comparison/9pl_18tp_encounters/swiftest_rmvs_vs_swifter_rmvs.ipynb
View file

Large diffs are not rendered by default.

73 changes: 38 additions & 35 deletions src/helio/helio_getacch.f90
View file
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
submodule (helio_classes) s_helio_getacch
use swiftest
contains
module subroutine helio_getacch_pl(self, system, param, t)
module subroutine helio_getacch_pl(self, system, param, t, lbeg)
!! author: David A. Minton
!!
!! Compute heliocentric accelerations of massive bodies
Expand All @@ -14,20 +14,19 @@ module subroutine helio_getacch_pl(self, system, param, t)
class(swiftest_nbody_system), intent(inout) :: system !! WHM nbody system object
class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters of
real(DP), intent(in) :: t !! Current simulation time
! Internals
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(:,:) = pl%ahi(:,:)
if (param%loblatecb) call pl%obl_acc(cb)
if (param%lextra_force) call pl%user_getacch(system, param, t)
!if (param%lgr) call pl%gr_get_accel(param)
if (param%loblatecb) call pl%accel_obl(system)
if (param%lextra_force) call pl%accel_user(system, param, t)
!if (param%lgr) call pl%gr_accel(param)
end associate

return
end subroutine helio_getacch_pl

module subroutine helio_getacch_tp(self, system, param, t, xhp)
module subroutine helio_getacch_tp(self, system, param, t, lbeg)
!! author: David A. Minton
!!
!! Compute heliocentric accelerations of test particles
Expand All @@ -40,22 +39,19 @@ module subroutine helio_getacch_tp(self, system, param, t, xhp)
class(swiftest_nbody_system), intent(inout) :: system !! WHM nbody system object
class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters of
real(DP), intent(in) :: t !! Current time
real(DP), dimension(:,:), intent(in) :: xhp !! Heliocentric positions of planets at the current substep
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
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)
select type(pl => system%pl)
class is (helio_pl)
call helio_getacch_int_tp(tp, pl, t, xhp)
tp%ah(:,:) = tp%ahi(:,:)
if (param%loblatecb) call tp%obl_acc(cb)
if (param%lextra_force) call tp%user_getacch(system, param, t)
!if (param%lgr) call tp%gr_get_accel(param)
end select
associate(tp => self, ntp => self%nbody, cb => system%cb, npl => system%pl%nbody)
if (present(lbeg)) system%lbeg = lbeg
call helio_getacch_int_tp(tp, system, param, t)
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
Expand All @@ -77,24 +73,24 @@ subroutine helio_getacch_int_pl(pl, t)
real(DP), dimension(NDIM) :: dx

associate(npl => pl%nbody)
pl%ahi(:,:) = 0.0_DP
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%ahi(:,i) = pl%ahi(:,i) + facj * dx(:)
pl%ahi(:,j) = pl%ahi(:,j) - faci * dx(:)
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, pl, t, xhp)
subroutine helio_getacch_int_tp(tp, system, param, t)
!! author: David A. Minton
!!
!! Compute direct cross term heliocentric accelerations of test particles
Expand All @@ -103,24 +99,31 @@ subroutine helio_getacch_int_tp(tp, pl, t, xhp)
!! Adapted from Hal Levison's Swift routine getacch_ah3_tp.f
implicit none
! Arguments
class(helio_tp), intent(inout) :: tp !! Helio test particle data structure
class(swiftest_pl), intent(inout) :: pl !! Helio massive body particle data structure
real(DP), intent(in) :: t !! Current time
real(DP), dimension(:,:), intent(in) :: xhp !! Heliocentric positions of planets
class(helio_tp), intent(inout) :: tp !! WHM test particle data structure
class(swiftest_nbody_system), intent(inout) :: system !! WHM nbody system object
class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters of
real(DP), intent(in) :: t !! Current times
! Internals
integer(I4B) :: i, j
real(DP) :: r2, fac
real(DP), dimension(NDIM) :: dx
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)
if (system%lbeg) then
allocate(xhp, source=pl%xbeg)
else
allocate(xhp, source=pl%xend)
end if

associate(ntp => tp%nbody, npl => pl%nbody)
tp%ahi(:,:) = 0.0_DP
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%ahi(:,i) = tp%ahi(:,i) - fac * dx(:)
tp%ah(:,i) = tp%ah(:,i) - fac * dx(:)
end do
end if
end do
Expand Down
8 changes: 4 additions & 4 deletions src/helio/helio_setup.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -27,8 +27,8 @@ end subroutine helio_setup_system
call setup_pl(self, n)
if (n <= 0) return

allocate(self%ahi(NDIM, n))
self%ahi(:,:) = 0.0_DP
allocate(self%ah(NDIM, n))
self%ah(:,:) = 0.0_DP
return
end procedure helio_setup_pl

Expand All @@ -44,8 +44,8 @@ end subroutine helio_setup_system
call setup_tp(self, n)
if (n <= 0) return

allocate(self%ahi(NDIM, n))
self%ahi(:,:) = 0.0_DP
allocate(self%ah(NDIM, n))
self%ah(:,:) = 0.0_DP

return
end procedure helio_setup_tp
Expand Down
15 changes: 9 additions & 6 deletions src/helio/helio_step.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -16,8 +16,8 @@ module subroutine helio_step_system(self, param, t, dt)
real(DP), intent(in) :: dt !! Current stepsize

associate(system => self, cb => self%cb, pl => self%pl, tp => self%tp)
call pl%set_rhill(cb)
tp%lfirst = pl%lfirst
call pl%set_rhill(cb)
call pl%step(system, param, t, dt)
call tp%step(system, param, t, dt)
end associate
Expand All @@ -41,7 +41,8 @@ module subroutine helio_step_pl(self, system, param, t, dt)
! Internals
integer(I4B) :: i
real(DP) :: dth, msys


if (self%nbody == 0) return
select type(system)
class is (helio_nbody_system)
associate(pl => self, cb => system%cb, ptb => system%ptb, pte => system%pte)
Expand All @@ -51,12 +52,12 @@ module subroutine helio_step_pl(self, system, param, t, dt)
pl%lfirst = .false.
end if
call pl%lindrift(system, dth, ptb)
call pl%get_accel(system, param, t)
call pl%accel(system, param, t)
call pl%kick(dth)
call pl%set_beg_end(xbeg = pl%xh)
call pl%drift(system, param, dt)
call pl%set_beg_end(xend = pl%xh)
call pl%get_accel(system, param, t + dt)
call pl%accel(system, param, t + dt)
call pl%kick(dth)
call pl%lindrift(system, dth, pte)
call pl%vb2vh(cb)
Expand Down Expand Up @@ -85,6 +86,8 @@ module subroutine helio_step_tp(self, system, param, t, dt)
! Internals
real(DP) :: dth !! Half step size

if (self%nbody == 0) return

select type(system)
class is (helio_nbody_system)
associate(tp => self, cb => system%cb, pl => system%pl, ptb => system%ptb, pte => system%pte)
Expand All @@ -94,10 +97,10 @@ module subroutine helio_step_tp(self, system, param, t, dt)
tp%lfirst = .false.
end if
call tp%lindrift(system, dth, ptb)
call tp%get_accel(system, param, t, pl%xbeg)
call tp%accel(system, param, t, lbeg=.true.)
call tp%kick(dth)
call tp%drift(system, param, dt)
call tp%get_accel(system, param, t + dt, pl%xend)
call tp%accel(system, param, t + dt, lbeg=.false.)
call tp%kick(dth)
call tp%lindrift(system, dth, pte)
call tp%vb2vh(vbcb = -pte)
Expand Down
8 changes: 4 additions & 4 deletions src/io/io.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1080,7 +1080,7 @@ module subroutine io_write_discard(self, param)
if (param%lgr) then
select type(discards)
class is (whm_tp)
call discards%gr_pv2vh(param)
call discards%pv2vh(param)
end select
end if
write(LUN, HDRFMT) t, nsp, param%lbig_discard
Expand All @@ -1097,7 +1097,7 @@ module subroutine io_write_discard(self, param)
allocate(pltemp, source = pl)
select type(pltemp)
class is (whm_pl)
call pltemp%gr_pv2vh(param)
call pltemp%pv2vh(param)
allocate(vh, source = pltemp%vh)
end select
deallocate(pltemp)
Expand Down Expand Up @@ -1316,11 +1316,11 @@ module subroutine io_write_frame_system(self, iu, param)
associate(vh => pl%vh, vht => tp%vh)
select type(pl)
class is (whm_pl)
call pl%gr_pv2vh(param)
call pl%pv2vh(param)
end select
select type(tp)
class is (whm_tp)
call tp%gr_pv2vh(param)
call tp%pv2vh(param)
end select
end associate
end if
Expand Down
40 changes: 20 additions & 20 deletions src/modules/helio_classes.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -37,14 +37,14 @@ module helio_classes
type, public, extends(swiftest_pl) :: helio_pl
real(DP), dimension(:,:), allocatable :: ahi !! heliocentric acceleration due to interactions
contains
procedure, public :: vh2vb => helio_coord_vh2vb_pl !! Convert massive bodies from heliocentric to barycentric coordinates (velocity only)
procedure, public :: vb2vh => helio_coord_vb2vh_pl !! Convert massive bodies from barycentric to heliocentric coordinates (velocity only)
procedure, public :: drift => helio_drift_pl !! Method for Danby drift in Democratic Heliocentric coordinates
procedure, public :: lindrift => helio_drift_linear_pl !! Method for linear drift of massive bodies due to barycentric momentum of Sun
procedure, public :: get_accel => helio_getacch_pl !! Compute heliocentric accelerations of massive bodies
procedure, public :: kick => helio_kickvb_pl !! Kicks the barycentric velocities
procedure, public :: setup => helio_setup_pl !! Constructor method - Allocates space for number of particles
procedure, public :: step => helio_step_pl !! Steps the body forward one stepsize
procedure, public :: vh2vb => helio_coord_vh2vb_pl !! Convert massive bodies from heliocentric to barycentric coordinates (velocity only)
procedure, public :: vb2vh => helio_coord_vb2vh_pl !! Convert massive bodies from barycentric to heliocentric coordinates (velocity only)
procedure, public :: drift => helio_drift_pl !! Method for Danby drift in Democratic Heliocentric coordinates
procedure, public :: lindrift => helio_drift_linear_pl !! Method for linear drift of massive bodies due to barycentric momentum of Sun
procedure, public :: accel => helio_getacch_pl !! Compute heliocentric accelerations of massive bodies
procedure, public :: kick => helio_kickvb_pl !! Kicks the barycentric velocities
procedure, public :: setup => helio_setup_pl !! Constructor method - Allocates space for number of particles
procedure, public :: step => helio_step_pl !! Steps the body forward one stepsize
end type helio_pl

!********************************************************************************************************************************
Expand All @@ -53,18 +53,17 @@ module helio_classes

!! Helio test particle class
type, public, extends(swiftest_tp) :: helio_tp
real(DP), dimension(:,:), allocatable :: ahi !! heliocentric acceleration due to interactions
real(DP), dimension(NDIM) :: ptbeg !! negative barycentric velocity of the Sun at beginning of time step
real(DP), dimension(NDIM) :: ptend !! negative barycentric velocity of the Sun at beginning of time step
contains
procedure, public :: vh2vb => helio_coord_vh2vb_tp !! Convert test particles from heliocentric to barycentric coordinates (velocity only)
procedure, public :: vb2vh => helio_coord_vb2vh_tp !! Convert test particles from barycentric to heliocentric coordinates (velocity only)
procedure, public :: drift => helio_drift_tp !! Method for Danby drift in Democratic Heliocentric coordinates
procedure, public :: lindrift => helio_drift_linear_tp !! Method for linear drift of massive bodies due to barycentric momentum of Sun
procedure, public :: get_accel => helio_getacch_tp !! Compute heliocentric accelerations of massive bodies
procedure, public :: kick => helio_kickvb_tp !! Kicks the barycentric velocities
procedure, public :: setup => helio_setup_tp !! Constructor method - Allocates space for number of particles
procedure, public :: step => helio_step_tp !! Steps the body forward one stepsize
procedure, public :: vh2vb => helio_coord_vh2vb_tp !! Convert test particles from heliocentric to barycentric coordinates (velocity only)
procedure, public :: vb2vh => helio_coord_vb2vh_tp !! Convert test particles from barycentric to heliocentric coordinates (velocity only)
procedure, public :: drift => helio_drift_tp !! Method for Danby drift in Democratic Heliocentric coordinates
procedure, public :: lindrift => helio_drift_linear_tp !! Method for linear drift of massive bodies due to barycentric momentum of Sun
procedure, public :: accel => helio_getacch_tp !! Compute heliocentric accelerations of massive bodies
procedure, public :: kick => helio_kickvb_tp !! Kicks the barycentric velocities
procedure, public :: setup => helio_setup_tp !! Constructor method - Allocates space for number of particles
procedure, public :: step => helio_step_tp !! Steps the body forward one stepsize
end type helio_tp

interface
Expand Down Expand Up @@ -128,23 +127,24 @@ module subroutine helio_drift_linear_tp(self, system, dt, pt)
real(DP), dimension(:), intent(in) :: pt !! negative barycentric velocity of the Sun
end subroutine helio_drift_linear_tp

module subroutine helio_getacch_pl(self, system, param, t)
module subroutine helio_getacch_pl(self, system, param, t, lbeg)
use swiftest_classes, only : swiftest_parameters, swiftest_nbody_system
implicit none
class(helio_pl), intent(inout) :: self !! Helio massive body object
class(swiftest_nbody_system), intent(inout) :: system !! WHM nbody system object
class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters of
real(DP), intent(in) :: t !! Current simulation time
logical, optional, intent(in) :: lbeg !! Optional argument that determines whether or not this is the beginning or end of the step
end subroutine helio_getacch_pl

module subroutine helio_getacch_tp(self, system, param, t, xhp)
module subroutine helio_getacch_tp(self, system, param, t, lbeg)
use swiftest_classes, only : swiftest_parameters, swiftest_nbody_system
implicit none
class(helio_tp), intent(inout) :: self !! Helio test particle object
class(swiftest_nbody_system), intent(inout) :: system !! WHM nbody system object
class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters
real(DP), intent(in) :: t !! Current time
real(DP), dimension(:,:), intent(in) :: xhp !! Heliocentric positions of planets
logical, optional, intent(in) :: lbeg !! Optional argument that determines whether or not this is the beginning or end of the step
end subroutine helio_getacch_tp

module subroutine helio_kickvb_pl(self, dt)
Expand Down
Loading

0 comments on commit 3b0041f

Please sign in to comment.