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

Commit

Permalink
Modified linear drift interfaces to match and now save the ptb and pt…
Browse files Browse the repository at this point in the history 
…e values inside the helio_nbody_system object
  • Loading branch information
@daminton
daminton committed Jul 7, 2021
1 parent b782c6c commit 886b63b
Show file tree
Hide file tree
Showing 3 changed files with 61 additions and 55 deletions.
43 changes: 22 additions & 21 deletions src/helio/helio_drift.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -56,7 +56,7 @@ module subroutine helio_drift_pl(self, system, param, dt)
return
end subroutine helio_drift_pl

module subroutine helio_drift_linear_pl(self, cb, dt, pt)
module subroutine helio_drift_linear_pl(self, system, dt, pt)
!! author: David A. Minton
!!
!! Perform linear drift of massive bodies due to barycentric momentum of Sun
Expand All @@ -65,23 +65,23 @@ module subroutine helio_drift_linear_pl(self, cb, dt, pt)
!! Adapted from Hal Levison's Swift routine helio_lindrift.f
implicit none
! Arguments
class(helio_pl), intent(inout) :: self !! Helio massive body object
class(swiftest_cb), intent(in) :: cb !! Helio central body object
real(DP), intent(in) :: dt !! Stepsize
real(DP), dimension(:), intent(out) :: pt !! negative barycentric velocity of the central body
class(helio_pl), intent(inout) :: self !! Helio massive body object
class(helio_nbody_system), intent(in) :: system !! Swiftest nbody system object
real(DP), intent(in) :: dt !! Stepsize
real(DP), dimension(:), intent(out) :: pt !! negative barycentric velocity of the central body
! Internals
integer(I4B) :: i
real(DP),dimension(NDIM) :: pttmp !intent(out) variables don't play nicely
!with openmp's reduction for some reason

associate(npl => self%nbody, xh => self%xh, vb => self%vb, GMpl => self%Gmass, GMcb => cb%Gmass)
associate(pl => self, npl => self%nbody, cb => system%cb)
pttmp(:) = 0.0_DP
do i = 2, npl
pttmp(:) = pttmp(:) + GMpl(i) * vb(:,i)
pttmp(:) = pttmp(:) + pl%Gmass(i) * pl%vb(:,i)
end do
pttmp(:) = pttmp(:) / GMcb
pttmp(:) = pttmp(:) / cb%Gmass
do i = 2, npl
xh(:,i) = xh(:,i) + pttmp(:) * dt
pl%xh(:,i) = pl%xh(:,i) + pttmp(:) * dt
end do
pt(:) = pttmp(:)
end associate
Expand All @@ -100,13 +100,13 @@ module subroutine helio_drift_tp(self, system, param, dt)
implicit none
! Arguments
class(helio_tp), intent(inout) :: self !! Helio test particle object
class(swiftest_nbody_system), intent(inout) :: system !! WHM nbody system object
class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object
class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters of
real(DP), intent(in) :: dt !! Stepsiz
real(DP), intent(in) :: dt !! Stepsize
! Internals
integer(I4B) :: i !! Loop counter
real(DP) :: rmag, vmag2, energy
real(DP), dimension(:), allocatable :: dtp
real(DP), dimension(:), allocatable :: dtp
integer(I4B), dimension(:),allocatable :: iflag !! Vectorized error code flag

associate(tp => self, ntp => self%nbody)
Expand Down Expand Up @@ -139,7 +139,7 @@ module subroutine helio_drift_tp(self, system, param, dt)
return
end subroutine helio_drift_tp

module subroutine helio_drift_linear_tp(self, dt, pt)
module subroutine helio_drift_linear_tp(self, system, dt, pt)
!! author: David A. Minton
!!
!! Perform linear drift of test particles due to barycentric momentum of Sun
Expand All @@ -149,15 +149,16 @@ module subroutine helio_drift_linear_tp(self, dt, pt)
!! Adapted from Hal Levison's Swift routine helio_lindrift_tp.f
implicit none
! Arguments
class(helio_tp), intent(inout) :: self !! Helio test particle data structure
real(DP), intent(in) :: dt !! Stepsize
real(DP), dimension(:), intent(in) :: pt !! negative barycentric velocity of the Sun
class(helio_tp), intent(inout) :: self !! Helio test particleb object
class(helio_nbody_system), intent(in) :: system !! Swiftest nbody system object
real(DP), intent(in) :: dt !! Stepsize
real(DP), dimension(:), intent(in) :: pt !! negative barycentric velocity of the central body

associate(ntp => self%nbody, xh => self%xh, status => self%status)
where (status(1:ntp) == ACTIVE)
xh(1, 1:ntp) = xh(1, 1:ntp) + pt(1) * dt
xh(2, 1:ntp) = xh(2, 1:ntp) + pt(2) * dt
xh(3, 1:ntp) = xh(3, 1:ntp) + pt(3) * dt
associate(tp => self, ntp => self%nbody)
where (tp%status(1:ntp) == ACTIVE)
tp%xh(1, 1:ntp) = tp%xh(1, 1:ntp) + pt(1) * dt
tp%xh(2, 1:ntp) = tp%xh(2, 1:ntp) + pt(2) * dt
tp%xh(3, 1:ntp) = tp%xh(3, 1:ntp) + pt(3) * dt
end where
end associate

