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

Commit

Permalink
Restructured kick methods for entire code with new masked interface. …
Browse files Browse the repository at this point in the history 
…WHM still works, but Helio does not.
  • Loading branch information
@daminton
daminton committed Jul 26, 2021
1 parent 178c1cc commit d85ee28
Show file tree
Hide file tree
Showing 10 changed files with 259 additions and 165 deletions.
43 changes: 25 additions & 18 deletions src/helio/helio_drift.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -71,7 +71,7 @@ module subroutine helio_drift_tp(self, system, param, dt, mask)
return
end subroutine helio_drift_tp

module subroutine helio_drift_linear_pl(self, cb, dt, lbeg)
module subroutine helio_drift_linear_pl(self, cb, dt, mask, lbeg)
!! author: David A. Minton
!!
!! Perform linear drift of massive bodies due to barycentric momentum of Sun
Expand All @@ -80,21 +80,26 @@ module subroutine helio_drift_linear_pl(self, cb, dt, lbeg)
!! Adapted from Hal Levison's Swift routine helio_lindrift.f
implicit none
! Arguments
class(helio_pl), intent(inout) :: self !! Helio massive body object
class(helio_cb), intent(inout) :: cb !! Helio central bod
real(DP), intent(in) :: dt !! Stepsize
logical, intent(in) :: lbeg !! Argument that determines whether or not this is the beginning or end of the step
class(helio_pl), intent(inout) :: self !! Helio massive body object
class(helio_cb), intent(inout) :: cb !! Helio central body
real(DP), intent(in) :: dt !! Stepsize
logical, dimension(:), intent(in) :: mask !! Mask that determines which bodies to kick
logical, intent(in) :: lbeg !! Argument that determines whether or not this is the beginning or end of the step
! Internals
real(DP), dimension(NDIM) :: pt !! negative barycentric velocity of the central body
real(DP), dimension(NDIM) :: pt !! negative barycentric velocity of the central body
integer(I4B) :: i

associate(pl => self, npl => self%nbody)
pt(1) = sum(pl%Gmass(1:npl) * pl%vb(1,1:npl))
pt(2) = sum(pl%Gmass(1:npl) * pl%vb(2,1:npl))
pt(3) = sum(pl%Gmass(1:npl) * pl%vb(3,1:npl))
if (npl == 0) return
pt(1) = sum(pl%Gmass(1:npl) * pl%vb(1,1:npl), mask)
pt(2) = sum(pl%Gmass(1:npl) * pl%vb(2,1:npl), mask)
pt(3) = sum(pl%Gmass(1:npl) * pl%vb(3,1:npl), mask)
pt(:) = pt(:) / cb%Gmass
pl%xh(1,1:npl) = pl%xh(1,1:npl) + pt(1) * dt
pl%xh(2,1:npl) = pl%xh(2,1:npl) + pt(2) * dt
pl%xh(3,1:npl) = pl%xh(3,1:npl) + pt(3) * dt
do concurrent(i = 1:npl, mask(i))
pl%xh(1,1:npl) = pl%xh(1,1:npl) + pt(1) * dt
pl%xh(2,1:npl) = pl%xh(2,1:npl) + pt(2) * dt
pl%xh(3,1:npl) = pl%xh(3,1:npl) + pt(3) * dt
end do

if (lbeg) then
cb%ptbeg = pt(:)
Expand All @@ -106,7 +111,7 @@ module subroutine helio_drift_linear_pl(self, cb, dt, lbeg)
return
end subroutine helio_drift_linear_pl

module subroutine helio_drift_linear_tp(self, cb, dt, lbeg)
module subroutine helio_drift_linear_tp(self, cb, dt, mask, lbeg)
!! author: David A. Minton
!!
!! Perform linear drift of test particles due to barycentric momentum of Sun
Expand All @@ -116,20 +121,22 @@ module subroutine helio_drift_linear_tp(self, cb, dt, lbeg)
!! Adapted from Hal Levison's Swift routine helio_lindrift_tp.f
implicit none
! Arguments
class(helio_tp), intent(inout) :: self !! Helio test particleb object
class(helio_cb), intent(in) :: cb !! Helio central body
real(DP), intent(in) :: dt !! Stepsize
logical, intent(in) :: lbeg !! Argument that determines whether or not this is the beginning or end of the step
class(helio_tp), intent(inout) :: self !! Helio test particleb object
class(helio_cb), intent(in) :: cb !! Helio central body
real(DP), intent(in) :: dt !! Stepsize
logical, dimension(:), intent(in) :: mask !! Mask that determines which bodies to kick
logical, intent(in) :: lbeg !! Argument that determines whether or not this is the beginning or end of the step
! Internals
real(DP), dimension(NDIM) :: pt !! negative barycentric velocity of the central body

