From 5b5d8e8891c2a663343a44e6bbce0b7d1e921f4d Mon Sep 17 00:00:00 2001 From: David A Minton Date: Wed, 28 Jul 2021 10:27:51 -0400 Subject: [PATCH 1/7] Cleaned up a few comments --- src/modules/symba_classes.f90 | 8 ++++---- src/symba/symba_util.f90 | 10 +++++----- 2 files changed, 9 insertions(+), 9 deletions(-) diff --git a/src/modules/symba_classes.f90 b/src/modules/symba_classes.f90 index fc3520079..7eab9e5c7 100644 --- a/src/modules/symba_classes.f90 +++ b/src/modules/symba_classes.f90 @@ -390,27 +390,27 @@ end subroutine symba_util_resize_pltpenc module subroutine symba_util_sort_pl(self, sortby, ascending) implicit none - class(symba_pl), intent(inout) :: self !! Symba massive body object + class(symba_pl), intent(inout) :: self !! SyMBA massive body object character(*), intent(in) :: sortby !! Sorting attribute logical, intent(in) :: ascending !! Logical flag indicating whether or not the sorting should be in ascending or descending order end subroutine symba_util_sort_pl module subroutine symba_util_sort_tp(self, sortby, ascending) implicit none - class(symba_tp), intent(inout) :: self !! Swiftest test particle object + class(symba_tp), intent(inout) :: self !! SyMBA test particle object character(*), intent(in) :: sortby !! Sorting attribute logical, intent(in) :: ascending !! Logical flag indicating whether or not the sorting should be in ascending or descending order end subroutine symba_util_sort_tp module subroutine symba_util_sort_rearrange_pl(self, ind) implicit none - class(symba_pl), intent(inout) :: self !! Symba massive body object + class(symba_pl), intent(inout) :: self !! SyMBA massive body object integer(I4B), dimension(:), intent(in) :: ind !! Index array used to restructure the body (should contain all 1:n index values in the desired order) end subroutine symba_util_sort_rearrange_pl module subroutine symba_util_sort_rearrange_tp(self, ind) implicit none - class(symba_tp), intent(inout) :: self !! Symba massive body object + class(symba_tp), intent(inout) :: self !! SyMBA massive body object integer(I4B), dimension(:), intent(in) :: ind !! Index array used to restructure the body (should contain all 1:n index values in the desired order) end subroutine symba_util_sort_rearrange_tp diff --git a/src/symba/symba_util.f90 b/src/symba/symba_util.f90 index a89da7b6a..be47c9b96 100644 --- a/src/symba/symba_util.f90 +++ b/src/symba/symba_util.f90 @@ -68,11 +68,11 @@ end subroutine symba_util_resize_pltpenc module subroutine symba_util_sort_pl(self, sortby, ascending) !! author: David A. Minton !! - !! Sort a Swiftest test particle object in-place. + !! Sort a SyMBA massive body object in-place. !! sortby is a string indicating which array component to sort. implicit none ! Arguments - class(symba_pl), intent(inout) :: self !! Symba massive body object + class(symba_pl), intent(inout) :: self !! SyMBA massive body object character(*), intent(in) :: sortby !! Sorting attribute logical, intent(in) :: ascending !! Logical flag indicating whether or not the sorting should be in ascending or descending order ! Internals @@ -128,11 +128,11 @@ end subroutine symba_util_sort_pl module subroutine symba_util_sort_tp(self, sortby, ascending) !! author: David A. Minton !! - !! Sort a Swiftest test particle object in-place. + !! Sort a SyMBA test particle object in-place. !! sortby is a string indicating which array component to sort. implicit none ! Arguments - class(symba_tp), intent(inout) :: self !! Swiftest test particle object + class(symba_tp), intent(inout) :: self !! SyMBA test particle object character(*), intent(in) :: sortby !! Sorting attribute logical, intent(in) :: ascending !! Logical flag indicating whether or not the sorting should be in ascending or descending order ! Internals @@ -207,7 +207,7 @@ end subroutine symba_util_sort_rearrange_pl module subroutine symba_util_sort_rearrange_tp(self, ind) !! author: David A. Minton !! - !! Rearrange SyMBA test particle object in-place from an index list. + !! Rearrange SyMBA test particle object in-place from an index list. !! This is a helper utility used to make polymorphic sorting work on Swiftest structures. implicit none ! Arguments From 570717efc6a58047d98b25880e6ac4b597d2da54 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Wed, 28 Jul 2021 10:39:09 -0400 Subject: [PATCH 2/7] Added sort and rearrange methods to WHM --- src/modules/whm_classes.f90 | 43 ++++++++++++++++--------- src/whm/whm_util.f90 | 62 +++++++++++++++++++++++++++++++++++++ 2 files changed, 91 insertions(+), 14 deletions(-) diff --git a/src/modules/whm_classes.f90 b/src/modules/whm_classes.f90 index a6ab59958..c242d2521 100644 --- a/src/modules/whm_classes.f90 +++ b/src/modules/whm_classes.f90 @@ -30,20 +30,22 @@ module whm_classes !! Note to developers: If you add componenets to this class, be sure to update methods and subroutines that traverse the !! component list, such as whm_setup_pl and whm_util_spill_pl contains - procedure :: h2j => whm_coord_h2j_pl !! Convert position and velcoity vectors from heliocentric to Jacobi coordinates - procedure :: j2h => whm_coord_j2h_pl !! Convert position and velcoity vectors from Jacobi to helliocentric coordinates - procedure :: vh2vj => whm_coord_vh2vj_pl !! Convert velocity vectors from heliocentric to Jacobi coordinates - procedure :: drift => whm_drift_pl !! Loop through massive bodies and call Danby drift routine to jacobi coordinates - procedure :: fill => whm_util_fill_pl !! "Fills" bodies from one object into another depending on the results of a mask (uses the MERGE intrinsic) - procedure :: accel => whm_kick_getacch_pl !! Compute heliocentric accelerations of massive bodies - procedure :: kick => whm_kick_vh_pl !! Kick heliocentric velocities of massive bodies - procedure :: accel_gr => whm_gr_kick_getacch_pl !! Acceleration term arising from the post-Newtonian correction - procedure :: gr_pos_kick => whm_gr_p4_pl !! Position kick due to p**4 term in the post-Newtonian correction - procedure :: setup => whm_setup_pl !! Constructor method - Allocates space for number of particles - procedure :: set_mu => whm_util_set_mu_eta_pl !! Sets the Jacobi mass value for all massive bodies. - procedure :: set_ir3 => whm_util_set_ir3j !! Sets both the heliocentric and jacobi inverse radius terms (1/rj**3 and 1/rh**3) - procedure :: step => whm_step_pl !! Steps the body forward one stepsize - procedure :: spill => whm_util_spill_pl !!"Spills" bodies from one object to another depending on the results of a mask (uses the PACK intrinsic) + procedure :: h2j => whm_coord_h2j_pl !! Convert position and velcoity vectors from heliocentric to Jacobi coordinates + procedure :: j2h => whm_coord_j2h_pl !! Convert position and velcoity vectors from Jacobi to helliocentric coordinates + procedure :: vh2vj => whm_coord_vh2vj_pl !! Convert velocity vectors from heliocentric to Jacobi coordinates + procedure :: drift => whm_drift_pl !! Loop through massive bodies and call Danby drift routine to jacobi coordinates + procedure :: fill => whm_util_fill_pl !! "Fills" bodies from one object into another depending on the results of a mask (uses the MERGE intrinsic) + procedure :: accel => whm_kick_getacch_pl !! Compute heliocentric accelerations of massive bodies + procedure :: kick => whm_kick_vh_pl !! Kick heliocentric velocities of massive bodies + procedure :: accel_gr => whm_gr_kick_getacch_pl !! Acceleration term arising from the post-Newtonian correction + procedure :: gr_pos_kick => whm_gr_p4_pl !! Position kick due to p**4 term in the post-Newtonian correction + procedure :: setup => whm_setup_pl !! Constructor method - Allocates space for number of particles + procedure :: set_mu => whm_util_set_mu_eta_pl !! Sets the Jacobi mass value for all massive bodies. + procedure :: set_ir3 => whm_util_set_ir3j !! Sets both the heliocentric and jacobi inverse radius terms (1/rj**3 and 1/rh**3) + procedure :: sort => whm_util_sort_pl !! Sort a WHM massive body object in-place. + procedure :: rearrange => whm_util_sort_rearrange_pl !! Rearranges the order of array elements of body based on an input index array. Used in sorting methods + procedure :: step => whm_step_pl !! Steps the body forward one stepsize + procedure :: spill => whm_util_spill_pl !!"Spills" bodies from one object to another depending on the results of a mask (uses the PACK intrinsic) end type whm_pl !******************************************************************************************************************************** @@ -209,6 +211,19 @@ module subroutine whm_util_set_mu_eta_pl(self, cb) class(swiftest_cb), intent(inout) :: cb !! Swiftest central body object end subroutine whm_util_set_mu_eta_pl + module subroutine whm_util_sort_pl(self, sortby, ascending) + implicit none + class(whm_pl), intent(inout) :: self !! WHM massive body object + character(*), intent(in) :: sortby !! Sorting attribute + logical, intent(in) :: ascending !! Logical flag indicating whether or not the sorting should be in ascending or descending order + end subroutine whm_util_sort_pl + + module subroutine whm_util_sort_rearrange_pl(self, ind) + implicit none + class(whm_pl), intent(inout) :: self !! WHM massive body object + integer(I4B), dimension(:), intent(in) :: ind !! Index array used to restructure the body (should contain all 1:n index values in the desired order) + end subroutine whm_util_sort_rearrange_pl + module subroutine whm_setup_initialize_system(self, param) use swiftest_classes, only : swiftest_parameters implicit none diff --git a/src/whm/whm_util.f90 b/src/whm/whm_util.f90 index 67c7ef4a1..791b5d651 100644 --- a/src/whm/whm_util.f90 +++ b/src/whm/whm_util.f90 @@ -112,4 +112,66 @@ module subroutine whm_util_set_ir3j(self) end if end subroutine whm_util_set_ir3j + module subroutine whm_util_sort_pl(self, sortby, ascending) + !! author: David A. Minton + !! + !! Sort a WHM massive body object in-place. + !! sortby is a string indicating which array component to sort. + implicit none + ! Arguments + class(whm_pl), intent(inout) :: self !! WHM massive body object + character(*), intent(in) :: sortby !! Sorting attribute + logical, intent(in) :: ascending !! Logical flag indicating whether or not the sorting should be in ascending or descending order + ! Internals + integer(I4B), dimension(self%nbody) :: ind + integer(I4B) :: direction + + if (ascending) then + direction = 1 + else + direction = -1 + end if + associate(pl => self, npl => self%nbody) + select case(sortby) + case("eta") + call util_sort(direction * pl%eta(1:npl), ind(1:npl)) + case("muj") + call util_sort(direction * pl%muj(1:npl), ind(1:npl)) + case("ir3j") + call util_sort(direction * pl%ir3j(1:npl), ind(1:npl)) + case default + call util_sort_pl(pl, sortby, ascending) + return + end select + call pl%rearrange(ind) + end associate + return + end subroutine whm_util_sort_pl + + module subroutine whm_util_sort_rearrange_pl(self, ind) + !! author: David A. Minton + !! + !! Rearrange WHM massive body structure in-place from an index list. + !! This is a helper utility used to make polymorphic sorting work on Swiftest structures. + implicit none + ! Arguments + class(whm_pl), intent(inout) :: self !! WHM massive body object + integer(I4B), dimension(:), intent(in) :: ind !! Index array used to restructure the body (should contain all 1:n index values in the desired order) + ! Internals + class(whm_pl), allocatable :: pl_sorted !! Temporary holder for sorted body + integer(I4B) :: i + + associate(pl => self, npl => self%nbody) + call util_sort_rearrange_pl(pl,ind) + allocate(pl_sorted, source=self) + pl%eta(1:npl) = pl_sorted%eta(ind(1:npl)) + pl%xj(:,1:npl) = pl_sorted%xj(:,ind(1:npl)) + pl%vj(:,1:npl) = pl_sorted%vj(:,ind(1:npl)) + pl%muj(1:npl) = pl_sorted%muj(ind(1:npl)) + pl%ir3j(1:npl) = pl_sorted%ir3j(ind(1:npl)) + deallocate(pl_sorted) + end associate + return + end subroutine whm_util_sort_rearrange_pl + end submodule s_whm_util From f3ef56c63ae28a76c4f262a9befea5fe92e4ab09 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Wed, 28 Jul 2021 10:45:06 -0400 Subject: [PATCH 3/7] Simplified sorting methods with direction variable --- src/symba/symba_util.f90 | 72 +++++++++--------------- src/util/util_sort.f90 | 118 ++++++++++++++------------------------- src/whm/whm_util.f90 | 3 + 3 files changed, 72 insertions(+), 121 deletions(-) diff --git a/src/symba/symba_util.f90 b/src/symba/symba_util.f90 index be47c9b96..397cd789e 100644 --- a/src/symba/symba_util.f90 +++ b/src/symba/symba_util.f90 @@ -77,50 +77,35 @@ module subroutine symba_util_sort_pl(self, sortby, ascending) logical, intent(in) :: ascending !! Logical flag indicating whether or not the sorting should be in ascending or descending order ! Internals integer(I4B), dimension(self%nbody) :: ind + integer(I4B) :: direction + + if (ascending) then + direction = 1 + else + direction = -1 + end if associate(pl => self, npl => self%nbody) select case(sortby) case("nplenc") - if (ascending) then - call util_sort(pl%nplenc(1:npl), ind(1:npl)) - else - call util_sort(-pl%nplenc(1:npl), ind(1:npl)) - end if + call util_sort(direction * pl%nplenc(1:npl), ind(1:npl)) case("ntpenc") - if (ascending) then - call util_sort(pl%ntpenc(1:npl), ind(1:npl)) - else - call util_sort(-pl%ntpenc(1:npl), ind(1:npl)) - end if + call util_sort(direction * pl%ntpenc(1:npl), ind(1:npl)) case("levelg") - if (ascending) then - call util_sort(pl%levelg(1:npl), ind(1:npl)) - else - call util_sort(-pl%levelg(1:npl), ind(1:npl)) - end if + call util_sort(direction * pl%levelg(1:npl), ind(1:npl)) case("levelm") - if (ascending) then - call util_sort(pl%levelm(1:npl), ind(1:npl)) - else - call util_sort(-pl%levelm(1:npl), ind(1:npl)) - end if + call util_sort(direction * pl%levelm(1:npl), ind(1:npl)) case("peri") - if (ascending) then - call util_sort(pl%peri(1:npl), ind(1:npl)) - else - call util_sort(-pl%peri(1:npl), ind(1:npl)) - end if + call util_sort(direction * pl%peri(1:npl), ind(1:npl)) case("atp") - if (ascending) then - call util_sort(pl%atp(1:npl), ind(1:npl)) - else - call util_sort(-pl%atp(1:npl), ind(1:npl)) - end if + call util_sort(direction * pl%atp(1:npl), ind(1:npl)) case default call util_sort_pl(pl, sortby, ascending) return end select + call pl%rearrange(ind) + end associate return end subroutine symba_util_sort_pl @@ -137,32 +122,29 @@ module subroutine symba_util_sort_tp(self, sortby, ascending) logical, intent(in) :: ascending !! Logical flag indicating whether or not the sorting should be in ascending or descending order ! Internals integer(I4B), dimension(self%nbody) :: ind + integer(I4B) :: direction + + if (ascending) then + direction = 1 + else + direction = -1 + end if associate(tp => self, ntp => self%nbody) select case(sortby) case("nplenc") - if (ascending) then - call util_sort(tp%nplenc(1:ntp), ind(1:ntp)) - else - call util_sort(-tp%nplenc(1:ntp), ind(1:ntp)) - end if + call util_sort(direction * tp%nplenc(1:ntp), ind(1:ntp)) case("levelg") - if (ascending) then - call util_sort(tp%levelg(1:ntp), ind(1:ntp)) - else - call util_sort(-tp%levelg(1:ntp), ind(1:ntp)) - end if + call util_sort(direction * tp%levelg(1:ntp), ind(1:ntp)) case("levelm") - if (ascending) then - call util_sort(tp%levelm(1:ntp), ind(1:ntp)) - else - call util_sort(-tp%levelm(1:ntp), ind(1:ntp)) - end if + call util_sort(direction * tp%levelm(1:ntp), ind(1:ntp)) case default call util_sort_tp(tp, sortby, ascending) return end select + call tp%rearrange(ind) + end associate return end subroutine symba_util_sort_tp diff --git a/src/util/util_sort.f90 b/src/util/util_sort.f90 index 34fe600ed..a3e2d37a4 100644 --- a/src/util/util_sort.f90 +++ b/src/util/util_sort.f90 @@ -13,50 +13,35 @@ module subroutine util_sort_body(self, sortby, ascending) logical, intent(in) :: ascending !! Logical flag indicating whether or not the sorting should be in ascending or descending order ! Internals integer(I4B), dimension(self%nbody) :: ind + integer(I4B) :: direction + + if (ascending) then + direction = 1 + else + direction = -1 + end if associate(body => self, n => self%nbody) select case(sortby) case("id") - if (ascending) then - call util_sort(body%id(1:n), ind(1:n)) - else - call util_sort(-body%id(1:n), ind(1:n)) - end if + call util_sort(direction * body%id(1:n), ind(1:n)) case("a") - if (ascending) then - call util_sort(body%a(1:n), ind(1:n)) - else - call util_sort(-body%a(1:n), ind(1:n)) - end if + call util_sort(direction * body%a(1:n), ind(1:n)) case("e") - if (ascending) then - call util_sort(body%e(1:n), ind(1:n)) - else - call util_sort(-body%e(1:n), ind(1:n)) - end if + call util_sort(direction * body%e(1:n), ind(1:n)) case("inc") - if (ascending) then - call util_sort(body%inc(1:n), ind(1:n)) - else - call util_sort(-body%inc(1:n), ind(1:n)) - end if + call util_sort(direction * body%inc(1:n), ind(1:n)) case("capom") - if (ascending) then - call util_sort(body%capom(1:n), ind(1:n)) - else - call util_sort(-body%capom(1:n), ind(1:n)) - end if + call util_sort(direction * body%capom(1:n), ind(1:n)) case("mu") - if (ascending) then - call util_sort(body%mu(1:n), ind(1:n)) - else - call util_sort(-body%mu(1:n), ind(1:n)) - end if + call util_sort(direction * body%mu(1:n), ind(1:n)) case default - write(*,*) 'Cannot sort structure by component '//trim(adjustl(sortby)) + write(*,*) 'Cannot sort structure by component ' // trim(adjustl(sortby)) return end select + call body%rearrange(ind) + end associate return @@ -74,56 +59,37 @@ module subroutine util_sort_pl(self, sortby, ascending) logical, intent(in) :: ascending !! Logical flag indicating whether or not the sorting should be in ascending or descending order ! Internals integer(I4B), dimension(self%nbody) :: ind + integer(I4B) :: direction + + if (ascending) then + direction = 1 + else + direction = -1 + end if associate(pl => self, npl => self%nbody) select case(sortby) case("Gmass","mass") - if (ascending) then - call util_sort(pl%Gmass(1:npl), ind(1:npl)) - else - call util_sort(-pl%Gmass(1:npl), ind(1:npl)) - end if + call util_sort(direction * pl%Gmass(1:npl), ind(1:npl)) case("rhill") - if (ascending) then - call util_sort(pl%rhill(1:npl), ind(1:npl)) - else - call util_sort(-pl%rhill(1:npl), ind(1:npl)) - end if + call util_sort(direction * pl%rhill(1:npl), ind(1:npl)) case("radius") - if (ascending) then - call util_sort(pl%radius(1:npl), ind(1:npl)) - else - call util_sort(-pl%radius(1:npl), ind(1:npl)) - end if + call util_sort(direction * pl%radius(1:npl), ind(1:npl)) case("density") - if (ascending) then - call util_sort(pl%density(1:npl), ind(1:npl)) - else - call util_sort(-pl%density(1:npl), ind(1:npl)) - end if + call util_sort(direction * pl%density(1:npl), ind(1:npl)) case("k2") - if (ascending) then - call util_sort(pl%k2(1:npl), ind(1:npl)) - else - call util_sort(-pl%k2(1:npl), ind(1:npl)) - end if + call util_sort(direction * pl%k2(1:npl), ind(1:npl)) case("Q") - if (ascending) then - call util_sort(pl%Q(1:npl), ind(1:npl)) - else - call util_sort(-pl%Q(1:npl), ind(1:npl)) - end if + call util_sort(direction * pl%Q(1:npl), ind(1:npl)) case("tlag") - if (ascending) then - call util_sort(pl%tlag(1:npl), ind(1:npl)) - else - call util_sort(-pl%tlag(1:npl), ind(1:npl)) - end if + call util_sort(direction * pl%tlag(1:npl), ind(1:npl)) case default call util_sort_body(pl, sortby, ascending) return end select + call pl%rearrange(ind) + end associate return @@ -141,27 +107,27 @@ module subroutine util_sort_tp(self, sortby, ascending) logical, intent(in) :: ascending !! Logical flag indicating whether or not the sorting should be in ascending or descending order ! Internals integer(I4B), dimension(self%nbody) :: ind + integer(I4B) :: direction + + if (ascending) then + direction = 1 + else + direction = -1 + end if associate(tp => self, ntp => self%nbody) select case(sortby) case("peri") - if (ascending) then - call util_sort(tp%peri(1:ntp), ind(1:ntp)) - else - call util_sort(-tp%peri(1:ntp), ind(1:ntp)) - end if + call util_sort(direction * tp%peri(1:ntp), ind(1:ntp)) case("atp") - if (ascending) then - call util_sort(tp%atp(1:ntp), ind(1:ntp)) - else - call util_sort(-tp%atp(1:ntp), ind(1:ntp)) - end if + call util_sort(direction * tp%atp(1:ntp), ind(1:ntp)) case default call util_sort_body(tp, sortby, ascending) return end select call tp%rearrange(ind) + end associate return diff --git a/src/whm/whm_util.f90 b/src/whm/whm_util.f90 index 791b5d651..9ebeb0b3f 100644 --- a/src/whm/whm_util.f90 +++ b/src/whm/whm_util.f90 @@ -131,6 +131,7 @@ module subroutine whm_util_sort_pl(self, sortby, ascending) else direction = -1 end if + associate(pl => self, npl => self%nbody) select case(sortby) case("eta") @@ -143,7 +144,9 @@ module subroutine whm_util_sort_pl(self, sortby, ascending) call util_sort_pl(pl, sortby, ascending) return end select + call pl%rearrange(ind) + end associate return end subroutine whm_util_sort_pl From 1bd12ac47a0d7b2185341586722d7b97155978de Mon Sep 17 00:00:00 2001 From: David A Minton Date: Wed, 28 Jul 2021 10:47:56 -0400 Subject: [PATCH 4/7] Rearranged rmvs_util implementations into alphabetical order --- src/rmvs/rmvs_util.f90 | 134 ++++++++++++++++++++--------------------- 1 file changed, 65 insertions(+), 69 deletions(-) diff --git a/src/rmvs/rmvs_util.f90 b/src/rmvs/rmvs_util.f90 index 0781c6429..379c1d4c0 100644 --- a/src/rmvs/rmvs_util.f90 +++ b/src/rmvs/rmvs_util.f90 @@ -1,6 +1,70 @@ submodule(rmvs_classes) s_rmvs_util use swiftest contains + module subroutine rmvs_util_fill_pl(self, inserts, lfill_list) + !! author: David A. Minton + !! + !! Insert new RMVS massive body structure into an old one. + !! This is the inverse of a fill operation. + !! + implicit none + ! Arguments + class(rmvs_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 + integer(I4B) :: i + + associate(keeps => self) + select type(inserts) + class is (rmvs_pl) + + keeps%nenc(:) = unpack(keeps%nenc(:), .not.lfill_list(:), keeps%nenc(:)) + keeps%nenc(:) = unpack(inserts%nenc(:), lfill_list(:), keeps%nenc(:)) + + call whm_util_fill_pl(keeps, inserts, lfill_list) + class default + write(*,*) 'Error! spill method called for incompatible return type on rmvs_pl' + end select + end associate + + return + end subroutine rmvs_util_fill_pl + + module subroutine rmvs_util_fill_tp(self, inserts, lfill_list) + !! author: David A. Minton + !! + !! Insert new RMVS test particle structure into an old one. + !! This is the inverse of a fill operation. + !! + implicit none + ! Arguments + class(rmvs_tp), 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 + + associate(keeps => self) + select type(inserts) + class is (rmvs_tp) + + keeps%lperi(:) = unpack(keeps%lperi(:), .not.lfill_list(:), keeps%lperi(:)) + keeps%lperi(:) = unpack(inserts%lperi(:), lfill_list(:), keeps%lperi(:)) + + keeps%plperP(:) = unpack(keeps%plperP(:), .not.lfill_list(:), keeps%plperP(:)) + keeps%plperP(:) = unpack(inserts%plperP(:), lfill_list(:), keeps%plperP(:)) + + keeps%plencP(:) = unpack(keeps%plencP(:), .not.lfill_list(:), keeps%plencP(:)) + keeps%plencP(:) = unpack(inserts%plencP(:), lfill_list(:), keeps%plencP(:)) + + call util_fill_tp(keeps, inserts, lfill_list) + class default + write(*,*) 'Error! fill method called for incompatible return type on rmvs_tp' + end select + end associate + + return + end subroutine rmvs_util_fill_tp + module subroutine rmvs_util_spill_pl(self, discards, lspill_list) !! author: David A. Minton !! @@ -30,39 +94,7 @@ module subroutine rmvs_util_spill_pl(self, discards, lspill_list) end associate return - - end subroutine rmvs_util_spill_pl - - module subroutine rmvs_util_fill_pl(self, inserts, lfill_list) - !! author: David A. Minton - !! - !! Insert new RMVS massive body structure into an old one. - !! This is the inverse of a fill operation. - !! - implicit none - ! Arguments - class(rmvs_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 - integer(I4B) :: i - - associate(keeps => self) - select type(inserts) - class is (rmvs_pl) - - keeps%nenc(:) = unpack(keeps%nenc(:), .not.lfill_list(:), keeps%nenc(:)) - keeps%nenc(:) = unpack(inserts%nenc(:), lfill_list(:), keeps%nenc(:)) - - call whm_util_fill_pl(keeps, inserts, lfill_list) - class default - write(*,*) 'Error! spill method called for incompatible return type on rmvs_pl' - end select - end associate - - return - - end subroutine rmvs_util_fill_pl + end subroutine rmvs_util_spill_pl module subroutine rmvs_util_spill_tp(self, discards, lspill_list) !! author: David A. Minton @@ -97,42 +129,6 @@ module subroutine rmvs_util_spill_tp(self, discards, lspill_list) end associate return - end subroutine rmvs_util_spill_tp - module subroutine rmvs_util_fill_tp(self, inserts, lfill_list) - !! author: David A. Minton - !! - !! Insert new RMVS test particle structure into an old one. - !! This is the inverse of a fill operation. - !! - implicit none - ! Arguments - class(rmvs_tp), 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 - - associate(keeps => self) - select type(inserts) - class is (rmvs_tp) - - keeps%lperi(:) = unpack(keeps%lperi(:), .not.lfill_list(:), keeps%lperi(:)) - keeps%lperi(:) = unpack(inserts%lperi(:), lfill_list(:), keeps%lperi(:)) - - keeps%plperP(:) = unpack(keeps%plperP(:), .not.lfill_list(:), keeps%plperP(:)) - keeps%plperP(:) = unpack(inserts%plperP(:), lfill_list(:), keeps%plperP(:)) - - keeps%plencP(:) = unpack(keeps%plencP(:), .not.lfill_list(:), keeps%plencP(:)) - keeps%plencP(:) = unpack(inserts%plencP(:), lfill_list(:), keeps%plencP(:)) - - call util_fill_tp(keeps, inserts, lfill_list) - class default - write(*,*) 'Error! fill method called for incompatible return type on rmvs_tp' - end select - end associate - - return - - end subroutine rmvs_util_fill_tp - end submodule s_rmvs_util From 2536428866532ec74cf5abee0005d4754e3a3bde Mon Sep 17 00:00:00 2001 From: David A Minton Date: Wed, 28 Jul 2021 10:48:57 -0400 Subject: [PATCH 5/7] Cleaned up formatting and corrected comment typos --- src/rmvs/rmvs_util.f90 | 20 ++++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/src/rmvs/rmvs_util.f90 b/src/rmvs/rmvs_util.f90 index 379c1d4c0..29d0b7b2a 100644 --- a/src/rmvs/rmvs_util.f90 +++ b/src/rmvs/rmvs_util.f90 @@ -13,7 +13,7 @@ module subroutine rmvs_util_fill_pl(self, inserts, lfill_list) class(swiftest_body), intent(inout) :: inserts !! Inserted object logical, dimension(:), intent(in) :: lfill_list !! Logical array of bodies to merge into the keeps ! Internals - integer(I4B) :: i + integer(I4B) :: i associate(keeps => self) select type(inserts) @@ -39,9 +39,9 @@ module subroutine rmvs_util_fill_tp(self, inserts, lfill_list) !! implicit none ! Arguments - class(rmvs_tp), 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 + class(rmvs_tp), intent(inout) :: self !! RMVS test particle object + class(swiftest_body), intent(inout) :: inserts !! Inserted object + logical, dimension(:), intent(in) :: lfill_list !! Logical array of bodies to merge into the keeps associate(keeps => self) select type(inserts) @@ -73,11 +73,11 @@ module subroutine rmvs_util_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_pl), intent(inout) :: self !! RMVS 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 - integer(I4B) :: i + integer(I4B) :: i associate(keeps => self) select type(discards) @@ -104,11 +104,11 @@ module subroutine rmvs_util_spill_tp(self, discards, lspill_list) !! Adapted from David E. Kaufmann's Swifter routine whm_discard_spill.f90 implicit none ! Arguments - class(rmvs_tp), intent(inout) :: self !! RMVS test particle object - class(swiftest_body), intent(inout) :: discards !! Discarded object - logical, dimension(:), intent(in) :: lspill_list !! Logical array of bodies to spill into the discards + class(rmvs_tp), intent(inout) :: self !! RMVS test particle 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 - integer(I4B) :: i + integer(I4B) :: i associate(keeps => self) select type(discards) From 952ffc8966e926a60fa2492d6f9c7978d0c53365 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Wed, 28 Jul 2021 11:13:12 -0400 Subject: [PATCH 6/7] Added RMVS sorting methods and cleaned up some comment typos and formatting --- src/modules/rmvs_classes.f90 | 84 +++++++++++++++------- src/modules/symba_classes.f90 | 4 +- src/rmvs/rmvs_util.f90 | 130 ++++++++++++++++++++++++++++++++++ src/symba/symba_setup.f90 | 4 +- src/symba/symba_util.f90 | 10 +-- src/util/util_sort.f90 | 12 +++- src/whm/whm_util.f90 | 2 + 7 files changed, 208 insertions(+), 38 deletions(-) diff --git a/src/modules/rmvs_classes.f90 b/src/modules/rmvs_classes.f90 index 37a88993c..de4cdec4c 100644 --- a/src/modules/rmvs_classes.f90 +++ b/src/modules/rmvs_classes.f90 @@ -66,13 +66,15 @@ module rmvs_classes integer(I4B) :: ipleP !! index value of encountering planet logical :: lplanetocentric = .false. !! Flag that indicates that the object is a planetocentric set of masive bodies used for close encounter calculations contains - procedure :: discard => rmvs_discard_tp !! Check to see if test particles should be discarded based on pericenter passage distances with respect to planets encountered - procedure :: encounter_check => rmvs_encounter_check_tp !! Checks if any test particles are undergoing a close encounter with a massive body - procedure :: fill => rmvs_util_fill_tp !! "Fills" bodies from one object into another depending on the results of a mask (uses the MERGE intrinsic) - procedure :: accel => rmvs_kick_getacch_tp !! Calculates either the standard or modified version of the acceleration depending if the - !! if the test particle is undergoing a close encounter or not - procedure :: setup => rmvs_setup_tp !! Constructor method - Allocates space for number of particles - procedure :: spill => rmvs_util_spill_tp !! "Spills" bodies from one object to another depending on the results of a mask (uses the PACK intrinsic) + procedure :: discard => rmvs_discard_tp !! Check to see if test particles should be discarded based on pericenter passage distances with respect to planets encountered + procedure :: encounter_check => rmvs_encounter_check_tp !! Checks if any test particles are undergoing a close encounter with a massive body + procedure :: accel => rmvs_kick_getacch_tp !! Calculates either the standard or modified version of the acceleration depending if the + !! if the test particle is undergoing a close encounter or not + procedure :: setup => rmvs_setup_tp !! Constructor method - Allocates space for number of particles + procedure :: fill => rmvs_util_fill_tp !! "Fills" bodies from one object into another depending on the results of a mask (uses the MERGE intrinsic) + procedure :: sort => rmvs_util_sort_tp !! Sorts body arrays by a sortable componen + procedure :: rearrange => rmvs_util_sort_rearrange_tp !! Rearranges the order of array elements of body based on an input index array. Used in sorting methods + procedure :: spill => rmvs_util_spill_tp !! "Spills" bodies from one object to another depending on the results of a mask (uses the PACK intrinsic) end type rmvs_tp !******************************************************************************************************************************** @@ -89,9 +91,11 @@ module rmvs_classes class(rmvs_nbody_system), dimension(:), allocatable :: planetocentric !! Planetocentric version of the massive body objects (one for each massive body) logical :: lplanetocentric = .false. !! Flag that indicates that the object is a planetocentric set of masive bodies used for close encounter calculations contains - procedure :: fill => rmvs_util_fill_pl !! "Fills" bodies from one object into another depending on the results of a mask (uses the MERGE intrinsic) - procedure :: setup => rmvs_setup_pl !! Constructor method - Allocates space for number of particles - procedure :: spill => rmvs_util_spill_pl !! "Spills" bodies from one object to another depending on the results of a mask (uses the PACK intrinsic) + procedure :: setup => rmvs_setup_pl !! Constructor method - Allocates space for number of particles + procedure :: sort => rmvs_util_sort_pl !! Sorts body arrays by a sortable componen + procedure :: rearrange => rmvs_util_sort_rearrange_pl !! Rearranges the order of array elements of body based on an input index array. Used in sorting methods + procedure :: fill => rmvs_util_fill_pl !! "Fills" bodies from one object into another depending on the results of a mask (uses the MERGE intrinsic) + procedure :: spill => rmvs_util_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 @@ -117,22 +121,6 @@ module function rmvs_encounter_check_tp(self, system, dt) result(lencounter) logical :: lencounter !! Returns true if there is at least one close encounter end function rmvs_encounter_check_tp - module subroutine rmvs_util_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(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_util_fill_pl - - module subroutine rmvs_util_fill_tp(self, inserts, lfill_list) - use swiftest_classes, only : swiftest_body - implicit none - class(rmvs_tp), 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_util_fill_tp - module subroutine rmvs_kick_getacch_tp(self, system, param, t, lbeg) use swiftest_classes, only : swiftest_nbody_system, swiftest_parameters implicit none @@ -162,10 +150,52 @@ module subroutine rmvs_setup_tp(self,n) integer, intent(in) :: n !! Number of test particles to allocate end subroutine rmvs_setup_tp + module subroutine rmvs_util_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(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_util_fill_pl + + module subroutine rmvs_util_fill_tp(self, inserts, lfill_list) + use swiftest_classes, only : swiftest_body + implicit none + class(rmvs_tp), 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_util_fill_tp + + module subroutine rmvs_util_sort_pl(self, sortby, ascending) + implicit none + class(rmvs_pl), intent(inout) :: self !! RMVS massive body object + character(*), intent(in) :: sortby !! Sorting attribute + logical, intent(in) :: ascending !! Logical flag indicating whether or not the sorting should be in ascending or descending order + end subroutine rmvs_util_sort_pl + + module subroutine rmvs_util_sort_tp(self, sortby, ascending) + implicit none + class(rmvs_tp), intent(inout) :: self !! RMVS test particle object + character(*), intent(in) :: sortby !! Sorting attribute + logical, intent(in) :: ascending !! Logical flag indicating whether or not the sorting should be in ascending or descending order + end subroutine rmvs_util_sort_tp + + module subroutine rmvs_util_sort_rearrange_pl(self, ind) + implicit none + class(rmvs_pl), intent(inout) :: self !! RMVS massive body object + integer(I4B), dimension(:), intent(in) :: ind !! Index array used to restructure the body (should contain all 1:n index values in the desired order) + end subroutine rmvs_util_sort_rearrange_pl + + module subroutine rmvs_util_sort_rearrange_tp(self, ind) + implicit none + class(rmvs_tp), intent(inout) :: self !! RMVS test particle object + integer(I4B), dimension(:), intent(in) :: ind !! Index array used to restructure the body (should contain all 1:n index values in the desired order) + end subroutine rmvs_util_sort_rearrange_tp + module subroutine rmvs_util_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_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_util_spill_pl diff --git a/src/modules/symba_classes.f90 b/src/modules/symba_classes.f90 index 7eab9e5c7..01fb7bbf4 100644 --- a/src/modules/symba_classes.f90 +++ b/src/modules/symba_classes.f90 @@ -316,13 +316,13 @@ end subroutine symba_setup_pl module subroutine symba_setup_pltpenc(self,n) implicit none - class(symba_pltpenc), intent(inout) :: self !! Symba pl-tp encounter structure + class(symba_pltpenc), intent(inout) :: self !! SyMBA pl-tp encounter structure integer, intent(in) :: n !! Number of encounters to allocate space for end subroutine symba_setup_pltpenc module subroutine symba_setup_plplenc(self,n) implicit none - class(symba_plplenc), intent(inout) :: self !! Symba pl-tp encounter structure + class(symba_plplenc), intent(inout) :: self !! SyMBA pl-tp encounter structure integer, intent(in) :: n !! Number of encounters to allocate space for end subroutine symba_setup_plplenc diff --git a/src/rmvs/rmvs_util.f90 b/src/rmvs/rmvs_util.f90 index 29d0b7b2a..745888a64 100644 --- a/src/rmvs/rmvs_util.f90 +++ b/src/rmvs/rmvs_util.f90 @@ -65,6 +65,136 @@ module subroutine rmvs_util_fill_tp(self, inserts, lfill_list) return end subroutine rmvs_util_fill_tp + module subroutine rmvs_util_sort_pl(self, sortby, ascending) + !! author: David A. Minton + !! + !! Sort a RMVS massive body object in-place. + !! sortby is a string indicating which array component to sort. + implicit none + ! Arguments + class(rmvs_pl), intent(inout) :: self !! RMVS massive body object + character(*), intent(in) :: sortby !! Sorting attribute + logical, intent(in) :: ascending !! Logical flag indicating whether or not the sorting should be in ascending or descending order + ! Internals + integer(I4B), dimension(self%nbody) :: ind + integer(I4B) :: direction + + if (ascending) then + direction = 1 + else + direction = -1 + end if + + associate(pl => self, npl => self%nbody) + select case(sortby) + case("nenc") + call util_sort(direction * pl%nenc(1:npl), ind(1:npl)) + case("tpenc1P") + call util_sort(direction * pl%tpenc1P(1:npl), ind(1:npl)) + case("plind") + call util_sort(direction * pl%plind(1:npl), ind(1:npl)) + case("outer", "inner", "planetocentric", "lplanetocentric") + write(*,*) 'Cannot sort by ' // trim(adjustl(sortby)) // '. Component not sortable!' + case default ! Look for components in the parent class + call whm_util_sort_pl(pl, sortby, ascending) + return + end select + + call pl%rearrange(ind) + + end associate + return + end subroutine rmvs_util_sort_pl + + module subroutine rmvs_util_sort_tp(self, sortby, ascending) + !! author: David A. Minton + !! + !! Sort a RMVS test particle object in-place. + !! sortby is a string indicating which array component to sort. + implicit none + ! Arguments + class(rmvs_tp), intent(inout) :: self !! RMVS test particle object + character(*), intent(in) :: sortby !! Sorting attribute + logical, intent(in) :: ascending !! Logical flag indicating whether or not the sorting should be in ascending or descending order + ! Internals + integer(I4B), dimension(self%nbody) :: ind + integer(I4B) :: direction + + if (ascending) then + direction = 1 + else + direction = -1 + end if + + associate(tp => self, ntp => self%nbody) + select case(sortby) + case("plperP") + call util_sort(direction * tp%plperP(1:ntp), ind(1:ntp)) + case("plencP") + call util_sort(direction * tp%plencP(1:ntp), ind(1:ntp)) + case("lperi", "cb_heliocentric", "xheliocentric", "index", "ipleP", "lplanetocentric") + write(*,*) 'Cannot sort by ' // trim(adjustl(sortby)) // '. Component not sortable!' + case default ! Look for components in the parent class (*NOTE whm_tp does not need its own sort method, so we go straight to the swiftest_tp method) + call util_sort_tp(tp, sortby, ascending) + return + end select + + call tp%rearrange(ind) + + end associate + return + end subroutine rmvs_util_sort_tp + + module subroutine rmvs_util_sort_rearrange_pl(self, ind) + !! author: David A. Minton + !! + !! Rearrange RMVS massive body structure in-place from an index list. + !! This is a helper utility used to make polymorphic sorting work on Swiftest structures. + implicit none + ! Arguments + class(rmvs_pl), intent(inout) :: self !! RMVS massive body object + integer(I4B), dimension(:), intent(in) :: ind !! Index array used to restructure the body (should contain all 1:n index values in the desired order) + ! Internals + class(rmvs_pl), allocatable :: pl_sorted !! Temporary holder for sorted body + integer(I4B) :: i + + associate(pl => self, npl => self%nbody) + call util_sort_rearrange_pl(pl,ind) + allocate(pl_sorted, source=self) + pl%eta(1:npl) = pl_sorted%eta(ind(1:npl)) + pl%xj(:,1:npl) = pl_sorted%xj(:,ind(1:npl)) + pl%vj(:,1:npl) = pl_sorted%vj(:,ind(1:npl)) + pl%muj(1:npl) = pl_sorted%muj(ind(1:npl)) + pl%ir3j(1:npl) = pl_sorted%ir3j(ind(1:npl)) + deallocate(pl_sorted) + end associate + return + end subroutine rmvs_util_sort_rearrange_pl + + module subroutine rmvs_util_sort_rearrange_tp(self, ind) + !! author: David A. Minton + !! + !! Rearrange RMVS test particle object in-place from an index list. + !! This is a helper utility used to make polymorphic sorting work on Swiftest structures. + implicit none + ! Arguments + class(rmvs_tp), intent(inout) :: self !! RMVS test particle object + integer(I4B), dimension(:), intent(in) :: ind !! Index array used to restructure the body (should contain all 1:n index values in the desired order) + ! Internals + class(rmvs_tp), allocatable :: tp_sorted !! Temporary holder for sorted body + + associate(tp => self, ntp => self%nbody) + call util_sort_rearrange_tp(tp,ind) + allocate(tp_sorted, source=self) + tp%lperi(1:ntp) = tp_sorted%lperi(ind(1:ntp)) + tp%plperP(1:ntp) = tp_sorted%plperP(ind(1:ntp)) + tp%plencP(1:ntp) = tp_sorted%plencP(ind(1:ntp)) + tp%xheliocentric(:,1:ntp) = tp_sorted%xheliocentric(:,ind(1:ntp)) + deallocate(tp_sorted) + end associate + return + end subroutine rmvs_util_sort_rearrange_tp + module subroutine rmvs_util_spill_pl(self, discards, lspill_list) !! author: David A. Minton !! diff --git a/src/symba/symba_setup.f90 b/src/symba/symba_setup.f90 index 51aaf69ba..8ae223228 100644 --- a/src/symba/symba_setup.f90 +++ b/src/symba/symba_setup.f90 @@ -50,7 +50,7 @@ module subroutine symba_setup_pltpenc(self, n) !! implicit none ! Arguments - class(symba_pltpenc), intent(inout) :: self !! Symba pl-tp encounter structure + class(symba_pltpenc), intent(inout) :: self !! SyMBA pl-tp encounter structure integer, intent(in) :: n !! Number of encounters to allocate space for self%nenc = n @@ -80,7 +80,7 @@ module subroutine symba_setup_plplenc(self,n) ! implicit none ! Arguments - class(symba_plplenc), intent(inout) :: self !! Symba pl-tp encounter structure + class(symba_plplenc), intent(inout) :: self !! SyMBA pl-tp encounter structure integer, intent(in) :: n !! Number of encounters to allocate space for call symba_setup_pltpenc(self, n) diff --git a/src/symba/symba_util.f90 b/src/symba/symba_util.f90 index 397cd789e..71db85cc3 100644 --- a/src/symba/symba_util.f90 +++ b/src/symba/symba_util.f90 @@ -99,7 +99,9 @@ module subroutine symba_util_sort_pl(self, sortby, ascending) call util_sort(direction * pl%peri(1:npl), ind(1:npl)) case("atp") call util_sort(direction * pl%atp(1:npl), ind(1:npl)) - case default + case("lcollision", "lencounter", "lmtiny", "nplm", "nplplm", "kin", "info") + write(*,*) 'Cannot sort by ' // trim(adjustl(sortby)) // '. Component not sortable!' + case default ! Look for components in the parent class call util_sort_pl(pl, sortby, ascending) return end select @@ -138,7 +140,7 @@ module subroutine symba_util_sort_tp(self, sortby, ascending) call util_sort(direction * tp%levelg(1:ntp), ind(1:ntp)) case("levelm") call util_sort(direction * tp%levelm(1:ntp), ind(1:ntp)) - case default + case default ! Look for components in the parent class call util_sort_tp(tp, sortby, ascending) return end select @@ -156,7 +158,7 @@ module subroutine symba_util_sort_rearrange_pl(self, ind) !! This is a helper utility used to make polymorphic sorting work on Swiftest structures. implicit none ! Arguments - class(symba_pl), intent(inout) :: self !! Symba massive body object + class(symba_pl), intent(inout) :: self !! SyMBA massive body object integer(I4B), dimension(:), intent(in) :: ind !! Index array used to restructure the body (should contain all 1:n index values in the desired order) ! Internals class(symba_pl), allocatable :: pl_sorted !! Temporary holder for sorted body @@ -193,7 +195,7 @@ module subroutine symba_util_sort_rearrange_tp(self, ind) !! This is a helper utility used to make polymorphic sorting work on Swiftest structures. implicit none ! Arguments - class(symba_tp), intent(inout) :: self !! Symba massive body object + class(symba_tp), intent(inout) :: self !! SyMBA test particle object integer(I4B), dimension(:), intent(in) :: ind !! Index array used to restructure the body (should contain all 1:n index values in the desired order) ! Internals class(symba_tp), allocatable :: tp_sorted !! Temporary holder for sorted body diff --git a/src/util/util_sort.f90 b/src/util/util_sort.f90 index a3e2d37a4..c81ced32d 100644 --- a/src/util/util_sort.f90 +++ b/src/util/util_sort.f90 @@ -35,8 +35,10 @@ module subroutine util_sort_body(self, sortby, ascending) call util_sort(direction * body%capom(1:n), ind(1:n)) case("mu") call util_sort(direction * body%mu(1:n), ind(1:n)) + case("lfirst", "nbody","xh", "vh", "xb", "vb", "ah", "aobl", "atide", "agr") + write(*,*) 'Cannot sort by ' // trim(adjustl(sortby)) // '. Component not sortable!' case default - write(*,*) 'Cannot sort structure by component ' // trim(adjustl(sortby)) + write(*,*) 'Cannot sort by ' // trim(adjustl(sortby)) // '. Component not found!' return end select @@ -83,7 +85,9 @@ module subroutine util_sort_pl(self, sortby, ascending) call util_sort(direction * pl%Q(1:npl), ind(1:npl)) case("tlag") call util_sort(direction * pl%tlag(1:npl), ind(1:npl)) - case default + case("xbeg", "xend", "vbeg", "Ip", "rot", "k_plpl", "nplpl") + write(*,*) 'Cannot sort by ' // trim(adjustl(sortby)) // '. Component not sortable!' + case default ! Look for components in the parent class call util_sort_body(pl, sortby, ascending) return end select @@ -121,7 +125,9 @@ module subroutine util_sort_tp(self, sortby, ascending) call util_sort(direction * tp%peri(1:ntp), ind(1:ntp)) case("atp") call util_sort(direction * tp%atp(1:ntp), ind(1:ntp)) - case default + case("isperi") + write(*,*) 'Cannot sort by ' // trim(adjustl(sortby)) // '. Component not sortable!' + case default ! Look for components in the parent class call util_sort_body(tp, sortby, ascending) return end select diff --git a/src/whm/whm_util.f90 b/src/whm/whm_util.f90 index 9ebeb0b3f..7e1b02f50 100644 --- a/src/whm/whm_util.f90 +++ b/src/whm/whm_util.f90 @@ -140,6 +140,8 @@ module subroutine whm_util_sort_pl(self, sortby, ascending) call util_sort(direction * pl%muj(1:npl), ind(1:npl)) case("ir3j") call util_sort(direction * pl%ir3j(1:npl), ind(1:npl)) + case("xj", "vj") + write(*,*) 'Cannot sort by ' // trim(adjustl(sortby)) // '. Component not sortable!' case default call util_sort_pl(pl, sortby, ascending) return From bc943d3c9102302782238e41eb4327ad9367143a Mon Sep 17 00:00:00 2001 From: David A Minton Date: Wed, 28 Jul 2021 11:23:11 -0400 Subject: [PATCH 7/7] Fixed swiftest_body sort method to accept all defined componenbts. Added sort method to initilization to ensure bodies are in correct heliocentric order before computing Jacobis --- src/util/util_sort.f90 | 6 +++++- src/whm/whm_setup.f90 | 7 +++++-- 2 files changed, 10 insertions(+), 3 deletions(-) diff --git a/src/util/util_sort.f90 b/src/util/util_sort.f90 index c81ced32d..c08343cee 100644 --- a/src/util/util_sort.f90 +++ b/src/util/util_sort.f90 @@ -25,6 +25,10 @@ module subroutine util_sort_body(self, sortby, ascending) select case(sortby) case("id") call util_sort(direction * body%id(1:n), ind(1:n)) + case("status") + call util_sort(direction * body%status(1:n), ind(1:n)) + case("ir3h") + call util_sort(direction * body%ir3h(1:n), ind(1:n)) case("a") call util_sort(direction * body%a(1:n), ind(1:n)) case("e") @@ -35,7 +39,7 @@ module subroutine util_sort_body(self, sortby, ascending) call util_sort(direction * body%capom(1:n), ind(1:n)) case("mu") call util_sort(direction * body%mu(1:n), ind(1:n)) - case("lfirst", "nbody","xh", "vh", "xb", "vb", "ah", "aobl", "atide", "agr") + case("lfirst", "nbody", "ldiscard", "xh", "vh", "xb", "vb", "ah", "aobl", "atide", "agr") write(*,*) 'Cannot sort by ' // trim(adjustl(sortby)) // '. Component not sortable!' case default write(*,*) 'Cannot sort by ' // trim(adjustl(sortby)) // '. Component not found!' diff --git a/src/whm/whm_setup.f90 b/src/whm/whm_setup.f90 index 940ba0b26..733adc871 100644 --- a/src/whm/whm_setup.f90 +++ b/src/whm/whm_setup.f90 @@ -53,8 +53,8 @@ module subroutine whm_util_set_mu_eta_pl(self, cb) !! Sets the Jacobi mass value eta for all massive bodies implicit none ! Arguments - class(whm_pl), intent(inout) :: self !! Swiftest system object - class(swiftest_cb), intent(inout) :: cb !! Swiftest central body particle data structure + class(whm_pl), intent(inout) :: self !! WHM system object + class(swiftest_cb), intent(inout) :: cb !! Swiftest central body object ! Internals integer(I4B) :: i @@ -82,6 +82,9 @@ module subroutine whm_setup_initialize_system(self, param) class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters call setup_initialize_system(self, param) + ! First we need to make sure that the massive bodies are sorted by heliocentric distance before computing jacobies + call util_set_ir3h(self%pl) + call self%pl%sort("ir3h", ascending=.false.) ! Make sure that the discard list gets allocated initially call self%tp_discards%setup(self%tp%nbody) call self%pl%set_mu(self%cb)