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

Commit

Permalink
Fixed bug in RMVS that accidentally advanced system by two steps if t…
Browse files Browse the repository at this point in the history
…here were no test particles
  • Loading branch information
daminton committed Oct 12, 2022
1 parent 1bd4d53 commit 4738e4b
Showing 1 changed file with 40 additions and 37 deletions.
77 changes: 40 additions & 37 deletions src/rmvs/rmvs_step.f90
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down

0 comments on commit 4738e4b

Please sign in to comment.