diff --git a/src/modules/swiftest_classes.f90 b/src/modules/swiftest_classes.f90 index af751fdcf..1aec7dc9d 100644 --- a/src/modules/swiftest_classes.f90 +++ b/src/modules/swiftest_classes.f90 @@ -289,6 +289,7 @@ module swiftest_classes logical :: lbeg !! True if this is the beginning of a step. This is used so that test particle steps can be calculated !! separately from massive bodies. Massive body variables are saved at half steps, and passed to !! the test particles + integer(I4B) :: maxid = -1 !! The current maximum particle id number contains !> Each integrator will have its own version of the step procedure(abstract_step_system), deferred :: step diff --git a/src/setup/setup.f90 b/src/setup/setup.f90 index 33b5c07c4..6cba6d27b 100644 --- a/src/setup/setup.f90 +++ b/src/setup/setup.f90 @@ -128,6 +128,7 @@ module subroutine setup_initialize_system(self, param) call self%pl%initialize(param) call self%tp%initialize(param) call util_valid(self%pl, self%tp) + self%maxid = maxval([self%pl%id(:), self%tp%id(:)]) call self%set_msys() call self%pl%set_mu(self%cb) call self%tp%set_mu(self%cb) diff --git a/src/symba/symba_fragmentation.f90 b/src/symba/symba_fragmentation.f90 index 3ad475e4d..fea7dea91 100644 --- a/src/symba/symba_fragmentation.f90 +++ b/src/symba/symba_fragmentation.f90 @@ -1,5 +1,8 @@ submodule (symba_classes) s_symba_fragmentation use swiftest + + integer(I4B), parameter :: NFRAG_DISRUPT = 12 + integer(I4B), parameter :: NFRAG_SUPERCAT = 20 contains module function symba_fragmentation_casedisruption(system, param, family, x, v, mass, radius, L_spin, Ip, mass_res, Qloss) result(status) @@ -19,8 +22,114 @@ module function symba_fragmentation_casedisruption(system, param, family, x, v, ! Result integer(I4B) :: status !! Status flag assigned to this outcome ! Internals + integer(I4B) :: i, istart, nfrag, ibiggest, nstart, nend + real(DP) :: mtot, avg_dens + real(DP), dimension(NDIM) :: xcom, vcom, Ip_new + real(DP), dimension(2) :: vol + real(DP), dimension(:, :), allocatable :: vb_frag, xb_frag, rot_frag, Ip_frag + real(DP), dimension(:), allocatable :: m_frag, rad_frag + logical :: lfailure + class(symba_pl), allocatable :: plnew + + select type(pl => system%pl) + class is (symba_pl) + associate(mergeadd_list => system%mergeadd_list, mergesub_list => system%mergesub_list, cb => system%cb) + ! Collisional fragments will be uniformly distributed around the pre-impact barycenter + nfrag = NFRAG_DISRUPT + allocate(m_frag(nfrag)) + allocate(rad_frag(nfrag)) + allocate(xb_frag(NDIM, nfrag)) + allocate(vb_frag(NDIM, nfrag)) + allocate(rot_frag(NDIM, nfrag)) + allocate(Ip_frag(NDIM, nfrag)) + + mtot = sum(mass(:)) + xcom(:) = (mass(1) * x(:,1) + mass(2) * x(:,2)) / mtot + vcom(:) = (mass(1) * v(:,1) + mass(2) * v(:,2)) / mtot + + ! Get mass weighted mean of Ip and average density + Ip_new(:) = (mass(1) * Ip(:,1) + mass(2) * Ip(:,2)) / mtot + vol(:) = 4._DP / 3._DP * PI * radius(:)**3 + avg_dens = mtot / sum(vol(:)) + + ! Distribute the mass among fragments, with a branch to check for the size of the second largest fragment + m_frag(1) = mass_res(1) + if (mass_res(2) > mass_res(1) / 3._DP) then + m_frag(2) = mass_res(2) + istart = 3 + else + istart = 2 + end if + ! Distribute remaining mass among the remaining bodies + do i = istart, nfrag + m_frag(i) = (mtot - sum(m_frag(1:istart - 1))) / (nfrag - istart + 1) + end do + + ! Distribute any residual mass if there is any and set the radius + m_frag(nfrag) = m_frag(nfrag) + (mtot - sum(m_frag(:))) + rad_frag(:) = (3 * m_frag(:) / (4 * PI * avg_dens))**(1.0_DP / 3.0_DP) + + do i = 1, nfrag + Ip_frag(:, i) = Ip_new(:) + end do + + !call fragmentation_initialize(pl, param, family, x, v, L_spin, Ip, mass, radius, & + ! nfrag, Ip_frag, m_frag, rad_frag, xb_frag, vb_frag, rot_frag, Qloss, lfailure) + + if (lfailure) then + write(*,*) 'No fragment solution found, so treat as a pure hit-and-run' + status = ACTIVE + nfrag = 0 + else + ! Populate the list of new bodies + write(*,'("Generating ",I2.0," fragments")') nfrag + allocate(plnew, mold=pl) + call plnew%setup(nfrag, param) + + plnew%id(:) = [(i, i = system%maxid + 1, system%maxid + nfrag)] + system%maxid = system%maxid + nfrag + plnew%status(:) = ACTIVE + plnew%lcollision(:) = .false. + plnew%ldiscard(:) = .false. + plnew%xb(:,:) = xb_frag(:, :) + plnew%vb(:,:) = vb_frag(:, :) + do i = 1, nfrag + plnew%xh(:,i) = xb_frag(:, i) - cb%xb(:) + plnew%vh(:,i) = vb_frag(:, i) - cb%vb(:) + end do + plnew%mass(:) = m_frag(:) + plnew%Gmass(:) = param%GU * m_frag(:) + plnew%density(:) = avg_dens + plnew%radius(:) = rad_frag(:) + plnew%info(:)%origin_type = "Disruption" + plnew%info(:)%origin_time = param%t + do i = 1, nfrag + plnew%info(i)%origin_xh(:) = plnew%xh(:,i) + plnew%info(i)%origin_vh(:) = plnew%vh(:,i) + end do + if (param%lrotation) then + plnew%Ip(:,:) = Ip_frag(:,:) + plnew%rot(:,:) = rot_frag(:,:) + end if + if (param%ltides) then + ibiggest = maxloc(pl%Gmass(family(:)), dim=1) + plnew%Q = pl%Q(ibiggest) + plnew%k2 = pl%k2(ibiggest) + plnew%tlag = pl%tlag(ibiggest) + end if + + ! Append the new merged body to the list and record how many we made + nstart = mergeadd_list%nbody + 1 + nend = mergeadd_list%nbody + plnew%nbody + call mergeadd_list%append(plnew) + mergeadd_list%ncomp(nstart:nend) = plnew%nbody + + call plnew%setup(0, param) + deallocate(plnew) + end if - status = DISRUPTION + end associate + end select return end function symba_fragmentation_casedisruption @@ -157,8 +266,8 @@ module function symba_fragmentation_casemerge(system, param, family, x, v, mass, plnew%radius(1) = radius_new plnew%info(1) = pl%info(family(ibiggest)) if (param%lrotation) then - pl%Ip(:,1) = Ip_new(:) - pl%rot(:,1) = rot_new(:) + plnew%Ip(:,1) = Ip_new(:) + plnew%rot(:,1) = rot_new(:) end if if (param%ltides) then plnew%Q = pl%Q(ibiggest)