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

Commit

Permalink
Browse files Browse the repository at this point in the history 
Fixed bugs in data structures for planetocentric encounters
  • Loading branch information
@daminton
daminton committed Apr 16, 2021
1 parent 8cfa099 commit b9838f0
Show file tree
Hide file tree
Showing 5 changed files with 62 additions and 85 deletions.
39 changes: 19 additions & 20 deletions src/modules/rmvs_classes.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -16,13 +16,6 @@ module rmvs_classes
real(DP), parameter :: RHPSCALE = 1.0_DP
real(DP), parameter :: FACQDT = 2.0_DP

type, private :: rmvs_interp
real(DP), dimension(:, :), allocatable :: x !! interpolated heliocentric planet position for outer encounter
real(DP), dimension(:, :), allocatable :: v !! interpolated heliocentric planet velocity for outer encounter
real(DP), dimension(:, :), allocatable :: aobl !! Encountering planet's oblateness acceleration value
end type rmvs_interp


!********************************************************************************************************************************
! rmvs_nbody_system class definitions and method interfaces
!********************************************************************************************************************************
Expand Down Expand Up @@ -81,10 +74,24 @@ module rmvs_classes
! rmvs_pl class definitions and method interfaces
!*******************************************************************************************************************************

type, private :: rmvs_interp
real(DP), dimension(:, :), allocatable :: x !! interpolated heliocentric planet position for outer encounter
real(DP), dimension(:, :), allocatable :: v !! interpolated heliocentric planet velocity for outer encounter
real(DP), dimension(:, :), allocatable :: aobl !! Encountering planet's oblateness acceleration value
end type rmvs_interp

!> RMVS massive body particle class
type, private, extends(whm_pl) :: rmvs_base_pl
logical :: lplanetocentric !! Flag that indicates that the object is a planetocentric set of masive bodies used for close encounter calculations
integer(I4B), dimension(:), allocatable :: plind ! Connects the planetocentric indices back to the heliocentric planet list
logical :: lplanetocentric !! Flag that indicates that the object is a planetocentric set of masive bodies used for close encounter calculations
integer(I4B), dimension(:), allocatable :: nenc !! number of test particles encountering planet this full rmvs time step
integer(I4B), dimension(:), allocatable :: tpenc1P !! index of first test particle encountering planet
integer(I4B), dimension(:), allocatable :: plind ! Connects the planetocentric indices back to the heliocentric planet list
type(rmvs_interp), dimension(:), allocatable :: outer !! interpolated heliocentric central body position for outer encounters
type(rmvs_interp), dimension(:), allocatable :: inner !! interpolated heliocentric central body position for inner encounters
contains
procedure, public :: fill => rmvs_fill_pl !! "Fills" bodies from one object into another depending on the results of a mask (uses the MERGE intrinsic)
procedure, public :: setup => rmvs_setup_pl !! Constructor method - Allocates space for number of particles
procedure, public :: spill => rmvs_spill_pl !! "Spills" bodies from one object to another depending on the results of a mask (uses the PACK intrinsic)
end type rmvs_base_pl

type, public :: rmvs_encounter_system
Expand All @@ -94,17 +101,9 @@ module rmvs_classes
end type rmvs_encounter_system

type, public, extends(rmvs_base_pl) :: rmvs_pl
integer(I4B), dimension(:), allocatable :: nenc !! number of test particles encountering planet this full rmvs time step
integer(I4B), dimension(:), allocatable :: tpenc1P !! index of first test particle encountering planet
type(rmvs_encounter_system), dimension(:), allocatable :: planetocentric
type(rmvs_interp), dimension(:), allocatable :: outer !! interpolated heliocentric central body position for outer encounters
type(rmvs_interp), dimension(:), allocatable :: inner !! interpolated heliocentric central body position for inner encounters
!! 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
procedure, public :: fill => rmvs_fill_pl !! "Fills" bodies from one object into another depending on the results of a mask (uses the MERGE intrinsic)
procedure, public :: setup => rmvs_setup_pl !! Constructor method - Allocates space for number of particles
procedure, public :: spill => rmvs_spill_pl !! "Spills" bodies from one object to another depending on the results of a mask (uses the PACK intrinsic)
end type rmvs_pl

