diff --git a/src/modules/symba_classes.f90 b/src/modules/symba_classes.f90 index b0daf9496..71df4b92e 100644 --- a/src/modules/symba_classes.f90 +++ b/src/modules/symba_classes.f90 @@ -95,6 +95,7 @@ module symba_classes procedure :: append => symba_util_append_pl !! Appends elements from one structure to another procedure :: fill => symba_util_fill_pl !! "Fills" bodies from one object into another depending on the results of a mask (uses the UNPACK intrinsic) procedure :: get_peri => symba_util_peri_pl !! Determine system pericenter passages for massive bodies + procedure :: rearray => symba_util_rearray_pl !! Clean up the massive body structures to remove discarded bodies and add new bodies procedure :: resize => symba_util_resize_pl !! Checks the current size of a SyMBA massive body against the requested size and resizes it if it is too small. procedure :: sort => symba_util_sort_pl !! Sorts body arrays by a sortable componen procedure :: rearrange => symba_util_sort_rearrange_pl !! Rearranges the order of array elements of body based on an input index array. Used in sorting methods @@ -533,6 +534,13 @@ module subroutine symba_util_peri_pl(self, system, param) class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters end subroutine symba_util_peri_pl + + module subroutine symba_util_rearray_pl(self, system, param) + implicit none + class(symba_pl), intent(inout) :: self !! SyMBA massive body object + class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object + class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters + end subroutine symba_util_rearray_pl end interface interface util_resize diff --git a/src/symba/symba_collision.f90 b/src/symba/symba_collision.f90 index 0f8b4d519..ad0e64079 100644 --- a/src/symba/symba_collision.f90 +++ b/src/symba/symba_collision.f90 @@ -208,8 +208,8 @@ function symba_collision_consolidate_familes(pl, param, idx_parent, family, x, v fam_size = 2 + sum(nchild(:)) allocate(family(fam_size)) family = [parent_child_index_array(1)%idx(:),parent_child_index_array(2)%idx(:)] - fam_size = count(pl%status(family(:)) == COLLISION) - family = pack(family(:), pl%status(family(:)) == COLLISION) + fam_size = count(pl%lcollision(family(:))) + family = pack(family(:), pl%lcollision(family(:))) L_spin(:,:) = 0.0_DP Ip(:,:) = 0.0_DP @@ -226,8 +226,8 @@ function symba_collision_consolidate_familes(pl, param, idx_parent, family, x, v if (nchild(j) > 0) then do i = 1, nchild(j) ! Loop over all children and take the mass weighted mean of the properties idx_child = parent_child_index_array(j)%idx(i + 1) - if ((idx_child) /= COLLISION) cycle - mchild = pl%Gmass(idx_child) + if (.not. pl%lcollision(idx_child)) cycle + mchild = pl%mass(idx_child) xchild(:) = pl%xb(:, idx_child) vchild(:) = pl%vb(:, idx_child) volchild = (4.0_DP / 3.0_DP) * PI * pl%radius(idx_child)**3 @@ -266,6 +266,7 @@ function symba_collision_consolidate_familes(pl, param, idx_parent, family, x, v return end function symba_collision_consolidate_familes + module subroutine symba_collision_encounter_scrub(self, system, param) !! author: David A. Minton !! diff --git a/src/symba/symba_discard.f90 b/src/symba/symba_discard.f90 index 1c5b4a732..f8ab11d04 100644 --- a/src/symba/symba_discard.f90 +++ b/src/symba/symba_discard.f90 @@ -58,6 +58,108 @@ subroutine symba_discard_cb_pl(pl, system, param) end subroutine symba_discard_cb_pl + subroutine symba_discard_conserve_mtm(pl, system, param, ipl, lescape_body) + !! author: David A. Minton + !! + !! Conserves system momentum when a body is lost from the system or collides with central body + implicit none + ! Arguments + class(symba_pl), intent(inout) :: pl + class(symba_nbody_system), intent(inout) :: system + class(symba_parameters), intent(inout) :: param + integer(I4B), intent(in) :: ipl + logical, intent(in) :: lescape_body + ! Internals + real(DP), dimension(NDIM) :: Lpl, Ltot, Lcb, xcom, vcom + real(DP) :: pe, ke_orbit, ke_spin + integer(I4B) :: i, oldstat + + select type(cb => system%cb) + class is (symba_cb) + + ! Add the potential and kinetic energy of the lost body to the records + pe = -cb%mass * pl%mass(ipl) / norm2(pl%xb(:, ipl) - pl%xb(:, 1)) + ke_orbit = 0.5_DP * pl%mass(ipl) * dot_product(pl%vb(:, ipl), pl%vb(:, ipl)) + if (param%lrotation) then + ke_spin = 0.5_DP * pl%mass(ipl) * pl%radius(ipl)**2 * pl%Ip(3, ipl) * dot_product(pl%rot(:, ipl), pl%rot(:, ipl)) + else + ke_spin = 0.0_DP + end if + + ! Add the pre-collision ke of the central body to the records + ! Add planet mass to central body accumulator + if (lescape_body) then + system%Mescape = system%Mescape + pl%mass(ipl) + do i = 1, npl + if (i == ipl) cycle + pe = pe - pl%mass(i) * pl%mass(ipl) / norm2(xb(:, ipl) - xb(:, i)) + end do + + Ltot(:) = 0.0_DP + do i = 1, npl + Lpl(:) = mass(i) * pl%xb(:,i) .cross. pl%vb(:, i) + Ltot(:) = Ltot(:) + Lpl(:) + end do + Ltot(:) = Ltot(:) + cb%mass * cb%xb(:) .cross. cb%vb(:) + call pl%b2h(cb) + oldstat = status(ipl) + pl%status(ipl) = INACTIVE + call pl%h2b(cb) + pl%status(ipl) = oldstat + do i = 1, npl + if (i == ipl) cycle + Lpl(:) = mass(i) * pl%xb(:,i) .cross. pl%vb(:, i) + Ltot(:) = Ltot(:) - Lpl(:) + end do + Ltot(:) = Ltot(:) - cb%mass * cb%xb(:) .cross. cb%vb(:) + system%Lescape(:) = system%Lescape(:) + system%Ltot(:) + if (param%lrotation) system%Lescape(:) = system%Lescape + pl%mass(ipl) * pl%radius(ipl)**2 * pl%Ip(3, ipl) * pl%rot(:, ipl) + + else + xcom(:) = (pl%mass(ipl) * pl%xb(:, ipl) + cb%mass * cb%xb(:)) / (cb%mass + pl%mass(ipl)) + vcom(:) = (pl%mass(ipl) * pl%vb(:, ipl) + cb%mass * cb%vb(:)) / (cb%mass + pl%mass(ipl)) + Lpl(:) = (pl%xb(:,ipl) - xcom(:)) .cross. pL%vb(:,ipl) - vcom(:) + if (param%lrotation) Lpl(:) = pl%mass(ipl) * (Lpl(:) + pl%radius(ipl)**2 * pl%Ip(3,ipl) * pl%rot(:, ipl)) + + Lcb(:) = cb%mass * (cb%xb(:) - xcom(:)) .cross. (cb%vb(:) - vcom(:)) + + ke_orbit = ke_orbit + 0.5_DP * cb%mass * dot_product(cb%vb(:), cb%vb(:)) + if (param%lrotation) ke_spin = ke_spin + 0.5_DP * cb%mass * cb%radius**2 * cb%Ip(3) * dot_product(cb%rot(:), cb%rot(:)) + ! Update mass of central body to be consistent with its total mass + cb%dM = cb%dM + pl%mass(ipl) + cb%dR = cb%dR + 1.0_DP / 3.0_DP * (pl%radius(ipl) / cb%radius)**3 - 2.0_DP / 9.0_DP * (pl%radius(ipl) / cb%radius)**6 + cb%mass = cb%M0 + cb%dM + cb%Gmass = param%GU * cb%mass + cb%radius = cb%R0 + cb%dR + param%rmin = cb%radius + ! Add planet angular momentum to central body accumulator + cb%dL(:) = Lpl(:) + Lcb(:) + cb%dL(:) + ! Update rotation of central body to by consistent with its angular momentum + if (param%lrotation) then + cb%rot(:) = (cb%L0(:) + cb%dL(:)) / (cb%Ip(3) * cb%mass * cb%radius**2) + ke_spin = ke_spin - 0.5_DP * cb%mass * cb%radius**2 * cb%Ip(3) * dot_product(cb%rot(:), cb%rot(:)) + end if + cb%xb(:) = xcom(:) + cb%vb(:) = vcom(:) + ke_orbit = ke_orbit - 0.5_DP * cb%mass * dot_product(cb%vb(:), cb%vb(:)) + end if + call pl%b2h(cb) + + ! We must do this for proper book-keeping, since we can no longer track this body's contribution to energy directly + if (lescape_body) then + system%Ecollisions = system%Ecollisions + ke_orbit + ke_spin + pe + system%Euntracked = system%Euntracked - (ke_orbit + ke_spin + pe) + else + system%Ecollisions = system%Ecollisions + pe + system%Euntracked = system%Euntracked - pe + end if + + end select + return + + end subroutine symba_discard_conserve_mtm + + subroutine symba_discard_nonplpl(pl, system, param) !! author: David A. Minton !! @@ -89,6 +191,45 @@ subroutine symba_discard_nonplpl(pl, system, param) end subroutine symba_discard_nonplpl + subroutine symba_discard_nonplpl_conservation(pl, system, param) + !! author: David A. Minton + !! + !! If there are any bodies that are removed due to either colliding with the central body or escaping the systme, + !! we need to track the conserved quantities with the system bookkeeping terms. + implicit none + ! Arguments + class(symba_pl), intent(inout) :: pl !! SyMBA test particle object + class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object + class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters + ! Internals + integer(I4B) :: i, ndiscard, dstat + logical :: lescape + logical, dimension(pl%nbody) :: discard_l_pl + integer(I4B), dimension(:), allocatable :: discard_index_list + + associate(npl => pl%nbody) + discard_l_pl(1:npl) = pl%ldiscard(1:npl) .and. .not. pl%lcollision(1:npl) ! These are bodies that are discarded but not flagged as pl-pl collision + ndiscard = count(discard_l_pl(:)) + allocate(discard_index_list(ndiscard)) + discard_index_list(:) = pack([(i, i = 1, npl)], discard_l_pl(1:npl)) + do i = 1, ndiscard + dstat = pl%status(discard_index_list(i)) + if ((dstat == DISCARDED_RMIN) .or. (dstat == DISCARDED_PERI)) then + lescape = .false. + else if ((dstat == DISCARDED_RMAX) .or. (dstat == DISCARDED_RMAXU)) then + lescape = .true. + else + cycle + end if + ! Conserve all the quantities + call symba_discard_conserve_mtm(pl, system, param, discard_index_list(i), lescape) + end do + end associate + + return + end subroutine symba_discard_nonplpl_conservation + + subroutine symba_discard_peri_pl(pl, system, param) !! author: David A. Minton !! @@ -150,17 +291,28 @@ module subroutine symba_discard_pl(self, system, param) select type(param) class is (symba_parameters) associate(pl => self, plplenc_list => system%plplenc_list) + call pl%h2b(system%cb) + + ! First deal with the non pl-pl collisions call symba_discard_nonplpl(self, system, param) + + ! Scrub the pl-pl encounter list of any encounters that did not lead to a collision call plplenc_list%scrub_non_collision(system, param) - if (plplenc_list%nenc == 0) return ! No collisions to resolve - write(*, *) "Collision detected at time t = ",param%t - call pl%h2b(system%cb) - if (param%lfragmentation) then - call plplenc_list%resolve_fragmentations(system, param) - else - call plplenc_list%resolve_mergers(system, param) + if ((plplenc_list%nenc > 0) .and. any(pl%lcollision(:))) then + write(*, *) "Collision between massive bodies detected at time t = ",param%t + if (param%lfragmentation) then + call plplenc_list%resolve_fragmentations(system, param) + else + call plplenc_list%resolve_mergers(system, param) + end if end if + + if (any(pl%ldiscard(:))) then + call symba_discard_nonplpl_conservation(self, system, param) + call pl%rearray(self, system, param) + end if + end associate end select end select diff --git a/src/symba/symba_fragmentation.f90 b/src/symba/symba_fragmentation.f90 index ccc79892f..cafaacfd7 100644 --- a/src/symba/symba_fragmentation.f90 +++ b/src/symba/symba_fragmentation.f90 @@ -73,7 +73,7 @@ module function symba_fragmentation_casemerge(system, param, family, x, v, mass, pe = 0.0_DP do j = 1, nfamily do i = j + 1, nfamily - pe = pe - pl%Gmass(i) * pl%Gmass(j) / norm2(pl%xb(:, i) - pl%xb(:, j)) + pe = pe - pl%mass(i) * pl%mass(j) / norm2(pl%xb(:, i) - pl%xb(:, j)) end do end do system%Ecollisions = system%Ecollisions + pe diff --git a/src/symba/symba_util.f90 b/src/symba/symba_util.f90 index 7063e0910..c0276291a 100644 --- a/src/symba/symba_util.f90 +++ b/src/symba/symba_util.f90 @@ -369,6 +369,23 @@ module subroutine symba_util_peri_pl(self, system, param) end subroutine symba_util_peri_pl + module subroutine symba_util_rearray_pl(self, system, param) + !! Author: the Purdue Swiftest Team - David A. Minton, Carlisle A. Wishard, Jennifer L.L. Pouplin, and Jacob R. Elliott + !! + !! Clean up the massive body structures to remove discarded bodies and add new bodies + implicit none + ! Arguments + class(symba_pl), intent(inout) :: self !! SyMBA massive body object + class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object + class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters + + ! First + + return + end subroutine symba_util_rearray_pl + + + module subroutine symba_util_resize_arr_info(arr, nnew) !! author: David A. Minton !!