From 4738e4be6ffb05e0821c5a22ad468dc6288ff21b Mon Sep 17 00:00:00 2001 From: David A Minton Date: Wed, 12 Oct 2022 07:15:15 -0400 Subject: [PATCH] Fixed bug in RMVS that accidentally advanced system by two steps if there were no test particles --- src/rmvs/rmvs_step.f90 | 77 ++++++++++++++++++++++-------------------- 1 file changed, 40 insertions(+), 37 deletions(-) diff --git a/src/rmvs/rmvs_step.f90 b/src/rmvs/rmvs_step.f90 index 76c5a802c..2b7e9045a 100644 --- a/src/rmvs/rmvs_step.f90 +++ b/src/rmvs/rmvs_step.f90 @@ -21,45 +21,48 @@ module subroutine rmvs_step_system(self, param, t, dt) real(DP), dimension(:,:), allocatable :: xbeg, xend, vbeg integer(I4B) :: i - if (self%tp%nbody == 0) call whm_step_system(self, param, t, dt) - - select type(cb => self%cb) - class is (rmvs_cb) - select type(pl => self%pl) - class is (rmvs_pl) - select type(tp => self%tp) - class is (rmvs_tp) - associate(system => self, ntp => tp%nbody, npl => pl%nbody) - allocate(xbeg, source=pl%xh) - allocate(vbeg, source=pl%vh) - call pl%set_beg_end(xbeg = xbeg, vbeg = vbeg) - ! ****** Check for close encounters ***** ! - call pl%set_renc(RHSCALE) - lencounter = tp%encounter_check(param, system, dt) - if (lencounter) then - lfirstpl = pl%lfirst - pl%outer(0)%x(:, 1:npl) = xbeg(:, 1:npl) - pl%outer(0)%v(:, 1:npl) = vbeg(:, 1:npl) - call pl%step(system, param, t, dt) - pl%outer(NTENC)%x(:, 1:npl) = pl%xh(:, 1:npl) - pl%outer(NTENC)%v(:, 1:npl) = pl%vh(:, 1:npl) - call rmvs_interp_out(cb, pl, dt) - call rmvs_step_out(cb, pl, tp, system, param, t, dt) - tp%lmask(1:ntp) = .not. tp%lmask(1:ntp) - call pl%set_beg_end(xbeg = xbeg, xend = xend) - tp%lfirst = .true. - call tp%step(system, param, t, dt) - tp%lmask(1:ntp) = .true. - pl%lfirst = lfirstpl - tp%lfirst = .true. - if (param%ltides) call system%step_spin(param, t, dt) - else - call whm_step_system(system, param, t, dt) - end if - end associate + if (self%tp%nbody == 0) then + call whm_step_system(self, param, t, dt) + else + select type(cb => self%cb) + class is (rmvs_cb) + select type(pl => self%pl) + class is (rmvs_pl) + select type(tp => self%tp) + class is (rmvs_tp) + associate(system => self, ntp => tp%nbody, npl => pl%nbody) + allocate(xbeg, source=pl%xh) + allocate(vbeg, source=pl%vh) + call pl%set_beg_end(xbeg = xbeg, vbeg = vbeg) + ! ****** Check for close encounters ***** ! + call pl%set_renc(RHSCALE) + lencounter = tp%encounter_check(param, system, dt) + if (lencounter) then + lfirstpl = pl%lfirst + pl%outer(0)%x(:, 1:npl) = xbeg(:, 1:npl) + pl%outer(0)%v(:, 1:npl) = vbeg(:, 1:npl) + call pl%step(system, param, t, dt) + pl%outer(NTENC)%x(:, 1:npl) = pl%xh(:, 1:npl) + pl%outer(NTENC)%v(:, 1:npl) = pl%vh(:, 1:npl) + call rmvs_interp_out(cb, pl, dt) + call rmvs_step_out(cb, pl, tp, system, param, t, dt) + tp%lmask(1:ntp) = .not. tp%lmask(1:ntp) + call pl%set_beg_end(xbeg = xbeg, xend = xend) + tp%lfirst = .true. + call tp%step(system, param, t, dt) + tp%lmask(1:ntp) = .true. + pl%lfirst = lfirstpl + tp%lfirst = .true. + if (param%ltides) call system%step_spin(param, t, dt) + else + call whm_step_system(system, param, t, dt) + end if + end associate + end select end select end select - end select + end if + return end subroutine rmvs_step_system