From b9838f0f443cc40da6bdc85e50b1bc3a60887589 Mon Sep 17 00:00:00 2001 From: David Minton Date: Fri, 16 Apr 2021 17:00:02 -0400 Subject: [PATCH] Fixed bugs in data structures for planetocentric encounters --- src/modules/rmvs_classes.f90 | 39 +++++++------- src/rmvs/rmvs_getacch.f90 | 8 --- src/rmvs/rmvs_setup.f90 | 6 ++- src/rmvs/rmvs_spill_and_fill.f90 | 4 +- src/rmvs/rmvs_step.f90 | 90 +++++++++++++------------------- 5 files changed, 62 insertions(+), 85 deletions(-) diff --git a/src/modules/rmvs_classes.f90 b/src/modules/rmvs_classes.f90 index 9a3f32bd0..d5fe4861d 100644 --- a/src/modules/rmvs_classes.f90 +++ b/src/modules/rmvs_classes.f90 @@ -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 !******************************************************************************************************************************** @@ -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 @@ -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 @@ -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 @@ -174,7 +173,7 @@ 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 @@ -182,7 +181,7 @@ 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 diff --git a/src/rmvs/rmvs_getacch.f90 b/src/rmvs/rmvs_getacch.f90 index 9366404be..50dd38f1e 100644 --- a/src/rmvs/rmvs_getacch.f90 +++ b/src/rmvs/rmvs_getacch.f90 @@ -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. diff --git a/src/rmvs/rmvs_setup.f90 b/src/rmvs/rmvs_setup.f90 index 2f41870a8..3dd22ba43 100644 --- a/src/rmvs/rmvs_setup.f90 +++ b/src/rmvs/rmvs_setup.f90 @@ -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 @@ -96,6 +96,9 @@ 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 @@ -103,7 +106,6 @@ module subroutine rmvs_setup_system(self, config) 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 diff --git a/src/rmvs/rmvs_spill_and_fill.f90 b/src/rmvs/rmvs_spill_and_fill.f90 index d072e0274..4eb1e81dc 100644 --- a/src/rmvs/rmvs_spill_and_fill.f90 +++ b/src/rmvs/rmvs_spill_and_fill.f90 @@ -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 @@ -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 diff --git a/src/rmvs/rmvs_step.f90 b/src/rmvs/rmvs_step.f90 index 0f0b605ed..13fdf2230 100644 --- a/src/rmvs/rmvs_step.f90 +++ b/src/rmvs/rmvs_step.f90 @@ -411,20 +411,18 @@ 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 @@ -432,18 +430,11 @@ subroutine rmvs_step_in(pl, cb, tp, config, outer_time, dto) ! 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) @@ -451,11 +442,6 @@ subroutine rmvs_step_in(pl, cb, tp, config, outer_time, dto) 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) @@ -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 @@ -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 @@ -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 @@ -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 @@ -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