From 813c736fb8dca38ec5d85a6d7f9f415210a01ee6 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Fri, 2 Dec 2022 15:06:30 -0500 Subject: [PATCH] Added mass, radius, and name variables to encounter list. Also added collision time to SyMBAs version of the list. Added methods needed for basic operations on the lists --- src/encounter/encounter_setup.f90 | 26 ++++++++++++++++++++++--- src/encounter/encounter_util.f90 | 30 +++++++++++++++++++++++++---- src/modules/encounter_classes.f90 | 30 +++++++++++++++++------------ src/modules/symba_classes.f90 | 17 ++++++++-------- src/symba/symba_collision.f90 | 2 +- src/symba/symba_encounter_check.f90 | 15 ++++++++++++--- src/symba/symba_setup.f90 | 3 +++ src/symba/symba_step.f90 | 3 ++- src/symba/symba_util.f90 | 26 ++++++++++++++++++++++--- 9 files changed, 117 insertions(+), 35 deletions(-) diff --git a/src/encounter/encounter_setup.f90 b/src/encounter/encounter_setup.f90 index 95b2680a0..c741cf0a7 100644 --- a/src/encounter/encounter_setup.f90 +++ b/src/encounter/encounter_setup.f90 @@ -62,6 +62,8 @@ module subroutine encounter_setup_list(self, n) ! 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 @@ -75,10 +77,16 @@ module subroutine encounter_setup_list(self, n) if (allocated(self%x2)) deallocate(self%x2) if (allocated(self%v1)) deallocate(self%v1) if (allocated(self%v2)) deallocate(self%v2) - if (allocated(self%t)) deallocate(self%t) + if (allocated(self%Gmass1)) deallocate(self%Gmass1) + if (allocated(self%Gmass2)) deallocate(self%Gmass2) + if (allocated(self%radius1)) deallocate(self%radius1) + if (allocated(self%radius2)) deallocate(self%radius2) + if (allocated(self%name1)) deallocate(self%name1) + if (allocated(self%name2)) deallocate(self%name2) self%nenc = n if (n == 0_I8B) return + self%t = 0.0_DP allocate(self%lvdotr(n)) allocate(self%status(n)) @@ -90,7 +98,12 @@ module subroutine encounter_setup_list(self, n) allocate(self%x2(NDIM,n)) allocate(self%v1(NDIM,n)) allocate(self%v2(NDIM,n)) - allocate(self%t(n)) + allocate(self%Gmass1(n)) + allocate(self%Gmass2(n)) + allocate(self%radius1(n)) + allocate(self%radius2(n)) + allocate(self%name1(n)) + allocate(self%name2(n)) self%lvdotr(:) = .false. self%status(:) = INACTIVE @@ -102,7 +115,14 @@ module subroutine encounter_setup_list(self, n) self%x2(:,:) = 0.0_DP self%v1(:,:) = 0.0_DP self%v2(:,:) = 0.0_DP - self%t(:) = 0.0_DP + self%Gmass1(:) = 0.0_DP + self%Gmass2(:) = 0.0_DP + self%radius1(:) = 0.0_DP + self%radius2(:) = 0.0_DP + do i = 1_I8B, n + self%name1(i) = "UNNAMED" + self%name2(i) = "UNNAMED" + end do return end subroutine encounter_setup_list diff --git a/src/encounter/encounter_util.f90 b/src/encounter/encounter_util.f90 index c2597377a..e6d70dd53 100644 --- a/src/encounter/encounter_util.f90 +++ b/src/encounter/encounter_util.f90 @@ -36,7 +36,12 @@ module subroutine encounter_util_append_list(self, source, lsource_mask) call util_append(self%x2, source%x2, nold, nsrc, lsource_mask) call util_append(self%v1, source%v1, nold, nsrc, lsource_mask) call util_append(self%v2, source%v2, nold, nsrc, lsource_mask) - call util_append(self%t, source%t, nold, nsrc, lsource_mask) + call util_append(self%Gmass1, source%Gmass1, nold, nsrc, lsource_mask) + call util_append(self%Gmass2, source%Gmass2, nold, nsrc, lsource_mask) + call util_append(self%radius1, source%radius1, nold, nsrc, lsource_mask) + call util_append(self%radius2, source%radius2, nold, nsrc, lsource_mask) + call util_append(self%name1, source%name1, nold, nsrc, lsource_mask) + call util_append(self%name2, source%name2, nold, nsrc, lsource_mask) self%nenc = nold + count(lsource_mask(1:nsrc)) return @@ -54,6 +59,7 @@ module subroutine encounter_util_copy_list(self, source) associate(n => source%nenc) self%nenc = n + self%t = source%t self%lvdotr(1:n) = source%lvdotr(1:n) self%status(1:n) = source%status(1:n) self%index1(1:n) = source%index1(1:n) @@ -64,7 +70,12 @@ module subroutine encounter_util_copy_list(self, source) self%x2(:,1:n) = source%x2(:,1:n) self%v1(:,1:n) = source%v1(:,1:n) self%v2(:,1:n) = source%v2(:,1:n) - self%t(1:n) = source%t(1:n) + self%Gmass1(1:n) = source%Gmass1(1:n) + self%Gmass2(1:n) = source%Gmass2(1:n) + self%radius1(1:n) = source%radius1(1:n) + self%radius2(1:n) = source%radius2(1:n) + self%name1(1:n) = source%name1(1:n) + self%name2(1:n) = source%name2(1:n) end associate return @@ -104,7 +115,12 @@ module subroutine encounter_util_dealloc_list(self) if (allocated(self%x2)) deallocate(self%x2) if (allocated(self%v1)) deallocate(self%v1) if (allocated(self%v2)) deallocate(self%v2) - if (allocated(self%t)) deallocate(self%t) + if (allocated(self%Gmass1)) deallocate(self%Gmass1) + if (allocated(self%Gmass2)) deallocate(self%Gmass2) + if (allocated(self%radius1)) deallocate(self%radius1) + if (allocated(self%radius2)) deallocate(self%radius2) + if (allocated(self%name1)) deallocate(self%name1) + if (allocated(self%name2)) deallocate(self%name2) return end subroutine encounter_util_dealloc_list @@ -214,6 +230,7 @@ module subroutine encounter_util_resize_storage(self, nnew) return end subroutine encounter_util_resize_storage + module subroutine encounter_util_spill_list(self, discards, lspill_list, ldestructive) !! author: David A. Minton !! @@ -238,7 +255,12 @@ module subroutine encounter_util_spill_list(self, discards, lspill_list, ldestru call util_spill(keeps%x2, discards%x2, lspill_list, ldestructive) call util_spill(keeps%v1, discards%v1, lspill_list, ldestructive) call util_spill(keeps%v2, discards%v2, lspill_list, ldestructive) - call util_spill(keeps%t, discards%t, lspill_list, ldestructive) + call util_spill(keeps%Gmass1, discards%Gmass1, lspill_list, ldestructive) + call util_spill(keeps%Gmass2, discards%Gmass2, lspill_list, ldestructive) + call util_spill(keeps%radius1, discards%radius1, lspill_list, ldestructive) + call util_spill(keeps%radius2, discards%radius2, lspill_list, ldestructive) + call util_spill(keeps%name1, discards%name1, lspill_list, ldestructive) + call util_spill(keeps%name2, discards%name2, lspill_list, ldestructive) nenc_old = keeps%nenc diff --git a/src/modules/encounter_classes.f90 b/src/modules/encounter_classes.f90 index 10c9cceec..b8f26ce23 100644 --- a/src/modules/encounter_classes.f90 +++ b/src/modules/encounter_classes.f90 @@ -19,18 +19,24 @@ module encounter_classes integer(I4B), parameter :: SWEEPDIM = 3 type :: encounter_list - integer(I8B) :: nenc = 0 !! Total number of encounters - logical, dimension(:), allocatable :: lvdotr !! relative vdotr flag - integer(I4B), dimension(:), allocatable :: status !! status of the interaction - integer(I4B), dimension(:), allocatable :: index1 !! position of the first body in the encounter - integer(I4B), dimension(:), allocatable :: index2 !! position of the second body in the encounter - integer(I4B), dimension(:), allocatable :: id1 !! id of the first body in the encounter - integer(I4B), dimension(:), allocatable :: id2 !! id of the second body in the encounter - real(DP), dimension(:,:), allocatable :: x1 !! the position of body 1 in the encounter - real(DP), dimension(:,:), allocatable :: x2 !! the position of body 2 in the encounter - real(DP), dimension(:,:), allocatable :: v1 !! the velocity of body 1 in the encounter - real(DP), dimension(:,:), allocatable :: v2 !! the velocity of body 2 in the encounter - real(DP), dimension(:), allocatable :: t !! Time of encounter + integer(I8B) :: nenc = 0 !! Total number of encounters + real(DP) :: t !! Time of encounter + logical, dimension(:), allocatable :: lvdotr !! relative vdotr flag + integer(I4B), dimension(:), allocatable :: status !! status of the interaction + integer(I4B), dimension(:), allocatable :: index1 !! position of the first body in the encounter + integer(I4B), dimension(:), allocatable :: index2 !! position of the second body in the encounter + integer(I4B), dimension(:), allocatable :: id1 !! id of the first body in the encounter + integer(I4B), dimension(:), allocatable :: id2 !! id of the second body in the encounter + real(DP), dimension(:,:), allocatable :: x1 !! the position of body 1 in the encounter + real(DP), dimension(:,:), allocatable :: x2 !! the position of body 2 in the encounter + real(DP), dimension(:,:), allocatable :: v1 !! the velocity of body 1 in the encounter + real(DP), dimension(:,:), allocatable :: v2 !! the velocity of body 2 in the encounter + real(DP), dimension(:), allocatable :: Gmass1 !! G*mass of body 1 in the encounter + real(DP), dimension(:), allocatable :: Gmass2 !! G*mass of body 2 in the encounter + real(DP), dimension(:), allocatable :: radius1 !! radius of body 1 in the encounter + real(DP), dimension(:), allocatable :: radius2 !! radius of body 2 in the encounter + character(NAMELEN), dimension(:), allocatable :: name1 !! name body 1 in the encounter + character(NAMELEN), dimension(:), allocatable :: name2 !! name of body 2 in the encounter contains procedure :: setup => encounter_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 diff --git a/src/modules/symba_classes.f90 b/src/modules/symba_classes.f90 index f303d6143..139b221e0 100644 --- a/src/modules/symba_classes.f90 +++ b/src/modules/symba_classes.f90 @@ -142,7 +142,8 @@ module symba_classes !******************************************************************************************************************************* !> SyMBA class for tracking close encounters in a step type, extends(encounter_list) :: symba_encounter - integer(I4B), dimension(:), allocatable :: level !! encounter recursion level + integer(I4B), dimension(:), allocatable :: level !! encounter recursion level + real(DP), dimension(:), allocatable :: tcollision !! Time of collision contains procedure :: collision_check => symba_collision_check_encounter !! Checks if a test particle is going to collide with a massive body procedure :: encounter_check => symba_encounter_check !! Checks if massive bodies are going through close encounters with each other @@ -180,13 +181,13 @@ module symba_classes ! symba_nbody_system class definitions and method interfaces !******************************************************************************************************************************** type, extends(helio_nbody_system) :: symba_nbody_system - class(symba_merger), allocatable :: pl_adds !! List of added bodies in mergers or collisions - class(symba_pltpenc), allocatable :: pltpenc_list !! List of massive body-test particle encounters in a single step - class(symba_plplenc), allocatable :: plplenc_list !! List of massive body-massive body encounters in a single step - class(symba_plplenc), allocatable :: plplcollision_list !! List of massive body-massive body collisions in a single step - integer(I4B) :: irec !! System recursion level - type(encounter_storage(nframes=:)), allocatable :: encounter_history - integer(I4B) :: iframe = 0 !! Encounter history frame number + class(symba_merger), allocatable :: pl_adds !! List of added bodies in mergers or collisions + class(symba_pltpenc), allocatable :: pltpenc_list !! List of massive body-test particle encounters in a single step + class(symba_plplenc), allocatable :: plplenc_list !! List of massive body-massive body encounters in a single step + class(symba_plplenc), allocatable :: plplcollision_list !! List of massive body-massive body collisions in a single step + integer(I4B) :: irec !! System recursion level + type(encounter_storage(nframes=:)), allocatable :: encounter_history !! Stores encounter history for later retrieval and saving to file + integer(I4B) :: ienc_frame = 0 !! Encounter history frame number contains procedure :: write_discard => symba_io_write_discard !! Write out information about discarded and merged planets and test particles in SyMBA procedure :: initialize => symba_setup_initialize_system !! Performs SyMBA-specific initilization steps diff --git a/src/symba/symba_collision.f90 b/src/symba/symba_collision.f90 index e839af1de..04ec18b46 100644 --- a/src/symba/symba_collision.f90 +++ b/src/symba/symba_collision.f90 @@ -339,7 +339,7 @@ module function symba_collision_check_encounter(self, system, param, t, dt, irec i = self%index1(k) j = self%index2(k) if (lcollision(k)) self%status(k) = COLLISION - self%t(k) = t + self%tcollision(k) = t self%x1(:,k) = pl%xh(:,i) + system%cb%xb(:) self%v1(:,k) = pl%vb(:,i) if (isplpl) then diff --git a/src/symba/symba_encounter_check.f90 b/src/symba/symba_encounter_check.f90 index 8983dfa97..1e900a94a 100644 --- a/src/symba/symba_encounter_check.f90 +++ b/src/symba/symba_encounter_check.f90 @@ -34,7 +34,7 @@ module function symba_encounter_check_pl(self, param, system, dt, irec) result(l lany_encounter = .false. if (self%nbody == 0) return - associate(pl => self, plplenc_list => system%plplenc_list, cb => system%cb, iframe => system%iframe, encounter_history => system%encounter_history) + associate(pl => self, plplenc_list => system%plplenc_list, cb => system%cb, ienc_frame => system%ienc_frame, encounter_history => system%encounter_history) npl = pl%nbody nplm = pl%nplm @@ -69,6 +69,15 @@ module function symba_encounter_check_pl(self, param, system, dt, irec) result(l plplenc_list%x2(:,k) = pl%xh(:,j) plplenc_list%v1(:,k) = pl%vb(:,i) - cb%vb(:) plplenc_list%v2(:,k) = pl%vb(:,j) - cb%vb(:) + plplenc_list%Gmass1(k) = pl%Gmass(i) + plplenc_list%Gmass2(k) = pl%Gmass(j) + if (param%lclose) then + plplenc_list%radius1(k) = pl%radius(i) + plplenc_list%radius2(k) = pl%radius(j) + end if + plplenc_list%name1(k) = pl%info(i)%name + plplenc_list%name2(k) = pl%info(j)%name + pl%lencounter(i) = .true. pl%lencounter(j) = .true. pl%levelg(i) = irec @@ -78,8 +87,8 @@ module function symba_encounter_check_pl(self, param, system, dt, irec) result(l pl%nplenc(i) = pl%nplenc(i) + 1 pl%nplenc(j) = pl%nplenc(j) + 1 end do - iframe = iframe + 1 - encounter_history%frame(iframe) = plplenc_list + ienc_frame = ienc_frame + 1 + encounter_history%frame(ienc_frame) = plplenc_list end if end associate diff --git a/src/symba/symba_setup.f90 b/src/symba/symba_setup.f90 index 9187e4457..9a4ace98f 100644 --- a/src/symba/symba_setup.f90 +++ b/src/symba/symba_setup.f90 @@ -120,9 +120,12 @@ module subroutine symba_setup_encounter_list(self, n) if (n <= 0_I8B) return if (allocated(self%level)) deallocate(self%level) + if (allocated(self%tcollision)) deallocate(self%tcollision) allocate(self%level(n)) + allocate(self%tcollision(n)) self%level(:) = -1 + self%tcollision(:) = 0.0_DP return end subroutine symba_setup_encounter_list diff --git a/src/symba/symba_step.f90 b/src/symba/symba_step.f90 index d936ebcea..29011331d 100644 --- a/src/symba/symba_step.f90 +++ b/src/symba/symba_step.f90 @@ -275,7 +275,8 @@ module subroutine symba_step_reset_system(self, param) nenc_old = system%plplenc_list%nenc call system%plplenc_list%setup(0_I8B) call system%plplcollision_list%setup(0_I8B) - system%iframe = 0 + system%ienc_frame = 0 + if (allocated(system%encounter_history)) deallocate(system%encounter_history) if (npl > 0) then pl%lcollision(1:npl) = .false. call pl%reset_kinship([(i, i=1, npl)]) diff --git a/src/symba/symba_util.f90 b/src/symba/symba_util.f90 index 6f53d6bbd..805addee8 100644 --- a/src/symba/symba_util.f90 +++ b/src/symba/symba_util.f90 @@ -183,6 +183,7 @@ module subroutine symba_util_copy_encounter_list(self, source) class is (symba_encounter) associate(n => source%nenc) self%level(1:n) = source%level(1:n) + self%tcollision(1:n) = source%tcollision(1:n) end associate end select @@ -201,6 +202,7 @@ module subroutine symba_util_dealloc_encounter_list(self) class(symba_encounter), intent(inout) :: self !! SyMBA encounter list if (allocated(self%level)) deallocate(self%level) + if (allocated(self%tcollision)) deallocate(self%tcollision) return end subroutine symba_util_dealloc_encounter_list @@ -724,21 +726,33 @@ module subroutine symba_util_rearray_pl(self, system, param) ! This is an encounter we already know about, so save the old information system%plplenc_list%lvdotr(k) = plplenc_old%lvdotr(k) system%plplenc_list%status(k) = plplenc_old%status(k) + system%plplenc_list%Gmass1(k) = plplenc_old%Gmass1(k) + system%plplenc_list%Gmass2(k) = plplenc_old%Gmass2(k) + system%plplenc_list%radius1(k) = plplenc_old%radius1(k) + system%plplenc_list%radius2(k) = plplenc_old%radius2(k) + system%plplenc_list%name1(k) = plplenc_old%name1(k) + system%plplenc_list%name2(k) = plplenc_old%name2(k) system%plplenc_list%x1(:,k) = plplenc_old%x1(:,k) system%plplenc_list%x2(:,k) = plplenc_old%x2(:,k) system%plplenc_list%v1(:,k) = plplenc_old%v1(:,k) system%plplenc_list%v2(:,k) = plplenc_old%v2(:,k) - system%plplenc_list%t(k) = plplenc_old%t(k) + system%plplenc_list%tcollision(k) = plplenc_old%tcollision(k) system%plplenc_list%level(k) = plplenc_old%level(k) else if (((idnew1 == idold2) .and. (idnew2 == idold1))) then ! This is an encounter we already know about, but with the order reversed, so save the old information system%plplenc_list%lvdotr(k) = plplenc_old%lvdotr(k) system%plplenc_list%status(k) = plplenc_old%status(k) + system%plplenc_list%Gmass1(k) = plplenc_old%Gmass2(k) + system%plplenc_list%Gmass2(k) = plplenc_old%Gmass1(k) + system%plplenc_list%radius1(k) = plplenc_old%radius2(k) + system%plplenc_list%radius2(k) = plplenc_old%radius1(k) + system%plplenc_list%name1(k) = plplenc_old%name2(k) + system%plplenc_list%name2(k) = plplenc_old%name1(k) system%plplenc_list%x1(:,k) = plplenc_old%x2(:,k) system%plplenc_list%x2(:,k) = plplenc_old%x1(:,k) system%plplenc_list%v1(:,k) = plplenc_old%v2(:,k) system%plplenc_list%v2(:,k) = plplenc_old%v1(:,k) - system%plplenc_list%t(k) = plplenc_old%t(k) + system%plplenc_list%tcollision(k) = plplenc_old%tcollision(k) system%plplenc_list%level(k) = plplenc_old%level(k) end if system%plplenc_list%index1(k) = findloc(pl%id(1:npl), system%plplenc_list%id1(k), dim=1) @@ -761,7 +775,13 @@ module subroutine symba_util_rearray_pl(self, system, param) system%plplenc_list%id2(1:nencmin) = pack(system%plplenc_list%id2(1:nenc_old), lmask(1:nenc_old)) system%plplenc_list%lvdotr(1:nencmin) = pack(system%plplenc_list%lvdotr(1:nenc_old), lmask(1:nenc_old)) system%plplenc_list%status(1:nencmin) = pack(system%plplenc_list%status(1:nenc_old), lmask(1:nenc_old)) - system%plplenc_list%t(1:nencmin) = pack(system%plplenc_list%t(1:nenc_old), lmask(1:nenc_old)) + system%plplenc_list%Gmass1(1:nencmin) = pack(system%plplenc_list%Gmass1(1:nenc_old), lmask(1:nenc_old)) + system%plplenc_list%Gmass2(1:nencmin) = pack(system%plplenc_list%Gmass2(1:nenc_old), lmask(1:nenc_old)) + system%plplenc_list%radius1(1:nencmin) = pack(system%plplenc_list%radius1(1:nenc_old), lmask(1:nenc_old)) + system%plplenc_list%radius2(1:nencmin) = pack(system%plplenc_list%radius2(1:nenc_old), lmask(1:nenc_old)) + system%plplenc_list%name1(1:nencmin) = pack(system%plplenc_list%name1(1:nenc_old), lmask(1:nenc_old)) + system%plplenc_list%name2(1:nencmin) = pack(system%plplenc_list%name2(1:nenc_old), lmask(1:nenc_old)) + system%plplenc_list%tcollision(1:nencmin) = pack(system%plplenc_list%tcollision(1:nenc_old), lmask(1:nenc_old)) system%plplenc_list%level(1:nencmin) = pack(system%plplenc_list%level(1:nenc_old), lmask(1:nenc_old)) do i = 1, NDIM system%plplenc_list%x1(i, 1:nencmin) = pack(system%plplenc_list%x1(i, 1:nenc_old), lmask(1:nenc_old))