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 the wrong planet position was fed into the tp acce…
Browse files Browse the repository at this point in the history 
…leration
  • Loading branch information
@daminton
daminton committed Apr 19, 2021
1 parent e80355f commit 7ee863e
Show file tree
Hide file tree
Showing 2 changed files with 126 additions and 124 deletions.
16 changes: 8 additions & 8 deletions src/rmvs/rmvs_getacch.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -23,7 +23,7 @@ module subroutine rmvs_getacch_tp(self, cb, pl, config, t, xh)


associate(tp => self, ntp => self%nbody, ipleP => self%ipleP, &
enc_index => self%index, cb_heliocentric => self%cb_heliocentric)
inner_index => self%index, cb_heliocentric => self%cb_heliocentric)

if (tp%lplanetocentric) then ! This is a close encounter step, so any accelerations requiring heliocentric position values
! must be handeled outside the normal WHM method call
Expand All @@ -39,39 +39,39 @@ module subroutine rmvs_getacch_tp(self, cb, pl, config, t, xh)
config_planetocen%lextra_force = .false.
config_planetocen%lgr = .false.
! Now compute the planetocentric values of acceleration
call whm_getacch_tp(tp, cb, pl, config_planetocen, t, pl%xh)
call whm_getacch_tp(tp, cb, pl, config_planetocen, t, xh)

! Now compute any heliocentric values of acceleration
if (tp%lfirst) then
do i = 1, ntp
tp%xheliocentric(:,i) = tp%xh(:,i) + cb%inner(enc_index - 1)%x(:,1)
tp%xheliocentric(:,i) = tp%xh(:,i) + cb%inner(inner_index - 1)%x(:,1)
end do
else
do i = 1, ntp
tp%xheliocentric(:,i) = tp%xh(:,i) + cb%inner(enc_index )%x(:,1)
tp%xheliocentric(:,i) = tp%xh(:,i) + cb%inner(inner_index )%x(:,1)
end do
end if
! Swap the planetocentric and heliocentric position vectors
tp%xh(:,:) = tp%xheliocentric(:,:)
if (config%loblatecb) then
! Put in the current encountering planet's oblateness acceleration as the central body's
if (tp%lfirst) then
cb_heliocentric%aobl(:) = cb%inner(enc_index - 1)%aobl(:,1)
cb_heliocentric%aobl(:) = cb%inner(inner_index - 1)%aobl(:,1)
else
cb_heliocentric%aobl(:) = cb%inner(enc_index )%aobl(:,1)
cb_heliocentric%aobl(:) = cb%inner(inner_index )%aobl(:,1)
end if
call tp%obl_acc(cb_heliocentric)
end if

if (config%lextra_force) call tp%user_getacch(cb_heliocentric, config, t)
if (config%lgr) call tp%gr_getacch(cb_heliocentric, config)

call move_alloc(xh_original, tp%xh)
tp%xh(:,:) = xh_original(:,:)
end associate
end select
end select
else ! Not a close encounter, so just proceded with the standard WHM method
call whm_getacch_tp(tp, cb, pl, config, t, pl%xh)
call whm_getacch_tp(tp, cb, pl, config, t, xh)
end if

end associate
Expand Down
234 changes: 118 additions & 116 deletions src/rmvs/rmvs_step.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -63,6 +63,124 @@ module subroutine rmvs_step_system(self, config)

end subroutine rmvs_step_system

subroutine rmvs_step_out(pl, cb, tp, dt, config)
!! author: David A. Minton
!!
!! Step ACTIVE test particles ahead in the outer encounter region, setting up and calling the inner region
!! integration if necessar
!!
!! Adapted from Hal Levison's Swift routines rmvs3_step_out.f and rmvs3_step_out2.f
!! Adapted from David E. Kaufmann's Swifter routines rmvs_step_out.f90 and rmvs_step_out2.f90
implicit none
! Arguments
class(rmvs_pl), intent(inout) :: pl !! RMVS massive body object
class(rmvs_cb), intent(inout) :: cb !! RMVS central body object
class(rmvs_tp), intent(inout) :: tp !! RMVS test particle object
real(DP), intent(in) :: dt !! Step size
class(swiftest_configuration), intent(in) :: config !! Input collection of configuration parameters
! Internals
integer(I4B) :: outer_index, j, k
real(DP) :: dto, outer_time, rts
logical :: lencounter, lfirsttp