interface
Expand Down Expand Up @@ -138,7 +137,7 @@ end subroutine rmvs_setup_system

module subroutine rmvs_setup_pl(self,n)
implicit none
class(rmvs_pl), intent(inout) :: self !! RMVS test particle object
class(rmvs_base_pl), intent(inout) :: self !! RMVS test particle object
integer, intent(in) :: n !! Number of test particles to allocate
end subroutine rmvs_setup_pl

Expand Down Expand Up @@ -174,15 +173,15 @@ end function rmvs_encounter_check_tp
module subroutine rmvs_spill_pl(self, discards, lspill_list)
use swiftest_classes, only : swiftest_body
implicit none
class(rmvs_pl), intent(inout) :: self !! RMVS massive body object
class(rmvs_base_pl), intent(inout) :: self !! RMVS massive body object
class(swiftest_body), intent(inout) :: discards !! Discarded object
logical, dimension(:), intent(in) :: lspill_list !! Logical array of bodies to spill into the discards
end subroutine rmvs_spill_pl

module subroutine rmvs_fill_pl(self, inserts, lfill_list)
use swiftest_classes, only : swiftest_body
implicit none
class(rmvs_pl), intent(inout) :: self !! RMVS massive body object
class(rmvs_base_pl), intent(inout) :: self !! RMVS massive body object
class(swiftest_body), intent(inout) :: inserts !! Inserted object
logical, dimension(:), intent(in) :: lfill_list !! Logical array of bodies to merge into the keeps
end subroutine rmvs_fill_pl
Expand Down
8 changes: 0 additions & 8 deletions src/rmvs/rmvs_getacch.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -30,14 +30,6 @@ module subroutine rmvs_getacch_tp(self, cb, pl, config, t, xh)
select type(pl)
class is (rmvs_pl)
allocate(xh_original, source=tp%xh)
write(*,*) 'ipleP = ',ipleP
write(*,*) 'enc_index = ',enc_index
write(*,*) 'xin-1: ',pl%inner(enc_index - 1)%x(:, ipleP)
write(*,*) 'xin : ',pl%inner(enc_index )%x(:, ipleP)
write(*,*) 'xh :'
do i = 1, pl%nbody
write(*,*) i,pl%xh(:,i)
end do
config_planetocen = config
! Temporarily turn off the heliocentric-dependent acceleration terms during an inner encounter
config_planetocen%loblatecb = .false.
Expand Down
6 changes: 4 additions & 2 deletions src/rmvs/rmvs_setup.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -9,7 +9,7 @@ module subroutine rmvs_setup_pl(self,n)
!! Equivalent in functionality to David E. Kaufmann's Swifter routine rmvs_setup.f90
implicit none
! Arguments
class(rmvs_pl), intent(inout) :: self !! RMVS test particle object
class(rmvs_base_pl), intent(inout) :: self !! RMVS test particle object
integer(I4B), intent(in) :: n !! Number of test particles to allocate
! Internals
integer(I4B) :: i,j
Expand Down Expand Up @@ -96,14 +96,16 @@ module subroutine rmvs_setup_system(self, config)
select type (tp => self%tp)
class is (rmvs_tp)
tp%cb_heliocentric = cb
pl%lplanetocentric = .false.
tp%lplanetocentric = .false.
cb%lplanetocentric = .false.
associate(npl => pl%nbody)
allocate(pl%planetocentric(npl))
do i = 1, npl
allocate(pl%planetocentric(i)%cb, source=cb)
allocate(rmvs_pl :: pl%planetocentric(i)%pl)
associate(plenci => pl%planetocentric(i)%pl, cbenci => pl%planetocentric(i)%cb)
cbenci%lplanetocentric = .true.

