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

Commit

Permalink
Re-wrote array operations on encountering particles
Browse files Browse the repository at this point in the history
  • Loading branch information
@daminton
daminton committed Apr 2, 2021
1 parent 893407a commit d890867
Show file tree
Hide file tree
Showing 5 changed files with 164 additions and 125 deletions.
3 changes: 2 additions & 1 deletion src/modules/rmvs_classes.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -74,10 +74,11 @@ module rmvs_classes
real(DP), dimension(:, :, :), allocatable :: xin !! interpolated heliocentric planet position for inner encounter
real(DP), dimension(:, :, :), allocatable :: vin !! interpolated heliocentric planet velocity for inner encounter
type(rmvs_tp), dimension(:), allocatable :: tpenc !! array of encountering test particles with this planet in planetocentric coordinates
type(whm_pl) , dimension(:, :), allocatable :: plenc !! array of massive bodies that includes the Sun, but not the encountering planet in planetocentric coordinates
type(whm_pl) , dimension(:), allocatable :: plenc !! array of massive bodies that includes the Sun, but not the encountering planet in planetocentric coordinates
type(rmvs_cb), dimension(:), allocatable :: cbenc !! The planet acting as a central body for close encounters
real(DP), dimension(:, :, :), allocatable :: aoblin !! barycentric acceleration on planets due to central body oblateness during inner encounter
logical, dimension(:, :), allocatable :: encmask !! logical mask indicating which test particles are encountering a particular planet
integer(I4B), dimension(:,:), allocatable :: plind
!! Note to developers: If you add componenets to this class, be sure to update methods and subroutines that traverse the
!! component list, such as rmvs_setup_pl and rmvs_spill_pl
contains
Expand Down
35 changes: 17 additions & 18 deletions src/rmvs/rmvs_encounter_check.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -19,12 +19,13 @@ module function rmvs_encounter_check_tp(self, cb, pl, dt, rts) result(lencounter
! Internals
integer(I4B) :: i, j, k, nenc
real(DP) :: r2crit
real(DP), dimension(NDIM) :: xht, vht, xr, vr
real(DP), dimension(NDIM) :: xr, vr
integer(I4B) :: tpencPindex
logical :: lflag
logical, save :: lfirst = .true.

associate(tp => self, ntp => self%nbody, npl => pl%nbody, rhill => pl%rhill)
associate(tp => self, ntp => self%nbody, npl => pl%nbody, rhill => pl%rhill, &
xht => self%xh, vht => self%vh, xbeg => self%xbeg, vbeg => self%vbeg)
if (.not.allocated(pl%encmask)) allocate(pl%encmask(ntp, npl))
pl%encmask(:,:) = .false.
lencounter = .false.
Expand All @@ -38,30 +39,28 @@ module function rmvs_encounter_check_tp(self, cb, pl, dt, rts) result(lencounter
do i = 1, ntp
if (tp%status(i) == ACTIVE) then
tp%plencP(i) = 0
tp%tpencP(i) = 0
!tp%tpencP(i) = 0
lflag = .false.
xht(:) = tp%xh(:, i)
vht(:) = tp%vh(:, i)

do j = 1, npl
r2crit = (rts * pl%rhill(j))**2
xr(:) = xht(:) - tp%xbeg(:, j)
vr(:) = vht(:) - tp%vbeg(:, j)
r2crit = (rts * rhill(j))**2
xr(:) = xht(:, i) - xbeg(:, j)
vr(:) = vht(:, i) - vbeg(:, j)
lflag = rmvs_chk_ind(xr(:), vr(:), dt, r2crit)
if (lflag) then
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
!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
Expand Down
34 changes: 25 additions & 9 deletions src/rmvs/rmvs_getacch.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -32,18 +32,34 @@ module subroutine rmvs_getacch_in_tp(self, cb, pl, config, t, xh)
config_planetocen%loblatecb = .false.
config_planetocen%lextra_force = .false.
config_planetocen%lgr = .false.
! Now compute the heliocentric values of acceleration
select type(pl)
class is (rmvs_pl)
if (tp%lfirst) then
do i = 1, NDIM
tp%xheliocen(i,:) = tp%xh(i,:) + pl%xin(i,ipleP,index - 1)
end do
else
do i = 1, NDIM
tp%xheliocen(i,:) = tp%xh(i,:) + pl%xin(i,ipleP,index)
end do
end if
end select
call whm_getacch_tp(tp, cb, pl, config_planetocen, t, xh)

! Now compute the heliocentric values of acceleration
if (tp%lfirst) then
do i = 1, NDIM
tp%xheliocen(i,:) = tp%xh(i,:) + tp%xh_pl(i, index - 1)
end do
else
do i = 1, NDIM
tp%xheliocen(i,:) = tp%xh(i,:) + tp%xh_pl(i, index)
end do
end if
select type(pl)
class is (rmvs_pl)
if (tp%lfirst) then
do i = 1, NDIM
tp%xheliocen(i,:) = tp%xh(i,:) + pl%xin(i,ipleP,index - 1)
end do
else
do i = 1, NDIM
tp%xheliocen(i,:) = tp%xh(i,:) + pl%xin(i,ipleP,index)
end do
end if
end select

tp%xh(:,:) = tp%xheliocen(:,:)
if (config%loblatecb) then
Expand Down
39 changes: 37 additions & 2 deletions src/rmvs/rmvs_setup.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -10,7 +10,9 @@ module subroutine rmvs_setup_pl(self,n)
implicit none
! Arguments
class(rmvs_pl), intent(inout) :: self !! RMVS test particle object
integer, intent(in) :: n !! Number of test particles to allocate
integer(I4B), intent(in) :: n !! Number of test particles to allocate
! Internals
integer(I4B) :: i,j

!> Call allocation method for parent class
call whm_setup_pl(self, n)
Expand All @@ -23,6 +25,11 @@ module subroutine rmvs_setup_pl(self,n)
allocate(self%xin(NDIM, n, 0:NTPHENC))
allocate(self%vin(NDIM, n, 0:NTPHENC))
allocate(self%aoblin(NDIM, n, 0:NTPHENC))
allocate(self%plind(n,n))
allocate(self%tpenc(n))
allocate(self%plenc(n))
allocate(self%cbenc(n))
self%plenc(:)%nbody = n

self%nenc = 0
self%tpenc1P(:) = 0
Expand All @@ -32,6 +39,11 @@ module subroutine rmvs_setup_pl(self,n)
self%vin(:,:,:) = 0.0_DP
self%aoblin(:,:,:) = 0.0_DP

do j = 1, n
self%plind(j,:) = [(i,i=1,n)]
self%plind(j,2:n) = pack(self%plind(j,1:n), self%plind(j,1:n) /= j)
end do

return
end subroutine rmvs_setup_pl

Expand Down Expand Up @@ -72,10 +84,33 @@ module subroutine rmvs_setup_system(self, config)
! Arguments
class(rmvs_nbody_system), intent(inout) :: self !! RMVS system object
class(swiftest_configuration), intent(inout) :: config !! Input collection of configuration parameters

! Internals
integer(I4B) :: i
! Call parent method
call whm_setup_system(self, config)

! Set up the tp-planet encounter structures
select type(pl => self%pl)
class is(rmvs_pl)
select type(cb => self%cb)
class is (rmvs_cb)
select type (tp => self%tp)
class is (rmvs_tp)
associate(npl => pl%nbody)
tp%cb = cb
do i = 1, npl
allocate(pl%plenc(i)%Gmass(npl))
allocate(pl%plenc(i)%xh(NDIM,npl))
allocate(pl%plenc(i)%vh(NDIM,npl))
pl%cbenc(i) = cb
pl%cbenc(i)%Gmass = pl%Gmass(i)
pl%plenc(i)%Gmass(1) = cb%Gmass
end do
end associate
end select
end select
end select

end subroutine rmvs_setup_system

module subroutine rmvs_setup_set_beg_end(self, xbeg, xend, vbeg)
Expand Down
Loading

0 comments on commit d890867

Please sign in to comment.