associate(npl => pl%nbody, ntp => tp%nbody, t => config%t)
dto = dt / NTENC
where(tp%plencP(:) == 0)
tp%status(:) = INACTIVE
elsewhere
tp%lperi(:) = .false.
end where
do outer_index = 1, NTENC
outer_time = t + (outer_index - 1) * dto
call tp%set_beg_end(xbeg = pl%outer(outer_index - 1)%x(:, :), &
vbeg = pl%outer(outer_index - 1)%v(:, :), &
xend = pl%outer(outer_index )%x(:, :))
rts = RHPSCALE
lencounter = tp%encounter_check(cb, pl, dt, rts)
if (lencounter) then
associate(xbeg => pl%outer(outer_index - 1)%x(:, :), xend => pl%outer(outer_index )%x(:, :), &
tpxbeg => tp%xbeg, tpxend => tp%xend)
! Interpolate planets in inner encounter region
call rmvs_interp_in(pl, cb, dto, outer_index, config)
! Step through the inner region
call rmvs_step_in(pl, cb, tp, config, outer_time, dto)
lfirsttp = tp%lfirst
tp%lfirst = .true.
call tp%step(cb, pl, config, outer_time, dto)
tp%lfirst = lfirsttp
end associate
else
call tp%step(cb, pl, config, outer_time, dto)
end if
do j = 1, npl
if (pl%nenc(j) == 0) cycle
where((tp%plencP(:) == j) .and. (tp%status(:) == INACTIVE)) tp%status(:) = ACTIVE
end do
end do
end associate
return

end subroutine rmvs_step_out

subroutine rmvs_step_in(pl, cb, tp, config, outer_time, dto)
!! author: David A. Minton
!!
!! Step active test particles ahead in the inner encounter region
!!
!! Adapted from Hal Levison's Swift routine rmvs3_step_in.f
!! Adapted from David E. Kaufmann's Swifter routine rmvs_step_in.f90
implicit none
! Arguments
class(rmvs_pl), intent(inout) :: pl !! RMVS massive body object
class(rmvs_cb), intent(inout) :: cb !! RMVS central body object
class(rmvs_tp), intent(inout) :: tp !! RMVS test particle object
class(swiftest_configuration), intent(in) :: config !! Input collection of configuration parameters
real(DP), intent(in) :: outer_time !! Current time
real(DP), intent(in) :: dto !! Step size
! Internals
logical :: lfirsttp
integer(I4B) :: i, j, ipleP
real(DP) :: dti, inner_time
real(DP), dimension(:, :), allocatable :: xbeg, xend, vbeg

associate(npl => pl%nbody, nenc => pl%nenc)
dti = dto / NTPHENC
allocate(xbeg, mold=pl%xh)
allocate(xend, mold=pl%xh)
allocate(vbeg, mold=pl%vh)
if (config%loblatecb) call pl%obl_acc(cb)
call rmvs_make_planetocentric(pl, cb, tp, config)
do i = 1, npl
if (nenc(i) == 0) cycle
associate(cbenci => pl%planetocentric(i)%cb, plenci => pl%planetocentric(i)%pl, &
tpenci => pl%planetocentric(i)%tp)
associate(inner_index => tpenci%index)
! There are inner encounters with this planet...switch to planetocentric coordinates to proceed
tpenci%lfirst = .true.
inner_time = outer_time
call rmvs_peri_tp(tpenci, pl, inner_time, dti, .true., 0, i, config)
! now step the encountering test particles fully through the inner encounter
lfirsttp = .true.
do inner_index = 1, NTPHENC ! Integrate over the encounter region, using the "substitute" planetocentric systems at each level
plenci%xh(:,:) = plenci%inner(inner_index - 1)%x(:,:)
call tpenci%set_beg_end(xbeg = plenci%inner(inner_index - 1)%x, &
xend = plenci%inner(inner_index)%x)
call tpenci%step(cbenci, plenci, config, inner_time, dti)
do j = 1, nenc(i)
tpenci%xheliocentric(:, j) = tpenci%xh(:, j) + pl%inner(inner_index)%x(:,i)
end do
inner_time = outer_time + j * dti
call rmvs_peri_tp(tpenci, pl, inner_time, dti, .false., inner_index, i, config)
end do
where(tpenci%status(:) == ACTIVE) tpenci%status(:) = INACTIVE
end associate
end associate
end do
call rmvs_end_planetocentric(pl, cb, tp)
end associate
return
end subroutine rmvs_step_in

subroutine rmvs_interp_in(pl, cb, dt, outer_index, config)
!! author: David A. Minton
!!
Expand Down Expand Up @@ -332,122 +450,6 @@ subroutine rmvs_peri_tp(tp, pl, t, dt, lfirst, inner_index, ipleP, config)

end subroutine rmvs_peri_tp

subroutine rmvs_step_out(pl, cb, tp, dt, config)
!! author: David A. Minton
!!
!! Step ACTIVE test particles ahead in the outer encounter region, setting up and calling the inner region
!! integration if necessar
!!
!! Adapted from Hal Levison's Swift routines rmvs3_step_out.f and rmvs3_step_out2.f
!! Adapted from David E. Kaufmann's Swifter routines rmvs_step_out.f90 and rmvs_step_out2.f90
implicit none
! Arguments
class(rmvs_pl), intent(inout) :: pl !! RMVS massive body object
class(rmvs_cb), intent(inout) :: cb !! RMVS central body object
class(rmvs_tp), intent(inout) :: tp !! RMVS test particle object
real(DP), intent(in) :: dt !! Step size
class(swiftest_configuration), intent(in) :: config !! Input collection of configuration parameters
! Internals
integer(I4B) :: outer_index, j, k
real(DP) :: dto, outer_time, rts
logical :: lencounter, lfirsttp