plenci%lplanetocentric = .true.
call plenci%setup(npl)
plenci%status(:) = ACTIVE
Expand Down
4 changes: 2 additions & 2 deletions src/rmvs/rmvs_spill_and_fill.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -9,7 +9,7 @@ module subroutine rmvs_spill_pl(self, discards, lspill_list)
!! Adapted from David E. Kaufmann's Swifter routine discard_discard_spill.f90
implicit none
! Arguments
class(rmvs_pl), intent(inout) :: self !! Swiftest massive body body object
class(rmvs_base_pl), intent(inout) :: self !! Swiftest massive body body object
class(swiftest_body), intent(inout) :: discards !! Discarded object
logical, dimension(:), intent(in) :: lspill_list !! Logical array of bodies to spill into the discards
! Internals
Expand Down Expand Up @@ -41,7 +41,7 @@ module subroutine rmvs_fill_pl(self, inserts, lfill_list)
!!
implicit none
! Arguments
class(rmvs_pl), intent(inout) :: self !! RMVS massive body object
class(rmvs_base_pl), intent(inout) :: self !! RMVS massive body object
class(swiftest_body), intent(inout) :: inserts !! Inserted object
logical, dimension(:), intent(in) :: lfill_list !! Logical array of bodies to merge into the keeps
! Internals
Expand Down
90 changes: 37 additions & 53 deletions src/rmvs/rmvs_step.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -411,51 +411,37 @@ subroutine rmvs_step_in(pl, cb, tp, config, outer_time, dto)
real(DP), dimension(:, :), allocatable :: xbeg, xend, vbeg