associate(tp => self, ntp => self%nbody)
if (ntp == 0) return
if (lbeg) then
pt(:) = cb%ptbeg
else
pt(:) = cb%ptend
end if
where (tp%status(1:ntp) == ACTIVE)
where (mask(1:ntp))
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
Expand Down
33 changes: 25 additions & 8 deletions src/helio/helio_kick.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -66,7 +66,7 @@ module subroutine helio_kick_getacch_tp(self, system, param, t, lbeg)
return
end subroutine helio_kick_getacch_tp

module subroutine helio_kick_vb_pl(self, dt)
module subroutine helio_kick_vb_pl(self, system, param, t, dt, mask, lbeg)
!! author: David A. Minton
!!
!! Kick barycentric velocities of bodies
Expand All @@ -75,14 +75,25 @@ module subroutine helio_kick_vb_pl(self, dt)
!! Adapted from David E. Kaufmann's Swifter routine helio_kick_vb.f90
implicit none
! Arguments
class(helio_pl), intent(inout) :: self !! Swiftest generic body object
real(DP), intent(in) :: dt !! Stepsize
class(helio_pl), intent(inout) :: self !! Swiftest generic body 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 time
real(DP), intent(in) :: dt !! Stepsize
logical, dimension(:), intent(in) :: mask !! Mask that determines which bodies to kick
logical, intent(in) :: lbeg !! Logical flag indicating whether this is the beginning of the half step or not.
! Internals
integer(I4B) :: i

associate(pl => self, npl => self%nbody)
if (npl ==0) return
do concurrent(i = 1:npl, pl%status(i) == ACTIVE)
call pl%accel(system, param, t)
if (lbeg) then
call pl%set_beg_end(xbeg = pl%xh)
else
call pl%set_beg_end(xend = pl%xh)
end if
do concurrent(i = 1:npl, mask(i))
pl%vb(:, i) = pl%vb(:, i) + pl%ah(:, i) * dt
end do
end associate
Expand All @@ -91,7 +102,7 @@ module subroutine helio_kick_vb_pl(self, dt)

end subroutine helio_kick_vb_pl

module subroutine helio_kick_vb_tp(self, dt)
module subroutine helio_kick_vb_tp(self, system, param, t, dt, mask, lbeg)
!! author: David A. Minton
!!
!! Kick barycentric velocities of bodies
Expand All @@ -100,14 +111,20 @@ module subroutine helio_kick_vb_tp(self, dt)
!! Adapted from David E. Kaufmann's Swifter routine helio_kick_vb_tp.f90
implicit none
! Arguments
class(helio_tp), intent(inout) :: self !! Swiftest generic body object
real(DP), intent(in) :: dt !! Stepsize
class(helio_tp), intent(inout) :: self !! Swiftest generic body 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 time
real(DP), intent(in) :: dt !! Stepsize
logical, dimension(:), intent(in) :: mask !! Mask that determines which bodies to kick
logical, intent(in) :: lbeg !! Logical flag indicating whether this is the beginning of the half step or not.
! Internals
integer(I4B) :: i

associate(tp => self, ntp => self%nbody)
if (ntp ==0) return
do concurrent(i = 1:ntp, tp%status(i) == ACTIVE)
call tp%accel(system, param, t, lbeg)
do concurrent(i = 1:ntp, mask(i))
tp%vb(:, i) = tp%vb(:, i) + tp%ah(:, i) * dt
end do
end associate
Expand Down
31 changes: 12 additions & 19 deletions src/helio/helio_step.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -37,8 +37,7 @@ module subroutine helio_step_pl(self, system, param, t, dt)
real(DP), intent(in) :: t !! Current simulation time
real(DP), intent(in) :: dt !! Stepsize
! Internals
integer(I4B) :: i
real(DP) :: dth, msys
real(DP) :: dth !! Half step size