associate(npl => pl%nbody, ntp => tp%nbody, t => config%t)
dto = dt / NTENC
where(tp%plencP(:) == 0)
tp%status(:) = INACTIVE
elsewhere
tp%lperi(:) = .false.
end where
do outer_index = 1, NTENC
outer_time = t + (outer_index - 1) * dto
call tp%set_beg_end(xbeg = pl%outer(outer_index - 1)%x(:, :), &
vbeg = pl%outer(outer_index - 1)%v(:, :), &
xend = pl%outer(outer_index )%x(:, :))
rts = RHPSCALE
lencounter = tp%encounter_check(cb, pl, dt, rts)
if (lencounter) then
! Interpolate planets in inner encounter region
call rmvs_interp_in(pl, cb, dto, outer_index, config)
! Step through the inner region
call rmvs_step_in(pl, cb, tp, config, outer_time, dto)
lfirsttp = tp%lfirst
tp%lfirst = .true.
call tp%step(cb, pl, config, outer_time, dto)
tp%lfirst = lfirsttp
else
call tp%step(cb, pl, config, outer_time, dto)
end if
do j = 1, npl
if (pl%nenc(j) == 0) cycle
where((tp%plencP(:) == j) .and. (tp%status(:) == INACTIVE)) tp%status(:) = ACTIVE
end do
end do
end associate
return

end subroutine rmvs_step_out

subroutine rmvs_step_in(pl, cb, tp, config, outer_time, dto)
!! author: David A. Minton
!!
!! Step active test particles ahead in the inner encounter region
!!
!! Adapted from Hal Levison's Swift routine rmvs3_step_in.f
!! Adapted from David E. Kaufmann's Swifter routine rmvs_step_in.f90
implicit none
! Arguments
class(rmvs_pl), intent(inout) :: pl !! RMVS massive body object
class(rmvs_cb), intent(inout) :: cb !! RMVS central body object
class(rmvs_tp), intent(inout) :: tp !! RMVS test particle object
class(swiftest_configuration), intent(in) :: config !! Input collection of configuration parameters
real(DP), intent(in) :: outer_time !! Current time
real(DP), intent(in) :: dto !! Step size
! Internals
logical :: lfirsttp
integer(I4B) :: i, j, ipleP
real(DP) :: dti, inner_time
real(DP), dimension(:, :), allocatable :: xbeg, xend, vbeg


associate(npl => pl%nbody, nenc => pl%nenc)
dti = dto / NTPHENC
allocate(xbeg, mold=pl%xh)
allocate(xend, mold=pl%xh)
allocate(vbeg, mold=pl%vh)
if (config%loblatecb) call pl%obl_acc(cb)
call rmvs_make_planetocentric(pl, cb, tp, config)
do i = 1, npl
if (nenc(i) == 0) cycle
associate(cbenci => pl%planetocentric(i)%cb, plenci => pl%planetocentric(i)%pl, &
tpenci => pl%planetocentric(i)%tp)
associate(inner_index => tpenci%index)
! There are inner encounters with this planet...switch to planetocentric coordinates to proceed
tpenci%lfirst = .true.
inner_time = outer_time
call rmvs_peri_tp(tpenci, pl, inner_time, dti, .true., 0, i, config)
! now step the encountering test particles fully through the inner encounter
lfirsttp = .true.
do inner_index = 1, NTPHENC ! Integrate over the encounter region, using the "substitute" planetocentric systems at each level
plenci%xh(:,:) = plenci%inner(inner_index - 1)%x(:,:)
call tpenci%set_beg_end(xbeg = plenci%inner(inner_index - 1)%x, &
xend = plenci%inner(inner_index)%x)
call tpenci%step(cbenci, plenci, config, inner_time, dti)
do j = 1, nenc(i)
tpenci%xheliocentric(:, j) = tpenci%xh(:, j) + pl%inner(inner_index)%x(:,i)
end do
inner_time = outer_time + j * dti
call rmvs_peri_tp(tpenci, pl, inner_time, dti, .false., inner_index, i, config)
end do
where(tpenci%status(:) == ACTIVE) tpenci%status(:) = INACTIVE
end associate
end associate
end do
call rmvs_end_planetocentric(pl, cb, tp)
end associate
return
end subroutine rmvs_step_in

subroutine rmvs_make_planetocentric(pl, cb, tp, config)
!! author: David A. Minton
!!
Expand Down

0 comments on commit 7ee863e

Please sign in to comment.