diff --git a/src/discard/discard.f90 b/src/discard/discard.f90 index 6c799844f..55ad97f65 100644 --- a/src/discard/discard.f90 +++ b/src/discard/discard.f90 @@ -135,22 +135,22 @@ subroutine discard_cb_tp(tp, system, param) if ((param%rmax >= 0.0_DP) .and. (rh2 > rmax2)) then tp%status(i) = DISCARDED_RMAX write(idstr, *) tp%id(i) - write(timestr, *) param%t + write(timestr, *) system%t write(*, *) "Particle " // trim(adjustl(tp%info(i)%name)) // " (" // trim(adjustl(idstr)) // ")" // & " too far from the central body at t = " // trim(adjustl(timestr)) tp%ldiscard(i) = .true. tp%lmask(i) = .false. - call tp%info(i)%set_value(status="DISCARDED_RMAX", discard_time=param%t, discard_xh=tp%xh(:,i), & + call tp%info(i)%set_value(status="DISCARDED_RMAX", discard_time=system%t, discard_xh=tp%xh(:,i), & discard_vh=tp%vh(:,i)) else if ((param%rmin >= 0.0_DP) .and. (rh2 < rmin2)) then tp%status(i) = DISCARDED_RMIN write(idstr, *) tp%id(i) - write(timestr, *) param%t + write(timestr, *) system%t write(*, *) "Particle " // trim(adjustl(tp%info(i)%name)) // " (" // trim(adjustl(idstr)) // ")" // & " too close to the central body at t = " // trim(adjustl(timestr)) tp%ldiscard(i) = .true. tp%lmask(i) = .false. - call tp%info(i)%set_value(status="DISCARDED_RMIN", discard_time=param%t, discard_xh=tp%xh(:,i), & + call tp%info(i)%set_value(status="DISCARDED_RMIN", discard_time=system%t, discard_xh=tp%xh(:,i), & discard_vh=tp%vh(:,i), discard_body_id=cb%id) else if (param%rmaxu >= 0.0_DP) then rb2 = dot_product(tp%xb(:, i), tp%xb(:, i)) @@ -159,12 +159,12 @@ subroutine discard_cb_tp(tp, system, param) if ((energy > 0.0_DP) .and. (rb2 > rmaxu2)) then tp%status(i) = DISCARDED_RMAXU write(idstr, *) tp%id(i) - write(timestr, *) param%t + write(timestr, *) system%t write(*, *) "Particle " // trim(adjustl(tp%info(i)%name)) // " (" // trim(adjustl(idstr)) // ")" // & " is unbound and too far from barycenter at t = " // trim(adjustl(timestr)) tp%ldiscard(i) = .true. tp%lmask(i) = .false. - call tp%info(i)%set_value(status="DISCARDED_RMAXU", discard_time=param%t, discard_xh=tp%xh(:,i), & + call tp%info(i)%set_value(status="DISCARDED_RMAXU", discard_time=system%t, discard_xh=tp%xh(:,i), & discard_vh=tp%vh(:,i)) end if end if @@ -194,7 +194,7 @@ subroutine discard_peri_tp(tp, system, param) real(DP), dimension(NDIM) :: dx character(len=STRMAX) :: idstr, timestr - associate(cb => system%cb, ntp => tp%nbody, pl => system%pl, npl => system%pl%nbody, t => param%t) + associate(cb => system%cb, ntp => tp%nbody, pl => system%pl, npl => system%pl%nbody, t => system%t) call tp%get_peri(system, param) do i = 1, ntp if (tp%status(i) == ACTIVE) then @@ -211,11 +211,11 @@ subroutine discard_peri_tp(tp, system, param) (tp%peri(i) <= param%qmin)) then tp%status(i) = DISCARDED_PERI write(idstr, *) tp%id(i) - write(timestr, *) param%t + write(timestr, *) system%t write(*, *) "Particle " // trim(adjustl(tp%info(i)%name)) // " (" // trim(adjustl(idstr)) // ")" // & " perihelion distance too small at t = " // trim(adjustl(timestr)) tp%ldiscard(i) = .true. - call tp%info(i)%set_value(status="DISCARDED_PERI", discard_time=param%t, discard_xh=tp%xh(:,i), & + call tp%info(i)%set_value(status="DISCARDED_PERI", discard_time=system%t, discard_xh=tp%xh(:,i), & discard_vh=tp%vh(:,i), discard_body_id=pl%id(j)) end if end if @@ -246,7 +246,7 @@ subroutine discard_pl_tp(tp, system, param) real(DP), dimension(NDIM) :: dx, dv character(len=STRMAX) :: idstri, idstrj, timestr - associate(ntp => tp%nbody, pl => system%pl, npl => system%pl%nbody, t => param%t, dt => param%dt) + associate(ntp => tp%nbody, pl => system%pl, npl => system%pl%nbody, t => system%t, dt => param%dt) do i = 1, ntp if (tp%status(i) == ACTIVE) then do j = 1, npl @@ -260,12 +260,12 @@ subroutine discard_pl_tp(tp, system, param) pl%ldiscard(j) = .true. write(idstri, *) tp%id(i) write(idstrj, *) pl%id(j) - write(timestr, *) param%t + write(timestr, *) system%t write(*, *) "Test particle " // trim(adjustl(tp%info(i)%name)) // " (" // trim(adjustl(idstri)) // ")" & // " too close to massive body " // trim(adjustl(pl%info(j)%name)) // " (" // trim(adjustl(idstrj)) // ")" & // " at t = " // trim(adjustl(timestr)) tp%ldiscard(i) = .true. - call tp%info(i)%set_value(status="DISCARDED_PLR", discard_time=param%t, discard_xh=tp%xh(:,i), & + call tp%info(i)%set_value(status="DISCARDED_PLR", discard_time=system%t, discard_xh=tp%xh(:,i), & discard_vh=tp%vh(:,i), discard_body_id=pl%id(j)) exit end if diff --git a/src/io/io.f90 b/src/io/io.f90 index 43489e4b0..f77278424 100644 --- a/src/io/io.f90 +++ b/src/io/io.f90 @@ -35,7 +35,7 @@ module subroutine io_compact_output(self, param, timer) select type(timer) class is (walltimer) - formatted_output = fmt("ILOOP",param%iloop) // fmt("T",param%t) // fmt("NPL",self%pl%nbody) // fmt("NTP",self%tp%nbody) + formatted_output = fmt("ILOOP",param%iloop) // fmt("T",self%t) // fmt("NPL",self%pl%nbody) // fmt("NTP",self%tp%nbody) select type(pl => self%pl) class is (symba_pl) formatted_output = formatted_output // fmt("NPLM",pl%nplm) @@ -250,7 +250,7 @@ module subroutine io_dump_system(self, param) dump_param%in_netcdf = trim(adjustl(DUMP_NC_FILE(idx))) dump_param%nciu%id_chunk = self%pl%nbody + self%tp%nbody dump_param%nciu%time_chunk = 1 - dump_param%tstart = param%t + dump_param%tstart = self%t call dump_param%dump(param_file_name) diff --git a/src/main/swiftest_driver.f90 b/src/main/swiftest_driver.f90 index 14b06485b..1ec5df880 100644 --- a/src/main/swiftest_driver.f90 +++ b/src/main/swiftest_driver.f90 @@ -66,7 +66,7 @@ program swiftest_driver call setup_construct_system(nbody_system, param) call param%read_in(param_file_name) - associate(t => param%t, & + associate(t => nbody_system%t, & t0 => param%t0, & tstart => param%tstart, & dt => param%dt, & @@ -137,13 +137,13 @@ program swiftest_driver idump = dump_cadence end if - tfrac = (param%t - param%t0) / (param%tstop - param%t0) + tfrac = (t - t0) / (tstop - t0) select type(pl => nbody_system%pl) class is (symba_pl) - write(display_unit, symbastatfmt) param%t, tfrac, pl%nplm, pl%nbody, nbody_system%tp%nbody + write(display_unit, symbastatfmt) t, tfrac, pl%nplm, pl%nbody, nbody_system%tp%nbody class default - write(display_unit, statusfmt) param%t, tfrac, pl%nbody, nbody_system%tp%nbody + write(display_unit, statusfmt) t, tfrac, pl%nbody, nbody_system%tp%nbody end select if (param%lenergy) call nbody_system%conservation_report(param, lterminal=.true.) call integration_timer%report(message="Integration steps:", unit=display_unit, nsubsteps=istep_out) diff --git a/src/modules/swiftest_classes.f90 b/src/modules/swiftest_classes.f90 index 14f0e3369..878b5d308 100644 --- a/src/modules/swiftest_classes.f90 +++ b/src/modules/swiftest_classes.f90 @@ -37,7 +37,6 @@ module swiftest_classes integer(I4B) :: maxid = -1 !! The current maximum particle id number integer(I4B) :: maxid_collision = 0 !! The current maximum collision id number real(DP) :: t0 = 0.0_DP !! Integration reference time - real(DP) :: t = -1.0_DP !! Integration current time real(DP) :: tstart = -1.0_DP !! Integration start time real(DP) :: tstop = -1.0_DP !! Integration stop time real(DP) :: dt = -1.0_DP !! Time step @@ -347,6 +346,7 @@ module swiftest_classes class(swiftest_tp), allocatable :: tp !! Test particle data structure class(swiftest_tp), allocatable :: tp_discards !! Discarded test particle data structure class(swiftest_pl), allocatable :: pl_discards !! Discarded massive body particle data structure + real(DP) :: t = -1.0_DP !! Integration current time real(DP) :: GMtot = 0.0_DP !! Total system mass - used for barycentric coordinate conversion real(DP) :: ke_orbit = 0.0_DP !! System orbital kinetic energy real(DP) :: ke_spin = 0.0_DP !! System spin kinetic energy diff --git a/src/netcdf/netcdf.f90 b/src/netcdf/netcdf.f90 index 89824e4ec..abd03f1d6 100644 --- a/src/netcdf/netcdf.f90 +++ b/src/netcdf/netcdf.f90 @@ -793,7 +793,7 @@ module subroutine netcdf_read_hdr_system(self, iu, param) tslot = int(param%ioutput, kind=I4B) + 1 call check( nf90_inquire_dimension(iu%ncid, iu%id_dimid, len=idmax), "netcdf_read_frame_system nf90_inquire_dimension id_dimid" ) - call check( nf90_get_var(iu%ncid, iu%time_varid, param%t, start=[tslot]), "netcdf_read_hdr_system nf90_getvar time_varid" ) + call check( nf90_get_var(iu%ncid, iu%time_varid, self%t, start=[tslot]), "netcdf_read_hdr_system nf90_getvar time_varid" ) allocate(gmtemp(idmax)) allocate(tpmask(idmax)) @@ -1428,7 +1428,7 @@ module subroutine netcdf_write_hdr_system(self, iu, param) tslot = int(param%ioutput, kind=I4B) + 1 - call check( nf90_put_var(iu%ncid, iu%time_varid, param%t, start=[tslot]), "netcdf_write_hdr_system nf90_put_var time_varid" ) + call check( nf90_put_var(iu%ncid, iu%time_varid, self%t, start=[tslot]), "netcdf_write_hdr_system nf90_put_var time_varid" ) call check( nf90_put_var(iu%ncid, iu%npl_varid, self%pl%nbody, start=[tslot]), "netcdf_write_hdr_system nf90_put_var npl_varid" ) call check( nf90_put_var(iu%ncid, iu%ntp_varid, self%tp%nbody, start=[tslot]), "netcdf_write_hdr_system nf90_put_var ntp_varid" ) select type(pl => self%pl) diff --git a/src/rmvs/rmvs_discard.f90 b/src/rmvs/rmvs_discard.f90 index 732cbdea0..60be2f6b0 100644 --- a/src/rmvs/rmvs_discard.f90 +++ b/src/rmvs/rmvs_discard.f90 @@ -29,7 +29,7 @@ module subroutine rmvs_discard_tp(self, system, param) if (self%nbody == 0) return - associate(tp => self, ntp => self%nbody, pl => system%pl, t => param%t) + associate(tp => self, ntp => self%nbody, pl => system%pl, t => system%t) do i = 1, ntp associate(iplperP => tp%plperP(i)) if ((tp%status(i) == ACTIVE) .and. (tp%lperi(i))) then diff --git a/src/setup/symba_collision.f90 b/src/setup/symba_collision.f90 new file mode 100644 index 000000000..e839af1de --- /dev/null +++ b/src/setup/symba_collision.f90 @@ -0,0 +1,1090 @@ +!! 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_classes) s_symba_collision + use swiftest + +contains + + module function symba_collision_casedisruption(system, param, colliders, frag) result(status) + !! author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton + !! + !! Create the fragments resulting from a non-catastrophic disruption collision + !! + implicit none + ! Arguments + class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object + class(symba_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA additions + class(fraggle_colliders), intent(inout) :: colliders !! Fraggle colliders object + class(fraggle_fragments), intent(inout) :: frag !! Fraggle fragmentation system object + ! Result + integer(I4B) :: status !! Status flag assigned to this outcome + ! Internals + integer(I4B) :: i, nfrag + logical :: lfailure + character(len=STRMAX) :: message + + select case(frag%regime) + case(COLLRESOLVE_REGIME_DISRUPTION) + message = "Disruption between" + case(COLLRESOLVE_REGIME_SUPERCATASTROPHIC) + message = "Supercatastrophic disruption between" + end select + call symba_collision_collider_message(system%pl, colliders%idx, message) + call io_log_one_message(FRAGGLE_LOG_OUT, message) + + ! Collisional fragments will be uniformly distributed around the pre-impact barycenter + call frag%set_mass_dist(colliders, param) + + ! Generate the position and velocity distributions of the fragments + call frag%generate_fragments(colliders, system, param, lfailure) + + if (lfailure) then + call io_log_one_message(FRAGGLE_LOG_OUT, "No fragment solution found, so treat as a pure hit-and-run") + status = ACTIVE + nfrag = 0 + select type(pl => system%pl) + class is (symba_pl) + pl%status(colliders%idx(:)) = status + pl%ldiscard(colliders%idx(:)) = .false. + pl%lcollision(colliders%idx(:)) = .false. + end select + else + ! Populate the list of new bodies + nfrag = frag%nbody + write(message, *) nfrag + call io_log_one_message(FRAGGLE_LOG_OUT, "Generating " // trim(adjustl(message)) // " fragments") + select case(frag%regime) + case(COLLRESOLVE_REGIME_DISRUPTION) + status = DISRUPTION + case(COLLRESOLVE_REGIME_SUPERCATASTROPHIC) + status = SUPERCATASTROPHIC + end select + frag%id(1:nfrag) = [(i, i = param%maxid + 1, param%maxid + nfrag)] + param%maxid = frag%id(nfrag) + call symba_collision_mergeaddsub(system, param, colliders, frag, status) + end if + + return + end function symba_collision_casedisruption + + + module function symba_collision_casehitandrun(system, param, colliders, frag) result(status) + !! author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton + !! + !! Create the fragments resulting from a non-catastrophic hit-and-run collision + !! + implicit none + ! Arguments + class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object + class(symba_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA additions + class(fraggle_colliders), intent(inout) :: colliders !! Fraggle colliders object + class(fraggle_fragments), intent(inout) :: frag !! Fraggle fragmentation system object + ! Result + integer(I4B) :: status !! Status flag assigned to this outcom + ! Internals + integer(I4B) :: i, ibiggest, nfrag, jtarg, jproj + logical :: lpure + character(len=STRMAX) :: message + + message = "Hit and run between" + call symba_collision_collider_message(system%pl, colliders%idx, message) + call io_log_one_message(FRAGGLE_LOG_OUT, trim(adjustl(message))) + + if (colliders%mass(1) > colliders%mass(2)) then + jtarg = 1 + jproj = 2 + else + jtarg = 2 + jproj = 1 + end if + + if (frag%mass_dist(2) > 0.9_DP * colliders%mass(jproj)) then ! Pure hit and run, so we'll just keep the two bodies untouched + call io_log_one_message(FRAGGLE_LOG_OUT, "Pure hit and run. No new fragments generated.") + nfrag = 0 + lpure = .true. + else ! Imperfect hit and run, so we'll keep the largest body and destroy the other + lpure = .false. + call frag%set_mass_dist(colliders, param) + + ! Generate the position and velocity distributions of the fragments + call frag%generate_fragments(colliders, system, param, lpure) + + if (lpure) then + call io_log_one_message(FRAGGLE_LOG_OUT, "Should have been a pure hit and run instead") + nfrag = 0 + else + nfrag = frag%nbody + write(message, *) nfrag + call io_log_one_message(FRAGGLE_LOG_OUT, "Generating " // trim(adjustl(message)) // " fragments") + end if + end if + if (lpure) then ! Reset these bodies back to being active so that nothing further is done to them + status = HIT_AND_RUN_PURE + select type(pl => system%pl) + class is (symba_pl) + pl%status(colliders%idx(:)) = ACTIVE + pl%ldiscard(colliders%idx(:)) = .false. + pl%lcollision(colliders%idx(:)) = .false. + end select + else + ibiggest = colliders%idx(maxloc(system%pl%Gmass(colliders%idx(:)), dim=1)) + frag%id(1) = system%pl%id(ibiggest) + frag%id(2:nfrag) = [(i, i = param%maxid + 1, param%maxid + nfrag - 1)] + param%maxid = frag%id(nfrag) + status = HIT_AND_RUN_DISRUPT + call symba_collision_mergeaddsub(system, param, colliders, frag, status) + end if + + return + end function symba_collision_casehitandrun + + + module function symba_collision_casemerge(system, param, colliders, frag) result(status) + !! author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton + !! + !! Merge massive bodies. + !! + !! Adapted from David E. Kaufmann's Swifter routines symba_merge_pl.f90 and symba_discard_merge_pl.f90 + !! + !! Adapted from Hal Levison's Swift routines symba5_merge.f and discard_mass_merge.f + implicit none + ! Arguments + class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object + class(symba_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA additions + class(fraggle_colliders), intent(inout) :: colliders !! Fraggle colliders object + class(fraggle_fragments), intent(inout) :: frag !! Fraggle fragmentation system object + ! Result + integer(I4B) :: status !! Status flag assigned to this outcome + ! Internals + integer(I4B) :: i, j, k, ibiggest + real(DP) :: pe + real(DP), dimension(NDIM) :: L_spin_new + character(len=STRMAX) :: message + + message = "Merging" + call symba_collision_collider_message(system%pl, colliders%idx, message) + call io_log_one_message(FRAGGLE_LOG_OUT, message) + + select type(pl => system%pl) + class is (symba_pl) + + call frag%set_mass_dist(colliders, param) + ibiggest = colliders%idx(maxloc(pl%Gmass(colliders%idx(:)), dim=1)) + frag%id(1) = pl%id(ibiggest) + frag%xb(:,1) = frag%xbcom(:) + frag%vb(:,1) = frag%vbcom(:) + + if (param%lrotation) then + ! Conserve angular momentum by putting pre-impact orbital momentum into spin of the new body + L_spin_new(:) = colliders%L_orbit(:,1) + colliders%L_orbit(:,2) + colliders%L_spin(:,1) + colliders%L_spin(:,2) + + ! Assume prinicpal axis rotation on 3rd Ip axis + frag%rot(:,1) = L_spin_new(:) / (frag%Ip(3,1) * frag%mass(1) * frag%radius(1)**2) + else ! If spin is not enabled, we will consider the lost pre-collision angular momentum as "escaped" and add it to our bookkeeping variable + param%Lescape(:) = param%Lescape(:) + colliders%L_orbit(:,1) + colliders%L_orbit(:,2) + end if + + ! Keep track of the component of potential energy due to the pre-impact colliders%idx for book-keeping + pe = 0.0_DP + do j = 1, colliders%ncoll + do i = j + 1, colliders%ncoll + pe = pe - pl%Gmass(i) * pl%mass(j) / norm2(pl%xb(:, i) - pl%xb(:, j)) + end do + end do + system%Ecollisions = system%Ecollisions + pe + system%Euntracked = system%Euntracked - pe + + ! Update any encounter lists that have the removed bodies in them so that they instead point to the new + do k = 1, system%plplenc_list%nenc + do j = 1, colliders%ncoll + i = colliders%idx(j) + if (i == ibiggest) cycle + if (system%plplenc_list%id1(k) == pl%id(i)) then + system%plplenc_list%id1(k) = pl%id(ibiggest) + system%plplenc_list%index1(k) = i + end if + if (system%plplenc_list%id2(k) == pl%id(i)) then + system%plplenc_list%id2(k) = pl%id(ibiggest) + system%plplenc_list%index2(k) = i + end if + if (system%plplenc_list%id1(k) == system%plplenc_list%id2(k)) system%plplenc_list%status(k) = INACTIVE + end do + end do + + status = MERGED + + call symba_collision_mergeaddsub(system, param, colliders, frag, status) + + end select + + return + end function symba_collision_casemerge + + + subroutine symba_collision_collider_message(pl, collidx, collider_message) + !! author: David A. Minton + !! + !! Prints a nicely formatted message about which bodies collided, including their names and ids. + !! This subroutine appends the body names and ids to an input message. + implicit none + ! Arguments + class(swiftest_pl), intent(in) :: pl !! Swiftest massive body object + integer(I4B), dimension(:), intent(in) :: collidx !! Index of collisional colliders%idx members + character(*), intent(inout) :: collider_message !! The message to print to the screen. + ! Internals + integer(I4B) :: i, n + character(len=STRMAX) :: idstr + + n = size(collidx) + if (n == 0) return + + do i = 1, n + if (i > 1) collider_message = trim(adjustl(collider_message)) // " and " + collider_message = " " // trim(adjustl(collider_message)) // " " // trim(adjustl(pl%info(collidx(i))%name)) + write(idstr, '(I10)') pl%id(collidx(i)) + collider_message = trim(adjustl(collider_message)) // " (" // trim(adjustl(idstr)) // ") " + end do + + return + end subroutine symba_collision_collider_message + + + module function symba_collision_check_encounter(self, system, param, t, dt, irec) result(lany_collision) + !! 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.f90 and symba_merge_tp.f90 + !! + !! Adapted from Hal Levison's Swift routine symba5_merge.f + implicit none + ! Arguments + class(symba_encounter), 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 + ! Result + logical :: lany_collision !! Returns true if cany pair of encounters resulted in a collision + ! Internals + logical, dimension(:), allocatable :: lcollision, lmask + real(DP), dimension(NDIM) :: xr, vr + integer(I4B) :: i, j, k, nenc + real(DP) :: rlim, Gmtot + logical :: isplpl + character(len=STRMAX) :: timestr, idstri, idstrj, message + class(symba_encounter), allocatable :: tmp + + lany_collision = .false. + 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) + nenc = self%nenc + allocate(lmask(nenc)) + lmask(:) = ((self%status(1:nenc) == ACTIVE) .and. (pl%levelg(self%index1(1:nenc)) >= irec)) + if (isplpl) then + lmask(:) = lmask(:) .and. (pl%levelg(self%index2(1:nenc)) >= irec) + else + lmask(:) = lmask(:) .and. (tp%levelg(self%index2(1:nenc)) >= irec) + end if + if (.not.any(lmask(:))) return + + allocate(lcollision(nenc)) + lcollision(:) = .false. + + if (isplpl) then + do concurrent(k = 1:nenc, lmask(k)) + i = self%index1(k) + j = self%index2(k) + xr(:) = pl%xh(:, i) - pl%xh(:, j) + vr(:) = pl%vb(:, i) - pl%vb(:, j) + rlim = pl%radius(i) + pl%radius(j) + Gmtot = pl%Gmass(i) + pl%Gmass(j) + lcollision(k) = symba_collision_check_one(xr(1), xr(2), xr(3), vr(1), vr(2), vr(3), & + Gmtot, rlim, dt, self%lvdotr(k)) + end do + else + do concurrent(k = 1:nenc, lmask(k)) + i = self%index1(k) + j = self%index2(k) + xr(:) = pl%xh(:, i) - tp%xh(:, j) + vr(:) = pl%vb(:, i) - tp%vb(:, j) + lcollision(k) = symba_collision_check_one(xr(1), xr(2), xr(3), vr(1), vr(2), vr(3), & + pl%Gmass(i), pl%radius(i), dt, self%lvdotr(k)) + end do + end if + + lany_collision = any(lcollision(:)) + if (lany_collision) then + call pl%xh2xb(system%cb) ! Update the central body barycenteric position vector to get us out of DH and into bary + do k = 1, nenc + i = self%index1(k) + j = self%index2(k) + if (lcollision(k)) self%status(k) = COLLISION + self%t(k) = t + self%x1(:,k) = pl%xh(:,i) + system%cb%xb(:) + self%v1(:,k) = pl%vb(:,i) + if (isplpl) then + self%x2(:,k) = pl%xh(:,j) + system%cb%xb(:) + self%v2(:,k) = pl%vb(:,j) + if (lcollision(k)) then + ! Check to see if either of these bodies has been involved with a collision before, and if so, make this a collisional colliders%idx + if (pl%lcollision(i) .or. pl%lcollision(j)) call pl%make_colliders([i,j]) + + ! Set the collision flag for these to bodies to true in case they become involved in another collision later in the step + pl%lcollision([i, j]) = .true. + pl%status([i, j]) = COLLISION + call pl%info(i)%set_value(status="COLLISION", discard_time=t, discard_xh=pl%xh(:,i), discard_vh=pl%vh(:,i)) + call pl%info(j)%set_value(status="COLLISION", discard_time=t, discard_xh=pl%xh(:,j), discard_vh=pl%vh(:,j)) + end if + else + self%x2(:,k) = tp%xh(:,j) + system%cb%xb(:) + self%v2(:,k) = tp%vb(:,j) + if (lcollision(k)) then + tp%status(j) = DISCARDED_PLR + tp%ldiscard(j) = .true. + write(idstri, *) pl%id(i) + write(idstrj, *) tp%id(j) + write(timestr, *) t + call tp%info(j)%set_value(status="DISCARDED_PLR", discard_time=t, discard_xh=tp%xh(:,j), discard_vh=tp%vh(:,j)) + write(message, *) "Particle " // trim(adjustl(tp%info(j)%name)) // " (" // trim(adjustl(idstrj)) // ")" & + // " collided with massive body " // trim(adjustl(pl%info(i)%name)) // " (" // trim(adjustl(idstri)) // ")" & + // " at t = " // trim(adjustl(timestr)) + call io_log_one_message(FRAGGLE_LOG_OUT, message) + end if + end if + end do + end if + end select + end select + + ! Extract the pl-pl or pl-tp encounter list and return the pl-pl or pl-tp collision_list + if (lany_collision) then + select type(self) + class is (symba_plplenc) + call self%extract_collisions(system, param) + class default + allocate(tmp, mold=self) + call self%spill(tmp, lcollision, ldestructive=.true.) ! Remove this encounter pair from the encounter list + end select + end if + + return + end function symba_collision_check_encounter + + + pure elemental function symba_collision_check_one(xr, yr, zr, vxr, vyr, vzr, Gmtot, rlim, dt, lvdotr) result(lcollision) + !! author: David A. Minton + !! + !! Check for a merger between a single pair of particles + !! + !! Adapted from David E. Kaufmann's Swifter routines symba_merge_tp.f90 and symba_merge_pl.f90 + !! + !! Adapted from Hal Levison's Swift routine symba5_merge.f + implicit none + ! Arguments + real(DP), intent(in) :: xr, yr, zr !! Relative position vector components + real(DP), intent(in) :: vxr, vyr, vzr !! Relative velocity vector components + real(DP), intent(in) :: Gmtot !! Sum of G*mass of colliding bodies + real(DP), intent(in) :: rlim !! Collision limit - Typically the sum of the radii of colliding bodies + real(DP), intent(in) :: dt !! Step size + logical, intent(in) :: lvdotr !! Logical flag indicating that these two bodies are approaching in the current substep + ! Result + logical :: lcollision !! Logical flag indicating whether these two bodies will collide or not + ! Internals + real(DP) :: r2, rlim2, a, e, q, vdotr, tcr2, dt2 + + r2 = xr**2 + yr**2 + zr**2 + rlim2 = rlim**2 + + if (r2 <= rlim2) then ! checks if bodies are actively colliding in this time step + lcollision = .true. + else ! if they are not actively colliding in this time step, checks if they are going to collide next time step based on velocities and q + lcollision = .false. + vdotr = xr * vxr + yr * vyr + zr * vzr + if (lvdotr .and. (vdotr > 0.0_DP)) then + tcr2 = r2 / (vxr**2 + vyr**2 + vzr**2) + dt2 = dt**2 + if (tcr2 <= dt2) then + call orbel_xv2aeq(Gmtot, xr, yr, zr, vxr, vyr, vzr, a, e, q) + lcollision = (q < rlim) + end if + end if + end if + + return + end function symba_collision_check_one + + + function symba_collision_consolidate_colliders(pl, cb, param, idx_parent, colliders) result(lflag) + !! author: David A. Minton + !! + !! Loops through the pl-pl collision list and groups families together by index. Outputs the indices of all colliders%idx members, + !! and pairs of quantities (x and v vectors, mass, radius, L_spin, and Ip) that can be used to resolve the collisional outcome. + implicit none + ! Arguments + class(symba_pl), intent(inout) :: pl !! SyMBA massive body object + class(symba_cb), intent(inout) :: cb !! SyMBA central body object + class(symba_parameters), intent(in) :: param !! Current run configuration parameters with SyMBA additions + integer(I4B), dimension(2), intent(inout) :: idx_parent !! Index of the two bodies considered the "parents" of the collision + class(fraggle_colliders), intent(out) :: colliders + ! Result + logical :: lflag !! Logical flag indicating whether a colliders%idx was successfully created or not + ! Internals + type collidx_array + integer(I4B), dimension(:), allocatable :: id + integer(I4B), dimension(:), allocatable :: idx + end type collidx_array + type(collidx_array), dimension(2) :: parent_child_index_array + integer(I4B), dimension(2) :: nchild + integer(I4B) :: i, j, ncolliders, idx_child + real(DP), dimension(2) :: volume, density + real(DP) :: mchild, volchild + real(DP), dimension(NDIM) :: xc, vc, xcom, vcom, xchild, vchild, xcrossv + real(DP), dimension(NDIM,2) :: mxc, vcc + + nchild(:) = pl%kin(idx_parent(:))%nchild + ! If all of these bodies share a parent, but this is still a unique collision, move the last child + ! out of the parent's position and make it the secondary body + if (idx_parent(1) == idx_parent(2)) then + if (nchild(1) == 0) then ! There is only one valid body recorded in this pair (this could happen due to restructuring of the kinship relationships, though it should be rare) + lflag = .false. + call pl%reset_kinship([idx_parent(1)]) + return + end if + idx_parent(2) = pl%kin(idx_parent(1))%child(nchild(1)) + nchild(1) = nchild(1) - 1 + nchild(2) = 0 + pl%kin(idx_parent(:))%nchild = nchild(:) + pl%kin(idx_parent(2))%parent = idx_parent(1) + end if + + colliders%mass(:) = pl%mass(idx_parent(:)) ! Note: This is meant to mass, not G*mass, as the collisional regime determination uses mass values that will be converted to Si + colliders%radius(:) = pl%radius(idx_parent(:)) + volume(:) = (4.0_DP / 3.0_DP) * PI * colliders%radius(:)**3 + + ! Group together the ids and indexes of each collisional parent and its children + do j = 1, 2 + allocate(parent_child_index_array(j)%idx(nchild(j)+ 1)) + allocate(parent_child_index_array(j)%id(nchild(j)+ 1)) + associate(idx_arr => parent_child_index_array(j)%idx, & + id_arr => parent_child_index_array(j)%id, & + ncj => nchild(j), & + pl => pl, & + plkinj => pl%kin(idx_parent(j))) + idx_arr(1) = idx_parent(j) + if (ncj > 0) idx_arr(2:ncj + 1) = plkinj%child(1:ncj) + id_arr(:) = pl%id(idx_arr(:)) + end associate + end do + + ! Consolidate the groups of collsional parents with any children they may have into a single "colliders%idx" index array + ncolliders = 2 + sum(nchild(:)) + allocate(colliders%idx(ncolliders)) + colliders%idx = [parent_child_index_array(1)%idx(:),parent_child_index_array(2)%idx(:)] + + colliders%ncoll = count(pl%lcollision(colliders%idx(:))) + colliders%idx = pack(colliders%idx(:), pl%lcollision(colliders%idx(:))) + colliders%L_spin(:,:) = 0.0_DP + colliders%Ip(:,:) = 0.0_DP + + ! Find the barycenter of each body along with its children, if it has any + do j = 1, 2 + colliders%xb(:, j) = pl%xh(:, idx_parent(j)) + cb%xb(:) + colliders%vb(:, j) = pl%vb(:, idx_parent(j)) + ! Assume principal axis rotation about axis corresponding to highest moment of inertia (3rd Ip) + if (param%lrotation) then + colliders%Ip(:, j) = colliders%mass(j) * pl%Ip(:, idx_parent(j)) + colliders%L_spin(:, j) = colliders%Ip(3, j) * colliders%radius(j)**2 * pl%rot(:, idx_parent(j)) + end if + + 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 (.not. pl%lcollision(idx_child)) cycle + mchild = pl%mass(idx_child) + xchild(:) = pl%xh(:, idx_child) + cb%xb(:) + vchild(:) = pl%vb(:, idx_child) + volchild = (4.0_DP / 3.0_DP) * PI * pl%radius(idx_child)**3 + volume(j) = volume(j) + volchild + ! Get angular momentum of the child-parent pair and add that to the spin + ! Add the child's spin + if (param%lrotation) then + xcom(:) = (colliders%mass(j) * colliders%xb(:,j) + mchild * xchild(:)) / (colliders%mass(j) + mchild) + vcom(:) = (colliders%mass(j) * colliders%vb(:,j) + mchild * vchild(:)) / (colliders%mass(j) + mchild) + xc(:) = colliders%xb(:, j) - xcom(:) + vc(:) = colliders%vb(:, j) - vcom(:) + xcrossv(:) = xc(:) .cross. vc(:) + colliders%L_spin(:, j) = colliders%L_spin(:, j) + colliders%mass(j) * xcrossv(:) + + xc(:) = xchild(:) - xcom(:) + vc(:) = vchild(:) - vcom(:) + xcrossv(:) = xc(:) .cross. vc(:) + colliders%L_spin(:, j) = colliders%L_spin(:, j) + mchild * xcrossv(:) + + colliders%L_spin(:, j) = colliders%L_spin(:, j) + mchild * pl%Ip(3, idx_child) & + * pl%radius(idx_child)**2 & + * pl%rot(:, idx_child) + colliders%Ip(:, j) = colliders%Ip(:, j) + mchild * pl%Ip(:, idx_child) + end if + + ! Merge the child and parent + colliders%mass(j) = colliders%mass(j) + mchild + colliders%xb(:, j) = xcom(:) + colliders%vb(:, j) = vcom(:) + end do + end if + density(j) = colliders%mass(j) / volume(j) + colliders%radius(j) = (3 * volume(j) / (4 * PI))**(1.0_DP / 3.0_DP) + if (param%lrotation) colliders%Ip(:, j) = colliders%Ip(:, j) / colliders%mass(j) + end do + lflag = .true. + + xcom(:) = (colliders%mass(1) * colliders%xb(:, 1) + colliders%mass(2) * colliders%xb(:, 2)) / sum(colliders%mass(:)) + vcom(:) = (colliders%mass(1) * colliders%vb(:, 1) + colliders%mass(2) * colliders%vb(:, 2)) / sum(colliders%mass(:)) + mxc(:, 1) = colliders%mass(1) * (colliders%xb(:, 1) - xcom(:)) + mxc(:, 2) = colliders%mass(2) * (colliders%xb(:, 2) - xcom(:)) + vcc(:, 1) = colliders%vb(:, 1) - vcom(:) + vcc(:, 2) = colliders%vb(:, 2) - vcom(:) + colliders%L_orbit(:,:) = mxc(:,:) .cross. vcc(:,:) + + ! Destroy the kinship relationships for all members of this colliders%idx + call pl%reset_kinship(colliders%idx(:)) + + return + end function symba_collision_consolidate_colliders + + + module subroutine symba_collision_encounter_extract_collisions(self, system, param) + !! author: David A. Minton + !! + !! Processes the pl-pl encounter list remove only those encounters that led to a collision + !! + implicit none + ! Arguments + class(symba_plplenc), intent(inout) :: self !! SyMBA pl-pl encounter list + class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object + class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters + ! Internals + logical, dimension(:), allocatable :: lplpl_collision + logical, dimension(:), allocatable :: lplpl_unique_parent + integer(I4B), dimension(:), pointer :: plparent + integer(I4B), dimension(:), allocatable :: collision_idx, unique_parent_idx + integer(I4B) :: i, index_coll, ncollisions, nunique_parent, nplplenc + + select type (pl => system%pl) + class is (symba_pl) + associate(plplenc_list => self, idx1 => self%index1, idx2 => self%index2, plparent => pl%kin%parent) + nplplenc = plplenc_list%nenc + allocate(lplpl_collision(nplplenc)) + lplpl_collision(:) = plplenc_list%status(1:nplplenc) == COLLISION + if (.not.any(lplpl_collision)) return + ! Collisions have been detected in this step. So we need to determine which of them are between unique bodies. + + ! Get the subset of pl-pl encounters that lead to a collision + ncollisions = count(lplpl_collision(:)) + allocate(collision_idx(ncollisions)) + collision_idx = pack([(i, i=1, nplplenc)], lplpl_collision) + + ! Get the subset of collisions that involve a unique pair of parents + allocate(lplpl_unique_parent(ncollisions)) + + lplpl_unique_parent(:) = plparent(idx1(collision_idx(:))) /= plparent(idx2(collision_idx(:))) + nunique_parent = count(lplpl_unique_parent(:)) + allocate(unique_parent_idx(nunique_parent)) + unique_parent_idx = pack(collision_idx(:), lplpl_unique_parent(:)) + + ! Scrub all pl-pl collisions involving unique pairs of parents, which will remove all duplicates and leave behind + ! all pairs that have themselves as parents but are not part of the unique parent list. This can hapepn in rare cases + ! due to restructuring of parent/child relationships when there are large numbers of multi-body collisions in a single + ! step + lplpl_unique_parent(:) = .true. + do index_coll = 1, ncollisions + associate(ip1 => plparent(idx1(collision_idx(index_coll))), ip2 => plparent(idx2(collision_idx(index_coll)))) + lplpl_unique_parent(:) = .not. ( any(plparent(idx1(unique_parent_idx(:))) == ip1) .or. & + any(plparent(idx2(unique_parent_idx(:))) == ip1) .or. & + any(plparent(idx1(unique_parent_idx(:))) == ip2) .or. & + any(plparent(idx2(unique_parent_idx(:))) == ip2) ) + end associate + end do + + ! Reassemble collision index list to include only those containing the unique pairs of parents, plus all the non-unique pairs that don't + ! contain a parent body on the unique parent list. + ncollisions = nunique_parent + count(lplpl_unique_parent) + collision_idx = [unique_parent_idx(:), pack(collision_idx(:), lplpl_unique_parent(:))] + + ! Create a mask that contains only the pl-pl encounters that did not result in a collision, and then discard them + lplpl_collision(:) = .false. + lplpl_collision(collision_idx(:)) = .true. + call plplenc_list%spill(system%plplcollision_list, lplpl_collision, ldestructive=.true.) ! Extract any encounters that are not collisions from the list. + end associate + end select + + return + end subroutine symba_collision_encounter_extract_collisions + + + module subroutine symba_collision_make_colliders_pl(self, idx) + !! author: Jennifer L.L. Pouplin, Carlisle A. wishard, and David A. Minton + !! + !! When a single body is involved in more than one collision in a single step, it becomes part of a colliders%idx. + !! The largest body involved in a multi-body collision is the "parent" and all bodies that collide with it are its "children," + !! including those that collide with the children. + !! + !! 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_pl), intent(inout) :: self !! SyMBA massive body object + integer(I4B), dimension(2), intent(in) :: idx !! Array holding the indices of the two bodies involved in the collision + ! Internals + integer(I4B) :: i, j, index_parent, index_child, p1, p2 + integer(I4B) :: nchild_inherit, nchild_orig, nchild_new + integer(I4B), dimension(:), allocatable :: temp + + associate(pl => self) + p1 = pl%kin(idx(1))%parent + p2 = pl%kin(idx(2))%parent + if (p1 == p2) return ! This is a collision between to children of a shared parent. We will ignore it. + + if (pl%mass(p1) > pl%mass(p2)) then + index_parent = p1 + index_child = p2 + else + index_parent = p2 + index_child = p1 + end if + + ! Expand the child array (or create it if necessary) and copy over the previous lists of children + nchild_orig = pl%kin(index_parent)%nchild + nchild_inherit = pl%kin(index_child)%nchild + nchild_new = nchild_orig + nchild_inherit + 1 + allocate(temp(nchild_new)) + + if (nchild_orig > 0) temp(1:nchild_orig) = pl%kin(index_parent)%child(1:nchild_orig) + ! Find out if the child body has any children of its own. The new parent wil inherit these children + if (nchild_inherit > 0) then + temp(nchild_orig+1:nchild_orig+nchild_inherit) = pl%kin(index_child)%child(1:nchild_inherit) + do i = 1, nchild_inherit + j = pl%kin(index_child)%child(i) + ! Set the childrens' parent to the new parent + pl%kin(j)%parent = index_parent + end do + end if + call pl%reset_kinship([index_child]) + ! Add the new child to its parent + pl%kin(index_child)%parent = index_parent + temp(nchild_new) = index_child + ! Save the new child array to the parent + pl%kin(index_parent)%nchild = nchild_new + call move_alloc(from=temp, to=pl%kin(index_parent)%child) + end associate + + return + end subroutine symba_collision_make_colliders_pl + + + subroutine symba_collision_mergeaddsub(system, param, colliders, frag, status) + !! author: David A. Minton + !! + !! Fills the pl_discards and pl_adds with removed and added bodies + !! + implicit none + ! Arguments + class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object + class(symba_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA additions + class(fraggle_colliders), intent(inout) :: colliders !! Fraggle colliders object + class(fraggle_fragments), intent(inout) :: frag !! Fraggle fragmentation system object + integer(I4B), intent(in) :: status !! Status flag to assign to adds + ! Internals + integer(I4B) :: i, ibiggest, ismallest, iother, nstart, nend, ncolliders, nfrag + logical, dimension(system%pl%nbody) :: lmask + class(symba_pl), allocatable :: plnew, plsub + character(*), parameter :: FRAGFMT = '("Newbody",I0.7)' + character(len=NAMELEN) :: newname + + select type(pl => system%pl) + class is (symba_pl) + select type(pl_discards => system%pl_discards) + class is (symba_merger) + associate(info => pl%info, pl_adds => system%pl_adds, cb => system%cb, npl => pl%nbody) + ! Add the colliders%idx bodies to the subtraction list + ncolliders = colliders%ncoll + nfrag = frag%nbody + + param%maxid_collision = max(param%maxid_collision, maxval(system%pl%info(:)%collision_id)) + param%maxid_collision = param%maxid_collision + 1 + + ! Setup new bodies + allocate(plnew, mold=pl) + call plnew%setup(nfrag, param) + ibiggest = colliders%idx(maxloc(pl%Gmass(colliders%idx(:)), dim=1)) + ismallest = colliders%idx(minloc(pl%Gmass(colliders%idx(:)), dim=1)) + + ! Copy over identification, information, and physical properties of the new bodies from the fragment list + plnew%id(1:nfrag) = frag%id(1:nfrag) + plnew%xb(:, 1:nfrag) = frag%xb(:, 1:nfrag) + plnew%vb(:, 1:nfrag) = frag%vb(:, 1:nfrag) + call pl%vb2vh(cb) + call pl%xh2xb(cb) + do i = 1, nfrag + plnew%xh(:,i) = frag%xb(:, i) - cb%xb(:) + plnew%vh(:,i) = frag%vb(:, i) - cb%vb(:) + end do + plnew%mass(1:nfrag) = frag%mass(1:nfrag) + plnew%Gmass(1:nfrag) = param%GU * frag%mass(1:nfrag) + plnew%radius(1:nfrag) = frag%radius(1:nfrag) + plnew%density(1:nfrag) = frag%mass(1:nfrag) / frag%radius(1:nfrag) + call plnew%set_rhill(cb) + + select case(status) + case(DISRUPTION) + plnew%status(1:nfrag) = NEW_PARTICLE + do i = 1, nfrag + write(newname, FRAGFMT) frag%id(i) + call plnew%info(i)%set_value(origin_type="Disruption", origin_time=system%t, name=newname, & + origin_xh=plnew%xh(:,i), & + origin_vh=plnew%vh(:,i), collision_id=param%maxid_collision) + end do + do i = 1, ncolliders + if (colliders%idx(i) == ibiggest) then + iother = ismallest + else + iother = ibiggest + end if + call pl%info(colliders%idx(i))%set_value(status="Disruption", discard_time=system%t, & + discard_xh=pl%xh(:,i), discard_vh=pl%vh(:,i), discard_body_id=iother) + end do + case(SUPERCATASTROPHIC) + plnew%status(1:nfrag) = NEW_PARTICLE + do i = 1, nfrag + write(newname, FRAGFMT) frag%id(i) + call plnew%info(i)%set_value(origin_type="Supercatastrophic", origin_time=system%t, name=newname, & + origin_xh=plnew%xh(:,i), origin_vh=plnew%vh(:,i), & + collision_id=param%maxid_collision) + end do + do i = 1, ncolliders + if (colliders%idx(i) == ibiggest) then + iother = ismallest + else + iother = ibiggest + end if + call pl%info(colliders%idx(i))%set_value(status="Supercatastrophic", discard_time=system%t, & + discard_xh=pl%xh(:,i), discard_vh=pl%vh(:,i), & + discard_body_id=iother) + end do + case(HIT_AND_RUN_DISRUPT) + call plnew%info(1)%copy(pl%info(ibiggest)) + plnew%status(1) = OLD_PARTICLE + do i = 2, nfrag + write(newname, FRAGFMT) frag%id(i) + call plnew%info(i)%set_value(origin_type="Hit and run fragment", origin_time=system%t, name=newname, & + origin_xh=plnew%xh(:,i), origin_vh=plnew%vh(:,i), & + collision_id=param%maxid_collision) + end do + do i = 1, ncolliders + if (colliders%idx(i) == ibiggest) cycle + iother = ibiggest + call pl%info(colliders%idx(i))%set_value(status="Hit and run fragmention", discard_time=system%t, & + discard_xh=pl%xh(:,i), discard_vh=pl%vh(:,i), & + discard_body_id=iother) + end do + case(MERGED) + call plnew%info(1)%copy(pl%info(ibiggest)) + plnew%status(1) = OLD_PARTICLE + do i = 1, ncolliders + if (colliders%idx(i) == ibiggest) cycle + + iother = ibiggest + call pl%info(colliders%idx(i))%set_value(status="MERGED", discard_time=system%t, discard_xh=pl%xh(:,i), & + discard_vh=pl%vh(:,i), discard_body_id=iother) + end do + end select + + if (param%lrotation) then + plnew%Ip(:, 1:nfrag) = frag%Ip(:, 1:nfrag) + plnew%rot(:, 1:nfrag) = frag%rot(:, 1:nfrag) + end if + + ! if (param%ltides) then + ! plnew%Q = pl%Q(ibiggest) + ! plnew%k2 = pl%k2(ibiggest) + ! plnew%tlag = pl%tlag(ibiggest) + ! end if + + !Copy over or set integration parameters for new bodies + plnew%lcollision(1:nfrag) = .false. + plnew%ldiscard(1:nfrag) = .false. + plnew%levelg(1:nfrag) = pl%levelg(ibiggest) + plnew%levelm(1:nfrag) = pl%levelm(ibiggest) + + ! Log the properties of the new bodies + call fraggle_io_log_pl(plnew, param) + + ! Append the new merged body to the list + nstart = pl_adds%nbody + 1 + nend = pl_adds%nbody + nfrag + call pl_adds%append(plnew, lsource_mask=[(.true., i=1, nfrag)]) + ! Record how many bodies were added in this event + pl_adds%ncomp(nstart:nend) = plnew%nbody + + ! Add the discarded bodies to the discard list + pl%status(colliders%idx(:)) = MERGED + pl%ldiscard(colliders%idx(:)) = .true. + pl%lcollision(colliders%idx(:)) = .true. + lmask(:) = .false. + lmask(colliders%idx(:)) = .true. + + call plnew%setup(0, param) + deallocate(plnew) + + allocate(plsub, mold=pl) + call pl%spill(plsub, lmask, ldestructive=.false.) + + nstart = pl_discards%nbody + 1 + nend = pl_discards%nbody + ncolliders + call pl_discards%append(plsub, lsource_mask=[(.true., i = 1, ncolliders)]) + + ! Record how many bodies were subtracted in this event + pl_discards%ncomp(nstart:nend) = ncolliders + + call plsub%setup(0, param) + deallocate(plsub) + end associate + end select + end select + + return + end subroutine symba_collision_mergeaddsub + + + module subroutine symba_collision_resolve_fragmentations(self, system, param) + !! author: David A. Minton + !! + !! Process list of collisions, determine the collisional regime, and then create fragments. + !! + implicit none + ! Arguments + class(symba_plplenc), intent(inout) :: self !! SyMBA pl-pl encounter list + class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object + class(symba_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA additions + ! Internals + ! Internals + integer(I4B), dimension(2) :: idx_parent !! Index of the two bodies considered the "parents" of the collision + logical :: lgoodcollision + integer(I4B) :: i + type(fraggle_colliders) :: colliders !! Fraggle colliders object + type(fraggle_fragments) :: frag !! Fraggle fragmentation system object + + associate(plplcollision_list => self, ncollisions => self%nenc, idx1 => self%index1, idx2 => self%index2) + select type(pl => system%pl) + class is (symba_pl) + select type (cb => system%cb) + class is (symba_cb) + do i = 1, ncollisions + idx_parent(1) = pl%kin(idx1(i))%parent + idx_parent(2) = pl%kin(idx2(i))%parent + lgoodcollision = symba_collision_consolidate_colliders(pl, cb, param, idx_parent, colliders) + if ((.not. lgoodcollision) .or. any(pl%status(idx_parent(:)) /= COLLISION)) cycle + + call colliders%regime(frag, system, param) + + select case (frag%regime) + case (COLLRESOLVE_REGIME_DISRUPTION, COLLRESOLVE_REGIME_SUPERCATASTROPHIC) + plplcollision_list%status(i) = symba_collision_casedisruption(system, param, colliders, frag) + case (COLLRESOLVE_REGIME_HIT_AND_RUN) + plplcollision_list%status(i) = symba_collision_casehitandrun(system, param, colliders, frag) + case (COLLRESOLVE_REGIME_MERGE, COLLRESOLVE_REGIME_GRAZE_AND_MERGE) + plplcollision_list%status(i) = symba_collision_casemerge(system, param, colliders, frag) + case default + write(*,*) "Error in symba_collision, unrecognized collision regime" + call util_exit(FAILURE) + end select + end do + end select + end select + end associate + + return + end subroutine symba_collision_resolve_fragmentations + + + module subroutine symba_collision_resolve_mergers(self, system, param) + !! author: David A. Minton + !! + !! Process list of collisions and merge colliding bodies together. + !! + implicit none + ! Arguments + class(symba_plplenc), intent(inout) :: self !! SyMBA pl-pl encounter list + class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object + class(symba_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA additions + ! Internals + integer(I4B), dimension(2) :: idx_parent !! Index of the two bodies considered the "parents" of the collision + logical :: lgoodcollision + integer(I4B) :: i + type(fraggle_colliders) :: colliders !! Fraggle colliders object + type(fraggle_fragments) :: frag !! Fraggle fragmentation system object + + associate(plplcollision_list => self, ncollisions => self%nenc, idx1 => self%index1, idx2 => self%index2) + select type(pl => system%pl) + class is (symba_pl) + select type(cb => system%cb) + class is (symba_cb) + do i = 1, ncollisions + idx_parent(1) = pl%kin(idx1(i))%parent + idx_parent(2) = pl%kin(idx2(i))%parent + lgoodcollision = symba_collision_consolidate_colliders(pl, cb, param, idx_parent, colliders) + if (.not. lgoodcollision) cycle + if (any(pl%status(idx_parent(:)) /= COLLISION)) cycle ! One of these two bodies has already been resolved + + frag%regime = COLLRESOLVE_REGIME_MERGE + frag%mtot = sum(colliders%mass(:)) + frag%mass_dist(1) = frag%mtot + frag%mass_dist(2) = 0.0_DP + frag%mass_dist(3) = 0.0_DP + frag%xbcom(:) = (colliders%mass(1) * colliders%xb(:,1) + colliders%mass(2) * colliders%xb(:,2)) / frag%mtot + frag%vbcom(:) = (colliders%mass(1) * colliders%vb(:,1) + colliders%mass(2) * colliders%vb(:,2)) / frag%mtot + plplcollision_list%status(i) = symba_collision_casemerge(system, param, colliders, frag) + end do + end select + end select + end associate + + return + end subroutine symba_collision_resolve_mergers + + + module subroutine symba_collision_resolve_plplenc(self, system, param, t, dt, irec) + !! author: David A. Minton + !! + !! Process the pl-pl collision list, then modifiy the massive bodies based on the outcome of the collision + !! + implicit none + ! Arguments + class(symba_plplenc), intent(inout) :: self !! SyMBA pl-pl encounter list + class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA additions + real(DP), intent(in) :: t !! Current simulation time + real(DP), intent(in) :: dt !! Current simulation step size + integer(I4B), intent(in) :: irec !! Current recursion level + ! Internals + real(DP) :: Eorbit_before, Eorbit_after + logical :: lplpl_collision + character(len=STRMAX) :: timestr + class(symba_parameters), allocatable :: tmp_param + + associate(plplenc_list => self, plplcollision_list => system%plplcollision_list) + select type(pl => system%pl) + class is (symba_pl) + select type(param) + class is (symba_parameters) + if (plplcollision_list%nenc == 0) return ! No collisions to resolve + ! Make sure that the heliocentric and barycentric coordinates are consistent with each other + call pl%vb2vh(system%cb) + call pl%xh2xb(system%cb) + + ! Get the energy before the collision is resolved + if (param%lenergy) then + call system%get_energy_and_momentum(param) + Eorbit_before = system%te + end if + + do + write(timestr,*) t + call io_log_one_message(FRAGGLE_LOG_OUT, "") + call io_log_one_message(FRAGGLE_LOG_OUT, "***********************************************************" // & + "***********************************************************") + call io_log_one_message(FRAGGLE_LOG_OUT, "Collision between massive bodies detected at time t = " // & + trim(adjustl(timestr))) + call io_log_one_message(FRAGGLE_LOG_OUT, "***********************************************************" // & + "***********************************************************") + allocate(tmp_param, source=param) + if (param%lfragmentation) then + call plplcollision_list%resolve_fragmentations(system, param) + else + call plplcollision_list%resolve_mergers(system, param) + end if + + ! Destroy the collision list now that the collisions are resolved + call plplcollision_list%setup(0_I8B) + + if ((system%pl_adds%nbody == 0) .and. (system%pl_discards%nbody == 0)) exit + + ! Save the add/discard information to file + call system%write_discard(tmp_param) + + ! Rearrange the arrays: Remove discarded bodies, add any new bodies, resort, and recompute all indices and encounter lists + call pl%rearray(system, tmp_param) + + ! Destroy the add/discard list so that we don't append the same body multiple times if another collision is detected + call system%pl_discards%setup(0, param) + call system%pl_adds%setup(0, param) + deallocate(tmp_param) + + ! Check whether or not any of the particles that were just added are themselves in a collision state. This will generate a new plplcollision_list + lplpl_collision = plplenc_list%collision_check(system, param, t, dt, irec) + + if (.not.lplpl_collision) exit + end do + + if (param%lenergy) then + call system%get_energy_and_momentum(param) + Eorbit_after = system%te + system%Ecollisions = system%Ecollisions + (Eorbit_after - Eorbit_before) + end if + + end select + end select + end associate + + return + end subroutine symba_collision_resolve_plplenc + + + module subroutine symba_collision_resolve_pltpenc(self, system, param, t, dt, irec) + !! author: David A. Minton + !! + !! Process the pl-tp collision list, then modifiy the massive bodies based on the outcome of the collision + !! + implicit none + ! Arguments + class(symba_pltpenc), intent(inout) :: self !! SyMBA pl-pl encounter list + class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA additions + real(DP), intent(in) :: t !! Current simulation tim + real(DP), intent(in) :: dt !! Current simulation step size + integer(I4B), intent(in) :: irec !! Current recursion level + + ! Make sure coordinate systems are all synced up due to being inside the recursion at this point + call system%pl%vb2vh(system%cb) + call system%tp%vb2vh(system%cb%vb) + call system%pl%b2h(system%cb) + call system%tp%b2h(system%cb) + + ! Discard the collider + call system%tp%discard(system, param) + + return + end subroutine symba_collision_resolve_pltpenc + +end submodule s_symba_collision \ No newline at end of file diff --git a/src/symba/symba_collision.f90 b/src/symba/symba_collision.f90 index 1423adc7e..e839af1de 100644 --- a/src/symba/symba_collision.f90 +++ b/src/symba/symba_collision.f90 @@ -761,7 +761,7 @@ subroutine symba_collision_mergeaddsub(system, param, colliders, frag, status) plnew%status(1:nfrag) = NEW_PARTICLE do i = 1, nfrag write(newname, FRAGFMT) frag%id(i) - call plnew%info(i)%set_value(origin_type="Disruption", origin_time=param%t, name=newname, & + call plnew%info(i)%set_value(origin_type="Disruption", origin_time=system%t, name=newname, & origin_xh=plnew%xh(:,i), & origin_vh=plnew%vh(:,i), collision_id=param%maxid_collision) end do @@ -771,14 +771,14 @@ subroutine symba_collision_mergeaddsub(system, param, colliders, frag, status) else iother = ibiggest end if - call pl%info(colliders%idx(i))%set_value(status="Disruption", discard_time=param%t, & + call pl%info(colliders%idx(i))%set_value(status="Disruption", discard_time=system%t, & discard_xh=pl%xh(:,i), discard_vh=pl%vh(:,i), discard_body_id=iother) end do case(SUPERCATASTROPHIC) plnew%status(1:nfrag) = NEW_PARTICLE do i = 1, nfrag write(newname, FRAGFMT) frag%id(i) - call plnew%info(i)%set_value(origin_type="Supercatastrophic", origin_time=param%t, name=newname, & + call plnew%info(i)%set_value(origin_type="Supercatastrophic", origin_time=system%t, name=newname, & origin_xh=plnew%xh(:,i), origin_vh=plnew%vh(:,i), & collision_id=param%maxid_collision) end do @@ -788,7 +788,7 @@ subroutine symba_collision_mergeaddsub(system, param, colliders, frag, status) else iother = ibiggest end if - call pl%info(colliders%idx(i))%set_value(status="Supercatastrophic", discard_time=param%t, & + call pl%info(colliders%idx(i))%set_value(status="Supercatastrophic", discard_time=system%t, & discard_xh=pl%xh(:,i), discard_vh=pl%vh(:,i), & discard_body_id=iother) end do @@ -797,14 +797,14 @@ subroutine symba_collision_mergeaddsub(system, param, colliders, frag, status) plnew%status(1) = OLD_PARTICLE do i = 2, nfrag write(newname, FRAGFMT) frag%id(i) - call plnew%info(i)%set_value(origin_type="Hit and run fragment", origin_time=param%t, name=newname, & + call plnew%info(i)%set_value(origin_type="Hit and run fragment", origin_time=system%t, name=newname, & origin_xh=plnew%xh(:,i), origin_vh=plnew%vh(:,i), & collision_id=param%maxid_collision) end do do i = 1, ncolliders if (colliders%idx(i) == ibiggest) cycle iother = ibiggest - call pl%info(colliders%idx(i))%set_value(status="Hit and run fragmention", discard_time=param%t, & + call pl%info(colliders%idx(i))%set_value(status="Hit and run fragmention", discard_time=system%t, & discard_xh=pl%xh(:,i), discard_vh=pl%vh(:,i), & discard_body_id=iother) end do @@ -815,7 +815,7 @@ subroutine symba_collision_mergeaddsub(system, param, colliders, frag, status) if (colliders%idx(i) == ibiggest) cycle iother = ibiggest - call pl%info(colliders%idx(i))%set_value(status="MERGED", discard_time=param%t, discard_xh=pl%xh(:,i), & + call pl%info(colliders%idx(i))%set_value(status="MERGED", discard_time=system%t, discard_xh=pl%xh(:,i), & discard_vh=pl%vh(:,i), discard_body_id=iother) end do end select @@ -1019,7 +1019,6 @@ module subroutine symba_collision_resolve_plplenc(self, system, param, t, dt, ir call io_log_one_message(FRAGGLE_LOG_OUT, "***********************************************************" // & "***********************************************************") allocate(tmp_param, source=param) - tmp_param%t = t if (param%lfragmentation) then call plplcollision_list%resolve_fragmentations(system, param) else diff --git a/src/symba/symba_discard.f90 b/src/symba/symba_discard.f90 index f00b222cb..76bcbaf41 100644 --- a/src/symba/symba_discard.f90 +++ b/src/symba/symba_discard.f90 @@ -44,7 +44,7 @@ subroutine symba_discard_cb_pl(pl, system, param) pl%lcollision(i) = .false. pl%status(i) = DISCARDED_RMAX write(idstr, *) pl%id(i) - write(timestr, *) param%t + write(timestr, *) system%t write(message, *) trim(adjustl(pl%info(i)%name)) // " (" // trim(adjustl(idstr)) // ")" // & " too far from the central body at t = " // trim(adjustl(timestr)) call io_log_one_message(FRAGGLE_LOG_OUT, "") @@ -54,14 +54,14 @@ subroutine symba_discard_cb_pl(pl, system, param) call io_log_one_message(FRAGGLE_LOG_OUT, "***********************************************************" // & "***********************************************************") call io_log_one_message(FRAGGLE_LOG_OUT, "") - call pl%info(i)%set_value(status="DISCARDED_RMAX", discard_time=param%t, discard_xh=pl%xh(:,i), & + call pl%info(i)%set_value(status="DISCARDED_RMAX", discard_time=system%t, discard_xh=pl%xh(:,i), & discard_vh=pl%vh(:,i)) else if ((param%rmin >= 0.0_DP) .and. (rh2 < rmin2)) then pl%ldiscard(i) = .true. pl%lcollision(i) = .false. pl%status(i) = DISCARDED_RMIN write(idstr, *) pl%id(i) - write(timestr, *) param%t + write(timestr, *) system%t write(message, *) trim(adjustl(pl%info(i)%name)) // " (" // trim(adjustl(idstr)) // ")" // & " too close to the central body at t = " // trim(adjustl(timestr)) call io_log_one_message(FRAGGLE_LOG_OUT, "") @@ -71,7 +71,7 @@ subroutine symba_discard_cb_pl(pl, system, param) call io_log_one_message(FRAGGLE_LOG_OUT, "************************************************************" // & "************************************************************") call io_log_one_message(FRAGGLE_LOG_OUT, "") - call pl%info(i)%set_value(status="DISCARDED_RMIN", discard_time=param%t, discard_xh=pl%xh(:,i), & + call pl%info(i)%set_value(status="DISCARDED_RMIN", discard_time=system%t, discard_xh=pl%xh(:,i), & discard_vh=pl%vh(:,i), discard_body_id=cb%id) else if (param%rmaxu >= 0.0_DP) then rb2 = dot_product(pl%xb(:,i), pl%xb(:,i)) @@ -82,7 +82,7 @@ subroutine symba_discard_cb_pl(pl, system, param) pl%lcollision(i) = .false. pl%status(i) = DISCARDED_RMAXU write(idstr, *) pl%id(i) - write(timestr, *) param%t + write(timestr, *) system%t write(message, *) trim(adjustl(pl%info(i)%name)) // " (" // trim(adjustl(idstr)) // ")" // & " is unbound and too far from barycenter at t = " // trim(adjustl(timestr)) call io_log_one_message(FRAGGLE_LOG_OUT, "") @@ -92,7 +92,7 @@ subroutine symba_discard_cb_pl(pl, system, param) call io_log_one_message(FRAGGLE_LOG_OUT, "************************************************************" // & "************************************************************") call io_log_one_message(FRAGGLE_LOG_OUT, "") - call pl%info(i)%set_value(status="DISCARDED_RMAXU", discard_time=param%t, discard_xh=pl%xh(:,i), & + call pl%info(i)%set_value(status="DISCARDED_RMAXU", discard_time=system%t, discard_xh=pl%xh(:,i), & discard_vh=pl%vh(:,i)) end if end if @@ -325,11 +325,11 @@ subroutine symba_discard_peri_pl(pl, system, param) pl%ldiscard(i) = .true. pl%lcollision(i) = .false. pl%status(i) = DISCARDED_PERI - write(timestr, *) param%t + write(timestr, *) system%t write(idstr, *) pl%id(i) write(*, *) trim(adjustl(pl%info(i)%name)) // " (" // trim(adjustl(idstr)) // & ") perihelion distance too small at t = " // trim(adjustl(timestr)) - call pl%info(i)%set_value(status="DISCARDED_PERI", discard_time=param%t, & + call pl%info(i)%set_value(status="DISCARDED_PERI", discard_time=system%t, & discard_xh=pl%xh(:,i), discard_vh=pl%vh(:,i), discard_body_id=system%cb%id) end if end if diff --git a/src/walltime/walltime.f90 b/src/walltime/walltime.f90 index 104b63d95..491a2e478 100644 --- a/src/walltime/walltime.f90 +++ b/src/walltime/walltime.f90 @@ -314,7 +314,6 @@ module subroutine walltime_interaction_time_this_loop(self, param, ninteractions end select self%is_on = .true. - write(tstr,*) param%t select case(self%stage) case(1) self%stage1_ninteractions = ninteractions