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

Commit

Permalink
Fixed bug in which acceleration was not initialized to 0 on on half o…
Browse files Browse the repository at this point in the history 
…f a WHM pl step
  • Loading branch information
@daminton
daminton committed Jul 26, 2021
1 parent 028bc69 commit ebc2796
Show file tree
Hide file tree
Showing 4 changed files with 25 additions and 20 deletions.
4 changes: 2 additions & 2 deletions src/helio/helio_kick.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -17,7 +17,6 @@ module subroutine helio_kick_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)
pl%ah(:,:) = 0.0_DP
call pl%accel_int()
if (param%loblatecb) then
cb%aoblbeg = cb%aobl
Expand Down Expand Up @@ -52,7 +51,6 @@ 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

associate(tp => self, cb => system%cb, pl => system%pl, npl => system%pl%nbody)
tp%ah(:,:) = 0.0_DP
if (present(lbeg)) system%lbeg = lbeg
if (system%lbeg) then
call tp%accel_int(pl%Gmass(:), pl%xbeg(:,:), npl)
Expand Down Expand Up @@ -87,6 +85,7 @@ module subroutine helio_kick_vb_pl(self, system, param, t, dt, mask, lbeg)

associate(pl => self, npl => self%nbody)
if (npl ==0) return
pl%ah(:,:) = 0.0_DP
call pl%accel(system, param, t)
if (lbeg) then
call pl%set_beg_end(xbeg = pl%xh)
Expand Down Expand Up @@ -123,6 +122,7 @@ module subroutine helio_kick_vb_tp(self, system, param, t, dt, mask, lbeg)

associate(tp => self, ntp => self%nbody)
if (ntp ==0) return
tp%ah(:,:) = 0.0_DP
call tp%accel(system, param, t, lbeg)
do concurrent(i = 1:ntp, mask(i))
tp%vb(:, i) = tp%vb(:, i) + tp%ah(:, i) * dt
Expand Down
14 changes: 7 additions & 7 deletions src/symba/symba_kick.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -12,16 +12,16 @@ module subroutine symba_kick_pltpenc(self, system, dt, irec, sgn)
!! Adapted from Hal Levison's Swift routine symba5_kick.f
implicit none
! Arguments
class(symba_pltpenc), intent(in) :: self !! SyMBA pl-tp encounter list object
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
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
logical :: isplpl, lgoodlevel

select type(self)
class is (symba_plplenc)
Expand Down
26 changes: 16 additions & 10 deletions src/whm/whm_kick.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -25,7 +25,7 @@ module subroutine whm_kick_getacch_pl(self, system, param, t, lbeg)

ah0(:) = whm_kick_getacch_ah0(pl%Gmass(2:npl), pl%xh(:,2:npl), npl-1)
do i = 1, npl
pl%ah(:, i) = ah0(:)
pl%ah(:, i) = pl%ah(:, i) + ah0(:)
end do

call whm_kick_getacch_ah1(cb, pl)
Expand Down Expand Up @@ -75,13 +75,13 @@ module subroutine whm_kick_getacch_tp(self, system, param, t, lbeg)
if (system%lbeg) then
ah0(:) = whm_kick_getacch_ah0(pl%Gmass(:), pl%xbeg(:,:), npl)
do i = 1, ntp
tp%ah(:, i) = ah0(:)
tp%ah(:, i) = tp%ah(:, i) + ah0(:)
end do
call tp%accel_int(pl%Gmass(:), pl%xbeg(:,:), npl)
else
ah0(:) = whm_kick_getacch_ah0(pl%Gmass(:), pl%xend(:,:), npl)
do i = 1, ntp
tp%ah(:, i) = ah0(:)
tp%ah(:, i) = tp%ah(:, i) + ah0(:)
end do
call tp%accel_int(pl%Gmass(:), pl%xend(:,:), npl)
end if
Expand Down Expand Up @@ -200,16 +200,18 @@ module subroutine whm_kick_vh_pl(self, system, param, t, dt, mask, lbeg)

associate(pl => self, npl => self%nbody, cb => system%cb)
if (npl == 0) return
if (pl%lfirst) then
call pl%h2j(cb)
call pl%accel(system, param, t)
pl%lfirst = .false.
end if
if (lbeg) then
if (pl%lfirst) then
call pl%h2j(cb)
pl%ah(:,:) = 0.0_DP
call pl%accel(system, param, t)
pl%lfirst = .false.
end if
call pl%set_beg_end(xbeg = pl%xh)
else
call pl%set_beg_end(xend = pl%xh)
pl%ah(:,:) = 0.0_DP
call pl%accel(system, param, t)
call pl%set_beg_end(xend = pl%xh)
end if
do concurrent(i = 1:npl, mask(i))
pl%vh(:, i) = pl%vh(:, i) + pl%ah(:, i) * dt
Expand Down Expand Up @@ -241,10 +243,14 @@ module subroutine whm_kick_vh_tp(self, system, param, t, dt, mask, lbeg)
associate(tp => self, ntp => self%nbody)
if (ntp == 0) return
if (tp%lfirst) then
tp%ah(:,:) = 0.0_DP
call tp%accel(system, param, t, lbeg=.true.)
tp%lfirst = .false.
end if
if (.not.lbeg) call tp%accel(system, param, t, lbeg)
if (.not.lbeg) then
tp%ah(:,:) = 0.0_DP
call tp%accel(system, param, t, lbeg)
end if
do concurrent(i = 1:ntp, mask(i))
tp%vh(:, i) = tp%vh(:, i) + tp%ah(:, i) * dt
end do
Expand Down
1 change: 0 additions & 1 deletion src/whm/whm_step.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -54,7 +54,6 @@ module subroutine whm_step_pl(self, system, param, t, dt)
if (param%lgr) call pl%gr_pos_kick(param, dth)
call pl%j2h(cb)
call pl%kick(system, param, t + dt, dth, mask=(pl%status(:) == ACTIVE), lbeg=.false.)
call pl%set_beg_end(xend = pl%xh)
end associate
return
end subroutine whm_step_pl
Expand Down

0 comments on commit ebc2796

Please sign in to comment.