From 7ee863ea338c6fbbf302ef8e95fff6168bf224fc Mon Sep 17 00:00:00 2001 From: David Minton Date: Mon, 19 Apr 2021 14:46:31 -0400 Subject: [PATCH] Fixed bug in which the wrong planet position was fed into the tp acceleration --- src/rmvs/rmvs_getacch.f90 | 16 +-- src/rmvs/rmvs_step.f90 | 234 +++++++++++++++++++------------------- 2 files changed, 126 insertions(+), 124 deletions(-) diff --git a/src/rmvs/rmvs_getacch.f90 b/src/rmvs/rmvs_getacch.f90 index 21e80a716..550046131 100644 --- a/src/rmvs/rmvs_getacch.f90 +++ b/src/rmvs/rmvs_getacch.f90 @@ -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 @@ -39,16 +39,16 @@ 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 @@ -56,9 +56,9 @@ module subroutine rmvs_getacch_tp(self, cb, pl, config, t, xh) 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 @@ -66,12 +66,12 @@ module subroutine rmvs_getacch_tp(self, cb, pl, config, t, xh) 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 diff --git a/src/rmvs/rmvs_step.f90 b/src/rmvs/rmvs_step.f90 index ae4712efb..d36c701fd 100644 --- a/src/rmvs/rmvs_step.f90 +++ b/src/rmvs/rmvs_step.f90 @@ -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 !! @@ -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 !!