From 5e4bff866912c41098e60d2a2d5dce95b5e75d50 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Fri, 23 Dec 2022 06:31:56 -0500 Subject: [PATCH] Consolidated setup and util submodules --- src/CMakeLists.txt | 18 +- src/collision/collision_module.f90 | 20 +- src/collision/collision_setup.f90 | 68 --- src/collision/collision_util.f90 | 55 +- src/encounter/encounter_module.f90 | 12 +- src/encounter/encounter_setup.f90 | 108 ---- src/encounter/encounter_util.f90 | 93 ++++ src/fraggle/fraggle_module.f90 | 66 +-- src/fraggle/fraggle_set.f90 | 155 ------ src/fraggle/fraggle_setup.f90 | 31 -- src/fraggle/fraggle_util.f90 | 203 +++++++- src/helio/helio_module.f90 | 6 +- src/helio/{helio_setup.f90 => helio_util.f90} | 18 +- src/rmvs/rmvs_module.f90 | 118 +++-- src/rmvs/rmvs_setup.f90 | 166 ------ src/rmvs/rmvs_util.f90 | 226 +++++--- src/swiftest/swiftest_driver.f90 | 2 +- src/swiftest/swiftest_module.f90 | 36 +- src/swiftest/swiftest_setup.f90 | 385 -------------- src/swiftest/swiftest_util.f90 | 490 ++++++++++++++++-- src/symba/symba_module.f90 | 18 +- src/symba/symba_setup.f90 | 98 ---- src/symba/symba_util.f90 | 84 +++ src/whm/whm_module.f90 | 14 +- src/whm/whm_setup.f90 | 100 ---- src/whm/whm_util.f90 | 87 ++++ 26 files changed, 1274 insertions(+), 1403 deletions(-) delete mode 100644 src/collision/collision_setup.f90 delete mode 100644 src/encounter/encounter_setup.f90 delete mode 100644 src/fraggle/fraggle_set.f90 delete mode 100644 src/fraggle/fraggle_setup.f90 rename src/helio/{helio_setup.f90 => helio_util.f90} (72%) delete mode 100644 src/rmvs/rmvs_setup.f90 delete mode 100644 src/swiftest/swiftest_setup.f90 delete mode 100644 src/symba/symba_setup.f90 delete mode 100644 src/whm/whm_setup.f90 diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 86ffe6ee3..adfde33c7 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -44,51 +44,43 @@ SET(FAST_MATH_FILES ${SRC}/collision/collision_io.f90 ${SRC}/collision/collision_regime.f90 ${SRC}/collision/collision_resolve.f90 - ${SRC}/collision/collision_setup.f90 ${SRC}/collision/collision_util.f90 ${SRC}/encounter/encounter_check.f90 ${SRC}/encounter/encounter_io.f90 - ${SRC}/encounter/encounter_setup.f90 ${SRC}/encounter/encounter_util.f90 ${SRC}/fraggle/fraggle_generate.f90 - ${SRC}/fraggle/fraggle_set.f90 - ${SRC}/fraggle/fraggle_setup.f90 ${SRC}/fraggle/fraggle_util.f90 ${SRC}/helio/helio_drift.f90 ${SRC}/helio/helio_gr.f90 - ${SRC}/helio/helio_setup.f90 ${SRC}/helio/helio_step.f90 + ${SRC}/helio/helio_util.f90 ${SRC}/netcdf_io/netcdf_io_implementations.f90 ${SRC}/operator/operator_cross.f90 ${SRC}/operator/operator_mag.f90 ${SRC}/operator/operator_unit.f90 ${SRC}/rmvs/rmvs_discard.f90 ${SRC}/rmvs/rmvs_encounter_check.f90 - ${SRC}/rmvs/rmvs_setup.f90 ${SRC}/rmvs/rmvs_step.f90 ${SRC}/rmvs/rmvs_util.f90 ${SRC}/swiftest/swiftest_discard.f90 - ${SRC}/swiftest/swiftest_io.f90 - ${SRC}/swiftest/swiftest_obl.f90 - ${SRC}/swiftest/swiftest_util.f90 ${SRC}/swiftest/swiftest_drift.f90 - ${SRC}/swiftest/swiftest_orbel.f90 ${SRC}/swiftest/swiftest_gr.f90 + ${SRC}/swiftest/swiftest_io.f90 ${SRC}/swiftest/swiftest_kick.f90 - ${SRC}/swiftest/swiftest_setup.f90 + ${SRC}/swiftest/swiftest_obl.f90 + ${SRC}/swiftest/swiftest_orbel.f90 + ${SRC}/swiftest/swiftest_util.f90 ${SRC}/symba/symba_discard.f90 ${SRC}/symba/symba_drift.f90 ${SRC}/symba/symba_encounter_check.f90 ${SRC}/symba/symba_gr.f90 ${SRC}/symba/symba_io.f90 - ${SRC}/symba/symba_setup.f90 ${SRC}/symba/symba_step.f90 ${SRC}/symba/symba_util.f90 ${SRC}/walltime/walltime_implementations.f90 ${SRC}/whm/whm_coord.f90 ${SRC}/whm/whm_drift.f90 ${SRC}/whm/whm_gr.f90 - ${SRC}/whm/whm_setup.f90 ${SRC}/whm/whm_step.f90 ${SRC}/whm/whm_util.f90 ${SRC}/swiftest/swiftest_driver.f90 diff --git a/src/collision/collision_module.f90 b/src/collision/collision_module.f90 index 98463fb93..40ef8b012 100644 --- a/src/collision/collision_module.f90 +++ b/src/collision/collision_module.f90 @@ -130,9 +130,9 @@ module collision real(DP), dimension(2) :: pe !! Before/after potential energy real(DP), dimension(2) :: Etot !! Before/after total nbody_system energy contains - procedure :: setup => collision_setup_system !! Initializer for the encounter collision system and the before/after snapshots - procedure :: setup_impactors => collision_setup_impactors_system !! Initializer for the impactors for the encounter collision system. Deallocates old impactors before creating new ones - procedure :: setup_fragments => collision_setup_fragments_system !! Initializer for the fragments of the collision system. + procedure :: setup => collision_util_setup_system !! Initializer for the encounter collision system and the before/after snapshots + procedure :: setup_impactors => collision_util_setup_impactors_system !! Initializer for the impactors for the encounter collision system. Deallocates old impactors before creating new ones + procedure :: setup_fragments => collision_util_setup_fragments_system !! Initializer for the fragments of the collision system. procedure :: add_fragments => collision_util_add_fragments_to_system !! Add fragments to nbody_system procedure :: construct_temporary_system => collision_util_construct_temporary_system !! Constructs temporary n-body nbody_system in order to compute pre- and post-impact energy and momentum procedure :: get_energy_and_momentum => collision_util_get_energy_momentum !! Calculates total nbody_system energy in either the pre-collision outcome state (lbefore = .true.) or the post-collision outcome state (lbefore = .false.) @@ -352,22 +352,22 @@ module subroutine collision_util_set_coordinate_system(self) class(collision_merge), intent(inout) :: self !! Collisional nbody_system end subroutine collision_util_set_coordinate_system - module subroutine collision_setup_system(self, nbody_system) + module subroutine collision_util_setup_system(self, nbody_system) implicit none class(collision_merge), intent(inout) :: self !! Encounter collision system object class(base_nbody_system), intent(in) :: nbody_system !! Current nbody system. Used as a mold for the before/after snapshots - end subroutine collision_setup_system + end subroutine collision_util_setup_system - module subroutine collision_setup_impactors_system(self) + module subroutine collision_util_setup_impactors_system(self) implicit none class(collision_merge), intent(inout) :: self !! Encounter collision system object - end subroutine collision_setup_impactors_system + end subroutine collision_util_setup_impactors_system - module subroutine collision_setup_fragments_system(self, nfrag) + module subroutine collision_util_setup_fragments_system(self, nfrag) implicit none class(collision_merge), intent(inout) :: self !! Encounter collision system object integer(I4B), intent(in) :: nfrag !! Number of fragments to create - end subroutine collision_setup_fragments_system + end subroutine collision_util_setup_fragments_system module subroutine collision_util_add_fragments_to_system(self, nbody_system, param) implicit none @@ -438,7 +438,6 @@ end subroutine collision_util_snapshot contains - subroutine collision_final_fragments(self) !! author: David A. Minton !! @@ -466,6 +465,7 @@ subroutine collision_final_impactors(self) return end subroutine collision_final_impactors + subroutine collision_final_netcdf_parameters(self) !! author: David A. Minton !! diff --git a/src/collision/collision_setup.f90 b/src/collision/collision_setup.f90 deleted file mode 100644 index 2c4880cb0..000000000 --- a/src/collision/collision_setup.f90 +++ /dev/null @@ -1,68 +0,0 @@ -!! Copyright 2022 - David Minton, Carlisle Wishard, Jennifer Pouplin, Jake Elliott, & Dana Singh -!! This file is part of Swiftest. -!! Swiftest is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License -!! as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version. -!! Swiftest is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty -!! of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. -!! You should have received a copy of the GNU General Public License along with Swiftest. -!! If not, see: https://www.gnu.org/licenses. - -submodule (collision) s_collision_setup - use swiftest -contains - - module subroutine collision_setup_system(self, nbody_system) - !! author: David A. Minton - !! - !! Initializer for the encounter collision system. Sets up impactors and the before/after snapshots, - !! but not fragments. Those are setup later when the number of fragments is known. - implicit none - ! Arguments - class(collision_merge), intent(inout) :: self !! Encounter collision system object - class(base_nbody_system), intent(in) :: nbody_system !! Current nbody system. Used as a mold for the before/after snapshots - - call self%setup_impactors() - if (allocated(self%before)) deallocate(self%before) - if (allocated(self%after)) deallocate(self%after) - - allocate(self%before, mold=nbody_system) - allocate(self%after, mold=nbody_system) - - return - end subroutine collision_setup_system - - - module subroutine collision_setup_impactors_system(self) - !! author: David A. Minton - !! - !! Initializer for the impactors for the encounter collision system. Deallocates old impactors before creating new ones - implicit none - ! Arguments - class(collision_merge), intent(inout) :: self !! Encounter collision system object - - if (allocated(self%impactors)) deallocate(self%impactors) - allocate(collision_impactors :: self%impactors) - - return - end subroutine collision_setup_impactors_system - - - module subroutine collision_setup_fragments_system(self, nfrag) - !! author: David A. Minton - !! - !! Initializer for the fragments of the collision system. - implicit none - ! Arguments - class(collision_merge), intent(inout) :: self !! Encounter collision system object - integer(I4B), intent(in) :: nfrag !! Number of fragments to create - - if (allocated(self%fragments)) deallocate(self%fragments) - allocate(collision_fragments(nfrag) :: self%fragments) - self%fragments%nbody = nfrag - - return - end subroutine collision_setup_fragments_system - -end submodule s_collision_setup - - \ No newline at end of file diff --git a/src/collision/collision_util.f90 b/src/collision/collision_util.f90 index 410090107..62e275d24 100644 --- a/src/collision/collision_util.f90 +++ b/src/collision/collision_util.f90 @@ -88,7 +88,7 @@ module subroutine collision_util_construct_temporary_system(self, nbody_system, if (allocated(tmpparam)) deallocate(tmpparam) if (allocated(tmpsys)) deallocate(tmpsys) allocate(tmpparam, source=param) - call swiftest_setup_construct_system(tmpsys_local, tmpparam_local) + call swiftest_util_setup_construct_system(tmpsys_local, tmpparam_local) ! No test particles necessary for energy/momentum calcs call tmpsys_local%tp%setup(0, tmpparam_local) @@ -571,6 +571,59 @@ module subroutine collision_util_set_mass_dist(self, param) end subroutine collision_util_set_mass_dist + module subroutine collision_util_setup_system(self, nbody_system) + !! author: David A. Minton + !! + !! Initializer for the encounter collision system. Sets up impactors and the before/after snapshots, + !! but not fragments. Those are setup later when the number of fragments is known. + implicit none + ! Arguments + class(collision_merge), intent(inout) :: self !! Encounter collision system object + class(base_nbody_system), intent(in) :: nbody_system !! Current nbody system. Used as a mold for the before/after snapshots + + call self%setup_impactors() + if (allocated(self%before)) deallocate(self%before) + if (allocated(self%after)) deallocate(self%after) + + allocate(self%before, mold=nbody_system) + allocate(self%after, mold=nbody_system) + + return + end subroutine collision_util_setup_system + + + module subroutine collision_util_setup_impactors_system(self) + !! author: David A. Minton + !! + !! Initializer for the impactors for the encounter collision system. Deallocates old impactors before creating new ones + implicit none + ! Arguments + class(collision_merge), intent(inout) :: self !! Encounter collision system object + + if (allocated(self%impactors)) deallocate(self%impactors) + allocate(collision_impactors :: self%impactors) + + return + end subroutine collision_util_setup_impactors_system + + + module subroutine collision_util_setup_fragments_system(self, nfrag) + !! author: David A. Minton + !! + !! Initializer for the fragments of the collision system. + implicit none + ! Arguments + class(collision_merge), intent(inout) :: self !! Encounter collision system object + integer(I4B), intent(in) :: nfrag !! Number of fragments to create + + if (allocated(self%fragments)) deallocate(self%fragments) + allocate(collision_fragments(nfrag) :: self%fragments) + self%fragments%nbody = nfrag + + return + end subroutine collision_util_setup_fragments_system + + subroutine collision_util_save_snapshot(collision_history, snapshot) !! author: David A. Minton !! diff --git a/src/encounter/encounter_module.f90 b/src/encounter/encounter_module.f90 index abb3c3422..0f7ee1055 100644 --- a/src/encounter/encounter_module.f90 +++ b/src/encounter/encounter_module.f90 @@ -37,7 +37,7 @@ module encounter real(DP), dimension(:,:), allocatable :: v2 !! the velocity of body 2 in the encounter integer(I4B), dimension(:), allocatable :: level !! Recursion level (used in SyMBA) contains - procedure :: setup => encounter_setup_list !! A constructor that sets the number of encounters and allocates and initializes all arrays + procedure :: setup => encounter_util_setup_list !! A constructor that sets the number of encounters and allocates and initializes all arrays procedure :: append => encounter_util_append_list !! Appends elements from one structure to another procedure :: copy => encounter_util_copy_list !! Copies elements from the source encounter list into self. procedure :: dealloc => encounter_util_dealloc_list !! Deallocates all allocatables @@ -98,7 +98,7 @@ module encounter type encounter_bounding_box type(encounter_bounding_box_1D), dimension(SWEEPDIM) :: aabb contains - procedure :: setup => encounter_setup_aabb !! Setup a new axis-aligned bounding box structure + procedure :: setup => encounter_util_setup_aabb !! Setup a new axis-aligned bounding box structure procedure :: sweep_single => encounter_check_sweep_aabb_single_list !! Sweeps the sorted bounding box extents and returns the encounter candidates procedure :: sweep_double => encounter_check_sweep_aabb_double_list !! Sweeps the sorted bounding box extents and returns the encounter candidates generic :: sweep => sweep_single, sweep_double @@ -236,18 +236,18 @@ module subroutine encounter_io_netcdf_write_frame_snapshot(self, history, param) class(base_parameters), intent(inout) :: param !! Current run configuration parameters end subroutine encounter_io_netcdf_write_frame_snapshot - module subroutine encounter_setup_aabb(self, n, n_last) + module subroutine encounter_util_setup_aabb(self, n, n_last) implicit none class(encounter_bounding_box), intent(inout) :: self !! Swiftest encounter structure integer(I4B), intent(in) :: n !! Number of objects with bounding box extents integer(I4B), intent(in) :: n_last !! Number of objects with bounding box extents the previous time this was called - end subroutine encounter_setup_aabb + end subroutine encounter_util_setup_aabb - module subroutine encounter_setup_list(self, n) + module subroutine encounter_util_setup_list(self, n) implicit none class(encounter_list), intent(inout) :: self !! Swiftest encounter structure integer(I8B), intent(in) :: n !! Number of encounters to allocate space for - end subroutine encounter_setup_list + end subroutine encounter_util_setup_list module subroutine encounter_util_append_list(self, source, lsource_mask) implicit none diff --git a/src/encounter/encounter_setup.f90 b/src/encounter/encounter_setup.f90 deleted file mode 100644 index 3cf165129..000000000 --- a/src/encounter/encounter_setup.f90 +++ /dev/null @@ -1,108 +0,0 @@ -!! Copyright 2022 - David Minton, Carlisle Wishard, Jennifer Pouplin, Jake Elliott, & Dana Singh -!! This file is part of Swiftest. -!! Swiftest is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License -!! as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version. -!! Swiftest is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty -!! of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. -!! You should have received a copy of the GNU General Public License along with Swiftest. -!! If not, see: https://www.gnu.org/licenses. - -submodule (encounter) s_encounter_setup - use swiftest -contains - - module subroutine encounter_setup_aabb(self, n, n_last) - !! author: David A. Minton - !! - !! Sets up or modifies an axis-aligned bounding box structure. - implicit none - ! Arguments - class(encounter_bounding_box), intent(inout) :: self !! Swiftest encounter structure - integer(I4B), intent(in) :: n !! Number of objects with bounding box extents - integer(I4B), intent(in) :: n_last !! Number of objects with bounding box extents the previous time this was called - ! Internals - integer(I4B) :: next, next_last, k, dim - integer(I4B), dimension(:), allocatable :: itmp - - next = 2 * n - next_last = 2 * n_last - - if (n > n_last) then ! The number of bodies has grown. Resize and append the new bodies - do dim = 1, SWEEPDIM - allocate(itmp(next)) - if (n_last > 0) itmp(1:next_last) = self%aabb(dim)%ind(1:next_last) - call move_alloc(itmp, self%aabb(dim)%ind) - self%aabb(dim)%ind(next_last+1:next) = [(k, k = next_last+1, next)] - end do - else ! The number of bodies has gone down. Resize and chop of the old indices - do dim = 1, SWEEPDIM - allocate(itmp(next)) - itmp(1:next) = pack(self%aabb(dim)%ind(1:next_last), self%aabb(dim)%ind(1:next_last) <= next) - call move_alloc(itmp, self%aabb(dim)%ind) - end do - end if - - do dim = 1, SWEEPDIM - if (allocated(self%aabb(dim)%ibeg)) deallocate(self%aabb(dim)%ibeg) - allocate(self%aabb(dim)%ibeg(n)) - if (allocated(self%aabb(dim)%iend)) deallocate(self%aabb(dim)%iend) - allocate(self%aabb(dim)%iend(n)) - end do - - return - end subroutine encounter_setup_aabb - - - module subroutine encounter_setup_list(self, n) - !! author: David A. Minton - !! - !! A constructor that sets the number of encounters and allocates and initializes all arrays - !! - implicit none - ! Arguments - class(encounter_list), intent(inout) :: self !! Swiftest encounter structure - integer(I8B), intent(in) :: n !! Number of encounters to allocate space for - ! Internals - integer(I8B) :: i - - if (n < 0) return - call self%dealloc() - - self%nenc = n - if (n == 0_I8B) return - self%t = 0.0_DP - - allocate(self%tcollision(n)) - allocate(self%lvdotr(n)) - allocate(self%lclosest(n)) - allocate(self%status(n)) - allocate(self%index1(n)) - allocate(self%index2(n)) - allocate(self%id1(n)) - allocate(self%id2(n)) - allocate(self%r1(NDIM,n)) - allocate(self%r2(NDIM,n)) - allocate(self%v1(NDIM,n)) - allocate(self%v2(NDIM,n)) - allocate(self%level(n)) - - self%tcollision(:) = 0.0_DP - self%lvdotr(:) = .false. - self%lclosest(:) = .false. - self%status(:) = INACTIVE - self%index1(:) = 0 - self%index2(:) = 0 - self%id1(:) = 0 - self%id2(:) = 0 - self%r1(:,:) = 0.0_DP - self%r2(:,:) = 0.0_DP - self%v1(:,:) = 0.0_DP - self%v2(:,:) = 0.0_DP - self%level(:) = 0 - - return - end subroutine encounter_setup_list - -end submodule s_encounter_setup - - \ No newline at end of file diff --git a/src/encounter/encounter_util.f90 b/src/encounter/encounter_util.f90 index 06648ca9f..45e6f48fb 100644 --- a/src/encounter/encounter_util.f90 +++ b/src/encounter/encounter_util.f90 @@ -281,6 +281,99 @@ module subroutine encounter_util_resize_list(self, nnew) end subroutine encounter_util_resize_list + module subroutine encounter_util_setup_aabb(self, n, n_last) + !! author: David A. Minton + !! + !! Sets up or modifies an axis-aligned bounding box structure. + implicit none + ! Arguments + class(encounter_bounding_box), intent(inout) :: self !! Swiftest encounter structure + integer(I4B), intent(in) :: n !! Number of objects with bounding box extents + integer(I4B), intent(in) :: n_last !! Number of objects with bounding box extents the previous time this was called + ! Internals + integer(I4B) :: next, next_last, k, dim + integer(I4B), dimension(:), allocatable :: itmp + + next = 2 * n + next_last = 2 * n_last + + if (n > n_last) then ! The number of bodies has grown. Resize and append the new bodies + do dim = 1, SWEEPDIM + allocate(itmp(next)) + if (n_last > 0) itmp(1:next_last) = self%aabb(dim)%ind(1:next_last) + call move_alloc(itmp, self%aabb(dim)%ind) + self%aabb(dim)%ind(next_last+1:next) = [(k, k = next_last+1, next)] + end do + else ! The number of bodies has gone down. Resize and chop of the old indices + do dim = 1, SWEEPDIM + allocate(itmp(next)) + itmp(1:next) = pack(self%aabb(dim)%ind(1:next_last), self%aabb(dim)%ind(1:next_last) <= next) + call move_alloc(itmp, self%aabb(dim)%ind) + end do + end if + + do dim = 1, SWEEPDIM + if (allocated(self%aabb(dim)%ibeg)) deallocate(self%aabb(dim)%ibeg) + allocate(self%aabb(dim)%ibeg(n)) + if (allocated(self%aabb(dim)%iend)) deallocate(self%aabb(dim)%iend) + allocate(self%aabb(dim)%iend(n)) + end do + + return + end subroutine encounter_util_setup_aabb + + + module subroutine encounter_util_setup_list(self, n) + !! author: David A. Minton + !! + !! A constructor that sets the number of encounters and allocates and initializes all arrays + !! + implicit none + ! Arguments + class(encounter_list), intent(inout) :: self !! Swiftest encounter structure + integer(I8B), intent(in) :: n !! Number of encounters to allocate space for + ! Internals + integer(I8B) :: i + + if (n < 0) return + call self%dealloc() + + self%nenc = n + if (n == 0_I8B) return + self%t = 0.0_DP + + allocate(self%tcollision(n)) + allocate(self%lvdotr(n)) + allocate(self%lclosest(n)) + allocate(self%status(n)) + allocate(self%index1(n)) + allocate(self%index2(n)) + allocate(self%id1(n)) + allocate(self%id2(n)) + allocate(self%r1(NDIM,n)) + allocate(self%r2(NDIM,n)) + allocate(self%v1(NDIM,n)) + allocate(self%v2(NDIM,n)) + allocate(self%level(n)) + + self%tcollision(:) = 0.0_DP + self%lvdotr(:) = .false. + self%lclosest(:) = .false. + self%status(:) = INACTIVE + self%index1(:) = 0 + self%index2(:) = 0 + self%id1(:) = 0 + self%id2(:) = 0 + self%r1(:,:) = 0.0_DP + self%r2(:,:) = 0.0_DP + self%v1(:,:) = 0.0_DP + self%v2(:,:) = 0.0_DP + self%level(:) = 0 + + return + end subroutine encounter_util_setup_list + + module subroutine encounter_util_spill_list(self, discards, lspill_list, ldestructive) !! author: David A. Minton !! diff --git a/src/fraggle/fraggle_module.f90 b/src/fraggle/fraggle_module.f90 index 2471b9975..72c0e5dc3 100644 --- a/src/fraggle/fraggle_module.f90 +++ b/src/fraggle/fraggle_module.f90 @@ -44,10 +44,10 @@ module fraggle real(DP) :: Lscale = 1.0_DP !! Angular momentum scale factor (a convenience unit that is derived from dscale, tscale, and mscale) contains procedure :: generate => fraggle_generate_system !! Generates a system of fragments in barycentric coordinates that conserves energy and momentum. - procedure :: set_budgets => fraggle_set_budgets !! Sets the energy and momentum budgets of the fragments based on the collider value - procedure :: set_natural_scale => fraggle_set_natural_scale_factors !! Scales dimenional quantities to ~O(1) with respect to the collisional system. - procedure :: set_original_scale => fraggle_set_original_scale_factors !! Restores dimenional quantities back to the original system units - procedure :: setup_fragments => fraggle_setup_fragments_system !! Initializer for the fragments of the collision system. + procedure :: set_budgets => fraggle_util_set_budgets !! Sets the energy and momentum budgets of the fragments based on the collider value + procedure :: set_natural_scale => fraggle_util_set_natural_scale_factors !! Scales dimenional quantities to ~O(1) with respect to the collisional system. + procedure :: set_original_scale => fraggle_util_set_original_scale_factors !! Restores dimenional quantities back to the original system units + procedure :: setup_fragments => fraggle_util_setup_fragments_system !! Initializer for the fragments of the collision system. procedure :: construct_temporary_system => fraggle_util_construct_temporary_system !! Constructs temporary n-body system in order to compute pre- and post-impact energy and momentum procedure :: reset => fraggle_util_reset_system !! Deallocates all allocatables final :: fraggle_final_system !! Finalizer will deallocate all allocatables @@ -89,31 +89,11 @@ module subroutine fraggle_generate_system(self, nbody_system, param, t) real(DP), intent(in) :: t !! Time of collision end subroutine fraggle_generate_system - module subroutine fraggle_set_budgets(self) - implicit none - class(collision_fraggle), intent(inout) :: self !! Fraggle collision system object - end subroutine fraggle_set_budgets - - module subroutine fraggle_set_natural_scale_factors(self) - implicit none - class(collision_fraggle), intent(inout) :: self !! Fraggle collision system object - end subroutine fraggle_set_natural_scale_factors - - module subroutine fraggle_set_original_scale_factors(self) - implicit none - class(collision_fraggle), intent(inout) :: self !! Fraggle collision system object - end subroutine fraggle_set_original_scale_factors - - module subroutine fraggle_setup_fragments_system(self, nfrag) + module subroutine fraggle_util_setup_fragments_system(self, nfrag) implicit none class(collision_fraggle), intent(inout) :: self !! Encounter collision system object integer(I4B), intent(in) :: nfrag !! Number of fragments to create - end subroutine fraggle_setup_fragments_system - - module subroutine fraggle_util_get_angular_momentum(self) - implicit none - class(fraggle_fragments(*)), intent(inout) :: self !! Fraggle fragment system object - end subroutine fraggle_util_get_angular_momentum + end subroutine fraggle_util_setup_fragments_system module subroutine fraggle_util_construct_temporary_system(self, nbody_system, param, tmpsys, tmpparam) implicit none @@ -124,10 +104,10 @@ module subroutine fraggle_util_construct_temporary_system(self, nbody_system, pa class(base_parameters), allocatable, intent(out) :: tmpparam !! Output temporary configuration run parameters end subroutine fraggle_util_construct_temporary_system - module subroutine fraggle_final_impactors(self) + module subroutine fraggle_util_get_angular_momentum(self) implicit none - type(collision_impactors), intent(inout) :: self !! Fraggle impactors object - end subroutine fraggle_final_impactors + class(fraggle_fragments(*)), intent(inout) :: self !! Fraggle fragment system object + end subroutine fraggle_util_get_angular_momentum module subroutine fraggle_util_reset_fragments(self) implicit none @@ -148,6 +128,22 @@ module subroutine fraggle_util_restructure(self, impactors, try, f_spin, r_max_s real(DP), intent(inout) :: r_max_start !! The maximum radial distance that the position calculation starts with. This increases after a failed attempt end subroutine fraggle_util_restructure + module subroutine fraggle_util_set_budgets(self) + implicit none + class(collision_fraggle), intent(inout) :: self !! Fraggle collision system object + end subroutine fraggle_util_set_budgets + + module subroutine fraggle_util_set_natural_scale_factors(self) + implicit none + class(collision_fraggle), intent(inout) :: self !! Fraggle collision system object + end subroutine fraggle_util_set_natural_scale_factors + + module subroutine fraggle_util_set_original_scale_factors(self) + implicit none + class(collision_fraggle), intent(inout) :: self !! Fraggle collision system object + end subroutine fraggle_util_set_original_scale_factors + + module subroutine fraggle_util_shift_vector_to_origin(m_frag, vec_frag) implicit none real(DP), dimension(:), intent(in) :: m_frag !! Fragment masses @@ -181,6 +177,18 @@ subroutine fraggle_final_fragments(self) end subroutine fraggle_final_fragments + subroutine fraggle_final_impactors(self) + !! author: David A. Minton + !! + !! Finalizer will deallocate all allocatables + implicit none + ! Arguments + type(collision_impactors), intent(inout) :: self !! Fraggle impactors object + call self%reset() + return + end subroutine fraggle_final_impactors + + subroutine fraggle_final_system(self) !! author: David A. Minton !! diff --git a/src/fraggle/fraggle_set.f90 b/src/fraggle/fraggle_set.f90 deleted file mode 100644 index 03f2e5e88..000000000 --- a/src/fraggle/fraggle_set.f90 +++ /dev/null @@ -1,155 +0,0 @@ -!! Copyright 2022 - David Minton, Carlisle Wishard, Jennifer Pouplin, Jake Elliott, & Dana Singh -!! This file is part of Swiftest. -!! Swiftest is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License -!! as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version. -!! Swiftest is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty -!! of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. -!! You should have received a copy of the GNU General Public License along with Swiftest. -!! If not, see: https://www.gnu.org/licenses. - -submodule(fraggle) s_fraggle_set - use swiftest - use symba -contains - - module subroutine fraggle_set_budgets(self) - !! author: David A. Minton - !! - !! Sets the energy and momentum budgets of the fragments based on the collider values and the before/after values of energy and momentum - implicit none - ! Arguments - class(collision_fraggle), intent(inout) :: self !! Fraggle collision system object - ! Internals - real(DP) :: dEtot - real(DP), dimension(NDIM) :: dL - - associate(impactors => self%impactors) - select type(fragments => self%fragments) - class is (fraggle_fragments(*)) - - dEtot = self%Etot(2) - self%Etot(1) - dL(:) = self%Ltot(:,2) - self%Ltot(:,1) - - fragments%L_budget(:) = -dL(:) - fragments%ke_budget = -(dEtot - 0.5_DP * fragments%mtot * dot_product(impactors%vbcom(:), impactors%vbcom(:))) - impactors%Qloss - - end select - end associate - - return - end subroutine fraggle_set_budgets - - - module subroutine fraggle_set_natural_scale_factors(self) - !! author: David A. Minton - !! - !! Scales dimenional quantities to ~O(1) with respect to the collisional system. - !! This scaling makes it easier for the non-linear minimization to converge on a solution - implicit none - ! Arguments - class(collision_fraggle), intent(inout) :: self !! Fraggle collision system object - ! Internals - integer(I4B) :: i - - associate(collision_merge => self, fragments => self%fragments, impactors => self%impactors) - ! Set scale factors - collision_merge%Escale = 0.5_DP * ( impactors%mass(1) * dot_product(impactors%vb(:,1), impactors%vb(:,1)) & - + impactors%mass(2) * dot_product(impactors%vb(:,2), impactors%vb(:,2))) - collision_merge%dscale = sum(impactors%radius(:)) - collision_merge%mscale = fragments%mtot - collision_merge%vscale = sqrt(collision_merge%Escale / collision_merge%mscale) - collision_merge%tscale = collision_merge%dscale / collision_merge%vscale - collision_merge%Lscale = collision_merge%mscale * collision_merge%dscale * collision_merge%vscale - - ! Scale all dimensioned quantities of impactors and fragments - impactors%rbcom(:) = impactors%rbcom(:) / collision_merge%dscale - impactors%vbcom(:) = impactors%vbcom(:) / collision_merge%vscale - impactors%rbimp(:) = impactors%rbimp(:) / collision_merge%dscale - impactors%rb(:,:) = impactors%rb(:,:) / collision_merge%dscale - impactors%vb(:,:) = impactors%vb(:,:) / collision_merge%vscale - impactors%mass(:) = impactors%mass(:) / collision_merge%mscale - impactors%radius(:) = impactors%radius(:) / collision_merge%dscale - impactors%Lspin(:,:) = impactors%Lspin(:,:) / collision_merge%Lscale - impactors%Lorbit(:,:) = impactors%Lorbit(:,:) / collision_merge%Lscale - - do i = 1, 2 - impactors%rot(:,i) = impactors%Lspin(:,i) / (impactors%mass(i) * impactors%radius(i)**2 * impactors%Ip(3, i)) - end do - - fragments%mtot = fragments%mtot / collision_merge%mscale - fragments%mass = fragments%mass / collision_merge%mscale - fragments%radius = fragments%radius / collision_merge%dscale - impactors%Qloss = impactors%Qloss / collision_merge%Escale - end associate - - return - end subroutine fraggle_set_natural_scale_factors - - - module subroutine fraggle_set_original_scale_factors(self) - !! author: David A. Minton - !! - !! Restores dimenional quantities back to the system units - use, intrinsic :: ieee_exceptions - implicit none - ! Arguments - class(collision_fraggle), intent(inout) :: self !! Fraggle fragment system object - ! Internals - integer(I4B) :: i - logical, dimension(size(IEEE_ALL)) :: fpe_halting_modes - - call ieee_get_halting_mode(IEEE_ALL,fpe_halting_modes) ! Save the current halting modes so we can turn them off temporarily - call ieee_set_halting_mode(IEEE_ALL,.false.) - - associate(collision_merge => self, fragments => self%fragments, impactors => self%impactors) - - ! Restore scale factors - impactors%rbcom(:) = impactors%rbcom(:) * collision_merge%dscale - impactors%vbcom(:) = impactors%vbcom(:) * collision_merge%vscale - impactors%rbimp(:) = impactors%rbimp(:) * collision_merge%dscale - - impactors%mass = impactors%mass * collision_merge%mscale - impactors%radius = impactors%radius * collision_merge%dscale - impactors%rb = impactors%rb * collision_merge%dscale - impactors%vb = impactors%vb * collision_merge%vscale - impactors%Lspin = impactors%Lspin * collision_merge%Lscale - do i = 1, 2 - impactors%rot(:,i) = impactors%Lspin(:,i) * (impactors%mass(i) * impactors%radius(i)**2 * impactors%Ip(3, i)) - end do - - fragments%mtot = fragments%mtot * collision_merge%mscale - fragments%mass = fragments%mass * collision_merge%mscale - fragments%radius = fragments%radius * collision_merge%dscale - fragments%rot = fragments%rot / collision_merge%tscale - fragments%rc = fragments%rc * collision_merge%dscale - fragments%vc = fragments%vc * collision_merge%vscale - - do i = 1, fragments%nbody - fragments%rb(:, i) = fragments%rc(:, i) + impactors%rbcom(:) - fragments%vb(:, i) = fragments%vc(:, i) + impactors%vbcom(:) - end do - - impactors%Qloss = impactors%Qloss * collision_merge%Escale - - collision_merge%Lorbit(:,:) = collision_merge%Lorbit(:,:) * collision_merge%Lscale - collision_merge%Lspin(:,:) = collision_merge%Lspin(:,:) * collision_merge%Lscale - collision_merge%Ltot(:,:) = collision_merge%Ltot(:,:) * collision_merge%Lscale - collision_merge%ke_orbit(:) = collision_merge%ke_orbit(:) * collision_merge%Escale - collision_merge%ke_spin(:) = collision_merge%ke_spin(:) * collision_merge%Escale - collision_merge%pe(:) = collision_merge%pe(:) * collision_merge%Escale - collision_merge%Etot(:) = collision_merge%Etot(:) * collision_merge%Escale - - collision_merge%mscale = 1.0_DP - collision_merge%dscale = 1.0_DP - collision_merge%vscale = 1.0_DP - collision_merge%tscale = 1.0_DP - collision_merge%Lscale = 1.0_DP - collision_merge%Escale = 1.0_DP - end associate - call ieee_set_halting_mode(IEEE_ALL,fpe_halting_modes) - - return - end subroutine fraggle_set_original_scale_factors - - -end submodule s_fraggle_set \ No newline at end of file diff --git a/src/fraggle/fraggle_setup.f90 b/src/fraggle/fraggle_setup.f90 deleted file mode 100644 index e54ca46fd..000000000 --- a/src/fraggle/fraggle_setup.f90 +++ /dev/null @@ -1,31 +0,0 @@ -!! Copyright 2022 - David Minton, Carlisle Wishard, Jennifer Pouplin, Jake Elliott, & Dana Singh -!! This file is part of Swiftest. -!! Swiftest is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License -!! as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version. -!! Swiftest is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty -!! of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. -!! You should have received a copy of the GNU General Public License along with Swiftest. -!! If not, see: https://www.gnu.org/licenses. - -submodule (fraggle) s_fraggle_setup - use swiftest - use symba -contains - - module subroutine fraggle_setup_fragments_system(self, nfrag) - !! author: David A. Minton - !! - !! Initializer for the fragments of the collision system. - implicit none - ! Arguments - class(collision_fraggle), intent(inout) :: self !! Encounter collision system object - integer(I4B), intent(in) :: nfrag !! Number of fragments to create - - if (allocated(self%fragments)) deallocate(self%fragments) - allocate(fraggle_fragments(nbody=nfrag) :: self%fragments) - self%fragments%nbody = nfrag - - return - end subroutine fraggle_setup_fragments_system - -end submodule s_fraggle_setup \ No newline at end of file diff --git a/src/fraggle/fraggle_util.f90 b/src/fraggle/fraggle_util.f90 index 07b5eef50..f3e6b87f6 100644 --- a/src/fraggle/fraggle_util.f90 +++ b/src/fraggle/fraggle_util.f90 @@ -12,28 +12,6 @@ use symba contains - module subroutine fraggle_util_get_angular_momentum(self) - !! Author: David A. Minton - !! - !! Calcualtes the current angular momentum of the fragments - implicit none - ! Arguments - class(fraggle_fragments(*)), intent(inout) :: self !! Fraggle fragment system object - ! Internals - integer(I4B) :: i - - associate(fragments => self, nfrag => self%nbody) - fragments%Lorbit(:) = 0.0_DP - fragments%Lspin(:) = 0.0_DP - - do i = 1, nfrag - fragments%Lorbit(:) = fragments%Lorbit(:) + fragments%mass(i) * (fragments%rc(:, i) .cross. fragments%vc(:, i)) - fragments%Lspin(:) = fragments%Lspin(:) + fragments%mass(i) * fragments%radius(i)**2 * fragments%Ip(:, i) * fragments%rot(:, i) - end do - end associate - - return - end subroutine fraggle_util_get_angular_momentum module subroutine fraggle_util_construct_temporary_system(self, nbody_system, param, tmpsys, tmpparam) @@ -62,6 +40,30 @@ module subroutine fraggle_util_construct_temporary_system(self, nbody_system, pa end subroutine fraggle_util_construct_temporary_system + module subroutine fraggle_util_get_angular_momentum(self) + !! Author: David A. Minton + !! + !! Calcualtes the current angular momentum of the fragments + implicit none + ! Arguments + class(fraggle_fragments(*)), intent(inout) :: self !! Fraggle fragment system object + ! Internals + integer(I4B) :: i + + associate(fragments => self, nfrag => self%nbody) + fragments%Lorbit(:) = 0.0_DP + fragments%Lspin(:) = 0.0_DP + + do i = 1, nfrag + fragments%Lorbit(:) = fragments%Lorbit(:) + fragments%mass(i) * (fragments%rc(:, i) .cross. fragments%vc(:, i)) + fragments%Lspin(:) = fragments%Lspin(:) + fragments%mass(i) * fragments%radius(i)**2 * fragments%Ip(:, i) * fragments%rot(:, i) + end do + end associate + + return + end subroutine fraggle_util_get_angular_momentum + + module subroutine fraggle_util_reset_fragments(self) !! author: David A. Minton !! @@ -159,6 +161,163 @@ module subroutine fraggle_util_restructure(self, impactors, try, f_spin, r_max_s end subroutine fraggle_util_restructure + module subroutine fraggle_util_set_budgets(self) + !! author: David A. Minton + !! + !! Sets the energy and momentum budgets of the fragments based on the collider values and the before/after values of energy and momentum + implicit none + ! Arguments + class(collision_fraggle), intent(inout) :: self !! Fraggle collision system object + ! Internals + real(DP) :: dEtot + real(DP), dimension(NDIM) :: dL + + associate(impactors => self%impactors) + select type(fragments => self%fragments) + class is (fraggle_fragments(*)) + + dEtot = self%Etot(2) - self%Etot(1) + dL(:) = self%Ltot(:,2) - self%Ltot(:,1) + + fragments%L_budget(:) = -dL(:) + fragments%ke_budget = -(dEtot - 0.5_DP * fragments%mtot * dot_product(impactors%vbcom(:), impactors%vbcom(:))) - impactors%Qloss + + end select + end associate + + return + end subroutine fraggle_util_set_budgets + + + module subroutine fraggle_util_set_natural_scale_factors(self) + !! author: David A. Minton + !! + !! Scales dimenional quantities to ~O(1) with respect to the collisional system. + !! This scaling makes it easier for the non-linear minimization to converge on a solution + implicit none + ! Arguments + class(collision_fraggle), intent(inout) :: self !! Fraggle collision system object + ! Internals + integer(I4B) :: i + + associate(collision_merge => self, fragments => self%fragments, impactors => self%impactors) + ! Set scale factors + collision_merge%Escale = 0.5_DP * ( impactors%mass(1) * dot_product(impactors%vb(:,1), impactors%vb(:,1)) & + + impactors%mass(2) * dot_product(impactors%vb(:,2), impactors%vb(:,2))) + collision_merge%dscale = sum(impactors%radius(:)) + collision_merge%mscale = fragments%mtot + collision_merge%vscale = sqrt(collision_merge%Escale / collision_merge%mscale) + collision_merge%tscale = collision_merge%dscale / collision_merge%vscale + collision_merge%Lscale = collision_merge%mscale * collision_merge%dscale * collision_merge%vscale + + ! Scale all dimensioned quantities of impactors and fragments + impactors%rbcom(:) = impactors%rbcom(:) / collision_merge%dscale + impactors%vbcom(:) = impactors%vbcom(:) / collision_merge%vscale + impactors%rbimp(:) = impactors%rbimp(:) / collision_merge%dscale + impactors%rb(:,:) = impactors%rb(:,:) / collision_merge%dscale + impactors%vb(:,:) = impactors%vb(:,:) / collision_merge%vscale + impactors%mass(:) = impactors%mass(:) / collision_merge%mscale + impactors%radius(:) = impactors%radius(:) / collision_merge%dscale + impactors%Lspin(:,:) = impactors%Lspin(:,:) / collision_merge%Lscale + impactors%Lorbit(:,:) = impactors%Lorbit(:,:) / collision_merge%Lscale + + do i = 1, 2 + impactors%rot(:,i) = impactors%Lspin(:,i) / (impactors%mass(i) * impactors%radius(i)**2 * impactors%Ip(3, i)) + end do + + fragments%mtot = fragments%mtot / collision_merge%mscale + fragments%mass = fragments%mass / collision_merge%mscale + fragments%radius = fragments%radius / collision_merge%dscale + impactors%Qloss = impactors%Qloss / collision_merge%Escale + end associate + + return + end subroutine fraggle_util_set_natural_scale_factors + + + module subroutine fraggle_util_set_original_scale_factors(self) + !! author: David A. Minton + !! + !! Restores dimenional quantities back to the system units + use, intrinsic :: ieee_exceptions + implicit none + ! Arguments + class(collision_fraggle), intent(inout) :: self !! Fraggle fragment system object + ! Internals + integer(I4B) :: i + logical, dimension(size(IEEE_ALL)) :: fpe_halting_modes + + call ieee_get_halting_mode(IEEE_ALL,fpe_halting_modes) ! Save the current halting modes so we can turn them off temporarily + call ieee_set_halting_mode(IEEE_ALL,.false.) + + associate(collision_merge => self, fragments => self%fragments, impactors => self%impactors) + + ! Restore scale factors + impactors%rbcom(:) = impactors%rbcom(:) * collision_merge%dscale + impactors%vbcom(:) = impactors%vbcom(:) * collision_merge%vscale + impactors%rbimp(:) = impactors%rbimp(:) * collision_merge%dscale + + impactors%mass = impactors%mass * collision_merge%mscale + impactors%radius = impactors%radius * collision_merge%dscale + impactors%rb = impactors%rb * collision_merge%dscale + impactors%vb = impactors%vb * collision_merge%vscale + impactors%Lspin = impactors%Lspin * collision_merge%Lscale + do i = 1, 2 + impactors%rot(:,i) = impactors%Lspin(:,i) * (impactors%mass(i) * impactors%radius(i)**2 * impactors%Ip(3, i)) + end do + + fragments%mtot = fragments%mtot * collision_merge%mscale + fragments%mass = fragments%mass * collision_merge%mscale + fragments%radius = fragments%radius * collision_merge%dscale + fragments%rot = fragments%rot / collision_merge%tscale + fragments%rc = fragments%rc * collision_merge%dscale + fragments%vc = fragments%vc * collision_merge%vscale + + do i = 1, fragments%nbody + fragments%rb(:, i) = fragments%rc(:, i) + impactors%rbcom(:) + fragments%vb(:, i) = fragments%vc(:, i) + impactors%vbcom(:) + end do + + impactors%Qloss = impactors%Qloss * collision_merge%Escale + + collision_merge%Lorbit(:,:) = collision_merge%Lorbit(:,:) * collision_merge%Lscale + collision_merge%Lspin(:,:) = collision_merge%Lspin(:,:) * collision_merge%Lscale + collision_merge%Ltot(:,:) = collision_merge%Ltot(:,:) * collision_merge%Lscale + collision_merge%ke_orbit(:) = collision_merge%ke_orbit(:) * collision_merge%Escale + collision_merge%ke_spin(:) = collision_merge%ke_spin(:) * collision_merge%Escale + collision_merge%pe(:) = collision_merge%pe(:) * collision_merge%Escale + collision_merge%Etot(:) = collision_merge%Etot(:) * collision_merge%Escale + + collision_merge%mscale = 1.0_DP + collision_merge%dscale = 1.0_DP + collision_merge%vscale = 1.0_DP + collision_merge%tscale = 1.0_DP + collision_merge%Lscale = 1.0_DP + collision_merge%Escale = 1.0_DP + end associate + call ieee_set_halting_mode(IEEE_ALL,fpe_halting_modes) + + return + end subroutine fraggle_util_set_original_scale_factors + + + module subroutine fraggle_util_setup_fragments_system(self, nfrag) + !! author: David A. Minton + !! + !! Initializer for the fragments of the collision system. + implicit none + ! Arguments + class(collision_fraggle), intent(inout) :: self !! Encounter collision system object + integer(I4B), intent(in) :: nfrag !! Number of fragments to create + + if (allocated(self%fragments)) deallocate(self%fragments) + allocate(fraggle_fragments(nbody=nfrag) :: self%fragments) + self%fragments%nbody = nfrag + + return + end subroutine fraggle_util_setup_fragments_system + + module subroutine fraggle_util_shift_vector_to_origin(m_frag, vec_frag) !! Author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton !! diff --git a/src/helio/helio_module.f90 b/src/helio/helio_module.f90 index b8d9d8b26..a9d41b1a9 100644 --- a/src/helio/helio_module.f90 +++ b/src/helio/helio_module.f90 @@ -21,7 +21,7 @@ module helio type, extends(whm_nbody_system) :: helio_nbody_system contains procedure :: step => helio_step_system !! Advance the Helio nbody system forward in time by one step - procedure :: initialize => helio_setup_initialize_system !! Performs Helio-specific initilization steps, including converting to DH coordinates + procedure :: initialize => helio_util_setup_initialize_system !! Performs Helio-specific initilization steps, including converting to DH coordinates final :: helio_final_system !! Finalizes the Helio nbody_system object - deallocates all allocatables end type helio_nbody_system @@ -168,11 +168,11 @@ module subroutine helio_kick_vb_tp(self, nbody_system, param, t, dt, lbeg) logical, intent(in) :: lbeg !! Logical flag indicating whether this is the beginning of the half step or not. end subroutine helio_kick_vb_tp - module subroutine helio_setup_initialize_system(self, param) + module subroutine helio_util_setup_initialize_system(self, param) implicit none class(helio_nbody_system), intent(inout) :: self !! Helio nbody system object class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters - end subroutine helio_setup_initialize_system + end subroutine helio_util_setup_initialize_system module subroutine helio_step_pl(self, nbody_system, param, t, dt) implicit none diff --git a/src/helio/helio_setup.f90 b/src/helio/helio_util.f90 similarity index 72% rename from src/helio/helio_setup.f90 rename to src/helio/helio_util.f90 index 80effc1a9..5de0564a0 100644 --- a/src/helio/helio_setup.f90 +++ b/src/helio/helio_util.f90 @@ -1,17 +1,17 @@ !! Copyright 2022 - David Minton, Carlisle Wishard, Jennifer Pouplin, Jake Elliott, & Dana Singh !! This file is part of Swiftest. -!! Swiftest is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License +!! Swiftest is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License !! as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version. -!! Swiftest is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty +!! Swiftest is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty !! of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. -!! You should have received a copy of the GNU General Public License along with Swiftest. -!! If not, see: https://www.gnu.org/licenses. +!! You should have received a copy of the GNU General Public License along with Swiftest. +!! If not, see: https://www.gnu.org/licenses. -submodule(helio) s_helio_setup +submodule(helio) s_helio_util use swiftest contains - module subroutine helio_setup_initialize_system(self, param) + module subroutine helio_util_setup_initialize_system(self, param) !! author: David A. Minton !! !! Initialize a Helio nbody system from files, converting all heliocentric quantities to barycentric. @@ -21,13 +21,13 @@ module subroutine helio_setup_initialize_system(self, param) class(helio_nbody_system), intent(inout) :: self !! Helio nbody system object class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters - call whm_setup_initialize_system(self, param) + call whm_util_setup_initialize_system(self, param) call self%pl%h2b(self%cb) call self%tp%h2b(self%cb) call self%pl%sort("mass", ascending=.false.) call self%pl%flatten(param) return - end subroutine helio_setup_initialize_system + end subroutine helio_util_setup_initialize_system -end submodule s_helio_setup \ No newline at end of file +end submodule s_helio_util diff --git a/src/rmvs/rmvs_module.f90 b/src/rmvs/rmvs_module.f90 index 51e824b8d..3417ea0e2 100644 --- a/src/rmvs/rmvs_module.f90 +++ b/src/rmvs/rmvs_module.f90 @@ -32,7 +32,7 @@ module rmvs real(DP), dimension(:,:), allocatable :: vbeg !! Planet velocities at beginning ot step contains !> Replace the abstract procedures with concrete ones - procedure :: initialize => rmvs_setup_initialize_system !! Performs RMVS-specific initilization steps, including generating the close encounter planetocentric structures + procedure :: initialize => rmvs_util_setup_initialize_system !! Performs RMVS-specific initilization steps, including generating the close encounter planetocentric structures procedure :: step => rmvs_step_system !! Advance the RMVS nbody system forward in time by one step final :: rmvs_final_system !! Finalizes the RMVS nbody system object - deallocates all allocatables end type rmvs_nbody_system @@ -62,7 +62,7 @@ module rmvs !! RMVS test particle class type, extends(whm_tp) :: rmvs_tp !! 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_tp and rmvs_util_spill_tp + !! component list, such as rmvs_util_setup_tp and rmvs_util_spill_tp ! encounter steps) logical, dimension(:), allocatable :: lperi !! planetocentric pericenter passage flag (persistent for a full rmvs time step) over a full RMVS time step) integer(I4B), dimension(:), allocatable :: plperP !! index of planet associated with pericenter distance peri (persistent over a full RMVS time step) @@ -79,7 +79,7 @@ module rmvs 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 the input number of bodiess + procedure :: setup => rmvs_util_setup_tp !! Constructor method - Allocates space for the input number of bodiess procedure :: append => rmvs_util_append_tp !! Appends elements from one structure to another procedure :: dealloc => rmvs_util_dealloc_tp !! Deallocates all allocatable arrays procedure :: fill => rmvs_util_fill_tp !! "Fills" bodies from one object into another depending on the results of a mask (uses the UNPACK intrinsic) @@ -101,7 +101,7 @@ module rmvs 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 :: setup => rmvs_setup_pl !! Constructor method - Allocates space for the input number of bodiess + procedure :: setup => rmvs_util_setup_pl !! Constructor method - Allocates space for the input number of bodiess procedure :: append => rmvs_util_append_pl !! Appends elements from one structure to another procedure :: dealloc => rmvs_util_dealloc_pl !! Deallocates all allocatable arrays procedure :: fill => rmvs_util_fill_pl !! "Fills" bodies from one object into another depending on the results of a mask (uses the UNPACK intrinsic) @@ -138,25 +138,25 @@ module subroutine rmvs_kick_getacch_tp(self, nbody_system, param, t, lbeg) logical, intent(in) :: lbeg !! Logical flag that determines whether or not this is the beginning or end of the step end subroutine rmvs_kick_getacch_tp - module subroutine rmvs_setup_pl(self, n, param) + module subroutine rmvs_util_setup_pl(self, n, param) implicit none class(rmvs_pl), intent(inout) :: self !! RMVS massive body object integer(I4B), intent(in) :: n !! Number of particles to allocate space for class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters - end subroutine rmvs_setup_pl + end subroutine rmvs_util_setup_pl - module subroutine rmvs_setup_initialize_system(self, param) + module subroutine rmvs_util_setup_initialize_system(self, param) implicit none class(rmvs_nbody_system), intent(inout) :: self !! RMVS system object class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters - end subroutine rmvs_setup_initialize_system + end subroutine rmvs_util_setup_initialize_system - module subroutine rmvs_setup_tp(self, n, param) + module subroutine rmvs_util_setup_tp(self, n, param) implicit none class(rmvs_tp), intent(inout) :: self !! RMVS test particle object integer(I4B), intent(in) :: n !! Number of particles to allocate space for class(swiftest_parameters), intent(in) :: param !! Current run configuration parametere - end subroutine rmvs_setup_tp + end subroutine rmvs_util_setup_tp module subroutine rmvs_util_append_pl(self, source, lsource_mask) implicit none @@ -206,31 +206,6 @@ module subroutine rmvs_util_fill_tp(self, inserts, lfill_list) logical, dimension(:), intent(in) :: lfill_list !! Logical array of bodies to merge into the keeps end subroutine rmvs_util_fill_tp - module subroutine rmvs_final_cb(self) - implicit none - type(rmvs_cb), intent(inout) :: self !! RMVS central body object - end subroutine rmvs_final_cb - - module subroutine rmvs_final_interp(self) - implicit none - type(rmvs_interp), intent(inout) :: self !! RMVS interpolated nbody_system variables object - end subroutine rmvs_final_interp - - module subroutine rmvs_final_pl(self) - implicit none - type(rmvs_pl), intent(inout) :: self !! RMVS massive body object - end subroutine rmvs_final_pl - - module subroutine rmvs_final_system(self) - implicit none - type(rmvs_nbody_system), intent(inout) :: self !! RMVS nbody system object - end subroutine rmvs_final_system - - module subroutine rmvs_final_tp(self) - implicit none - type(rmvs_tp), intent(inout) :: self !! RMVS test particle object - end subroutine rmvs_final_tp - module subroutine rmvs_util_resize_pl(self, nnew) implicit none class(rmvs_pl), intent(inout) :: self !! RMVS massive body object @@ -295,4 +270,77 @@ end subroutine rmvs_step_system end interface + contains + + subroutine rmvs_final_cb(self) + !! author: David A. Minton + !! + !! Finalize the RMVS massive body object - deallocates all allocatables + implicit none + ! Arguments + type(rmvs_cb), intent(inout) :: self !! RMVS central body object + + call self%dealloc() + + return + end subroutine rmvs_final_cb + + + subroutine rmvs_final_interp(self) + !! author: David A. Minton + !! + !! Finalize the RMVS nbody system object - deallocates all allocatables + implicit none + ! Arguments + type(rmvs_interp), intent(inout) :: self !! RMVS nbody system object + + call self%dealloc() + + return + end subroutine rmvs_final_interp + + + subroutine rmvs_final_pl(self) + !! author: David A. Minton + !! + !! Finalize the RMVS massive body object - deallocates all allocatables + implicit none + ! Arguments + type(rmvs_pl), intent(inout) :: self !! RMVS massive body object + + call self%dealloc() + + return + end subroutine rmvs_final_pl + + + subroutine rmvs_final_system(self) + !! author: David A. Minton + !! + !! Finalize the RMVS nbody system object - deallocates all allocatables + implicit none + ! Arguments + type(rmvs_nbody_system), intent(inout) :: self !! RMVS nbody system object + + if (allocated(self%vbeg)) deallocate(self%vbeg) + call whm_final_system(self%whm_nbody_system) + + return + end subroutine rmvs_final_system + + + subroutine rmvs_final_tp(self) + !! author: David A. Minton + !! + !! Finalize the RMVS test particle object - deallocates all allocatables + implicit none + ! Arguments + type(rmvs_tp), intent(inout) :: self !! RMVS test particle object + + call self%dealloc() + + return + end subroutine rmvs_final_tp + + end module rmvs diff --git a/src/rmvs/rmvs_setup.f90 b/src/rmvs/rmvs_setup.f90 deleted file mode 100644 index fe7d15b6d..000000000 --- a/src/rmvs/rmvs_setup.f90 +++ /dev/null @@ -1,166 +0,0 @@ -!! Copyright 2022 - David Minton, Carlisle Wishard, Jennifer Pouplin, Jake Elliott, & Dana Singh -!! This file is part of Swiftest. -!! Swiftest is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License -!! as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version. -!! Swiftest is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty -!! of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. -!! You should have received a copy of the GNU General Public License along with Swiftest. -!! If not, see: https://www.gnu.org/licenses. - -submodule(rmvs) s_rmvs_setup - use swiftest -contains - - module subroutine rmvs_setup_pl(self, n, param) - !! author: David A. Minton - !! - !! Allocate RMVS test particle structure - !! - !! 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 - integer(I4B), intent(in) :: n !! Number of particles to allocate space for - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameter - ! Internals - integer(I4B) :: i - - !> Call allocation method for parent class - associate(pl => self) - call whm_setup_pl(pl, n, param) - if (n == 0) return - - allocate(pl%outer(0:NTENC)) - allocate(pl%inner(0:NTPHENC)) - if (.not.pl%lplanetocentric) then - allocate(pl%nenc(n)) - pl%nenc(:) = 0 - ! Set up inner and outer planet interpolation vector storage containers - do i = 0, NTENC - allocate(pl%outer(i)%x(NDIM, n)) - allocate(pl%outer(i)%v(NDIM, n)) - pl%outer(i)%x(:,:) = 0.0_DP - pl%outer(i)%v(:,:) = 0.0_DP - end do - do i = 0, NTPHENC - allocate(pl%inner(i)%x(NDIM, n)) - allocate(pl%inner(i)%v(NDIM, n)) - allocate(pl%inner(i)%aobl(NDIM, n)) - pl%inner(i)%x(:,:) = 0.0_DP - pl%inner(i)%v(:,:) = 0.0_DP - pl%inner(i)%aobl(:,:) = 0.0_DP - end do - ! if (param%ltides) then - ! do i = 0, NTPHENC - ! allocate(pl%inner(i)%atide(NDIM, n)) - ! pl%inner(i)%atide(:,:) = 0.0_DP - ! end do - ! end if - end if - end associate - return - end subroutine rmvs_setup_pl - - - module subroutine rmvs_setup_initialize_system(self, param) - !! author: David A. Minton - !! - !! Initialize an RMVS nbody system from files and sets up the planetocentric structures. - !! - !! We currently rearrange the pl order to keep it consistent with the way Swifter does it - !! In Swifter, the central body occupies the first position in the pl list, and during - !! encounters, the encountering planet is skipped in loops. In Swiftest, we instantiate an - !! RMVS nbody system object attached to each pl to store planetocentric versions of the nbody_system - !! to use during close encounters. - implicit none - ! Arguments - class(rmvs_nbody_system), intent(inout) :: self !! RMVS system object - class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters - ! Internals - integer(I4B) :: i, j - - ! Call parent method - call whm_setup_initialize_system(self, param) - - ! Set up the pl-tp planetocentric encounter structures for pl and cb. The planetocentric tp structures are - ! generated as necessary during close encounter steps. - 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) - tp%cb_heliocentric = cb - pl%lplanetocentric = .false. - tp%lplanetocentric = .false. - cb%lplanetocentric = .false. - associate(npl => pl%nbody) - allocate(pl%planetocentric(npl)) - pl%planetocentric(:)%lplanetocentric = .true. - do i = 1, npl - allocate(pl%planetocentric(i)%cb, source=cb) - allocate(rmvs_pl :: pl%planetocentric(i)%pl) - select type(cbenci => pl%planetocentric(i)%cb) - class is (rmvs_cb) - select type(plenci => pl%planetocentric(i)%pl) - class is (rmvs_pl) - cbenci%lplanetocentric = .true. - plenci%lplanetocentric = .true. - call plenci%setup(npl, param) - plenci%status(:) = ACTIVE - plenci%lmask(:) = .true. - ! plind stores the heliocentric index value of a planetocentric planet - ! e.g. Consider an encounter with planet 3. - ! Then the following will be the values of plind: - ! pl%planetocentric(3)%pl%plind(1) = 0 (central body - never used) - ! pl%planetocentric(3)%pl%plind(2) = 1 - ! pl%planetocentric(3)%pl%plind(3) = 2 - ! pl%planetocentric(3)%pl%plind(4) = 4 - ! pl%planetocentric(3)%pl%plind(5) = 5 - ! etc. - allocate(plenci%plind(npl)) - plenci%plind(1:npl) = [(j,j=1,npl)] - plenci%plind(2:npl) = pack(plenci%plind(1:npl), plenci%plind(1:npl) /= i) - plenci%plind(1) = 0 - plenci%Gmass(1) = cb%Gmass - plenci%Gmass(2:npl) = pl%Gmass(plenci%plind(2:npl)) - cbenci%Gmass = pl%Gmass(i) - end select - end select - end do - end associate - end select - end select - end select - return - end subroutine rmvs_setup_initialize_system - - - module subroutine rmvs_setup_tp(self, n, param) - !! author: David A. Minton - !! - !! Allocate WHM test particle structure - !! - !! Equivalent in functionality to David E. Kaufmann's Swifter routine whm_setup.f90 - implicit none - ! Arguments - class(rmvs_tp), intent(inout) :: self !! RMVS test particle object - integer(I4B), intent(in) :: n !! Number of particles to allocate space for - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameter - - !> Call allocation method for parent class. - call self%whm_tp%setup(n, param) - if (n <= 0) return - - allocate(self%lperi(n)) - allocate(self%plperP(n)) - allocate(self%plencP(n)) - - if (self%lplanetocentric) allocate(self%rheliocentric(NDIM, n)) - - self%lperi(:) = .false. - - return - end subroutine rmvs_setup_tp - -end submodule s_rmvs_setup diff --git a/src/rmvs/rmvs_util.f90 b/src/rmvs/rmvs_util.f90 index 26306f064..07e4f9a51 100644 --- a/src/rmvs/rmvs_util.f90 +++ b/src/rmvs/rmvs_util.f90 @@ -181,78 +181,7 @@ module subroutine rmvs_util_fill_pl(self, inserts, lfill_list) return end subroutine rmvs_util_fill_pl - - module subroutine rmvs_final_cb(self) - !! author: David A. Minton - !! - !! Finalize the RMVS massive body object - deallocates all allocatables - implicit none - ! Arguments - type(rmvs_cb), intent(inout) :: self !! RMVS central body object - - call self%dealloc() - - return - end subroutine rmvs_final_cb - - - module subroutine rmvs_final_interp(self) - !! author: David A. Minton - !! - !! Finalize the RMVS nbody system object - deallocates all allocatables - implicit none - ! Arguments - type(rmvs_interp), intent(inout) :: self !! RMVS nbody system object - - call self%dealloc() - - return - end subroutine rmvs_final_interp - - - module subroutine rmvs_final_pl(self) - !! author: David A. Minton - !! - !! Finalize the RMVS massive body object - deallocates all allocatables - implicit none - ! Arguments - type(rmvs_pl), intent(inout) :: self !! RMVS massive body object - - call self%dealloc() - - return - end subroutine rmvs_final_pl - - - module subroutine rmvs_final_system(self) - !! author: David A. Minton - !! - !! Finalize the RMVS nbody system object - deallocates all allocatables - implicit none - ! Arguments - type(rmvs_nbody_system), intent(inout) :: self !! RMVS nbody system object - - if (allocated(self%vbeg)) deallocate(self%vbeg) - call whm_final_system(self%whm_nbody_system) - - return - end subroutine rmvs_final_system - - - module subroutine rmvs_final_tp(self) - !! author: David A. Minton - !! - !! Finalize the RMVS test particle object - deallocates all allocatables - implicit none - ! Arguments - type(rmvs_tp), intent(inout) :: self !! RMVS test particle object - - call self%dealloc() - - return - end subroutine rmvs_final_tp - - + module subroutine rmvs_util_fill_tp(self, inserts, lfill_list) !! author: David A. Minton !! @@ -327,6 +256,159 @@ module subroutine rmvs_util_resize_tp(self, nnew) end subroutine rmvs_util_resize_tp + module subroutine rmvs_util_setup_pl(self, n, param) + !! author: David A. Minton + !! + !! Allocate RMVS test particle structure + !! + !! Equivalent in functionality to David E. Kaufmann's Swifter routine rmvs_util_setup.f90 + implicit none + ! Arguments + class(rmvs_pl), intent(inout) :: self !! RMVS test particle object + integer(I4B), intent(in) :: n !! Number of particles to allocate space for + class(swiftest_parameters), intent(in) :: param !! Current run configuration parameter + ! Internals + integer(I4B) :: i + + !> Call allocation method for parent class + associate(pl => self) + call whm_util_setup_pl(pl, n, param) + if (n == 0) return + + allocate(pl%outer(0:NTENC)) + allocate(pl%inner(0:NTPHENC)) + if (.not.pl%lplanetocentric) then + allocate(pl%nenc(n)) + pl%nenc(:) = 0 + ! Set up inner and outer planet interpolation vector storage containers + do i = 0, NTENC + allocate(pl%outer(i)%x(NDIM, n)) + allocate(pl%outer(i)%v(NDIM, n)) + pl%outer(i)%x(:,:) = 0.0_DP + pl%outer(i)%v(:,:) = 0.0_DP + end do + do i = 0, NTPHENC + allocate(pl%inner(i)%x(NDIM, n)) + allocate(pl%inner(i)%v(NDIM, n)) + allocate(pl%inner(i)%aobl(NDIM, n)) + pl%inner(i)%x(:,:) = 0.0_DP + pl%inner(i)%v(:,:) = 0.0_DP + pl%inner(i)%aobl(:,:) = 0.0_DP + end do + ! if (param%ltides) then + ! do i = 0, NTPHENC + ! allocate(pl%inner(i)%atide(NDIM, n)) + ! pl%inner(i)%atide(:,:) = 0.0_DP + ! end do + ! end if + end if + end associate + return + end subroutine rmvs_util_setup_pl + + + module subroutine rmvs_util_setup_initialize_system(self, param) + !! author: David A. Minton + !! + !! Initialize an RMVS nbody system from files and sets up the planetocentric structures. + !! + !! We currently rearrange the pl order to keep it consistent with the way Swifter does it + !! In Swifter, the central body occupies the first position in the pl list, and during + !! encounters, the encountering planet is skipped in loops. In Swiftest, we instantiate an + !! RMVS nbody system object attached to each pl to store planetocentric versions of the nbody_system + !! to use during close encounters. + implicit none + ! Arguments + class(rmvs_nbody_system), intent(inout) :: self !! RMVS system object + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters + ! Internals + integer(I4B) :: i, j + + ! Call parent method + call whm_util_setup_initialize_system(self, param) + + ! Set up the pl-tp planetocentric encounter structures for pl and cb. The planetocentric tp structures are + ! generated as necessary during close encounter steps. + 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) + tp%cb_heliocentric = cb + pl%lplanetocentric = .false. + tp%lplanetocentric = .false. + cb%lplanetocentric = .false. + associate(npl => pl%nbody) + allocate(pl%planetocentric(npl)) + pl%planetocentric(:)%lplanetocentric = .true. + do i = 1, npl + allocate(pl%planetocentric(i)%cb, source=cb) + allocate(rmvs_pl :: pl%planetocentric(i)%pl) + select type(cbenci => pl%planetocentric(i)%cb) + class is (rmvs_cb) + select type(plenci => pl%planetocentric(i)%pl) + class is (rmvs_pl) + cbenci%lplanetocentric = .true. + plenci%lplanetocentric = .true. + call plenci%setup(npl, param) + plenci%status(:) = ACTIVE + plenci%lmask(:) = .true. + ! plind stores the heliocentric index value of a planetocentric planet + ! e.g. Consider an encounter with planet 3. + ! Then the following will be the values of plind: + ! pl%planetocentric(3)%pl%plind(1) = 0 (central body - never used) + ! pl%planetocentric(3)%pl%plind(2) = 1 + ! pl%planetocentric(3)%pl%plind(3) = 2 + ! pl%planetocentric(3)%pl%plind(4) = 4 + ! pl%planetocentric(3)%pl%plind(5) = 5 + ! etc. + allocate(plenci%plind(npl)) + plenci%plind(1:npl) = [(j,j=1,npl)] + plenci%plind(2:npl) = pack(plenci%plind(1:npl), plenci%plind(1:npl) /= i) + plenci%plind(1) = 0 + plenci%Gmass(1) = cb%Gmass + plenci%Gmass(2:npl) = pl%Gmass(plenci%plind(2:npl)) + cbenci%Gmass = pl%Gmass(i) + end select + end select + end do + end associate + end select + end select + end select + return + end subroutine rmvs_util_setup_initialize_system + + + module subroutine rmvs_util_setup_tp(self, n, param) + !! author: David A. Minton + !! + !! Allocate WHM test particle structure + !! + !! Equivalent in functionality to David E. Kaufmann's Swifter routine whm_util_setup.f90 + implicit none + ! Arguments + class(rmvs_tp), intent(inout) :: self !! RMVS test particle object + integer(I4B), intent(in) :: n !! Number of particles to allocate space for + class(swiftest_parameters), intent(in) :: param !! Current run configuration parameter + + !> Call allocation method for parent class. + call self%whm_tp%setup(n, param) + if (n <= 0) return + + allocate(self%lperi(n)) + allocate(self%plperP(n)) + allocate(self%plencP(n)) + + if (self%lplanetocentric) allocate(self%rheliocentric(NDIM, n)) + + self%lperi(:) = .false. + + return + end subroutine rmvs_util_setup_tp + + module subroutine rmvs_util_sort_pl(self, sortby, ascending) !! author: David A. Minton !! diff --git a/src/swiftest/swiftest_driver.f90 b/src/swiftest/swiftest_driver.f90 index 4c0b18ce1..e2ee1c054 100644 --- a/src/swiftest/swiftest_driver.f90 +++ b/src/swiftest/swiftest_driver.f90 @@ -71,7 +71,7 @@ program swiftest_driver if (dump_cadence == 0) dump_cadence = ceiling(nloops / (1.0_DP * istep_out), kind=I8B) ! Construct the main n-body nbody_system using the user-input integrator to choose the type of nbody_system - call swiftest_setup_construct_system(nbody_system, param) + call swiftest_util_setup_construct_system(nbody_system, param) !> Define the maximum number of threads nthreads = 1 ! In the *serial* case diff --git a/src/swiftest/swiftest_module.f90 b/src/swiftest/swiftest_module.f90 index ff94cb85c..d09f88374 100644 --- a/src/swiftest/swiftest_module.f90 +++ b/src/swiftest/swiftest_module.f90 @@ -144,7 +144,7 @@ module swiftest procedure :: accel_obl => swiftest_obl_acc_body !! Compute the barycentric accelerations of bodies due to the oblateness of the central body procedure :: el2xv => swiftest_orbel_el2xv_vec !! Convert orbital elements to position and velocity vectors procedure :: xv2el => swiftest_orbel_xv2el_vec !! Convert position and velocity vectors to orbital elements - procedure :: setup => swiftest_setup_body !! A constructor that sets the number of bodies and allocates all allocatable arrays + procedure :: setup => swiftest_util_setup_body !! A constructor that sets the number of bodies and allocates all allocatable arrays procedure :: accel_user => swiftest_user_kick_getacch_body !! Add user-supplied heliocentric accelerations to planets procedure :: append => swiftest_util_append_body !! Appends elements from one structure to another procedure :: dealloc => swiftest_util_dealloc_body !! Deallocates all allocatable arrays @@ -249,7 +249,7 @@ module swiftest procedure :: discard => swiftest_discard_pl !! Placeholder method for discarding massive bodies procedure :: accel_int => swiftest_kick_getacch_int_pl !! Compute direct cross (third) term heliocentric accelerations of massive bodies procedure :: accel_obl => swiftest_obl_acc_pl !! Compute the barycentric accelerations of bodies due to the oblateness of the central body - procedure :: setup => swiftest_setup_pl !! A base constructor that sets the number of bodies and allocates and initializes all arrays + procedure :: setup => swiftest_util_setup_pl !! A base constructor that sets the number of bodies and allocates and initializes all arrays ! procedure :: accel_tides => tides_kick_getacch_pl !! Compute the accelerations of bodies due to tidal interactions with the central body procedure :: append => swiftest_util_append_pl !! Appends elements from one structure to another procedure :: h2b => swiftest_util_coord_h2b_pl !! Convert massive bodies from heliocentric to barycentric coordinates (position and velocity) @@ -281,14 +281,14 @@ module swiftest integer(I8B) :: npltp !! Number of pl-tp comparisons in the flattened upper triangular matrix integer(I4B), dimension(:), allocatable :: nplenc !! number of encounters with planets this time step !! Note to developers: If you add components to this class, be sure to update methods and subroutines that traverse the - !! component list, such as swiftest_setup_tp and util_spill_tp + !! component list, such as swiftest_util_setup_tp and util_spill_tp contains ! Test particle-specific concrete methods ! These are concrete because they are the same implemenation for all integrators procedure :: discard => swiftest_discard_tp !! Check to see if test particles should be discarded based on their positions relative to the massive bodies procedure :: accel_int => swiftest_kick_getacch_int_tp !! Compute direct cross (third) term heliocentric accelerations of test particles by massive bodies procedure :: accel_obl => swiftest_obl_acc_tp !! Compute the barycentric accelerations of bodies due to the oblateness of the central body - procedure :: setup => swiftest_setup_tp !! A base constructor that sets the number of bodies and + procedure :: setup => swiftest_util_setup_tp !! A base constructor that sets the number of bodies and procedure :: append => swiftest_util_append_tp !! Appends elements from one structure to another procedure :: h2b => swiftest_util_coord_h2b_tp !! Convert test particles from heliocentric to barycentric coordinates (position and velocity) procedure :: b2h => swiftest_util_coord_b2h_tp !! Convert test particles from barycentric to heliocentric coordinates (position and velocity) @@ -385,8 +385,8 @@ module swiftest procedure :: read_in => swiftest_io_read_in_system !! Reads the initial conditions for an nbody system procedure :: read_particle_info => swiftest_io_netcdf_read_particle_info_system !! Read in particle metadata from file procedure :: obl_pot => swiftest_obl_pot_system !! Compute the contribution to the total gravitational potential due solely to the oblateness of the central body - procedure :: initialize => swiftest_setup_initialize_system !! Initialize the nbody_system from input files - procedure :: init_particle_info => swiftest_setup_initialize_particle_info_system !! Initialize the nbody_system from input files + procedure :: initialize => swiftest_util_setup_initialize_system !! Initialize the nbody_system from input files + procedure :: init_particle_info => swiftest_util_setup_initialize_particle_info_system !! Initialize the nbody_system from input files ! procedure :: step_spin => tides_step_spin_system !! Steps the spins of the massive & central bodies due to tides. procedure :: set_msys => swiftest_util_set_msys !! Sets the value of msys from the masses of nbody_system bodies. procedure :: get_energy_and_momentum => swiftest_util_get_energy_momentum_system !! Calculates the total nbody_system energy and momentum @@ -1007,44 +1007,44 @@ module subroutine swiftest_orbel_xv2el_vec(self, cb) class(swiftest_cb), intent(inout) :: cb !! Swiftest central body object end subroutine swiftest_orbel_xv2el_vec - module subroutine swiftest_setup_body(self, n, param) + module subroutine swiftest_util_setup_body(self, n, param) implicit none class(swiftest_body), intent(inout) :: self !! Swiftest body object integer(I4B), intent(in) :: n !! Number of particles to allocate space for class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters - end subroutine swiftest_setup_body + end subroutine swiftest_util_setup_body - module subroutine swiftest_setup_construct_system(nbody_system, param) + module subroutine swiftest_util_setup_construct_system(nbody_system, param) implicit none class(swiftest_nbody_system), allocatable, intent(inout) :: nbody_system !! Swiftest nbody_system object class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters - end subroutine swiftest_setup_construct_system + end subroutine swiftest_util_setup_construct_system - module subroutine swiftest_setup_initialize_particle_info_system(self, param) + module subroutine swiftest_util_setup_initialize_particle_info_system(self, param) implicit none class(swiftest_nbody_system), intent(inout) :: self !! Swiftest nbody system object class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters - end subroutine swiftest_setup_initialize_particle_info_system + end subroutine swiftest_util_setup_initialize_particle_info_system - module subroutine swiftest_setup_initialize_system(self, param) + module subroutine swiftest_util_setup_initialize_system(self, param) implicit none class(swiftest_nbody_system), intent(inout) :: self !! Swiftest nbody_system object class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters - end subroutine swiftest_setup_initialize_system + end subroutine swiftest_util_setup_initialize_system - module subroutine swiftest_setup_pl(self, n, param) + module subroutine swiftest_util_setup_pl(self, n, param) implicit none class(swiftest_pl), intent(inout) :: self !! Swiftest massive body object integer(I4B), intent(in) :: n !! Number of particles to allocate space for class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters - end subroutine swiftest_setup_pl + end subroutine swiftest_util_setup_pl - module subroutine swiftest_setup_tp(self, n, param) + module subroutine swiftest_util_setup_tp(self, n, param) implicit none class(swiftest_tp), intent(inout) :: self !! Swiftest test particle object integer(I4B), intent(in) :: n !! Number of particles to allocate space for class(swiftest_parameters), intent(in) :: param !! Current run configuration parametersr - end subroutine swiftest_setup_tp + end subroutine swiftest_util_setup_tp module subroutine swiftest_user_kick_getacch_body(self, nbody_system, param, t, lbeg) implicit none diff --git a/src/swiftest/swiftest_setup.f90 b/src/swiftest/swiftest_setup.f90 deleted file mode 100644 index c693c96ed..000000000 --- a/src/swiftest/swiftest_setup.f90 +++ /dev/null @@ -1,385 +0,0 @@ -!! Copyright 2022 - David Minton, Carlisle Wishard, Jennifer Pouplin, Jake Elliott, & Dana Singh -!! This file is part of Swiftest. -!! Swiftest is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License -!! as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version. -!! Swiftest is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty -!! of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. -!! You should have received a copy of the GNU General Public License along with Swiftest. -!! If not, see: https://www.gnu.org/licenses. - -submodule (swiftest) s_setup - use whm - use rmvs - use helio - use symba - use fraggle -contains - - module subroutine swiftest_setup_construct_system(nbody_system, param) - !! author: David A. Minton - !! - !! Constructor for a Swiftest nbody system. Creates the nbody system object based on the user-input integrator - !! - implicit none - ! Arguments - class(swiftest_nbody_system), allocatable, intent(inout) :: nbody_system !! Swiftest nbody_system object - class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters - ! Internals - type(encounter_storage) :: encounter_history - type(collision_storage) :: collision_history - - select type(param) - class is (swiftest_parameters) - allocate(swiftest_storage(param%dump_cadence) :: param%system_history) - allocate(swiftest_netcdf_parameters :: param%system_history%nc) - call param%system_history%reset() - - select case(param%integrator) - case (INT_BS) - write(*,*) 'Bulirsch-Stoer integrator not yet enabled' - case (INT_HELIO) - allocate(helio_nbody_system :: nbody_system) - select type(nbody_system) - class is (helio_nbody_system) - allocate(helio_cb :: nbody_system%cb) - allocate(helio_pl :: nbody_system%pl) - allocate(helio_tp :: nbody_system%tp) - allocate(helio_tp :: nbody_system%tp_discards) - end select - param%collision_model = "MERGE" - case (INT_RA15) - write(*,*) 'Radau integrator not yet enabled' - case (INT_TU4) - write(*,*) 'INT_TU4 integrator not yet enabled' - case (INT_WHM) - allocate(whm_nbody_system :: nbody_system) - select type(nbody_system) - class is (whm_nbody_system) - allocate(whm_cb :: nbody_system%cb) - allocate(whm_pl :: nbody_system%pl) - allocate(whm_tp :: nbody_system%tp) - allocate(whm_tp :: nbody_system%tp_discards) - end select - param%collision_model = "MERGE" - case (INT_RMVS) - allocate(rmvs_nbody_system :: nbody_system) - select type(nbody_system) - class is (rmvs_nbody_system) - allocate(rmvs_cb :: nbody_system%cb) - allocate(rmvs_pl :: nbody_system%pl) - allocate(rmvs_tp :: nbody_system%tp) - allocate(rmvs_tp :: nbody_system%tp_discards) - end select - param%collision_model = "MERGE" - case (INT_SYMBA) - allocate(symba_nbody_system :: nbody_system) - select type(nbody_system) - class is (symba_nbody_system) - allocate(symba_cb :: nbody_system%cb) - allocate(symba_pl :: nbody_system%pl) - allocate(symba_tp :: nbody_system%tp) - - allocate(symba_tp :: nbody_system%tp_discards) - allocate(symba_pl :: nbody_system%pl_adds) - allocate(symba_pl :: nbody_system%pl_discards) - - allocate(symba_list_pltp :: nbody_system%pltp_encounter) - allocate(symba_list_plpl :: nbody_system%plpl_encounter) - allocate(collision_list_plpl :: nbody_system%plpl_collision) - - if (param%lenc_save_trajectory .or. param%lenc_save_closest) then - allocate(encounter_netcdf_parameters :: encounter_history%nc) - call encounter_history%reset() - select type(nc => encounter_history%nc) - class is (encounter_netcdf_parameters) - nc%file_number = param%iloop / param%dump_cadence - end select - allocate(nbody_system%encounter_history, source=encounter_history) - end if - - allocate(collision_netcdf_parameters :: collision_history%nc) - call collision_history%reset() - select type(nc => collision_history%nc) - class is (collision_netcdf_parameters) - nc%file_number = param%iloop / param%dump_cadence - end select - allocate(nbody_system%collision_history, source=collision_history) - - end select - case (INT_RINGMOONS) - write(*,*) 'RINGMOONS-SyMBA integrator not yet enabled' - case default - write(*,*) 'Unkown integrator',param%integrator - call util_exit(FAILURE) - end select - - allocate(swiftest_particle_info :: nbody_system%cb%info) - - select case(param%collision_model) - case("MERGE") - allocate(collision_merge :: nbody_system%collider) - case("BOUNCE") - allocate(collision_bounce :: nbody_system%collider) - case("SIMPLE") - allocate(collision_simple_disruption :: nbody_system%collider) - case("FRAGGLE") - allocate(collision_fraggle :: nbody_system%collider) - end select - call nbody_system%collider%setup(nbody_system) - - end select - - return - end subroutine swiftest_setup_construct_system - - - module subroutine swiftest_setup_initialize_particle_info_system(self, param) - !! author: David A. Minton - !! - !! Setup up particle information metadata from initial conditions - ! - implicit none - ! Arguments - class(swiftest_nbody_system), intent(inout) :: self !! Swiftest nbody system object - class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters - ! Internals - integer(I4B) :: i - - associate(pl => self%pl, npl => self%pl%nbody, tp => self%tp, ntp => self%tp%nbody) - - if (.not. allocated(self%cb%info)) allocate(swiftest_particle_info :: self%cb%info) - - call self%cb%info%set_value(particle_type=CB_TYPE_NAME, status="ACTIVE", origin_type="Initial conditions", & - origin_time=param%t0, origin_rh=[0.0_DP, 0.0_DP, 0.0_DP], origin_vh=[0.0_DP, 0.0_DP, 0.0_DP]) - do i = 1, self%pl%nbody - call pl%info(i)%set_value(particle_type=PL_TYPE_NAME, status="ACTIVE", origin_type="Initial conditions", & - origin_time=param%t0, origin_rh=self%pl%rh(:,i), origin_vh=self%pl%vh(:,i)) - end do - do i = 1, self%tp%nbody - call tp%info(i)%set_value(particle_type=TP_TYPE_NAME, status="ACTIVE", origin_type="Initial conditions", & - origin_time=param%t0, origin_rh=self%tp%rh(:,i), origin_vh=self%tp%vh(:,i)) - end do - - end associate - - return - end subroutine swiftest_setup_initialize_particle_info_system - - - module subroutine swiftest_setup_initialize_system(self, param) - !! author: David A. Minton - !! - !! Wrapper method to initialize a basic Swiftest nbody system from files - !! - implicit none - ! Arguments - class(swiftest_nbody_system), intent(inout) :: self !! Swiftest nbody_system object - class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters - - associate(nbody_system => self, cb => self%cb, pl => self%pl, tp => self%tp) - - call nbody_system%read_in(param) - call nbody_system%validate_ids(param) - call nbody_system%set_msys() - call pl%set_mu(cb) - call tp%set_mu(cb) - if (param%in_form == "EL") then - call pl%el2xv(cb) - call tp%el2xv(cb) - end if - call pl%flatten(param) - if (.not.param%lrhill_present) call pl%set_rhill(cb) - pl%lfirst = param%lfirstkick - tp%lfirst = param%lfirstkick - - if (.not.param%lrestart) then - call nbody_system%init_particle_info(param) - end if - end associate - - return - end subroutine swiftest_setup_initialize_system - - - module subroutine swiftest_setup_body(self, n, param) - !! author: David A. Minton - !! - !! Constructor for base Swiftest particle class. Allocates space for all particles and - !! initializes all components with a value. - !! Note: Timing tests indicate that (NDIM, n) is more efficient than (NDIM, n) - implicit none - ! Arguments - class(swiftest_body), intent(inout) :: self !! Swiftest generic body object - integer(I4B), intent(in) :: n !! Number of particles to allocate space for - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameter - ! Internals - integer(I4B) :: i - - if (n < 0) return - - self%lfirst = .true. - - call self%dealloc() - - self%nbody = n - if (n == 0) return - - allocate(swiftest_particle_info :: self%info(n)) - allocate(self%id(n)) - allocate(self%status(n)) - allocate(self%ldiscard(n)) - allocate(self%lmask(n)) - allocate(self%mu(n)) - allocate(self%rh(NDIM, n)) - allocate(self%vh(NDIM, n)) - allocate(self%rb(NDIM, n)) - allocate(self%vb(NDIM, n)) - allocate(self%ah(NDIM, n)) - allocate(self%ir3h(n)) - allocate(self%aobl(NDIM, n)) - if (param%lclose) then - allocate(self%lcollision(n)) - allocate(self%lencounter(n)) - self%lcollision(:) = .false. - self%lencounter(:) = .false. - end if - - self%id(:) = 0 - do i = 1, n - call self%info(i)%set_value(& - name = "UNNAMED", & - particle_type = "UNKNOWN", & - status = "INACTIVE", & - origin_type = "UNKNOWN", & - collision_id = 0, & - origin_time = -huge(1.0_DP), & - origin_rh = [0.0_DP, 0.0_DP, 0.0_DP], & - origin_vh = [0.0_DP, 0.0_DP, 0.0_DP], & - discard_time = huge(1.0_DP), & - discard_rh = [0.0_DP, 0.0_DP, 0.0_DP], & - discard_vh = [0.0_DP, 0.0_DP, 0.0_DP], & - discard_body_id = -1 & - ) - end do - - self%status(:) = INACTIVE - self%ldiscard(:) = .false. - self%lmask(:) = .false. - self%mu(:) = 0.0_DP - self%rh(:,:) = 0.0_DP - self%vh(:,:) = 0.0_DP - self%rb(:,:) = 0.0_DP - self%vb(:,:) = 0.0_DP - self%ah(:,:) = 0.0_DP - self%ir3h(:) = 0.0_DP - self%aobl(:,:) = 0.0_DP - - if (param%ltides) then - allocate(self%atide(NDIM, n)) - self%atide(:,:) = 0.0_DP - end if - if (param%lgr) then - allocate(self%agr(NDIM, n)) - self%agr(:,:) = 0.0_DP - end if - - return - end subroutine swiftest_setup_body - - - module subroutine swiftest_setup_pl(self, n, param) - !! author: David A. Minton - !! - !! Constructor for base Swiftest massive body class. Allocates space for all particles and - !! initializes all components with a value. - implicit none - ! Arguments - class(swiftest_pl), intent(inout) :: self !! Swiftest massive body object - integer(I4B), intent(in) :: n !! Number of particles to allocate space for - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameter - - !> Call allocation method for parent class - !> The parent class here is the abstract swiftest_body class, so we can't use the type-bound procedure - call swiftest_setup_body(self, n, param) - if (n == 0) return - - allocate(self%mass(n)) - allocate(self%Gmass(n)) - allocate(self%rhill(n)) - allocate(self%renc(n)) - - self%mass(:) = 0.0_DP - self%Gmass(:) = 0.0_DP - self%rhill(:) = 0.0_DP - self%renc(:) = 0.0_DP - - self%nplpl = 0 - - if (param%lclose) then - allocate(self%nplenc(n)) - allocate(self%ntpenc(n)) - allocate(self%radius(n)) - allocate(self%density(n)) - - self%nplenc(:) = 0 - self%ntpenc(:) = 0 - self%radius(:) = 0.0_DP - self%density(:) = 1.0_DP - - end if - - if (param%lmtiny_pl) then - allocate(self%lmtiny(n)) - self%lmtiny(:) = .false. - end if - - if (param%lrotation) then - allocate(self%rot(NDIM, n)) - allocate(self%Ip(NDIM, n)) - self%rot(:,:) = 0.0_DP - self%Ip(:,:) = 0.0_DP - end if - - if (param%ltides) then - allocate(self%k2(n)) - allocate(self%Q(n)) - allocate(self%tlag(n)) - self%k2(:) = 0.0_DP - self%Q(:) = 0.0_DP - self%tlag(:) = 0.0_DP - end if - - return - end subroutine swiftest_setup_pl - - - module subroutine swiftest_setup_tp(self, n, param) - !! author: David A. Minton - !! - !! Constructor for base Swiftest test particle particle class. Allocates space for - !! all particles and initializes all components with a value. - implicit none - ! Arguments - class(swiftest_tp), intent(inout) :: self !! Swiftest test particle object - integer(I4B), intent(in) :: n !! Number of particles to allocate space for - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameter - - !> Call allocation method for parent class - !> The parent class here is the abstract swiftest_body class, so we can't use the type-bound procedure - call swiftest_setup_body(self, n, param) - if (n == 0) return - - allocate(self%isperi(n)) - allocate(self%peri(n)) - allocate(self%atp(n)) - allocate(self%nplenc(n)) - - self%isperi(:) = 0 - self%peri(:) = 0.0_DP - self%atp(:) = 0.0_DP - self%nplenc(:) = 0 - - return - end subroutine swiftest_setup_tp - -end submodule s_setup diff --git a/src/swiftest/swiftest_util.f90 b/src/swiftest/swiftest_util.f90 index b94dbbd07..56f2259f6 100644 --- a/src/swiftest/swiftest_util.f90 +++ b/src/swiftest/swiftest_util.f90 @@ -8,6 +8,11 @@ !! If not, see: https://www.gnu.org/licenses. submodule (swiftest) s_util + use whm + use rmvs + use helio + use symba + use fraggle contains module subroutine swiftest_util_append_arr_char_string(arr, source, nold, nsrc, lsource_mask) @@ -2541,6 +2546,374 @@ module subroutine swiftest_util_set_rhill_approximate(self,cb) end subroutine swiftest_util_set_rhill_approximate + module subroutine swiftest_util_setup_construct_system(nbody_system, param) + !! author: David A. Minton + !! + !! Constructor for a Swiftest nbody system. Creates the nbody system object based on the user-input integrator + !! + implicit none + ! Arguments + class(swiftest_nbody_system), allocatable, intent(inout) :: nbody_system !! Swiftest nbody_system object + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters + ! Internals + type(encounter_storage) :: encounter_history + type(collision_storage) :: collision_history + + select type(param) + class is (swiftest_parameters) + allocate(swiftest_storage(param%dump_cadence) :: param%system_history) + allocate(swiftest_netcdf_parameters :: param%system_history%nc) + call param%system_history%reset() + + select case(param%integrator) + case (INT_BS) + write(*,*) 'Bulirsch-Stoer integrator not yet enabled' + case (INT_HELIO) + allocate(helio_nbody_system :: nbody_system) + select type(nbody_system) + class is (helio_nbody_system) + allocate(helio_cb :: nbody_system%cb) + allocate(helio_pl :: nbody_system%pl) + allocate(helio_tp :: nbody_system%tp) + allocate(helio_tp :: nbody_system%tp_discards) + end select + param%collision_model = "MERGE" + case (INT_RA15) + write(*,*) 'Radau integrator not yet enabled' + case (INT_TU4) + write(*,*) 'INT_TU4 integrator not yet enabled' + case (INT_WHM) + allocate(whm_nbody_system :: nbody_system) + select type(nbody_system) + class is (whm_nbody_system) + allocate(whm_cb :: nbody_system%cb) + allocate(whm_pl :: nbody_system%pl) + allocate(whm_tp :: nbody_system%tp) + allocate(whm_tp :: nbody_system%tp_discards) + end select + param%collision_model = "MERGE" + case (INT_RMVS) + allocate(rmvs_nbody_system :: nbody_system) + select type(nbody_system) + class is (rmvs_nbody_system) + allocate(rmvs_cb :: nbody_system%cb) + allocate(rmvs_pl :: nbody_system%pl) + allocate(rmvs_tp :: nbody_system%tp) + allocate(rmvs_tp :: nbody_system%tp_discards) + end select + param%collision_model = "MERGE" + case (INT_SYMBA) + allocate(symba_nbody_system :: nbody_system) + select type(nbody_system) + class is (symba_nbody_system) + allocate(symba_cb :: nbody_system%cb) + allocate(symba_pl :: nbody_system%pl) + allocate(symba_tp :: nbody_system%tp) + + allocate(symba_tp :: nbody_system%tp_discards) + allocate(symba_pl :: nbody_system%pl_adds) + allocate(symba_pl :: nbody_system%pl_discards) + + allocate(symba_list_pltp :: nbody_system%pltp_encounter) + allocate(symba_list_plpl :: nbody_system%plpl_encounter) + allocate(collision_list_plpl :: nbody_system%plpl_collision) + + if (param%lenc_save_trajectory .or. param%lenc_save_closest) then + allocate(encounter_netcdf_parameters :: encounter_history%nc) + call encounter_history%reset() + select type(nc => encounter_history%nc) + class is (encounter_netcdf_parameters) + nc%file_number = param%iloop / param%dump_cadence + end select + allocate(nbody_system%encounter_history, source=encounter_history) + end if + + allocate(collision_netcdf_parameters :: collision_history%nc) + call collision_history%reset() + select type(nc => collision_history%nc) + class is (collision_netcdf_parameters) + nc%file_number = param%iloop / param%dump_cadence + end select + allocate(nbody_system%collision_history, source=collision_history) + + end select + case (INT_RINGMOONS) + write(*,*) 'RINGMOONS-SyMBA integrator not yet enabled' + case default + write(*,*) 'Unkown integrator',param%integrator + call util_exit(FAILURE) + end select + + allocate(swiftest_particle_info :: nbody_system%cb%info) + + select case(param%collision_model) + case("MERGE") + allocate(collision_merge :: nbody_system%collider) + case("BOUNCE") + allocate(collision_bounce :: nbody_system%collider) + case("SIMPLE") + allocate(collision_simple_disruption :: nbody_system%collider) + case("FRAGGLE") + allocate(collision_fraggle :: nbody_system%collider) + end select + call nbody_system%collider%setup(nbody_system) + + end select + + return + end subroutine swiftest_util_setup_construct_system + + + module subroutine swiftest_util_setup_initialize_particle_info_system(self, param) + !! author: David A. Minton + !! + !! Setup up particle information metadata from initial conditions + ! + implicit none + ! Arguments + class(swiftest_nbody_system), intent(inout) :: self !! Swiftest nbody system object + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters + ! Internals + integer(I4B) :: i + + associate(pl => self%pl, npl => self%pl%nbody, tp => self%tp, ntp => self%tp%nbody) + + if (.not. allocated(self%cb%info)) allocate(swiftest_particle_info :: self%cb%info) + + call self%cb%info%set_value(particle_type=CB_TYPE_NAME, status="ACTIVE", origin_type="Initial conditions", & + origin_time=param%t0, origin_rh=[0.0_DP, 0.0_DP, 0.0_DP], origin_vh=[0.0_DP, 0.0_DP, 0.0_DP]) + do i = 1, self%pl%nbody + call pl%info(i)%set_value(particle_type=PL_TYPE_NAME, status="ACTIVE", origin_type="Initial conditions", & + origin_time=param%t0, origin_rh=self%pl%rh(:,i), origin_vh=self%pl%vh(:,i)) + end do + do i = 1, self%tp%nbody + call tp%info(i)%set_value(particle_type=TP_TYPE_NAME, status="ACTIVE", origin_type="Initial conditions", & + origin_time=param%t0, origin_rh=self%tp%rh(:,i), origin_vh=self%tp%vh(:,i)) + end do + + end associate + + return + end subroutine swiftest_util_setup_initialize_particle_info_system + + + module subroutine swiftest_util_setup_initialize_system(self, param) + !! author: David A. Minton + !! + !! Wrapper method to initialize a basic Swiftest nbody system from files + !! + implicit none + ! Arguments + class(swiftest_nbody_system), intent(inout) :: self !! Swiftest nbody_system object + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters + + associate(nbody_system => self, cb => self%cb, pl => self%pl, tp => self%tp) + + call nbody_system%read_in(param) + call nbody_system%validate_ids(param) + call nbody_system%set_msys() + call pl%set_mu(cb) + call tp%set_mu(cb) + if (param%in_form == "EL") then + call pl%el2xv(cb) + call tp%el2xv(cb) + end if + call pl%flatten(param) + if (.not.param%lrhill_present) call pl%set_rhill(cb) + pl%lfirst = param%lfirstkick + tp%lfirst = param%lfirstkick + + if (.not.param%lrestart) then + call nbody_system%init_particle_info(param) + end if + end associate + + return + end subroutine swiftest_util_setup_initialize_system + + + module subroutine swiftest_util_setup_body(self, n, param) + !! author: David A. Minton + !! + !! Constructor for base Swiftest particle class. Allocates space for all particles and + !! initializes all components with a value. + !! Note: Timing tests indicate that (NDIM, n) is more efficient than (NDIM, n) + implicit none + ! Arguments + class(swiftest_body), intent(inout) :: self !! Swiftest generic body object + integer(I4B), intent(in) :: n !! Number of particles to allocate space for + class(swiftest_parameters), intent(in) :: param !! Current run configuration parameter + ! Internals + integer(I4B) :: i + + if (n < 0) return + + self%lfirst = .true. + + call self%dealloc() + + self%nbody = n + if (n == 0) return + + allocate(swiftest_particle_info :: self%info(n)) + allocate(self%id(n)) + allocate(self%status(n)) + allocate(self%ldiscard(n)) + allocate(self%lmask(n)) + allocate(self%mu(n)) + allocate(self%rh(NDIM, n)) + allocate(self%vh(NDIM, n)) + allocate(self%rb(NDIM, n)) + allocate(self%vb(NDIM, n)) + allocate(self%ah(NDIM, n)) + allocate(self%ir3h(n)) + allocate(self%aobl(NDIM, n)) + if (param%lclose) then + allocate(self%lcollision(n)) + allocate(self%lencounter(n)) + self%lcollision(:) = .false. + self%lencounter(:) = .false. + end if + + self%id(:) = 0 + do i = 1, n + call self%info(i)%set_value(& + name = "UNNAMED", & + particle_type = "UNKNOWN", & + status = "INACTIVE", & + origin_type = "UNKNOWN", & + collision_id = 0, & + origin_time = -huge(1.0_DP), & + origin_rh = [0.0_DP, 0.0_DP, 0.0_DP], & + origin_vh = [0.0_DP, 0.0_DP, 0.0_DP], & + discard_time = huge(1.0_DP), & + discard_rh = [0.0_DP, 0.0_DP, 0.0_DP], & + discard_vh = [0.0_DP, 0.0_DP, 0.0_DP], & + discard_body_id = -1 & + ) + end do + + self%status(:) = INACTIVE + self%ldiscard(:) = .false. + self%lmask(:) = .false. + self%mu(:) = 0.0_DP + self%rh(:,:) = 0.0_DP + self%vh(:,:) = 0.0_DP + self%rb(:,:) = 0.0_DP + self%vb(:,:) = 0.0_DP + self%ah(:,:) = 0.0_DP + self%ir3h(:) = 0.0_DP + self%aobl(:,:) = 0.0_DP + + if (param%ltides) then + allocate(self%atide(NDIM, n)) + self%atide(:,:) = 0.0_DP + end if + if (param%lgr) then + allocate(self%agr(NDIM, n)) + self%agr(:,:) = 0.0_DP + end if + + return + end subroutine swiftest_util_setup_body + + + module subroutine swiftest_util_setup_pl(self, n, param) + !! author: David A. Minton + !! + !! Constructor for base Swiftest massive body class. Allocates space for all particles and + !! initializes all components with a value. + implicit none + ! Arguments + class(swiftest_pl), intent(inout) :: self !! Swiftest massive body object + integer(I4B), intent(in) :: n !! Number of particles to allocate space for + class(swiftest_parameters), intent(in) :: param !! Current run configuration parameter + + !> Call allocation method for parent class + !> The parent class here is the abstract swiftest_body class, so we can't use the type-bound procedure + call swiftest_util_setup_body(self, n, param) + if (n == 0) return + + allocate(self%mass(n)) + allocate(self%Gmass(n)) + allocate(self%rhill(n)) + allocate(self%renc(n)) + + self%mass(:) = 0.0_DP + self%Gmass(:) = 0.0_DP + self%rhill(:) = 0.0_DP + self%renc(:) = 0.0_DP + + self%nplpl = 0 + + if (param%lclose) then + allocate(self%nplenc(n)) + allocate(self%ntpenc(n)) + allocate(self%radius(n)) + allocate(self%density(n)) + + self%nplenc(:) = 0 + self%ntpenc(:) = 0 + self%radius(:) = 0.0_DP + self%density(:) = 1.0_DP + + end if + + if (param%lmtiny_pl) then + allocate(self%lmtiny(n)) + self%lmtiny(:) = .false. + end if + + if (param%lrotation) then + allocate(self%rot(NDIM, n)) + allocate(self%Ip(NDIM, n)) + self%rot(:,:) = 0.0_DP + self%Ip(:,:) = 0.0_DP + end if + + if (param%ltides) then + allocate(self%k2(n)) + allocate(self%Q(n)) + allocate(self%tlag(n)) + self%k2(:) = 0.0_DP + self%Q(:) = 0.0_DP + self%tlag(:) = 0.0_DP + end if + + return + end subroutine swiftest_util_setup_pl + + + module subroutine swiftest_util_setup_tp(self, n, param) + !! author: David A. Minton + !! + !! Constructor for base Swiftest test particle particle class. Allocates space for + !! all particles and initializes all components with a value. + implicit none + ! Arguments + class(swiftest_tp), intent(inout) :: self !! Swiftest test particle object + integer(I4B), intent(in) :: n !! Number of particles to allocate space for + class(swiftest_parameters), intent(in) :: param !! Current run configuration parameter + + !> Call allocation method for parent class + !> The parent class here is the abstract swiftest_body class, so we can't use the type-bound procedure + call swiftest_util_setup_body(self, n, param) + if (n == 0) return + + allocate(self%isperi(n)) + allocate(self%peri(n)) + allocate(self%atp(n)) + allocate(self%nplenc(n)) + + self%isperi(:) = 0 + self%peri(:) = 0.0_DP + self%atp(:) = 0.0_DP + self%nplenc(:) = 0 + + return + end subroutine swiftest_util_setup_tp + + module subroutine swiftest_util_snapshot_system(self, param, nbody_system, t, arg) !! author: David A. Minton !! @@ -2628,7 +3001,7 @@ pure module subroutine swiftest_util_sort_dp(arr) ! Arguments real(DP), dimension(:), intent(inout) :: arr - call swiftest_qsort_DP(arr) + call swiftest_util_sort_qsort_DP(arr) return end subroutine swiftest_util_sort_dp @@ -2657,13 +3030,13 @@ pure module subroutine swiftest_util_sort_index_dp(arr, ind) end if allocate(tmparr, mold=arr) tmparr(:) = arr(ind(:)) - call swiftest_qsort_DP(tmparr, ind) + call swiftest_util_sort_qsort_DP(tmparr, ind) return end subroutine swiftest_util_sort_index_dp - recursive pure subroutine swiftest_qsort_DP(arr, ind) + recursive pure subroutine swiftest_util_sort_qsort_DP(arr, ind) !! author: David A. Minton !! !! Sort input DP precision array by index in ascending numerical order using quicksort sort. @@ -2677,21 +3050,21 @@ recursive pure subroutine swiftest_qsort_DP(arr, ind) if (size(arr) > 1) then if (present(ind)) then - call swiftest_partition_DP(arr, iq, ind) - call swiftest_qsort_DP(arr(:iq-1),ind(:iq-1)) - call swiftest_qsort_DP(arr(iq:), ind(iq:)) + call swiftest_util_sort_partition_DP(arr, iq, ind) + call swiftest_util_sort_qsort_DP(arr(:iq-1),ind(:iq-1)) + call swiftest_util_sort_qsort_DP(arr(iq:), ind(iq:)) else - call swiftest_partition_DP(arr, iq) - call swiftest_qsort_DP(arr(:iq-1)) - call swiftest_qsort_DP(arr(iq:)) + call swiftest_util_sort_partition_DP(arr, iq) + call swiftest_util_sort_qsort_DP(arr(:iq-1)) + call swiftest_util_sort_qsort_DP(arr(iq:)) end if end if return - end subroutine swiftest_qsort_DP + end subroutine swiftest_util_sort_qsort_DP - pure subroutine swiftest_partition_DP(arr, marker, ind) + pure subroutine swiftest_util_sort_partition_DP(arr, marker, ind) !! author: David A. Minton !! !! Partition function for quicksort on DP type @@ -2745,7 +3118,7 @@ pure subroutine swiftest_partition_DP(arr, marker, ind) end do return - end subroutine swiftest_partition_DP + end subroutine swiftest_util_sort_partition_DP pure module subroutine swiftest_util_sort_i4b(arr) @@ -2758,7 +3131,7 @@ pure module subroutine swiftest_util_sort_i4b(arr) ! Arguments integer(I4B), dimension(:), intent(inout) :: arr - call swiftest_qsort_I4B(arr) + call swiftest_util_sort_qsort_I4B(arr) return end subroutine swiftest_util_sort_i4b @@ -2786,7 +3159,7 @@ pure module subroutine swiftest_util_sort_index_I4B(arr, ind) end if allocate(tmparr, mold=arr) tmparr(:) = arr(ind(:)) - call swiftest_qsort_I4B(tmparr, ind) + call swiftest_util_sort_qsort_I4B(tmparr, ind) return end subroutine swiftest_util_sort_index_I4B @@ -2814,13 +3187,13 @@ pure module subroutine swiftest_util_sort_index_I4B_I8Bind(arr, ind) end if allocate(tmparr, mold=arr) tmparr(:) = arr(ind(:)) - call swiftest_qsort_I4B_I8Bind(tmparr, ind) + call swiftest_util_sort_qsort_I4B_I8Bind(tmparr, ind) return end subroutine swiftest_util_sort_index_I4B_I8Bind - recursive pure subroutine swiftest_qsort_I4B(arr, ind) + recursive pure subroutine swiftest_util_sort_qsort_I4B(arr, ind) !! author: David A. Minton !! !! Sort input I4B array by index in ascending numerical order using quicksort. @@ -2834,20 +3207,21 @@ recursive pure subroutine swiftest_qsort_I4B(arr, ind) if (size(arr) > 1) then if (present(ind)) then - call swiftest_partition_I4B(arr, iq, ind) - call swiftest_qsort_I4B(arr(:iq-1),ind(:iq-1)) - call swiftest_qsort_I4B(arr(iq:), ind(iq:)) + call swiftest_util_sort_partition_I4B(arr, iq, ind) + call swiftest_util_sort_qsort_I4B(arr(:iq-1),ind(:iq-1)) + call swiftest_util_sort_qsort_I4B(arr(iq:), ind(iq:)) else - call swiftest_partition_I4B(arr, iq) - call swiftest_qsort_I4B(arr(:iq-1)) - call swiftest_qsort_I4B(arr(iq:)) + call swiftest_util_sort_partition_I4B(arr, iq) + call swiftest_util_sort_qsort_I4B(arr(:iq-1)) + call swiftest_util_sort_qsort_I4B(arr(iq:)) end if end if return - end subroutine swiftest_qsort_I4B + end subroutine swiftest_util_sort_qsort_I4B - recursive pure subroutine swiftest_qsort_I4B_I8Bind(arr, ind) + + recursive pure subroutine swiftest_util_sort_qsort_I4B_I8Bind(arr, ind) !! author: David A. Minton !! !! Sort input I4B array by index in ascending numerical order using quicksort. @@ -2861,21 +3235,21 @@ recursive pure subroutine swiftest_qsort_I4B_I8Bind(arr, ind) if (size(arr) > 1_I8B) then if (present(ind)) then - call swiftest_partition_I4B_I8Bind(arr, iq, ind) - call swiftest_qsort_I4B_I8Bind(arr(:iq-1_I8B),ind(:iq-1_I8B)) - call swiftest_qsort_I4B_I8Bind(arr(iq:), ind(iq:)) + call swiftest_util_sort_partition_I4B_I8Bind(arr, iq, ind) + call swiftest_util_sort_qsort_I4B_I8Bind(arr(:iq-1_I8B),ind(:iq-1_I8B)) + call swiftest_util_sort_qsort_I4B_I8Bind(arr(iq:), ind(iq:)) else - call swiftest_partition_I4B_I8Bind(arr, iq) - call swiftest_qsort_I4B_I8Bind(arr(:iq-1_I8B)) - call swiftest_qsort_I4B_I8Bind(arr(iq:)) + call swiftest_util_sort_partition_I4B_I8Bind(arr, iq) + call swiftest_util_sort_qsort_I4B_I8Bind(arr(:iq-1_I8B)) + call swiftest_util_sort_qsort_I4B_I8Bind(arr(iq:)) end if end if return - end subroutine swiftest_qsort_I4B_I8Bind + end subroutine swiftest_util_sort_qsort_I4B_I8Bind - recursive pure subroutine swiftest_qsort_I8B_I8Bind(arr, ind) + recursive pure subroutine swiftest_util_sort_qsort_I8B_I8Bind(arr, ind) !! author: David A. Minton !! !! Sort input I8B array by index in ascending numerical order using quicksort. @@ -2889,21 +3263,21 @@ recursive pure subroutine swiftest_qsort_I8B_I8Bind(arr, ind) if (size(arr) > 1_I8B) then if (present(ind)) then - call swiftest_partition_I8B_I8Bind(arr, iq, ind) - call swiftest_qsort_I8B_I8Bind(arr(:iq-1_I8B),ind(:iq-1_I8B)) - call swiftest_qsort_I8B_I8Bind(arr(iq:), ind(iq:)) + call swiftest_util_sort_partition_I8B_I8Bind(arr, iq, ind) + call swiftest_util_sort_qsort_I8B_I8Bind(arr(:iq-1_I8B),ind(:iq-1_I8B)) + call swiftest_util_sort_qsort_I8B_I8Bind(arr(iq:), ind(iq:)) else - call swiftest_partition_I8B_I8Bind(arr, iq) - call swiftest_qsort_I8B_I8Bind(arr(:iq-1_I8B)) - call swiftest_qsort_I8B_I8Bind(arr(iq:)) + call swiftest_util_sort_partition_I8B_I8Bind(arr, iq) + call swiftest_util_sort_qsort_I8B_I8Bind(arr(:iq-1_I8B)) + call swiftest_util_sort_qsort_I8B_I8Bind(arr(iq:)) end if end if return - end subroutine swiftest_qsort_I8B_I8Bind + end subroutine swiftest_util_sort_qsort_I8B_I8Bind - pure subroutine swiftest_partition_I4B(arr, marker, ind) + pure subroutine swiftest_util_sort_partition_I4B(arr, marker, ind) !! author: David A. Minton !! !! Partition function for quicksort on I4B type @@ -2957,9 +3331,10 @@ pure subroutine swiftest_partition_I4B(arr, marker, ind) end do return - end subroutine swiftest_partition_I4B + end subroutine swiftest_util_sort_partition_I4B + - pure subroutine swiftest_partition_I4B_I8Bind(arr, marker, ind) + pure subroutine swiftest_util_sort_partition_I4B_I8Bind(arr, marker, ind) !! author: David A. Minton !! !! Partition function for quicksort on I4B type @@ -3013,9 +3388,10 @@ pure subroutine swiftest_partition_I4B_I8Bind(arr, marker, ind) end do return - end subroutine swiftest_partition_I4B_I8Bind + end subroutine swiftest_util_sort_partition_I4B_I8Bind + - pure subroutine swiftest_partition_I8B_I8Bind(arr, marker, ind) + pure subroutine swiftest_util_sort_partition_I8B_I8Bind(arr, marker, ind) !! author: David A. Minton !! !! Partition function for quicksort on I8B type with I8B index @@ -3069,7 +3445,7 @@ pure subroutine swiftest_partition_I8B_I8Bind(arr, marker, ind) end do return - end subroutine swiftest_partition_I8B_I8Bind + end subroutine swiftest_util_sort_partition_I8B_I8Bind pure module subroutine swiftest_util_sort_sp(arr) @@ -3081,7 +3457,7 @@ pure module subroutine swiftest_util_sort_sp(arr) ! Arguments real(SP), dimension(:), intent(inout) :: arr - call swiftest_qsort_SP(arr) + call swiftest_util_sort_qsort_SP(arr) return end subroutine swiftest_util_sort_sp @@ -3109,13 +3485,13 @@ pure module subroutine swiftest_util_sort_index_sp(arr, ind) end if allocate(tmparr, mold=arr) tmparr(:) = arr(ind(:)) - call swiftest_qsort_SP(tmparr, ind) + call swiftest_util_sort_qsort_SP(tmparr, ind) return end subroutine swiftest_util_sort_index_sp - recursive pure subroutine swiftest_qsort_SP(arr, ind) + recursive pure subroutine swiftest_util_sort_qsort_SP(arr, ind) !! author: David A. Minton !! !! Sort input DP precision array by index in ascending numerical order using quicksort. @@ -3129,21 +3505,21 @@ recursive pure subroutine swiftest_qsort_SP(arr, ind) if (size(arr) > 1) then if (present(ind)) then - call swiftest_partition_SP(arr, iq, ind) - call swiftest_qsort_SP(arr(:iq-1),ind(:iq-1)) - call swiftest_qsort_SP(arr(iq:), ind(iq:)) + call swiftest_util_sort_partition_SP(arr, iq, ind) + call swiftest_util_sort_qsort_SP(arr(:iq-1),ind(:iq-1)) + call swiftest_util_sort_qsort_SP(arr(iq:), ind(iq:)) else - call swiftest_partition_SP(arr, iq) - call swiftest_qsort_SP(arr(:iq-1)) - call swiftest_qsort_SP(arr(iq:)) + call swiftest_util_sort_partition_SP(arr, iq) + call swiftest_util_sort_qsort_SP(arr(:iq-1)) + call swiftest_util_sort_qsort_SP(arr(iq:)) end if end if return - end subroutine swiftest_qsort_SP + end subroutine swiftest_util_sort_qsort_SP - pure subroutine swiftest_partition_SP(arr, marker, ind) + pure subroutine swiftest_util_sort_partition_SP(arr, marker, ind) !! author: David A. Minton !! !! Partition function for quicksort on SP type @@ -3197,7 +3573,7 @@ pure subroutine swiftest_partition_SP(arr, marker, ind) end do return - end subroutine swiftest_partition_SP + end subroutine swiftest_util_sort_partition_SP module subroutine swiftest_util_sort_pl(self, sortby, ascending) diff --git a/src/symba/symba_module.f90 b/src/symba/symba_module.f90 index 280fe2af8..f120919fe 100644 --- a/src/symba/symba_module.f90 +++ b/src/symba/symba_module.f90 @@ -40,7 +40,7 @@ module symba procedure :: gr_pos_kick => symba_gr_p4_pl !! Position kick due to p**4 term in the post-Newtonian correction procedure :: accel_int => symba_kick_getacch_int_pl !! Compute direct cross (third) term heliocentric accelerations of massive bodiess, with no mutual interactions between bodies below GMTINY procedure :: accel => symba_kick_getacch_pl !! Compute heliocentric accelerations of massive bodies - procedure :: setup => symba_setup_pl !! Constructor method - Allocates space for the input number of bodies + procedure :: setup => symba_util_setup_pl !! Constructor method - Allocates space for the input number of bodies procedure :: append => symba_util_append_pl !! Appends elements from one structure to another procedure :: dealloc => symba_util_dealloc_pl !! Deallocates all allocatable arrays procedure :: fill => symba_util_fill_pl !! "Fills" bodies from one object into another depending on the results of a mask (uses the UNPACK intrinsic) @@ -62,7 +62,7 @@ module symba procedure :: encounter_check => symba_encounter_check_tp !! Checks if any test particles are undergoing a close encounter with a massive body procedure :: gr_pos_kick => symba_gr_p4_tp !! Position kick due to p**4 term in the post-Newtonian correction procedure :: accel => symba_kick_getacch_tp !! Compute heliocentric accelerations of test particles - procedure :: setup => symba_setup_tp !! Constructor method - Allocates space for the input number of bodies + procedure :: setup => symba_util_setup_tp !! Constructor method - Allocates space for the input number of bodies procedure :: append => symba_util_append_tp !! Appends elements from one structure to another procedure :: dealloc => symba_util_dealloc_tp !! Deallocates all allocatable arrays procedure :: fill => symba_util_fill_tp !! "Fills" bodies from one object into another depending on the results of a mask (uses the UNPACK intrinsic) @@ -95,7 +95,7 @@ module symba type, extends(helio_nbody_system) :: symba_nbody_system integer(I4B) :: irec !! nbody_system recursion level contains - procedure :: initialize => symba_setup_initialize_system !! Performs SyMBA-specific initilization steps + procedure :: initialize => symba_util_setup_initialize_system !! Performs SyMBA-specific initilization steps procedure :: step => symba_step_system !! Advance the SyMBA nbody system forward in time by one step procedure :: interp => symba_step_interp_system !! Perform an interpolation step on the SymBA nbody system procedure :: set_recur_levels => symba_step_set_recur_levels_system !! Sets recursion levels of bodies and encounter lists to the current nbody_system level @@ -243,25 +243,25 @@ module subroutine symba_kick_list_pltp(self, nbody_system, dt, irec, sgn) integer(I4B), intent(in) :: sgn !! sign to be applied to acceleration end subroutine symba_kick_list_pltp - module subroutine symba_setup_initialize_system(self, param) + module subroutine symba_util_setup_initialize_system(self, param) implicit none class(symba_nbody_system), intent(inout) :: self !! SyMBA nbody_system object class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters - end subroutine symba_setup_initialize_system + end subroutine symba_util_setup_initialize_system - module subroutine symba_setup_pl(self, n, param) + module subroutine symba_util_setup_pl(self, n, param) implicit none class(symba_pl), intent(inout) :: self !! SyMBA massive body object integer(I4B), intent(in) :: n !! Number of particles to allocate space for class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters - end subroutine symba_setup_pl + end subroutine symba_util_setup_pl - module subroutine symba_setup_tp(self, n, param) + module subroutine symba_util_setup_tp(self, n, param) implicit none class(symba_tp), intent(inout) :: self !! SyMBA test particle object integer(I4B), intent(in) :: n !! Number of particles to allocate space for class(swiftest_parameters), intent(in) :: param !! Current run configuration parameter - end subroutine symba_setup_tp + end subroutine symba_util_setup_tp module subroutine symba_step_system(self, param, t, dt) implicit none diff --git a/src/symba/symba_setup.f90 b/src/symba/symba_setup.f90 deleted file mode 100644 index 20e713c5c..000000000 --- a/src/symba/symba_setup.f90 +++ /dev/null @@ -1,98 +0,0 @@ -!! Copyright 2022 - David Minton, Carlisle Wishard, Jennifer Pouplin, Jake Elliott, & Dana Singh -!! This file is part of Swiftest. -!! Swiftest is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License -!! as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version. -!! Swiftest is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty -!! of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. -!! You should have received a copy of the GNU General Public License along with Swiftest. -!! If not, see: https://www.gnu.org/licenses. - -submodule(symba) s_symba_setup - use swiftest -contains - - module subroutine symba_setup_initialize_system(self, param) - !! author: David A. Minton - !! - !! Initialize an SyMBA nbody system from files and sets up the planetocentric structures. - !! This subroutine will also sort the massive bodies in descending order by mass - !! - implicit none - ! Arguments - class(symba_nbody_system), intent(inout) :: self !! SyMBA nbody_system object - class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters - - ! Call parent method - associate(nbody_system => self) - call helio_setup_initialize_system(nbody_system, param) - call nbody_system%pltp_encounter%setup(0_I8B) - call nbody_system%plpl_encounter%setup(0_I8B) - call nbody_system%plpl_collision%setup(0_I8B) - end associate - - return - end subroutine symba_setup_initialize_system - - - module subroutine symba_setup_pl(self, n, param) - !! author: David A. Minton - !! - !! Allocate SyMBA test particle structure - !! - !! Equivalent in functionality to David E. Kaufmann's Swifter routine symba_setup.f90 - implicit none - ! Arguments - class(symba_pl), intent(inout) :: self !! SyMBA massive body object - integer(I4B), intent(in) :: n !! Number of particles to allocate space for - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameter - ! Internals - integer(I4B) :: i - - !> Call allocation method for parent class. - call self%helio_pl%setup(n, param) - if (n == 0) return - - allocate(self%levelg(n)) - allocate(self%levelm(n)) - allocate(self%isperi(n)) - allocate(self%peri(n)) - allocate(self%atp(n)) - allocate(self%kin(n)) - - - self%levelg(:) = -1 - self%levelm(:) = -1 - self%isperi(:) = 0 - self%peri(:) = 0.0_DP - self%atp(:) = 0.0_DP - call self%reset_kinship([(i, i=1, n)]) - return - end subroutine symba_setup_pl - - - module subroutine symba_setup_tp(self, n, param) - !! author: David A. Minton - !! - !! Allocate WHM test particle structure - !! - !! Equivalent in functionality to David E. Kaufmann's Swifter routine whm_setup.f90 - implicit none - ! Arguments - class(symba_tp), intent(inout) :: self !! SyMBA test particle object - integer(I4B), intent(in) :: n !! Number of particles to allocate space for - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameter - - !> Call allocation method for parent class. - call self%helio_tp%setup(n, param) - if (n == 0) return - - allocate(self%levelg(n)) - allocate(self%levelm(n)) - - self%levelg(:) = -1 - self%levelm(:) = -1 - - return - end subroutine symba_setup_tp - -end submodule s_symba_setup diff --git a/src/symba/symba_util.f90 b/src/symba/symba_util.f90 index 18bbbdfdb..87ba8dacf 100644 --- a/src/symba/symba_util.f90 +++ b/src/symba/symba_util.f90 @@ -272,6 +272,90 @@ module subroutine symba_util_set_renc(self, scale) return end subroutine symba_util_set_renc + module subroutine symba_util_setup_initialize_system(self, param) + !! author: David A. Minton + !! + !! Initialize an SyMBA nbody system from files and sets up the planetocentric structures. + !! This subroutine will also sort the massive bodies in descending order by mass + !! + implicit none + ! Arguments + class(symba_nbody_system), intent(inout) :: self !! SyMBA nbody_system object + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters + + ! Call parent method + associate(nbody_system => self) + call helio_util_setup_initialize_system(nbody_system, param) + call nbody_system%pltp_encounter%setup(0_I8B) + call nbody_system%plpl_encounter%setup(0_I8B) + call nbody_system%plpl_collision%setup(0_I8B) + end associate + + return + end subroutine symba_util_setup_initialize_system + + + module subroutine symba_util_setup_pl(self, n, param) + !! author: David A. Minton + !! + !! Allocate SyMBA test particle structure + !! + !! Equivalent in functionality to David E. Kaufmann's Swifter routine symba_util_setup.f90 + implicit none + ! Arguments + class(symba_pl), intent(inout) :: self !! SyMBA massive body object + integer(I4B), intent(in) :: n !! Number of particles to allocate space for + class(swiftest_parameters), intent(in) :: param !! Current run configuration parameter + ! Internals + integer(I4B) :: i + + !> Call allocation method for parent class. + call self%helio_pl%setup(n, param) + if (n == 0) return + + allocate(self%levelg(n)) + allocate(self%levelm(n)) + allocate(self%isperi(n)) + allocate(self%peri(n)) + allocate(self%atp(n)) + allocate(self%kin(n)) + + + self%levelg(:) = -1 + self%levelm(:) = -1 + self%isperi(:) = 0 + self%peri(:) = 0.0_DP + self%atp(:) = 0.0_DP + call self%reset_kinship([(i, i=1, n)]) + return + end subroutine symba_util_setup_pl + + + module subroutine symba_util_setup_tp(self, n, param) + !! author: David A. Minton + !! + !! Allocate WHM test particle structure + !! + !! Equivalent in functionality to David E. Kaufmann's Swifter routine whm_util_setup.f90 + implicit none + ! Arguments + class(symba_tp), intent(inout) :: self !! SyMBA test particle object + integer(I4B), intent(in) :: n !! Number of particles to allocate space for + class(swiftest_parameters), intent(in) :: param !! Current run configuration parameter + + !> Call allocation method for parent class. + call self%helio_tp%setup(n, param) + if (n == 0) return + + allocate(self%levelg(n)) + allocate(self%levelm(n)) + + self%levelg(:) = -1 + self%levelm(:) = -1 + + return + end subroutine symba_util_setup_tp + module subroutine symba_util_sort_pl(self, sortby, ascending) !! author: David A. Minton diff --git a/src/whm/whm_module.f90 b/src/whm/whm_module.f90 index fad0d8a55..d902065f6 100644 --- a/src/whm/whm_module.f90 +++ b/src/whm/whm_module.f90 @@ -30,7 +30,7 @@ module whm real(DP), dimension(:), allocatable :: muj !! Jacobi mu: GMcb * eta(i) / eta(i - 1) real(DP), dimension(:), allocatable :: ir3j !! Third term of heliocentric acceleration !! 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 + !! component list, such as whm_util_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 @@ -49,7 +49,7 @@ module whm 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 :: spill => whm_util_spill_pl !!"Spills" bodies from one object to another depending on the results of a mask (uses the PACK intrinsic) - procedure :: setup => whm_setup_pl !! Constructor method - Allocates space for the input number of bodiess + procedure :: setup => whm_util_setup_pl !! Constructor method - Allocates space for the input number of bodiess procedure :: step => whm_step_pl !! Steps the body forward one stepsize final :: whm_final_pl !! Finalizes the WHM massive body object - deallocates all allocatables end type whm_pl @@ -72,7 +72,7 @@ module whm type, extends(swiftest_nbody_system) :: whm_nbody_system contains !> Replace the abstract procedures with concrete ones - procedure :: initialize => whm_setup_initialize_system !! Performs WHM-specific initilization steps, like calculating the Jacobi masses + procedure :: initialize => whm_util_setup_initialize_system !! Performs WHM-specific initilization steps, like calculating the Jacobi masses procedure :: step => whm_step_system !! Advance the WHM nbody system forward in time by one step final :: whm_final_system !! Finalizes the WHM nbody_system object - deallocates all allocatables end type whm_nbody_system @@ -173,18 +173,18 @@ pure module subroutine whm_gr_p4_tp(self, nbody_system, param, dt) end subroutine whm_gr_p4_tp !> Reads WHM massive body object in from file - module subroutine whm_setup_pl(self, n, param) + module subroutine whm_util_setup_pl(self, n, param) implicit none class(whm_pl), intent(inout) :: self !! WHM massive body objectobject integer(I4B), intent(in) :: n !! Number of particles to allocate space for class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters - end subroutine whm_setup_pl + end subroutine whm_util_setup_pl - module subroutine whm_setup_initialize_system(self, param) + module subroutine whm_util_setup_initialize_system(self, param) implicit none class(whm_nbody_system), intent(inout) :: self !! WHM nbody system object class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters - end subroutine whm_setup_initialize_system + end subroutine whm_util_setup_initialize_system module subroutine whm_step_pl(self, nbody_system, param, t, dt) implicit none diff --git a/src/whm/whm_setup.f90 b/src/whm/whm_setup.f90 deleted file mode 100644 index fe4745777..000000000 --- a/src/whm/whm_setup.f90 +++ /dev/null @@ -1,100 +0,0 @@ -!! Copyright 2022 - David Minton, Carlisle Wishard, Jennifer Pouplin, Jake Elliott, & Dana Singh -!! This file is part of Swiftest. -!! Swiftest is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License -!! as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version. -!! Swiftest is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty -!! of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. -!! You should have received a copy of the GNU General Public License along with Swiftest. -!! If not, see: https://www.gnu.org/licenses. - -submodule(whm) s_whm_setup - use swiftest -contains - - module subroutine whm_setup_pl(self, n, param) - !! author: David A. Minton - !! - !! Allocate WHM planet structure - !! - !! Equivalent in functionality to David E. Kaufmann's Swifter routine whm_setup.f90 - implicit none - ! Arguments - class(whm_pl), intent(inout) :: self !! Swiftest test particle object - integer(I4B), intent(in) :: n !! Number of particles to allocate space for - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameter - - !> Call allocation method for parent class - call swiftest_setup_pl(self, n, param) - if (n == 0) return - - allocate(self%eta(n)) - allocate(self%muj(n)) - allocate(self%xj(NDIM, n)) - allocate(self%vj(NDIM, n)) - allocate(self%ir3j(n)) - - self%eta(:) = 0.0_DP - self%muj(:) = 0.0_DP - self%xj(:,:) = 0.0_DP - self%vj(:,:) = 0.0_DP - self%ir3j(:) = 0.0_DP - - return - end subroutine whm_setup_pl - - - module subroutine whm_util_set_mu_eta_pl(self, cb) - !! author: David A. Minton - !! - !! Sets the Jacobi mass value eta for all massive bodies - implicit none - ! Arguments - class(whm_pl), intent(inout) :: self !! WHM nbody_system object - class(swiftest_cb), intent(inout) :: cb !! Swiftest central body object - ! Internals - integer(I4B) :: i - - associate(pl => self, npl => self%nbody) - if (npl == 0) return - call swiftest_util_set_mu_pl(pl, cb) - pl%eta(1) = cb%Gmass + pl%Gmass(1) - pl%muj(1) = pl%eta(1) - do i = 2, npl - pl%eta(i) = pl%eta(i - 1) + pl%Gmass(i) - pl%muj(i) = cb%Gmass * pl%eta(i) / pl%eta(i - 1) - end do - end associate - - return - end subroutine whm_util_set_mu_eta_pl - - - module subroutine whm_setup_initialize_system(self, param) - !! author: David A. Minton - !! - !! Initialize a WHM nbody system from files - !! - implicit none - ! Arguments - class(whm_nbody_system), intent(inout) :: self !! WHM nbody system object - class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters - - call swiftest_setup_initialize_system(self, param) - ! First we need to make sure that the massive bodies are sorted by heliocentric distance before computing jacobies - call swiftest_util_set_ir3h(self%pl) - call self%pl%sort("ir3h", ascending=.false.) - call self%pl%flatten(param) - - ! Make sure that the discard list gets allocated initially - call self%tp_discards%setup(0, param) - call self%pl%set_mu(self%cb) - call self%tp%set_mu(self%cb) - if (param%lgr .and. param%in_type == "ASCII") then !! pseudovelocity conversion for NetCDF input files is handled by NetCDF routines - call self%pl%v2pv(param) - call self%tp%v2pv(param) - end if - - return - end subroutine whm_setup_initialize_system - -end submodule s_whm_setup \ No newline at end of file diff --git a/src/whm/whm_util.f90 b/src/whm/whm_util.f90 index 6ce3ce8df..a9377536d 100644 --- a/src/whm/whm_util.f90 +++ b/src/whm/whm_util.f90 @@ -142,6 +142,93 @@ module subroutine whm_util_set_ir3j(self) end subroutine whm_util_set_ir3j + module subroutine whm_util_setup_pl(self, n, param) + !! author: David A. Minton + !! + !! Allocate WHM planet structure + !! + !! Equivalent in functionality to David E. Kaufmann's Swifter routine whm_util_setup.f90 + implicit none + ! Arguments + class(whm_pl), intent(inout) :: self !! Swiftest test particle object + integer(I4B), intent(in) :: n !! Number of particles to allocate space for + class(swiftest_parameters), intent(in) :: param !! Current run configuration parameter + + !> Call allocation method for parent class + call swiftest_util_setup_pl(self, n, param) + if (n == 0) return + + allocate(self%eta(n)) + allocate(self%muj(n)) + allocate(self%xj(NDIM, n)) + allocate(self%vj(NDIM, n)) + allocate(self%ir3j(n)) + + self%eta(:) = 0.0_DP + self%muj(:) = 0.0_DP + self%xj(:,:) = 0.0_DP + self%vj(:,:) = 0.0_DP + self%ir3j(:) = 0.0_DP + + return + end subroutine whm_util_setup_pl + + + module subroutine whm_util_set_mu_eta_pl(self, cb) + !! author: David A. Minton + !! + !! Sets the Jacobi mass value eta for all massive bodies + implicit none + ! Arguments + class(whm_pl), intent(inout) :: self !! WHM nbody_system object + class(swiftest_cb), intent(inout) :: cb !! Swiftest central body object + ! Internals + integer(I4B) :: i + + associate(pl => self, npl => self%nbody) + if (npl == 0) return + call swiftest_util_set_mu_pl(pl, cb) + pl%eta(1) = cb%Gmass + pl%Gmass(1) + pl%muj(1) = pl%eta(1) + do i = 2, npl + pl%eta(i) = pl%eta(i - 1) + pl%Gmass(i) + pl%muj(i) = cb%Gmass * pl%eta(i) / pl%eta(i - 1) + end do + end associate + + return + end subroutine whm_util_set_mu_eta_pl + + + module subroutine whm_util_setup_initialize_system(self, param) + !! author: David A. Minton + !! + !! Initialize a WHM nbody system from files + !! + implicit none + ! Arguments + class(whm_nbody_system), intent(inout) :: self !! WHM nbody system object + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters + + call swiftest_util_setup_initialize_system(self, param) + ! First we need to make sure that the massive bodies are sorted by heliocentric distance before computing jacobies + call swiftest_util_set_ir3h(self%pl) + call self%pl%sort("ir3h", ascending=.false.) + call self%pl%flatten(param) + + ! Make sure that the discard list gets allocated initially + call self%tp_discards%setup(0, param) + call self%pl%set_mu(self%cb) + call self%tp%set_mu(self%cb) + if (param%lgr .and. param%in_type == "ASCII") then !! pseudovelocity conversion for NetCDF input files is handled by NetCDF routines + call self%pl%v2pv(param) + call self%tp%v2pv(param) + end if + + return + end subroutine whm_util_setup_initialize_system + + module subroutine whm_util_sort_pl(self, sortby, ascending) !! author: David A. Minton !!