diff --git a/src/modules/rmvs_classes.f90 b/src/modules/rmvs_classes.f90 index 5a7895442..70330f1b7 100644 --- a/src/modules/rmvs_classes.f90 +++ b/src/modules/rmvs_classes.f90 @@ -254,7 +254,7 @@ module function rmvs_encounter_check_tp(self, cb, pl, dt, rts) result(lencounter end function rmvs_encounter_check_tp - module subroutine rmvs_peri_tp(self, cb, pl, t, dt, lfirst, index, nenc, ipleP, config) + module subroutine rmvs_peri_tp(self, cb, pl, t, dt, lfirst, index, ipleP, config) use swiftest_classes, only : swiftest_configuration class(rmvs_tp), intent(inout) :: self !! RMVS test particle object class(rmvs_cb), intent(inout) :: cb !! RMVS central body object @@ -263,7 +263,6 @@ module subroutine rmvs_peri_tp(self, cb, pl, t, dt, lfirst, index, nenc, ipleP, real(DP), intent(in) :: dt !! step size logical, intent(in) :: lfirst !! Logical flag indicating whether current invocation is the first integer(I4B), intent(in) :: index !! outer substep number within current set - integer(I4B), intent(in) :: nenc !! number of test particles encountering current planet integer(I4B), intent(in) :: ipleP !!index of RMVS planet being closely encountered class(swiftest_configuration), intent(in) :: config !! Input collection of configuration parameters end subroutine rmvs_peri_tp diff --git a/src/rmvs/rmvs_encounter_check.f90 b/src/rmvs/rmvs_encounter_check.f90 index c43415e16..ebfc2095c 100644 --- a/src/rmvs/rmvs_encounter_check.f90 +++ b/src/rmvs/rmvs_encounter_check.f90 @@ -30,16 +30,9 @@ module function rmvs_encounter_check_tp(self, cb, pl, dt, rts) result(lencounter pl%encmask(:,:) = .false. lencounter = .false. pl%nenc(:) = 0 - pl%tpenc1P(:) = 0 - ! if first time through, calc hill's sphere for the planets - !if (lfirst) then - - ! lfirst = .false. - !end if do i = 1, ntp if (tp%status(i) == ACTIVE) then tp%plencP(i) = 0 - !tp%tpencP(i) = 0 lflag = .false. do j = 1, npl @@ -51,16 +44,6 @@ module function rmvs_encounter_check_tp(self, cb, pl, dt, rts) result(lencounter lencounter = .true. pl%encmask(i,j) = .true. pl%nenc(j) = pl%nenc(j) + 1 - !nenc = pl%nenc(j) - !if (nenc == 1) then - ! pl%tpenc1P(j) = i - !else - ! tpencPindex = pl%tpenc1P(j) - ! do k = 2, nenc - 1 - ! tpencPindex = tp%tpencP(tpencPindex) - ! end do - ! tp%tpencP(tpencPindex) = i - !end if tp%plencP(i) = j exit end if diff --git a/src/rmvs/rmvs_peri.f90 b/src/rmvs/rmvs_peri.f90 index 7265c4693..4f8cf3612 100644 --- a/src/rmvs/rmvs_peri.f90 +++ b/src/rmvs/rmvs_peri.f90 @@ -1,7 +1,7 @@ submodule(rmvs_classes) s_rmvs_peri use swiftest contains - module subroutine rmvs_peri_tp(self, cb, pl, t, dt, lfirst, index, nenc, ipleP, config) + module subroutine rmvs_peri_tp(self, cb, pl, t, dt, lfirst, index, ipleP, config) !! author: David A. Minton !! !! Determine planetocentric pericenter passages for test particles in close encounters with a planet @@ -17,7 +17,6 @@ module subroutine rmvs_peri_tp(self, cb, pl, t, dt, lfirst, index, nenc, ipleP, real(DP), intent(in) :: dt !! step size logical, intent(in) :: lfirst !! Logical flag indicating whether current invocation is the first integer(I4B), intent(in) :: index !! outer substep number within current set - integer(I4B), intent(in) :: nenc !! number of test particles encountering current planet integer(I4B), intent(in) :: ipleP !!index of RMVS planet being closely encountered class(swiftest_configuration), intent(in) :: config !! Input collection of configuration parameters ! Internals diff --git a/src/rmvs/rmvs_setup.f90 b/src/rmvs/rmvs_setup.f90 index 48d9bcb87..2de07cc78 100644 --- a/src/rmvs/rmvs_setup.f90 +++ b/src/rmvs/rmvs_setup.f90 @@ -19,7 +19,6 @@ module subroutine rmvs_setup_pl(self,n) if (n <= 0) return allocate(self%nenc(n)) - allocate(self%tpenc1P(n)) allocate(self%xout(NDIM, n, 0:NTENC)) allocate(self%vout(NDIM, n, 0:NTENC)) allocate(self%xin(NDIM, n, 0:NTPHENC)) @@ -32,7 +31,6 @@ module subroutine rmvs_setup_pl(self,n) self%plenc(:)%nbody = n self%nenc = 0 - self%tpenc1P(:) = 0 self%xout(:,:,:) = 0.0_DP self%vout(:,:,:) = 0.0_DP self%xin(:,:,:) = 0.0_DP @@ -65,12 +63,10 @@ module subroutine rmvs_setup_tp(self,n) allocate(self%lperi(n)) allocate(self%plperP(n)) allocate(self%plencP(n)) - allocate(self%tpencP(n)) self%lperi(:) = .false. self%plperP(:) = 0 self%plencP(:) = 0 - self%tpencP(:) = 0 return end subroutine rmvs_setup_tp diff --git a/src/rmvs/rmvs_spill_and_fill.f90 b/src/rmvs/rmvs_spill_and_fill.f90 index 6369375e0..d072e0274 100644 --- a/src/rmvs/rmvs_spill_and_fill.f90 +++ b/src/rmvs/rmvs_spill_and_fill.f90 @@ -20,10 +20,8 @@ module subroutine rmvs_spill_pl(self, discards, lspill_list) class is (rmvs_pl) discards%nenc(:) = pack(keeps%nenc(:), lspill_list(:)) - discards%tpenc1P(:) = pack(keeps%tpenc1P(:), lspill_list(:)) if (count(.not.lspill_list(:)) > 0) then keeps%nenc(:) = pack(keeps%nenc(:), .not. lspill_list(:)) - keeps%tpenc1P(:) = pack(keeps%tpenc1P(:), .not. lspill_list(:)) end if call whm_spill_pl(keeps, discards, lspill_list) class default @@ -56,9 +54,6 @@ module subroutine rmvs_fill_pl(self, inserts, lfill_list) keeps%nenc(:) = unpack(keeps%nenc(:), .not.lfill_list(:), keeps%nenc(:)) keeps%nenc(:) = unpack(inserts%nenc(:), lfill_list(:), keeps%nenc(:)) - keeps%tpenc1P(:) = unpack(keeps%tpenc1P(:), .not.lfill_list(:), keeps%tpenc1P(:)) - keeps%tpenc1P(:) = unpack(inserts%tpenc1P(:), lfill_list(:), keeps%tpenc1P(:)) - call whm_fill_pl(keeps, inserts, lfill_list) class default write(*,*) 'Error! spill method called for incompatible return type on rmvs_pl' @@ -89,12 +84,10 @@ module subroutine rmvs_spill_tp(self, discards, lspill_list) discards%lperi(:) = pack(keeps%lperi(:), lspill_list(:)) discards%plperP(:) = pack(keeps%plperP(:), lspill_list(:)) discards%plencP(:) = pack(keeps%plencP(:), lspill_list(:)) - discards%tpencP(:) = pack(keeps%tpencP(:), lspill_list(:)) if (count(.not.lspill_list(:)) > 0) then keeps%lperi(:) = pack(keeps%lperi(:), .not. lspill_list(:)) keeps%plperP(:) = pack(keeps%plperP(:), .not. lspill_list(:)) keeps%plencP(:) = pack(keeps%plencP(:), .not. lspill_list(:)) - keeps%tpencP(:) = pack(keeps%tpencP(:), .not. lspill_list(:)) end if call util_spill_tp(keeps, discards, lspill_list) @@ -132,9 +125,6 @@ module subroutine rmvs_fill_tp(self, inserts, lfill_list) keeps%plencP(:) = unpack(keeps%plencP(:), .not.lfill_list(:), keeps%plencP(:)) keeps%plencP(:) = unpack(inserts%plencP(:), lfill_list(:), keeps%plencP(:)) - keeps%tpencP(:) = unpack(keeps%tpencP(:), .not.lfill_list(:), keeps%tpencP(:)) - keeps%tpencP(:) = unpack(inserts%tpencP(:), lfill_list(:), keeps%tpencP(:)) - call util_fill_tp(keeps, inserts, lfill_list) class default write(*,*) 'Error! fill method called for incompatible return type on rmvs_tp' diff --git a/src/rmvs/rmvs_step.f90 b/src/rmvs/rmvs_step.f90 index 25465606d..fb5fb3954 100644 --- a/src/rmvs/rmvs_step.f90 +++ b/src/rmvs/rmvs_step.f90 @@ -70,7 +70,7 @@ module subroutine rmvs_step_out(self, cb, tp, dt, config) real(DP), intent(in) :: dt !! Step size class(swiftest_configuration), intent(in) :: config !! Input collection of configuration parameters ! Internals - integer(I4B) :: i, j, k, nenc, itpp + integer(I4B) :: i, j, k real(DP) :: dto, time ! executable code @@ -161,48 +161,46 @@ module subroutine rmvs_step_in_pl(self, cb, tp, config, t, dt) real(DP), intent(in) :: dt !! Step size ! Internals logical :: lfirsttp - integer(I4B) :: i, j, nenc + integer(I4B) :: i, j real(DP) :: dti, time real(DP), dimension(NDIM, self%nbody) :: xbeg, xend, vbeg dti = dt / NTPHENC - associate(pl => self, npl => self%nbody, xht => tp%xh, vht => tp%vh, plind => self%plind) + associate(pl => self, npl => self%nbody, xht => tp%xh, vht => tp%vh, plind => self%plind, nenc => pl%nenc) if (config%loblatecb) call pl%obl_acc_in(cb) call pl%make_planetocentric(cb, tp, config) do i = 1, npl - nenc = pl%nenc(i) - if (nenc > 0) then + if (nenc(i) == 0) cycle ! There are inner encounters with this planet...switch to planetocentric coordinates to proceed - time = t - call pl%tpenc(i)%peri_pass(cb, pl, time, dti, .true., 0, nenc, i, config) - ! now step the encountering test particles fully through the inner encounter - lfirsttp = .true. - associate(index => pl%tpenc(i)%index, & - xpc => self%tpenc(i)%xh, vpc => self%tpenc(i)%vh, apc => self%tpenc(i)%ah) - pl%tpenc(i)%lfirst = .true. - do index = 1, NTPHENC ! Integrate over the encounter region, using the "substitute" planetocentric systems at each level - xbeg(:,1) = cb%xin(:) - pl%xin(:, i, index - 1) - xend(:,1) = cb%xin(:) - pl%xin(:, i, index) - vbeg(:,1) = cb%vin(:) - pl%vin(:, i, index - 1) - do j = 1, NDIM - xbeg(j,2:npl) = pl%xin(j, plind(i,2:npl), index - 1) - pl%xin(j, i, index - 1) - xend(j,2:npl) = pl%xin(j, plind(i,2:npl), index) - pl%xin(j, i, index) - vbeg(j,2:npl) = pl%vin(j, plind(i,2:npl), index - 1) - pl%vin(j, i, index - 1) - end do - pl%plenc(i)%xh(:,:) = xbeg(:,:) - pl%plenc(i)%vh(:,:) = vbeg(:,:) - call pl%tpenc(i)%set_beg_end(xbeg = xbeg, xend = xend) - call pl%tpenc(i)%step(pl%cbenc(i), pl%plenc(i), config, time, dti) - do j = 1, NDIM - pl%tpenc(i)%xheliocen(j, :) = pl%tpenc(i)%xh(j, :) + pl%xin(j, i, index) - end do - time = config%t + j * dti - call pl%tpenc(i)%peri_pass(cb, pl, time, dti, .false., index, nenc, i, config) + time = t + call pl%tpenc(i)%peri_pass(cb, pl, time, dti, .true., 0, i, config) + ! now step the encountering test particles fully through the inner encounter + lfirsttp = .true. + associate(index => pl%tpenc(i)%index, & + xpc => self%tpenc(i)%xh, vpc => self%tpenc(i)%vh, apc => self%tpenc(i)%ah) + pl%tpenc(i)%lfirst = .true. + do index = 1, NTPHENC ! Integrate over the encounter region, using the "substitute" planetocentric systems at each level + xbeg(:,1) = cb%xin(:) - pl%xin(:, i, index - 1) + xend(:,1) = cb%xin(:) - pl%xin(:, i, index) + vbeg(:,1) = cb%vin(:) - pl%vin(:, i, index - 1) + do j = 1, NDIM + xbeg(j,2:npl) = pl%xin(j, plind(i,2:npl), index - 1) - pl%xin(j, i, index - 1) + xend(j,2:npl) = pl%xin(j, plind(i,2:npl), index) - pl%xin(j, i, index) + vbeg(j,2:npl) = pl%vin(j, plind(i,2:npl), index - 1) - pl%vin(j, i, index - 1) end do - where(pl%tpenc(i)%status(:) == ACTIVE) pl%tpenc(i)%status(:) = INACTIVE - end associate - end if + pl%plenc(i)%xh(:,:) = xbeg(:,:) + pl%plenc(i)%vh(:,:) = vbeg(:,:) + call pl%tpenc(i)%set_beg_end(xbeg = xbeg, xend = xend) + call pl%tpenc(i)%step(pl%cbenc(i), pl%plenc(i), config, time, dti) + do j = 1, NDIM + pl%tpenc(i)%xheliocen(j, :) = pl%tpenc(i)%xh(j, :) + pl%xin(j, i, index) + end do + time = config%t + j * dti + call pl%tpenc(i)%peri_pass(cb, pl, time, dti, .false., index, i, config) + end do + where(pl%tpenc(i)%status(:) == ACTIVE) pl%tpenc(i)%status(:) = INACTIVE + end associate end do call pl%end_planetocentric(cb,tp) @@ -225,12 +223,8 @@ module subroutine rmvs_step_make_planetocentric(self, cb, tp, config) class(rmvs_tp), intent(inout) :: tp !! RMVS test particle object class(swiftest_configuration), intent(in) :: config !! Input collection of configuration parameters ! Internals - integer(I4B) :: i, j, k - type(rmvs_pl) :: cb_as_pl - logical, dimension(:), allocatable :: copyflag - + integer(I4B) :: i, j - !allocate(copyflag(self%nbody)) associate(pl => self, npl => self%nbody, nenc => self%nenc, tpenc => self%tpenc, cbenc => self%cbenc, & plenc => self%plenc, GMpl => self%Gmass) @@ -257,6 +251,8 @@ module subroutine rmvs_step_make_planetocentric(self, cb, tp, config) allocate(tpenc(i)%peri(nenc(i))) tpenc(i)%status(:) = ACTIVE tpenc(i)%lperi(:) = .false. + tpenc(i)%plperP(:) = 0 + tpenc(i)%name(:) = pack(tp%name(:), pl%encmask(:,i)) !call tp%spill(tpenc(i), pl%encmask(:,i)) ! Grab all the encountering test particles and convert them to a planetocentric frame do j = 1, NDIM @@ -277,6 +273,7 @@ module subroutine rmvs_step_make_planetocentric(self, cb, tp, config) end if end do end associate + return end subroutine rmvs_step_make_planetocentric module subroutine rmvs_step_end_planetocentric(self, cb, tp) @@ -324,9 +321,7 @@ module subroutine rmvs_step_end_planetocentric(self, cb, tp) end associate deallocate(self%encmask) - return - end subroutine rmvs_step_end_planetocentric end submodule s_rmvs_step \ No newline at end of file