Expand Down
51 changes: 27 additions & 24 deletions src/helio/helio_step.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -50,26 +50,29 @@ module subroutine helio_step_pl(self, system, param, t, dt)
! Internals
integer(I4B) :: i
real(DP) :: dth, msys
real(DP), dimension(NDIM) :: ptbeg, ptend !! TODO: Incorporate these into the tp structure
!real(DP), dimension(NDIM) :: ptbeg, ptend !! TODO: Incorporate these into the tp structure
logical, save :: lfirst = .true.

associate(pl => self, cb => system%cb)
dth = 0.5_DP * dt
if (lfirst) then
call pl%vh2vb(cb)
lfirst = .false.
end if
call pl%lindrift(cb, dth, ptbeg)
call pl%get_accel(system, param, t)
call pl%kick(dth)
call system%set_beg_end(xbeg = pl%xh)
call pl%drift(system, param, dt)
call system%set_beg_end(xend = pl%xh)
call pl%get_accel(system, param, t + dt)
call pl%kick(dth)
call pl%lindrift(cb, dth, ptend)
call pl%vb2vh(cb)
end associate

select type(system)
class is (helio_nbody_system)
associate(pl => self, cb => system%cb, ptb => system%ptb, pte => system%pte)
dth = 0.5_DP * dt
if (lfirst) then
call pl%vh2vb(cb)
lfirst = .false.
end if
call pl%lindrift(system, dth, ptb)
call pl%get_accel(system, param, t)
call pl%kick(dth)
call system%set_beg_end(xbeg = pl%xh)
call pl%drift(system, param, dt)
call system%set_beg_end(xend = pl%xh)
call pl%get_accel(system, param, t + dt)
call pl%kick(dth)
call pl%lindrift(system, dth, pte)
call pl%vb2vh(cb)
end associate
end select

return

Expand All @@ -96,20 +99,20 @@ module subroutine helio_step_tp(self, system, param, t, dt)

select type(system)
class is (helio_nbody_system)
associate(tp => self, cb => system%cb, pl => system%pl, xbeg => system%xbeg, xend => system%xend)
associate(tp => self, cb => system%cb, pl => system%pl, xbeg => system%xbeg, xend => system%xend, ptb => system%ptb, pte => system%pte)
dth = 0.5_DP * dt
if (lfirst) then
call tp%vh2vb(vbcb = -tp%ptbeg)
call tp%vh2vb(vbcb = -ptb)
lfirst = .false.
end if
call tp%lindrift(dth, tp%ptbeg)
call tp%lindrift(system, dth, ptb)
call tp%get_accel(system, param, t, xbeg)
call tp%kick(dth)
call tp%drift(system, param, dt)
call tp%get_accel(system, param, t + dt, xend)
call tp%kick(dth)
call tp%lindrift(dth, tp%ptend)
call tp%vb2vh(vbcb = -tp%ptend)
call tp%lindrift(system, dth, pte)
call tp%vb2vh(vbcb = -pte)
end associate
end select

Expand Down
22 changes: 12 additions & 10 deletions src/modules/helio_classes.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -13,6 +13,8 @@ module helio_classes
! helio_nbody_system class definitions and method interfaces
!********************************************************************************************************************************
type, public, extends(whm_nbody_system) :: helio_nbody_system
real(DP), dimension(NDIM) :: ptb !! negative barycentric velocity of the central body at the beginning of time step
real(DP), dimension(NDIM) :: pte !! negative barycentric velocity of the central body at the end of time step
contains
private
procedure, public :: initialize => helio_setup_system !! Performs Helio-specific initilization steps,
Expand Down Expand Up @@ -110,20 +112,20 @@ module subroutine helio_drift_tp(self, system, param, dt)
real(DP), intent(in) :: dt !! Stepsize
end subroutine helio_drift_tp

module subroutine helio_drift_linear_pl(self, cb, dt, pt)
use swiftest_classes, only : swiftest_cb
module subroutine helio_drift_linear_pl(self, system, dt, pt)
implicit none
class(helio_pl), intent(inout) :: self !! Helio test particle object
class(swiftest_cb), intent(in) :: cb !! Helio central body object
real(DP), intent(in) :: dt !! Stepsize
real(DP), dimension(:), intent(out) :: pt !! negative barycentric velocity of the central body
class(helio_pl), intent(inout) :: self !! Helio massive body object
class(helio_nbody_system), intent(in) :: system !! Helio nbody system object
real(DP), intent(in) :: dt !! Stepsize
real(DP), dimension(:), intent(out) :: pt !! negative barycentric velocity of the central body
end subroutine helio_drift_linear_pl

module subroutine helio_drift_linear_tp(self, dt, pt)
module subroutine helio_drift_linear_tp(self, system, dt, pt)
implicit none
class(helio_tp), intent(inout) :: self !! Helio test particle object
real(DP), intent(in) :: dt !! Stepsize
real(DP), dimension(:), intent(in) :: pt !! negative barycentric velocity of the Sun
class(helio_tp), intent(inout) :: self !! Helio test particle object
class(helio_nbody_system), intent(in) :: system !! Helio nbody system object
real(DP), intent(in) :: dt !! Stepsize
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)
Expand Down

0 comments on commit 886b63b

Please sign in to comment.