if (self%nbody == 0) return
associate(pl => self)
Expand All @@ -49,15 +48,11 @@ module subroutine helio_step_pl(self, system, param, t, dt)
call pl%vh2vb(cb)
pl%lfirst = .false.
end if
call pl%lindrift(cb, dth, lbeg=.true.)
call pl%accel(system, param, t)
call pl%kick(dth)
call pl%set_beg_end(xbeg = pl%xh)
call pl%drift(system, param, dt, pl%status(:) == ACTIVE)
call pl%set_beg_end(xend = pl%xh)
call pl%accel(system, param, t + dt)
call pl%kick(dth)
call pl%lindrift(cb, dth, lbeg=.false.)
call pl%lindrift(cb, dth, mask=(pl%status(:) == ACTIVE), lbeg=.true.)
call pl%kick(system, param, t, dth, mask=(pl%status(:) == ACTIVE), lbeg=.true.)
call pl%drift(system, param, dt, mask=(pl%status(:) == ACTIVE))
call pl%kick(system, param, t + dt, dth, mask=(pl%status(:) == ACTIVE), lbeg=.false.)
call pl%lindrift(cb, dth, mask=(pl%status(:) == ACTIVE), lbeg=.false.)
call pl%vb2vh(cb)
end select
end associate
Expand All @@ -80,9 +75,9 @@ module subroutine helio_step_tp(self, system, param, t, dt)
class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nboody system
class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters
real(DP), intent(in) :: t !! Current simulation time
real(DP), intent(in) :: dt !! Stepsiz
real(DP), intent(in) :: dt !! Stepsize
! Internals
real(DP) :: dth !! Half step size
real(DP) :: dth !! Half step size

if (self%nbody == 0) return

Expand All @@ -94,13 +89,11 @@ module subroutine helio_step_tp(self, system, param, t, dt)
call tp%vh2vb(vbcb = -cb%ptbeg)
tp%lfirst = .false.
end if
call tp%lindrift(cb, dth, lbeg=.true.)
call tp%accel(system, param, t, lbeg=.true.)
call tp%kick(dth)
call tp%lindrift(cb, dth, mask=(tp%status(:) == ACTIVE), lbeg=.true.)
call tp%kick(system, param, t, dth, mask=(tp%status(:) == ACTIVE), lbeg=.true.)
call tp%drift(system, param, dt, tp%status(:) == ACTIVE)
call tp%accel(system, param, t + dt, lbeg=.false.)
call tp%kick(dth)
call tp%lindrift(cb, dth, lbeg=.false.)
call tp%kick(system, param, t + dt, dth, mask=(tp%status(:) == ACTIVE), lbeg=.false.)
call tp%lindrift(cb, dth, mask=(tp%status(:) == ACTIVE), lbeg=.false.)
call tp%vb2vh(vbcb = -cb%ptend)
end select
end associate
Expand Down
26 changes: 0 additions & 26 deletions src/kick/kick.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -64,30 +64,4 @@ module pure subroutine kick_getacch_int_tp(self, GMpl, xhp, npl)
return
end subroutine kick_getacch_int_tp

module subroutine kick_vh_body(self, dt)
!! author: David A. Minton
!!
!! Kick heliocentric velocities of bodies
!!
!! Adapted from Martin Duncan and Hal Levison's Swift routine kickvh.f and kickvh_tp.f
!! Adapted from David E. Kaufmann's Swifter routine whm_kickvh.f90 and whm_kickvh_tp.f90
implicit none
! Arguments
class(swiftest_body), intent(inout) :: self !! Swiftest generic body object
real(DP), intent(in) :: dt !! Stepsize
! Internals
integer(I4B) :: i

associate(n => self%nbody, vh => self%vh, ah => self%ah, status => self%status)
if (n == 0) return
do i = 1, n
if (status(i) == ACTIVE) vh(:, i) = vh(:, i) + ah(:, i) * dt
end do
end associate

return
end subroutine kick_vh_body