associate(npl => pl%nbody, nenc => pl%nenc, &
cbenci => pl%planetocentric(i)%cb, &
plenci => pl%planetocentric(i)%pl, &
tpenci => pl%planetocentric(i)%tp)
associate(enc_index => tpenci%index, plind => plenci%plind)

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(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(enc_index => tpenci%index, plind => plenci%plind)
! 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 enc_index = 1, NTPHENC ! Integrate over the encounter region, using the "substitute" planetocentric systems at each level
write(*,*) 'enc_index: ',enc_index
write(*,*) 'xin-1 ',pl%inner(enc_index - 1)%x(:, i)
write(*,*) 'xin ',pl%inner(enc_index )%x(:, i)
xbeg(:,1) = cb%xin(:) - pl%inner(enc_index - 1)%x(:, i)
xend(:,1) = cb%xin(:) - pl%inner(enc_index )%x(:, i)
vbeg(:,1) = cb%vin(:) - pl%inner(enc_index - 1)%v(:, i)
do j = 2, npl
ipleP = plind(j)
write(*,*) 'sub planet ', j
write(*,*) 'plind(i,j) ',ipleP
write(*,*) 'xin-1 ',pl%inner(enc_index - 1)%x(:,ipleP)
write(*,*) 'xin ',pl%inner(enc_index )%x(:,ipleP)
xbeg(:,j) = pl%inner(enc_index - 1)%x(:,ipleP) - pl%inner(enc_index - 1)%x(:,i)
xend(:,j) = pl%inner(enc_index )%x(:,ipleP) - pl%inner(enc_index )%x(:,i)
vbeg(:,j) = pl%inner(enc_index - 1)%v(:,ipleP) - pl%inner(enc_index - 1)%v(:,i)
end do
plenci%xh(:,:) = xbeg(:,:)
plenci%vh(:,:) = vbeg(:,:)
call tpenci%set_beg_end(xbeg = xbeg, xend = xend)
write(*,*) 'xend check'
do j = 1, npl
write(*,*) 'xend: ',j,xend(:,j)
write(*,*) 'xenc: ',j,tpenci%xend(:,j)
end do
call tpenci%step(cbenci, plenci, config, inner_time, dti)
do j = 1, nenc(i)
tpenci%xheliocentric(:, j) = tpenci%xh(:, j) + pl%inner(enc_index )%x(:,i)
Expand All @@ -464,10 +450,10 @@ subroutine rmvs_step_in(pl, cb, tp, config, outer_time, dto)
call rmvs_peri_tp(tpenci, pl, inner_time, dti, .false., enc_index, i, config)
end do
where(tpenci%status(:) == ACTIVE) tpenci%status(:) = INACTIVE
end do
call rmvs_end_planetocentric(pl, cb, tp)
end associate

end associate
end associate
end do
call rmvs_end_planetocentric(pl, cb, tp)
end associate
return
end subroutine rmvs_step_in
Expand All @@ -489,8 +475,7 @@ subroutine rmvs_make_planetocentric(pl, cb, tp, config)
integer(I4B) :: i, j
logical, dimension(:), allocatable :: encmask

associate(npl => pl%nbody, ntp => tp%nbody, GMpl => pl%Gmass, nenc => pl%nenc, &
cbenci => pl%planetocentric(i)%cb)
associate(npl => pl%nbody, ntp => tp%nbody, GMpl => pl%Gmass, nenc => pl%nenc)
do i = 1, npl
if (nenc(i) == 0) cycle
! There are inner encounters with this planet
Expand All @@ -500,7 +485,7 @@ subroutine rmvs_make_planetocentric(pl, cb, tp, config)

! Create encountering test particle structure
allocate(rmvs_tp :: pl%planetocentric(i)%tp)
associate(tpenci => pl%planetocentric(i)%tp)
associate(tpenci => pl%planetocentric(i)%tp, cbenci => pl%planetocentric(i)%cb, plenci => pl%planetocentric(i)%pl)
tpenci%lplanetocentric = .true.
call tpenci%setup(nenc(i))
tpenci%ipleP = i
Expand All @@ -512,10 +497,8 @@ subroutine rmvs_make_planetocentric(pl, cb, tp, config)
tpenci%xh(j, :) = tpenci%xheliocentric(j, :) - pl%inner(0)%x(j, i)
tpenci%vh(j, :) = pack(tp%vh(j,:), encmask(:)) - pl%inner(0)%v(j, i)
end do

allocate(plenci%inner, source = pl%inner)
! Make sure that the test particles get the planetocentric value of mu
write(*,*) i,'nenc: ',nenc(i)
write(*,*) cbenci%Gmass
call tpenci%set_mu(cbenci)
end associate
end do
Expand All @@ -538,27 +521,28 @@ subroutine rmvs_end_planetocentric(pl, cb, tp)
integer(I4B), dimension(:), allocatable :: tpind
logical, dimension(:), allocatable :: encmask

associate(nenc => pl%nenc, npl => pl%nbody, ntp => tp%nbody, tpenci => pl%planetocentric(i)%tp)

associate(nenc => pl%nenc, npl => pl%nbody, ntp => tp%nbody)
do i = 1, npl
if (nenc(i) == 0) cycle
allocate(tpind(nenc(i)))
! Index array of encountering test particles
if (allocated(encmask)) deallocate(encmask)
allocate(encmask(ntp))
encmask(:) = tp%plencP(:) == i
tpind(:) = pack([(j,j=1,ntp)], encmask(:))

! Copy the results of the integration back over and shift back to heliocentric reference
tp%status(tpind(1:nenc(i))) = tpenci%status(1:nenc(i))
do j = 1, NDIM
tp%xh(j, tpind(1:nenc(i))) = tpenci%xh(j,1:nenc(i)) + pl%inner(NTPHENC)%x(j, i)
tp%vh(j, tpind(1:nenc(i))) = tpenci%vh(j,1:nenc(i)) + pl%inner(NTPHENC)%v(j, i)
end do
deallocate(pl%planetocentric(i)%tp)
associate(tpenci => pl%planetocentric(i)%tp)
allocate(tpind(nenc(i)))
! Index array of encountering test particles
if (allocated(encmask)) deallocate(encmask)
allocate(encmask(ntp))
encmask(:) = tp%plencP(:) == i
tpind(:) = pack([(j,j=1,ntp)], encmask(:))

! Copy the results of the integration back over and shift back to heliocentric reference
tp%status(tpind(1:nenc(i))) = tpenci%status(1:nenc(i))
do j = 1, NDIM
tp%xh(j, tpind(1:nenc(i))) = tpenci%xh(j,1:nenc(i)) + pl%inner(NTPHENC)%x(j, i)
tp%vh(j, tpind(1:nenc(i))) = tpenci%vh(j,1:nenc(i)) + pl%inner(NTPHENC)%v(j, i)
end do
deallocate(pl%planetocentric(i)%tp)
deallocate(pl%planetocentric(i)%pl%inner)
end associate
end do
end associate

return
end subroutine rmvs_end_planetocentric

Expand Down

0 comments on commit b9838f0

Please sign in to comment.