diff --git a/src/modules/swiftest_classes.f90 b/src/modules/swiftest_classes.f90 index b55dfb5e2..fdb2360ee 100644 --- a/src/modules/swiftest_classes.f90 +++ b/src/modules/swiftest_classes.f90 @@ -292,8 +292,14 @@ module swiftest_classes 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 + 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 contains - procedure :: copy => util_copy_encounter + procedure :: copy => util_copy_encounter + procedure :: resize => util_resize_encounter + procedure :: setup => setup_encounter end type swiftest_encounter abstract interface @@ -707,6 +713,12 @@ module subroutine setup_construct_system(system, param) class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters end subroutine setup_construct_system + module subroutine setup_encounter(self, n) + implicit none + class(swiftest_encounter), intent(inout) :: self !! Swiftest encounter structure + integer(I4B), intent(in) :: n !! Number of encounters to allocate space for + end subroutine setup_encounter + module subroutine setup_initialize_system(self, param) implicit none class(swiftest_nbody_system), intent(inout) :: self !! Swiftest system object @@ -952,6 +964,12 @@ module subroutine util_resize_body(self, nnew) integer(I4B), intent(in) :: nnew !! New size neded end subroutine util_resize_body + module subroutine util_resize_encounter(self, nnew) + implicit none + class(swiftest_encounter), intent(inout) :: self !! Swiftest encounter list + integer(I4B), intent(in) :: nnew !! New size of list needed + end subroutine util_resize_encounter + module subroutine util_resize_pl(self, nnew) implicit none class(swiftest_pl), intent(inout) :: self !! Swiftest massive body object diff --git a/src/modules/symba_classes.f90 b/src/modules/symba_classes.f90 index a8c6ff716..92caa0bb3 100644 --- a/src/modules/symba_classes.f90 +++ b/src/modules/symba_classes.f90 @@ -133,7 +133,6 @@ module symba_classes procedure :: kick => symba_kick_pltpenc !! Kick barycentric velocities of active test particles within SyMBA recursion procedure :: setup => symba_setup_pltpenc !! A constructor that sets the number of encounters and allocates and initializes all arrays procedure :: copy => symba_util_copy_pltpenc !! Copies all elements of one pltpenc list to another - procedure :: resize => symba_util_resize_pltpenc !! Checks the current size of the pltpenc_list against the required size and extends it by a factor of 2 more than requested if it is too small end type symba_pltpenc !******************************************************************************************************************************** @@ -141,14 +140,7 @@ module symba_classes !******************************************************************************************************************************* !> SyMBA class for tracking pl-pl close encounters in a step type, extends(symba_pltpenc) :: symba_plplenc - real(DP), dimension(:,:), allocatable :: xh1 !! the heliocentric position of parent 1 in encounter - real(DP), dimension(:,:), allocatable :: xh2 !! the heliocentric position of parent 2 in encounter - real(DP), dimension(:,:), allocatable :: vb1 !! the barycentric velocity of parent 1 in encounter - real(DP), dimension(:,:), allocatable :: vb2 !! the barycentric velocity of parent 2 in encounter contains - procedure :: collision_check => symba_collision_check_plplenc !! Checks if two massive bodies are going to collide - procedure :: setup => symba_setup_plplenc !! A constructor that sets the number of encounters and allocates and initializes all arrays - procedure :: copy => symba_util_copy_plplenc !! Copies all elements of one plplenc list to another end type symba_plplenc !******************************************************************************************************************************** @@ -182,17 +174,6 @@ module subroutine symba_collision_check_pltpenc(self, system, param, t, dt, irec integer(I4B), intent(in) :: irec !! Current recursion level end subroutine symba_collision_check_pltpenc - module subroutine symba_collision_check_plplenc(self, system, param, t, dt, irec) - use swiftest_classes, only : swiftest_parameters - implicit none - class(symba_plplenc), intent(inout) :: self !! SyMBA pl-tp encounter list object - class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters - real(DP), intent(in) :: t !! current time - real(DP), intent(in) :: dt !! step size - integer(I4B), intent(in) :: irec !! Current recursion level - end subroutine symba_collision_check_plplenc - module subroutine symba_collision_make_family_pl(self,idx) implicit none class(symba_pl), intent(inout) :: self !! SyMBA massive body object @@ -343,7 +324,6 @@ module subroutine symba_kick_pltpenc(self, system, dt, irec, sgn) integer(I4B), intent(in) :: sgn !! sign to be applied to acceleration end subroutine symba_kick_pltpenc - module subroutine symba_io_write_frame_info(self, iu, param) use swiftest_classes, only : swiftest_parameters implicit none @@ -465,13 +445,6 @@ module subroutine symba_util_copy_pltpenc(self, source) class(symba_pltpenc), intent(inout) :: self !! SyMBA pl-tp encounter list class(swiftest_encounter), intent(in) :: source !! Source object to copy into end subroutine symba_util_copy_pltpenc - - module subroutine symba_util_copy_plplenc(self, source) - use swiftest_classes, only : swifest_encounter - implicit none - class(symba_plplenc), intent(inout) :: self !! SyMBA pl-pl encounter list - class(swiftest_encounter), intent(in) :: source !! Source object to copy into - end subroutine symba_util_copy_plplenc end interface interface util_fill @@ -529,12 +502,6 @@ module subroutine symba_util_resize_pl(self, nnew) integer(I4B), intent(in) :: nnew !! New size neded end subroutine symba_util_resize_pl - module subroutine symba_util_resize_pltpenc(self, nnew) - implicit none - class(symba_pltpenc), intent(inout) :: self !! SyMBA pl-tp encounter list - integer(I4B), intent(in) :: nnew !! New size of list needed - end subroutine symba_util_resize_pltpenc - module subroutine symba_util_resize_tp(self, nnew) implicit none class(symba_tp), intent(inout) :: self !! SyMBA massive body object diff --git a/src/setup/setup.f90 b/src/setup/setup.f90 index 7de8ce5c3..ca5f38c6e 100644 --- a/src/setup/setup.f90 +++ b/src/setup/setup.f90 @@ -70,6 +70,50 @@ module subroutine setup_construct_system(system, param) end subroutine setup_construct_system + module subroutine setup_encounter(self, n) + !! author: David A. Minton + !! + !! A constructor that sets the number of encounters and allocates and initializes all arrays + !! + implicit none + ! Arguments + class(swiftest_encounter), intent(inout) :: self !! Swiftest encounter structure + integer(I4B), intent(in) :: n !! Number of encounters to allocate space for + + self%nenc = n + if (n == 0) return + + if (allocated(self%lvdotr)) deallocate(self%lvdotr) + if (allocated(self%status)) deallocate(self%status) + if (allocated(self%index1)) deallocate(self%index1) + if (allocated(self%index2)) deallocate(self%index2) + if (allocated(self%x1)) deallocate(self%x1) + if (allocated(self%x2)) deallocate(self%x2) + if (allocated(self%v1)) deallocate(self%v1) + if (allocated(self%v2)) deallocate(self%v2) + + allocate(self%lvdotr(n)) + allocate(self%status(n)) + allocate(self%index1(n)) + allocate(self%index2(n)) + allocate(self%x1(NDIM,n)) + allocate(self%x2(NDIM,n)) + allocate(self%v1(NDIM,n)) + allocate(self%v2(NDIM,n)) + + self%lvdotr(:) = .false. + self%status(:) = INACTIVE + self%index1(:) = 0 + self%index2(:) = 0 + self%x1(:,:) = 0.0_DP + self%x2(:,:) = 0.0_DP + self%v1(:,:) = 0.0_DP + self%v2(:,:) = 0.0_DP + + return + end subroutine setup_encounter + + module subroutine setup_initialize_system(self, param) !! author: David A. Minton !! diff --git a/src/symba/symba_collision.f90 b/src/symba/symba_collision.f90 index 2614a486a..bb67508fd 100644 --- a/src/symba/symba_collision.f90 +++ b/src/symba/symba_collision.f90 @@ -2,89 +2,12 @@ use swiftest contains - module subroutine symba_collision_check_plplenc(self, system, param, t, dt, irec) - !! author: Jennifer L.L. Pouplin, Carlisle A. wishard, and David A. Minton - !! - !! Check for merger between massive bodies in SyMBA. If the user has turned on the FRAGMENTATION feature, it will call the - !! symba_regime subroutine to determine what kind of collision will occur. - !! - !! Adapted from David E. Kaufmann's Swifter routine symba_merge_pl.f90 - !! - !! Adapted from Hal Levison's Swift routine symba5_merge.f - implicit none - ! Arguments - class(symba_plplenc), intent(inout) :: self !! SyMBA pl-tp encounter list object - class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters - real(DP), intent(in) :: t !! current time - real(DP), intent(in) :: dt !! step size - integer(I4B), intent(in) :: irec !! Current recursion level - ! Internals - logical, dimension(:), allocatable :: lcollision, lmask - real(DP), dimension(NDIM) :: xr, vr - integer(I4B) :: k - real(DP) :: rlim, mtot - - if (self%nenc == 0) return - - select type(pl => system%pl) - class is (symba_pl) - associate(plplenc_list => self, nplplenc => self%nenc, ind1 => self%index1, ind2 => self%index2) - allocate(lmask(nplplenc)) - lmask(:) = ((plplenc_list%status(1:nplplenc) == ACTIVE) & - .and. (pl%levelg(ind1(1:nplplenc)) >= irec) & - .and. (pl%levelg(ind2(1:nplplenc)) >= irec)) - if (.not.any(lmask(:))) return - - allocate(lcollision(nplplenc)) - lcollision(:) = .false. - - do concurrent(k = 1:nplplenc, lmask(k)) - xr(:) = pl%xh(:, ind1(k)) - pl%xh(:, ind2(k)) - vr(:) = pl%vb(:, ind1(k)) - pl%vb(:, ind2(k)) - rlim = pl%radius(ind1(k)) + pl%radius(ind2(k)) - mtot = pl%Gmass(ind1(k)) + pl%Gmass(ind2(k)) - lcollision(k) = symba_collision_check_one(xr(1), xr(2), xr(3), vr(1), vr(2), vr(3), mtot, rlim, dt, plplenc_list%lvdotr(k)) - end do - deallocate(lmask) - - if (any(lcollision(:))) then - allocate(lmask(pl%nbody)) - do k = 1, nplplenc - if (.not.lcollision(k)) cycle - - ! Set this encounter as a collision and save the position and velocity vectors at the time of the collision - plplenc_list%status(k) = COLLISION - plplenc_list%xh1(:,k) = pl%xh(:,ind1(k)) - plplenc_list%vb1(:,k) = pl%vb(:,ind1(k)) - plplenc_list%xh2(:,k) = pl%xh(:,ind2(k)) - plplenc_list%vb2(:,k) = pl%vb(:,ind2(k)) - - ! Check to see if either of these bodies has been involved with a collision before, and if so, make this a collisional family - if (pl%lcollision(ind1(k)) .or. pl%lcollision(ind2(k))) call pl%make_family([ind1(k),ind2(k)]) - - ! Set the collision flag for these to bodies to true in case they become involved in another collision later in the step - pl%lcollision([ind1(k), ind2(k)]) = .true. - pl%lcollision([ind1(k), ind2(k)]) = .true. - pl%status([ind1(k), ind2(k)]) = COLLISION - end do - - end if - end associate - end select - - return - - return - end subroutine symba_collision_check_plplenc - - module subroutine symba_collision_check_pltpenc(self, system, param, t, dt, irec) !! author: David A. Minton !! !! Check for merger between massive bodies and test particles in SyMBA !! - !! Adapted from David E. Kaufmann's Swifter routine symba_merge_tp.f90 + !! Adapted from David E. Kaufmann's Swifter routine symba_merge.f90 and symba_merge_tp.f90 !! !! Adapted from Hal Levison's Swift routine symba5_merge.f implicit none @@ -99,39 +22,74 @@ module subroutine symba_collision_check_pltpenc(self, system, param, t, dt, irec logical, dimension(:), allocatable :: lcollision, lmask real(DP), dimension(NDIM) :: xr, vr integer(I4B) :: k + real(DP) :: rlim, mtot + logical :: isplpl if (self%nenc == 0) return + select type(self) + class is (symba_plplenc) + isplpl = .true. + class default + isplpl = .false. + end select select type(pl => system%pl) class is (symba_pl) select type(tp => system%tp) class is (symba_tp) - associate(pltpenc_list => self, npltpenc => self%nenc, plind => self%index1, tpind => self%index2) - allocate(lmask(npltpenc)) - lmask(:) = ((pltpenc_list%status(1:npltpenc) == ACTIVE) & - .and. (pl%levelg(plind(1:npltpenc)) >= irec) & - .and. (tp%levelg(tpind(1:npltpenc)) >= irec)) + associate(nenc => self%nenc, ind1 => self%index1, ind2 => self%index2) + allocate(lmask(nenc)) + lmask(:) = ((self%status(1:nenc) == ACTIVE) .and. (pl%levelg(ind1(1:nenc)) >= irec)) + if (isplpl) then + lmask(:) = lmask(:) .and. (pl%levelg(ind2(1:nenc)) >= irec) + else + lmask(:) = lmask(:) .and. (tp%levelg(ind2(1:nenc)) >= irec) + end if if (.not.any(lmask(:))) return - allocate(lcollision(npltpenc)) + allocate(lcollision(nenc)) lcollision(:) = .false. - do concurrent(k = 1:npltpenc, lmask(k)) - xr(:) = pl%xh(:, plind(k)) - tp%xh(:, tpind(k)) - vr(:) = pl%vb(:, plind(k)) - tp%vb(:, tpind(k)) - lcollision(k) = symba_collision_check_one(xr(1), xr(2), xr(3), vr(1), vr(2), vr(3), pl%Gmass(plind(k)), pl%radius(plind(k)), dt, pltpenc_list%lvdotr(k)) - end do + if (isplpl) then + do concurrent(k = 1:nenc, lmask(k)) + xr(:) = pl%xh(:, ind1(k)) - pl%xh(:, ind2(k)) + vr(:) = pl%vb(:, ind1(k)) - pl%vb(:, ind2(k)) + rlim = pl%radius(ind1(k)) + pl%radius(ind2(k)) + mtot = pl%Gmass(ind1(k)) + pl%Gmass(ind2(k)) + lcollision(k) = symba_collision_check_one(xr(1), xr(2), xr(3), vr(1), vr(2), vr(3), mtot, rlim, dt, self%lvdotr(k)) + end do + else + do concurrent(k = 1:nenc, lmask(k)) + xr(:) = pl%xh(:, ind1(k)) - tp%xh(:, ind2(k)) + vr(:) = pl%vb(:, ind1(k)) - tp%vb(:, ind2(k)) + lcollision(k) = symba_collision_check_one(xr(1), xr(2), xr(3), vr(1), vr(2), vr(3), pl%Gmass(ind1(k)), pl%radius(ind1(k)), dt, self%lvdotr(k)) + end do + end if if (any(lcollision(:))) then - where(lcollision(1:npltpenc)) - pltpenc_list%status(1:npltpenc) = COLLISION - tp%status(tpind(1:npltpenc)) = DISCARDED_PLR - tp%ldiscard(tpind(1:npltpenc)) = .true. - end where - - do k = 1, npltpenc - if (pltpenc_list%status(k) /= COLLISION) cycle - write(*,*) 'Test particle ',tp%id(tpind(k)), ' collided with massive body ',pl%id(plind(k)), ' at time ',t + do k = 1, nenc + if (.not.lcollision(k)) cycle + self%status(k) = COLLISION + self%x1(:,k) = pl%xh(:,ind1(k)) + self%v1(:,k) = pl%vb(:,ind1(k)) + if (isplpl) then + self%x2(:,k) = pl%xh(:,ind2(k)) + self%v2(:,k) = pl%vb(:,ind2(k)) + + ! Check to see if either of these bodies has been involved with a collision before, and if so, make this a collisional family + if (pl%lcollision(ind1(k)) .or. pl%lcollision(ind2(k))) call pl%make_family([ind1(k),ind2(k)]) + + ! Set the collision flag for these to bodies to true in case they become involved in another collision later in the step + pl%lcollision([ind1(k), ind2(k)]) = .true. + pl%ldiscard([ind1(k), ind2(k)]) = .true. + pl%status([ind1(k), ind2(k)]) = COLLISION + else + self%x2(:,k) = tp%xh(:,ind2(k)) + self%v2(:,k) = tp%vb(:,ind2(k)) + tp%status(ind2(k)) = DISCARDED_PLR + tp%ldiscard(ind2(k)) = .true. + write(*,*) 'Test particle ',tp%id(ind2(k)), ' collided with massive body ',pl%id(ind1(k)), ' at time ',t + end if end do end if end associate diff --git a/src/symba/symba_setup.f90 b/src/symba/symba_setup.f90 index f2c8e63dd..524420609 100644 --- a/src/symba/symba_setup.f90 +++ b/src/symba/symba_setup.f90 @@ -105,63 +105,18 @@ module subroutine symba_setup_pltpenc(self, n) class(symba_pltpenc), intent(inout) :: self !! SyMBA pl-tp encounter structure integer(I4B), intent(in) :: n !! Number of encounters to allocate space for - self%nenc = n + call setup_encounter(self, n) if (n == 0) return - if (allocated(self%lvdotr)) deallocate(self%lvdotr) - if (allocated(self%status)) deallocate(self%status) if (allocated(self%level)) deallocate(self%level) - if (allocated(self%index1)) deallocate(self%index1) - if (allocated(self%index2)) deallocate(self%index2) - - allocate(self%lvdotr(n)) - allocate(self%status(n)) allocate(self%level(n)) - allocate(self%index1(n)) - allocate(self%index2(n)) - self%lvdotr(:) = .false. - self%status(:) = INACTIVE self%level(:) = -1 - self%index1(:) = 0 - self%index2(:) = 0 return end subroutine symba_setup_pltpenc - module subroutine symba_setup_plplenc(self, n) - !! author: David A. Minton - !! - !! A constructor that sets the number of encounters and allocates and initializes all arrays - ! - implicit none - ! Arguments - class(symba_plplenc), intent(inout) :: self !! SyMBA pl-tp encounter structure - integer(I4B), intent(in) :: n !! Number of encounters to allocate space for - - call symba_setup_pltpenc(self, n) - if (n == 0) return - - if (allocated(self%xh1)) deallocate(self%xh1) - if (allocated(self%xh2)) deallocate(self%xh2) - if (allocated(self%vb1)) deallocate(self%vb1) - if (allocated(self%vb2)) deallocate(self%vb2) - - allocate(self%xh1(NDIM,n)) - allocate(self%xh2(NDIM,n)) - allocate(self%vb1(NDIM,n)) - allocate(self%vb2(NDIM,n)) - - self%xh1(:,:) = 0.0_DP - self%xh2(:,:) = 0.0_DP - self%vb1(:,:) = 0.0_DP - self%vb2(:,:) = 0.0_DP - - return - end subroutine symba_setup_plplenc - - module subroutine symba_setup_tp(self, n, param) !! author: David A. Minton !! diff --git a/src/symba/symba_util.f90 b/src/symba/symba_util.f90 index e76c2faff..feda27a07 100644 --- a/src/symba/symba_util.f90 +++ b/src/symba/symba_util.f90 @@ -162,31 +162,6 @@ module subroutine symba_util_copy_pltpenc(self, source) return end subroutine symba_util_copy_pltpenc - - module subroutine symba_util_copy_plplenc(self, source) - !! author: David A. Minton - !! - !! Copies elements from the source encounter list into self. - implicit none - ! Arguments - class(symba_plplenc), intent(inout) :: self !! SyMBA pl-pl encounter list - class(swiftest_encounter), intent(in) :: source !! Source object to copy into - - call symba_util_copy_pltpenc(self, source) - associate(n => source%nenc) - select type(source) - class is (symba_plplenc) - self%xh1(:,1:n) = source%xh1(:,1:n) - self%xh2(:,1:n) = source%xh2(:,1:n) - self%vb1(:,1:n) = source%vb1(:,1:n) - self%vb2(:,1:n) = source%vb2(:,1:n) - end select - end associate - - return - end subroutine symba_util_copy_plplenc - - module subroutine symba_util_fill_arr_info(keeps, inserts, lfill_list) !! author: David A. Minton !! @@ -410,43 +385,6 @@ module subroutine symba_util_resize_tp(self, nnew) return end subroutine symba_util_resize_tp - - module subroutine symba_util_resize_pltpenc(self, nnew) - !! author: David A. Minton - !! - !! Checks the current size of the encounter list against the required size and extends it by a factor of 2 more than requested if it is too small. - !! Polymorphic method works on both symba_pltpenc and symba_plplenc types - implicit none - ! Arguments - class(symba_pltpenc), intent(inout) :: self !! SyMBA pl-tp encounter list - integer(I4B), intent(in) :: nnew !! New size of list needed - ! Internals - class(symba_pltpenc), allocatable :: enc_temp - integer(I4B) :: nold - logical :: lmalloc - - lmalloc = allocated(self%status) - if (lmalloc) then - nold = size(self%status) - else - nold = 0 - end if - if (nnew > nold) then - if (lmalloc) allocate(enc_temp, source=self) - call self%setup(2 * nnew) - if (lmalloc) then - call self%copy(enc_temp) - deallocate(enc_temp) - end if - else - self%status(nnew+1:nold) = INACTIVE - end if - self%nenc = nnew - - return - end subroutine symba_util_resize_pltpenc - - module subroutine symba_util_sort_pl(self, sortby, ascending) !! author: David A. Minton !! diff --git a/src/util/util_copy.f90 b/src/util/util_copy.f90 index 704e10d64..f44777eec 100644 --- a/src/util/util_copy.f90 +++ b/src/util/util_copy.f90 @@ -17,6 +17,10 @@ module subroutine util_copy_encounter(self, source) self%status(1:n) = source%status(1:n) self%index1(1:n) = source%index1(1:n) self%index2(1:n) = source%index2(1:n) + self%x1(:,1:n) = source%x1(:,1:n) + self%x2(:,1:n) = source%x2(:,1:n) + self%v1(:,1:n) = source%v1(:,1:n) + self%v2(:,1:n) = source%v2(:,1:n) end associate return diff --git a/src/util/util_resize.f90 b/src/util/util_resize.f90 index 9abd66189..3772a8207 100644 --- a/src/util/util_resize.f90 +++ b/src/util/util_resize.f90 @@ -207,6 +207,41 @@ module subroutine util_resize_body(self, nnew) end subroutine util_resize_body + module subroutine util_resize_encounter(self, nnew) + !! author: David A. Minton + !! + !! Checks the current size of the encounter list against the required size and extends it by a factor of 2 more than requested if it is too small. + implicit none + ! Arguments + class(swiftest_encounter), intent(inout) :: self !! Swiftest encounter list + integer(I4B), intent(in) :: nnew !! New size of list needed + ! Internals + class(swiftest_encounter), allocatable :: enc_temp + integer(I4B) :: nold + logical :: lmalloc + + lmalloc = allocated(self%status) + if (lmalloc) then + nold = size(self%status) + else + nold = 0 + end if + if (nnew > nold) then + if (lmalloc) allocate(enc_temp, source=self) + call self%setup(2 * nnew) + if (lmalloc) then + call self%copy(enc_temp) + deallocate(enc_temp) + end if + else + self%status(nnew+1:nold) = INACTIVE + end if + self%nenc = nnew + + return + end subroutine util_resize_encounter + + module subroutine util_resize_pl(self, nnew) !! author: David A. Minton !!