diff --git a/src/collision/collision_regime.f90 b/src/collision/collision_regime.f90 index d49202fa6..29f4f60c4 100644 --- a/src/collision/collision_regime.f90 +++ b/src/collision/collision_regime.f90 @@ -279,6 +279,8 @@ subroutine collision_regime_LS12_SI(Mcb, m1, m2, rad1, rad2, rh1, rh2, vb1, vb2, if (regime == COLLRESOLVE_REGIME_SUPERCATASTROPHIC) then Mlr = max(Mtot * 0.1_DP * (Qr / (Qrd_pstar * SUPERCAT_QRATIO))**(ETA), min_mfrag) !LS12 eq (44) + else if (regime == COLLRESOLVE_REGIME_HIT_AND_RUN) then + Mlr = m1 else Mlr = max((1.0_DP - Qr / Qrd_pstar / 2.0_DP) * Mtot, min_mfrag) ! [kg] # LS12 eq (5) end if diff --git a/src/fraggle/fraggle_generate.f90 b/src/fraggle/fraggle_generate.f90 index 7de02ef09..b4c6b6f57 100644 --- a/src/fraggle/fraggle_generate.f90 +++ b/src/fraggle/fraggle_generate.f90 @@ -52,6 +52,8 @@ module subroutine fraggle_generate(self, nbody_system, param, t) write(*,*) "Error in swiftest_collision, unrecognized collision regime" call base_util_exit(FAILURE) end select + call collision_io_collider_message(pl, impactors%id, message) + call swiftest_io_log_one_message(COLLISION_LOG_OUT, trim(adjustl(message))) call self%set_mass_dist(param) call self%disrupt(nbody_system, param, t, lfailure) if (lfailure) then @@ -266,12 +268,15 @@ module subroutine fraggle_generate_hitandrun(self, nbody_system, param, t) end if call self%set_mass_dist(param) message = "Hit and run between" + call collision_io_collider_message(nbody_system%pl, impactors%id, message) call swiftest_io_log_one_message(COLLISION_LOG_OUT, trim(adjustl(message))) - - ! Generate the position and velocity distributions of the fragments - call self%disrupt(nbody_system, param, t, lpure) + if (self%fragments%nbody > 2) then ! Hit and run with disruption + call self%disrupt(nbody_system, param, t, lpure) + else + lpure = .true. + end if if (lpure) then ! Disruption failed to find a solution. Convert to pure hit and run - call swiftest_io_log_one_message(COLLISION_LOG_OUT, "Disruptive hit and run failed to converge on a solution. Converting to pure hit and run. No new fragments generated.") + call swiftest_io_log_one_message(COLLISION_LOG_OUT, "Pure hit and run. No new fragments generated.") call self%collision_basic%hitandrun(nbody_system, param, t) return end if @@ -311,12 +316,13 @@ module subroutine fraggle_generate_pos_vec(collider, nbody_system, param, lfailu ! Internals real(DP) :: dis, direction, rdistance real(DP), dimension(NDIM,2) :: fragment_cloud_center + real(DP), dimension(NDIM) :: rwalk real(DP), dimension(2) :: fragment_cloud_radius logical, dimension(collider%fragments%nbody) :: loverlap real(DP), dimension(collider%fragments%nbody) :: mass_rscale, phi, theta, u integer(I4B) :: i, j, loop, istart logical :: lsupercat, lhitandrun - integer(I4B), parameter :: MAXLOOP = 200 + integer(I4B), parameter :: MAXLOOP = 100 real(DP), parameter :: rdistance_scale_factor = 1.0_DP ! Scale factor to apply to distance scaling of cloud centers in the event of overlap ! The distance is chosen to be close to the original locations of the impactors ! but far enough apart to prevent a collisional cascade between fragments @@ -335,33 +341,44 @@ module subroutine fraggle_generate_pos_vec(collider, nbody_system, param, lfailu ! Place the first two bodies at the centers of the two fragment clouds, but be sure they are sufficiently far apart to avoid overlap loverlap(:) = .true. if (lhitandrun) then - rdistance = impactors%radius(2) + rdistance = 1.0_DP else if (lsupercat) then rdistance = 0.5_DP * sum(impactors%radius(:)) else rdistance = 2 * impactors%radius(2) end if - ! Give the fragment positions a random value that is scaled with fragment mass so that the more massive bodies tend to be closer to the impact point - ! Later, velocities will be scaled such that the farther away a fragment is placed from the impact point, the higher will its velocity be. - call random_number(mass_rscale) - mass_rscale(:) = (mass_rscale(:) + 1.0_DP) / 2 - mass_rscale(:) = mass_rscale(:) * (fragments%mtot / fragments%mass(1:nfrag))**(0.125_DP) ! The power is arbitrary. It just gives the velocity a small mass dependence - mass_rscale(:) = mass_rscale(:) / maxval(mass_rscale(:)) + + if (lhitandrun) then + mass_rscale(:) = 1.0_DP + else + ! Give the fragment positions a random value that is scaled with fragment mass so that the more massive bodies tend to be closer to the impact point + ! Later, velocities will be scaled such that the farther away a fragment is placed from the impact point, the higher will its velocity be. + call random_number(mass_rscale) + mass_rscale(:) = (mass_rscale(:) + 1.0_DP) / 2 + mass_rscale(:) = mass_rscale(:) * (fragments%mtot / fragments%mass(1:nfrag))**(0.125_DP) ! The power is arbitrary. It just gives the velocity a small mass dependence + mass_rscale(:) = mass_rscale(:) / maxval(mass_rscale(:)) + end if istart = 3 do loop = 1, MAXLOOP if (.not.any(loverlap(:))) exit if (lhitandrun) then ! Keep the target unchanged and place the largest fragment at rdistance away from the projectile along its trajectory - fragment_cloud_radius(:) = impactors%radius(:) - fragment_cloud_center(:,1) = impactors%rc(:,1) - fragment_cloud_center(:,2) = impactors%rc(:,2) + rdistance * impactors%bounce_unit(:) + fragment_cloud_radius(1) = rbuffer * max(fragments%radius(1), impactors%radius(1)) + fragment_cloud_radius(2) = rbuffer * max(fragments%radius(2), impactors%radius(2)) + ! Initialize the largest body at the original target body position + fragments%rc(:,1) = impactors%rc(:,1) + ! Ensure that the second largest body does not overlap (including the buffer). Otherwise, shift it downrange + dis = max(1.00001_DP * sum(fragment_cloud_radius(1:2)) - .mag.(impactors%rc(:,2) - impactors%rc(:,1)), 0.0_DP) + fragments%rc(:,2) = impactors%rc(:,2) + dis * impactors%bounce_unit(:) + fragment_cloud_center(:,1) = fragments%rc(:,1) + fragment_cloud_center(:,2) = fragments%rc(:,2) + sum(fragment_cloud_radius(1:2)) * rdistance * impactors%bounce_unit(:) else ! Keep the first and second bodies at approximately their original location, but so as not to be overlapping fragment_cloud_center(:,1) = impactors%rcimp(:) - rbuffer * max(fragments%radius(1),impactors%radius(1)) * impactors%y_unit(:) fragment_cloud_center(:,2) = impactors%rcimp(:) + rbuffer * max(fragments%radius(2),impactors%radius(2)) * impactors%y_unit(:) fragment_cloud_radius(:) = cloud_size_scale_factor * rdistance + fragments%rc(:,1) = fragment_cloud_center(:,1) + fragments%rc(:,2) = fragment_cloud_center(:,2) end if - fragments%rc(:,1) = fragment_cloud_center(:,1) - fragments%rc(:,2) = fragment_cloud_center(:,2) do i = 1, nfrag if (loverlap(i)) then @@ -400,11 +417,11 @@ module subroutine fraggle_generate_pos_vec(collider, nbody_system, param, lfailu ! when the rest are not, we will randomly walk their position in space so as not to move them too far from their starting position if (all(.not.loverlap(istart:nfrag)) .and. any(loverlap(1:istart-1))) then do concurrent(i = 1:istart-1,loverlap(i)) - dis = 0.1_DP * fragments%radius(i) - fragments%rmag(i) = dis * u(i)**(THIRD) - fragments%rc(1,i) = fragments%rmag(i) * sin(theta(i)) * cos(phi(i)) - fragments%rc(2,i) = fragments%rmag(i) * sin(theta(i)) * sin(phi(i)) - fragments%rc(3,i) = fragments%rmag(i) * cos(theta(i)) + dis = 0.1_DP * fragments%radius(i) * u(i)**(THIRD) + rwalk(1) = fragments%rmag(i) * sin(theta(i)) * cos(phi(i)) + rwalk(2) = fragments%rmag(i) * sin(theta(i)) * sin(phi(i)) + rwalk(3) = fragments%rmag(i) * cos(theta(i)) + fragments%rc(:,i) = fragments%rc(:,i) + rwalk(:) end do end if @@ -417,7 +434,6 @@ module subroutine fraggle_generate_pos_vec(collider, nbody_system, param, lfailu do i = j + 1, nfrag dis = .mag.(fragments%rc(:,j) - fragments%rc(:,i)) loverlap(i) = loverlap(i) .or. (dis <= rbuffer * (fragments%radius(i) + fragments%radius(j))) - loverlap(j) = loverlap(j) .or. (dis <= rbuffer * (fragments%radius(i) + fragments%radius(j))) end do ! Check for overlaps with existing bodies that are not involved in the collision do i = 1, npl @@ -455,14 +471,41 @@ module subroutine fraggle_generate_rot_vec(collider, nbody_system, param) class(collision_fraggle), intent(inout) :: collider !! Fraggle collision system object class(swiftest_nbody_system), intent(inout) :: nbody_system !! Swiftest nbody system object class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters + ! Internals + integer(I4B) :: i + real(DP), parameter :: frag_rot_fac = 0.1_DP ! Fraction of projectile rotation magnitude to add as random noise to fragment rotation + real(DP) :: mass_fac + real(DP), dimension(NDIM) :: drot, dL + integer(I4B), parameter :: MAXLOOP = 10 associate(fragments => collider%fragments, impactors => collider%impactors, nfrag => collider%fragments%nbody) ! Initialize fragment rotations and velocities to be pre-impact rotation for body 1, and randomized for bodies >1 and scaled to the original rotation. ! This will get updated later when conserving angular momentum - fragments%rot(:,1) = impactors%rot(:,1) + mass_fac = fragments%mass(1) / impactors%mass(1) + fragments%rot(:,1) = mass_fac**(5.0_DP/3.0_DP) * impactors%rot(:,1) + + ! If mass was added, also add spin angular momentum + if (mass_fac > 1.0_DP) then + dL(:) = (fragments%mass(1) - impactors%mass(1)) * (impactors%rc(:,2) - impactors%rc(:,1)) .cross. (impactors%vc(:,2) - impactors%vc(:,1)) + drot(:) = dL(:) / (fragments%mass(1) * fragments%radius(1)**2 * fragments%Ip(3,1)) + ! Check to make sure we haven't broken the spin barrier. Reduce the rotation change if so + do i = 1, MAXLOOP + if (.mag.(fragments%rot(:,1) + drot(:)) < collider%max_rot) exit + if (i == MAXLOOP) drot(:) = 0.0_DP + where(drot(:) > TINY(1.0_DP)) + drot(:) = drot(:) / 2 + elsewhere + drot(:) = 0.0_DP + endwhere + end do + fragments%rot(:,1) = fragments%rot(:,1) + drot(:) + end if call random_number(fragments%rot(:,2:nfrag)) - fragments%rot(:,2:nfrag) = fragments%rot(:,2:nfrag) * .mag.impactors%rot(:,2) + do concurrent (i = 2:nfrag) + mass_fac = fragments%mass(i) / impactors%mass(2) + fragments%rot(:,i) = mass_fac**(5.0_DP/3.0_DP) * impactors%rot(:,2) + 2 * (fragments%rot(:,i) - 1.0_DP) * frag_rot_fac * .mag.impactors%rot(:,2) + end do fragments%rotmag(:) = .mag.fragments%rot(:,:) end associate @@ -489,11 +532,13 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu integer(I4B) :: i, j, loop, try, istart, nfrag, nsteps, nsteps_best logical :: lhitandrun, lsupercat real(DP), dimension(NDIM) :: vimp_unit, rimp, vrot, L_residual, L_residual_unit, L_residual_best, dL, drot, rot_new, dL_metric - real(DP) :: vimp, vmag, vesc, dE, E_residual, E_residual_best, E_residual_last, ke_min, ke_avail, ke_remove, dE_best, fscale, dE_metric, mfrag, rn + real(DP) :: vimp, vmag, vesc, dE, E_residual, E_residual_best, E_residual_last, ke_avail, ke_remove, dE_best, fscale, dE_metric, mfrag, rn, dL1_mag integer(I4B), dimension(collider%fragments%nbody) :: vsign real(DP), dimension(collider%fragments%nbody) :: vscale - real(DP), parameter :: L_ROT_VEL_RATIO = 0.2_DP ! Ratio of angular momentum to put into rotation relative to velocity shear of fragments + real(DP), dimension(collider%fragments%nbody) :: dLi_mag + real(DP), parameter :: L_ROT_VEL_RATIO = 0.9_DP ! Ratio of angular momentum to put into rotation relative to velocity shear of fragments ! For the initial "guess" of fragment velocities, this is the minimum and maximum velocity relative to escape velocity that the fragments will have + real(DP), parameter :: hitandrun_vscale = 0.25_DP real(DP) :: vmin_guess real(DP) :: vmax_guess real(DP) :: delta_v, GC @@ -529,16 +574,17 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu end if ! The minimum fragment velocity will be set by the escape velocity + vimp = .mag. (impactors%vc(:,2) - impactors%vc(:,1)) + vesc = sqrt(2 * sum(fragments%Gmass(istart:nfrag)) / sum(fragments%radius(istart:nfrag))) if (lhitandrun) then - vesc = sqrt(2 * impactors%Gmass(2) / impactors%radius(2)) + vmin_guess = .mag.impactors%vc(:,2) - vimp * hitandrun_vscale + vmax_guess = .mag.impactors%vc(:,2) + vimp * hitandrun_vscale else vesc = sqrt(2 * sum(impactors%Gmass(:)) / sum(impactors%radius(:))) + vmin_guess = 1.001_DP * vesc + vmax_guess = vimp end if - vimp = .mag. (impactors%vc(:,2) - impactors%vc(:,1)) - vmax_guess = 1.1_DP * vimp - vmin_guess = 1.001_DP * vesc - E_residual_best = huge(1.0_DP) lfailure = .false. dE_metric = huge(1.0_DP) @@ -552,23 +598,23 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu rimp(:) = fragments%rc(:,i) - impactors%rcimp(:) vscale(i) = .mag. rimp(:) / (.mag. (impactors%rb(:,2) - impactors%rb(:,1))) end do - vscale(1) = 0.9_DP * minval(vscale(2:nfrag)) ! Set the velocity scale factor to span from vmin/vesc to vmax/vesc - vscale(:) = vscale(:)/minval(vscale(:)) - fscale = log(vmax_guess - vmin_guess + 1.0_DP) / log(maxval(vscale(:))) - vscale(:) = vscale(:)**fscale + vmin_guess - 1.0_DP + if (nfrag == 2) then + vscale(:) = 1.0_DP + else + vscale(2:nfrag) = vscale(2:nfrag)/minval(vscale(2:nfrag)) + fscale = log(vmax_guess - vmin_guess + 1.0_DP) / log(maxval(vscale(2:nfrag))) + vscale(2:nfrag) = vscale(2:nfrag)**fscale + vmin_guess - 1.0_DP + end if ! Set the velocities of all fragments using all of the scale factors determined above - do concurrent(i = 1:nfrag) + fragments%vc(:,1) = impactors%vc(:,1) * impactors%mass(1) / fragments%mass(1) + do concurrent(i = 2:nfrag) j = fragments%origin_body(i) vrot(:) = impactors%rot(:,j) .cross. (fragments%rc(:,i) - impactors%rc(:,j)) if (lhitandrun) then - if (i ==1) then - fragments%vc(:,i) = impactors%vc(:,1) - else - fragments%vc(:,i) = vsign(i) * impactors%bounce_unit(:) * vscale(i) - end if + fragments%vc(:,i) = vsign(i) * impactors%bounce_unit(:) * vscale(i) + vrot(:) else vmag = vscale(i) rimp(:) = fragments%rc(:,i) - impactors%rcimp(:) @@ -581,48 +627,59 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu ! Every time the collision-frame velocities are altered, we need to be sure to shift everything back to the center-of-mass frame call collision_util_shift_vector_to_origin(fragments%mass, fragments%vc) call fragments%set_coordinate_system() - ke_min = 0.5_DP * fragments%mtot * vesc**2 E_residual = huge(1.0_DP) inner: do loop = 1, MAXINNER nsteps = nsteps + 1 + mfrag = sum(fragments%mass(istart:nfrag)) ! Try to put residual angular momentum into the spin, but if this would go past the spin barrier, then put it into velocity shear instead + call collider_local%get_energy_and_momentum(nbody_system, param, phase="after") + L_residual(:) = (collider_local%L_total(:,2) - collider_local%L_total(:,1)) + L_residual_unit(:) = .unit. L_residual(:) + + ! Use equipartition of spin kinetic energy to distribution spin angular momentum + do concurrent(i = istart:nfrag) + dLi_mag(i) = ((fragments%mass(i) / fragments%mass(istart)) * & + (fragments%radius(i) / fragments%radius(istart))**2 * & + (fragments%Ip(3,i) / fragments%Ip(3,istart)))**(1.5_DP) + end do + dL1_mag = .mag.L_residual(:) / sum(dLi_mag(istart:nfrag)) + + do i = istart,nfrag + dL(:) = -dL1_mag * dLi_mag(i) * L_residual_unit(:) + drot(:) = L_ROT_VEL_RATIO * dL(:) / (fragments%mass(i) * fragments%Ip(3,i) * fragments%radius(i)**2) + rot_new(:) = fragments%rot(:,i) + drot(:) + if (.mag.rot_new(:) < collider_local%max_rot) then + fragments%rot(:,i) = rot_new(:) + fragments%rotmag(i) = .mag.fragments%rot(:,i) + else ! We would break the spin barrier here. Put less into spin and more into velocity shear. + call random_number(drot) + call random_number(rn) + drot(:) = (rn * collider_local%max_rot - fragments%rotmag(i)) * 2 * (drot(:) - 0.5_DP) + fragments%rot(:,i) = fragments%rot(:,i) + drot(:) + fragments%rotmag(i) = .mag.fragments%rot(:,i) + if (fragments%rotmag(i) > collider%max_rot) then + fragments%rotmag(i) = 0.5_DP * collider%max_rot + fragments%rot(:,i) = fragments%rotmag(i) * .unit. fragments%rot(:,i) + end if + end if + L_residual(:) = L_residual(:) - drot(:) * fragments%Ip(3,i) * fragments%mass(i) * fragments%radius(i)**2 + end do + + ! Put any remaining residual into velocity shear angmtm: do j = 1, MAXANGMTM call collider_local%get_energy_and_momentum(nbody_system, param, phase="after") L_residual(:) = (collider_local%L_total(:,2) - collider_local%L_total(:,1)) dL_metric(:) = abs(L_residual(:)) / .mag.(collider_local%L_total(:,1)) / MOMENTUM_SUCCESS_METRIC if (all(dL_metric(:) <= 1.0_DP)) exit angmtm - L_residual_unit(:) = .unit. L_residual(:) - do i = fragments%nbody,1,-1 - mfrag = sum(fragments%mass(i:fragments%nbody)) - drot(:) = -L_ROT_VEL_RATIO * L_residual(:) / (fragments%mtot * fragments%Ip(3,i) * fragments%radius(i)**2) - rot_new(:) = fragments%rot(:,i) + drot(:) - if (.mag.rot_new(:) < collider_local%max_rot) then - fragments%rot(:,i) = rot_new(:) - fragments%rotmag(i) = .mag.fragments%rot(:,i) - else ! We would break the spin barrier here. Put less into spin and more into velocity shear. - if (i >= istart) then - call random_number(drot) - call random_number(rn) - drot(:) = (rn * collider_local%max_rot - fragments%rotmag(i)) * 2 * (drot(:) - 0.5_DP) - fragments%rot(:,i) = fragments%rot(:,i) + drot(:) - fragments%rotmag(i) = .mag.fragments%rot(:,i) - if (fragments%rotmag(i) > collider%max_rot) then - fragments%rotmag(i) = 0.5_DP * collider%max_rot - fragments%rot(:,i) = fragments%rotmag(i) * .unit. fragments%rot(:,i) - end if - else - drot(:) = 0.0_DP - end if - end if - - dL(:) = -L_residual(:) * fragments%mass(i) / fragments%mtot - drot(:) * fragments%Ip(3,i) * fragments%mass(i) * fragments%radius(i)**2 + + do i = istart, nfrag + dL(:) = -L_residual(:) * fragments%mass(i) / sum(fragments%mass(istart:nfrag)) call collision_util_velocity_torque(dL, fragments%mass(i), fragments%rc(:,i), fragments%vc(:,i)) call collision_util_shift_vector_to_origin(fragments%mass, fragments%vc) fragments%vmag(i) = .mag.fragments%vc(:,i) - end do end do angmtm diff --git a/src/fraggle/fraggle_util.f90 b/src/fraggle/fraggle_util.f90 index e5956933f..e8647d589 100644 --- a/src/fraggle/fraggle_util.f90 +++ b/src/fraggle/fraggle_util.f90 @@ -30,7 +30,7 @@ module subroutine fraggle_util_set_mass_dist(self, param) real(DP), parameter :: BETA = 2.85_DP integer(I4B), parameter :: MASS_NOISE_FACTOR = 4 !! The number of digits of random noise that get added to the minimum mass value to prevent identical masses from being generated in a single run integer(I4B), parameter :: NFRAGMAX = 100 !! Maximum number of fragments that can be generated - integer(I4B), parameter :: NFRAGMIN = 7 !! Minimum number of fragments that can be generated (set by the fraggle_generate algorithm for constraining momentum and energy) + integer(I4B), parameter :: NFRAGMIN = 1 !! Minimum number of fragments that can be generated (set by the fraggle_generate algorithm for constraining momentum and energy) integer(I4B), parameter :: NFRAG_SIZE_MULTIPLIER = 3 !! Log-space scale factor that scales the number of fragments by the collisional system mass integer(I4B), parameter :: iMlr = 1 integer(I4B), parameter :: iMslr = 2 @@ -76,15 +76,19 @@ module subroutine fraggle_util_set_mass_dist(self, param) nfrag = NFRAGMAX end select - i = iMrem mremaining = impactors%mass_dist(iMrem) - do while (i <= nfrag) - mfrag = (1 + i - iMslr)**(-3._DP / BETA) * impactors%mass_dist(iMslr) - if (mremaining - mfrag <= 0.0_DP) exit - mremaining = mremaining - mfrag - i = i + 1 - end do - if (i < nfrag) nfrag = max(i, NFRAGMIN) ! The sfd would actually give us fewer fragments than our maximum + if (mremaining > 0.0_DP) then + i = iMrem + do while (i <= nfrag) + mfrag = max((1 + i - iMslr)**(-3._DP / BETA) * impactors%mass_dist(iMslr), min_mfrag) + if (mremaining - mfrag <= 0.0_DP) exit + mremaining = mremaining - mfrag + i = i + 1 + end do + else + i = iMrem - 1 + end if + nfrag = i call self%setup_fragments(nfrag) case (COLLRESOLVE_REGIME_MERGE, COLLRESOLVE_REGIME_GRAZE_AND_MERGE) @@ -113,7 +117,9 @@ module subroutine fraggle_util_set_mass_dist(self, param) ! Distribute the remaining mass the 3:nfrag bodies following the model SFD given by slope BETA mremaining = impactors%mass_dist(iMrem) do i = iMrem, nfrag - mfrag = (1 + i - iMslr)**(-3._DP / BETA) * impactors%mass_dist(iMslr) + call random_number(mass_noise) + mass_noise = 1.0_DP + mass_noise * epsilon(1.0_DP) * 10**(MASS_NOISE_FACTOR) + mfrag = max((1 + i - iMslr)**(-3._DP / BETA) * impactors%mass_dist(iMslr), min_mfrag) * mass_noise mass(i) = mfrag mremaining = mremaining - mfrag end do @@ -134,8 +140,7 @@ module subroutine fraggle_util_set_mass_dist(self, param) ! Sort the distribution in descending order by mass so that the largest fragment is always the first call swiftest_util_sort(-mass, ind) call swiftest_util_sort_rearrange(mass, ind, nfrag) - fragments%mass(:) = mass(:) - deallocate(mass) + call move_alloc(mass, fragments%mass) fragments%Gmass(:) = G * fragments%mass(:)