end submodule s_kick
48 changes: 31 additions & 17 deletions src/modules/helio_classes.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -116,20 +116,22 @@ module subroutine helio_drift_tp(self, system, param, dt, mask)
logical, dimension(:), intent(in) :: mask !! Logical mask of size self%nbody that determines which bodies to drift
end subroutine helio_drift_tp

module subroutine helio_drift_linear_pl(self, cb, dt, lbeg)
module subroutine helio_drift_linear_pl(self, cb, dt, mask, lbeg)
implicit none
class(helio_pl), intent(inout) :: self !! Helio massive body object
class(helio_cb), intent(inout) :: cb !! Helio central body object
real(DP), intent(in) :: dt !! Stepsize
logical, intent(in) :: lbeg !! Argument that determines whether or not this is the beginning or end of the step
class(helio_pl), intent(inout) :: self !! Helio massive body object
class(helio_cb), intent(inout) :: cb !! Helio central body
real(DP), intent(in) :: dt !! Stepsize
logical, dimension(:), intent(in) :: mask !! Mask that determines which bodies to kick
logical, intent(in) :: lbeg !! Argument that determines whether or not this is the beginning or end of the step
end subroutine helio_drift_linear_pl

module subroutine helio_drift_linear_tp(self, cb, dt, lbeg)
module subroutine helio_drift_linear_tp(self, cb, dt, mask, lbeg)
implicit none
class(helio_tp), intent(inout) :: self !! Helio test particle object
class(helio_cb), intent(in) :: cb !! Helio nbody system object
real(DP), intent(in) :: dt !! Stepsize
logical, intent(in) :: lbeg !! Argument that determines whether or not this is the beginning or end of the step
class(helio_tp), intent(inout) :: self !! Helio test particle object
class(helio_cb), intent(in) :: cb !! Helio central body
real(DP), intent(in) :: dt !! Stepsize
logical, dimension(:), intent(in) :: mask !! Mask that determines which bodies to kick
logical, intent(in) :: lbeg !! Argument that determines whether or not this is the beginning or end of the step
end subroutine helio_drift_linear_tp

module subroutine helio_kick_getacch_pl(self, system, param, t, lbeg)
Expand All @@ -143,7 +145,7 @@ module subroutine helio_kick_getacch_pl(self, system, param, t, lbeg)
end subroutine helio_kick_getacch_pl

module subroutine helio_kick_getacch_tp(self, system, param, t, lbeg)
use swiftest_classes, only : swiftest_parameters, swiftest_nbody_system
use swiftest_classes, only : swiftest_nbody_system, swiftest_parameters
implicit none
class(helio_tp), intent(inout) :: self !! Helio test particle object
class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object
Expand All @@ -152,16 +154,28 @@ module subroutine helio_kick_getacch_tp(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
end subroutine helio_kick_getacch_tp

module subroutine helio_kick_vb_pl(self, dt)
module subroutine helio_kick_vb_pl(self, system, param, t, dt, mask, lbeg)
use swiftest_classes, only : swiftest_nbody_system, swiftest_parameters
implicit none
class(helio_pl), intent(inout) :: self !! Helio massive body object
real(DP), intent(in) :: dt !! Stepsize
class(helio_pl), intent(inout) :: self !! Helio massive body 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 time
real(DP), intent(in) :: dt !! Stepsize
logical, dimension(:), intent(in) :: mask !! Mask that determines which bodies to kick
logical, intent(in) :: lbeg !! Logical flag indicating whether this is the beginning of the half step or not.
end subroutine helio_kick_vb_pl

module subroutine helio_kick_vb_tp(self, dt)
module subroutine helio_kick_vb_tp(self, system, param, t, dt, mask, lbeg)
use swiftest_classes, only : swiftest_nbody_system, swiftest_parameters
implicit none
class(helio_tp), intent(inout) :: self !! Helio test particle object
real(DP), intent(in) :: dt !! Stepsize
class(helio_tp), intent(inout) :: self !! 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 time
real(DP), intent(in) :: dt !! Stepsize
logical, dimension(:), intent(in) :: mask !! Mask that determines which bodies to kick
logical, intent(in) :: lbeg !! Logical flag indicating whether this is the beginning of the half step or not.
end subroutine helio_kick_vb_tp

module subroutine helio_step_pl(self, system, param, t, dt)
Expand Down
Loading

0 comments on commit d85ee28

Please sign in to comment.