diff --git a/src/collision/collision_resolve.f90 b/src/collision/collision_resolve.f90 index b50b21334..525b53969 100644 --- a/src/collision/collision_resolve.f90 +++ b/src/collision/collision_resolve.f90 @@ -152,9 +152,6 @@ module subroutine collision_resolve_consolidate_impactors(self, nbody_system, pa ! Destroy the kinship relationships for all members of this impactors%id call pl%reset_kinship(impactors%id(:)) - ! Set the coordinate system of the impactors - call impactors%set_coordinate_system() - end associate end select return diff --git a/src/fraggle/fraggle_generate.f90 b/src/fraggle/fraggle_generate.f90 index 3b70fe70c..bcf7e611e 100644 --- a/src/fraggle/fraggle_generate.f90 +++ b/src/fraggle/fraggle_generate.f90 @@ -35,6 +35,9 @@ module subroutine fraggle_generate(self, nbody_system, param, t) select type(pl => nbody_system%pl) class is (swiftest_pl) associate(impactors => self%impactors, status => self%status, maxid => nbody_system%maxid) + ! Set the coordinate system of the impactors + call impactors%set_coordinate_system() + select case (impactors%regime) case (COLLRESOLVE_REGIME_HIT_AND_RUN) call self%hitandrun(nbody_system, param, t) @@ -407,7 +410,7 @@ module subroutine fraggle_generate_rot_vec(collider, nbody_system, param) ! Internals real(DP), dimension(NDIM) :: Lbefore, Lafter, L_spin, rotdir real(DP) :: v_init, v_final, mass_init, mass_final, rotmag, dKE, KE_init, KE_final - real(DP), parameter :: random_scale_factor = 0.001_DP !! The relative scale factor to apply to the random component of the rotation vector + real(DP), parameter :: random_scale_factor = 0.01_DP !! The relative scale factor to apply to the random component of the rotation vector integer(I4B) :: i associate(fragments => collider%fragments, impactors => collider%impactors, nfrag => collider%fragments%nbody) @@ -438,8 +441,6 @@ module subroutine fraggle_generate_rot_vec(collider, nbody_system, param) Lafter(:) = mass_final * (impactors%rb(:,2) - impactors%rb(:,1)) .cross. (v_final * impactors%bounce_unit(:)) L_spin(:) = impactors%L_spin(:,1) + random_scale_factor * (Lbefore(:) - Lafter(:)) - !fragments%rot(:,1) = L_spin(:) / (fragments%mass(1) * fragments%radius(1)**2 * fragments%Ip(3,1)) - !fragments%rotmag(1) = .mag.fragments%rot(:,1) ! Add in some random spin noise. The magnitude will be scaled by the before-after amount and the direction will be random do i = 2,nfrag call random_number(rotdir) @@ -447,7 +448,7 @@ module subroutine fraggle_generate_rot_vec(collider, nbody_system, param) rotdir = rotdir - 0.5_DP rotdir = .unit. rotdir fragments%rotmag(i) = random_scale_factor * rotmag * .mag.L_spin(:) / ((nfrag - 1) * fragments%mass(i) * fragments%radius(i)**2 * fragments%Ip(3,i)) - fragments%rot(:,i) = fragments%rotmag(i) * rotdir + fragments%rot(:,i) = fragments%rot(:,i) + fragments%rotmag(i) * rotdir end do end associate @@ -470,16 +471,16 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu ! Internals integer(I4B) :: i, j, loop, try, istart, nfrag, nsteps logical :: lhitandrun, lsupercat - real(DP), dimension(NDIM) :: vimp_unit, rimp, vrot, L_residual - real(DP) :: vmag, vesc, dE, E_residual, ke_min, ke_avail, ke_remove, dE_best, E_residual_best, fscale, f_spin, f_orbit, dE_metric, dM + real(DP), dimension(NDIM) :: vimp_unit, rimp, vrot, L_residual, vdir, L_residual_unit, Li, Lrat, r_lever + real(DP) :: rmag, vmag, vesc, dE, E_residual, ke_min, ke_avail, ke_remove, dE_best, E_residual_best, fscale, f_spin, f_orbit, dE_metric, dM, mfrag integer(I4B), dimension(collider%fragments%nbody) :: vsign real(DP), dimension(collider%fragments%nbody) :: vscale, ke_rot_remove, volume ! 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) :: vmin_guess = 1.01_DP real(DP) :: vmax_guess = 8.0_DP real(DP) :: delta_v, GC - integer(I4B), parameter :: MAXLOOP = 100 - integer(I4B), parameter :: MAXTRY = 20 + integer(I4B), parameter :: MAXLOOP = 50 + integer(I4B), parameter :: MAXTRY = 50 real(DP), parameter :: mass_reduction_ratio = 0.1_DP ! Ratio of difference between first and second fragment mass to remove from the largest fragment in case of a failure real(DP), parameter :: SUCCESS_METRIC = 1.0e-2_DP class(collision_fraggle), allocatable :: collider_local @@ -506,7 +507,7 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu vsign(:) = 1 end where - ! Hit and run collisions should only affect the runner + ! Hit and run collisions should only affect the runners, not the target if (lhitandrun) then istart = 2 else @@ -546,7 +547,7 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu if (i ==1) then fragments%vc(:,i) = impactors%vc(:,1) else - fragments%vc(:,i) = impactors%vc(:,2) + vsign(i) * impactors%bounce_unit(:) * (vscale(i) - 0.5_DP *(maxval(vscale(:)) - minval(vscale(:)))) + fragments%vc(:,i) = impactors%vc(:,2) + vsign(i) * impactors%bounce_unit(:) * vscale(i) end if else vmag = vesc * vscale(i) @@ -569,17 +570,44 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu ke_avail = ke_avail + 0.5_DP * fragments%mass(i) * max(fragments%vmag(i) - vesc,0.0_DP)**2 end do - ! Check for any residual angular momentum, and if there is any, put it into spin of the largest body + ! Check for any residual angular momentum, and if there is any, put it into spin/shear depending on the type of collision (hit and runs prioritize velocity shear over spin) L_residual(:) = collider_local%L_total(:,2) - collider_local%L_total(:,1) - if (ke_avail < epsilon(1.0_DP)) then + if (lhitandrun .or. (ke_avail < epsilon(1.0_DP))) then + ! Start by putting residual angular momentum into velocity shear + mfrag = sum(fragments%mass(istart:fragments%nbody)) + L_residual_unit(:) = .unit. L_residual(:) + fragments%r_unit(:,:) = .unit.fragments%rc(:,:) do i = 1, fragments%nbody + r_lever(:) = (L_residual_unit(:) .cross. fragments%r_unit(:,i)) + rmag = .mag.r_lever(:) + if (rmag > epsilon(1.0_DP)) then + vdir(:) = -.unit. r_lever(:) + vmag = .mag.L_residual(:) / (fragments%mtot * .mag.r_lever(:) * fragments%rmag(i)) + Li(:) = fragments%mass(i) * fragments%rc(:,i) .cross. (vmag * vdir(:)) + Lrat(:) = L_residual(:) / Li(:) + fragments%vc(:,i) = fragments%vc(:,i) + vmag * vdir(:) + end if + end do + fragments%vmag(:) = .mag.fragments%vc(:,:) + + ! Update the coordinate system now that the velocities have changed + call collision_util_shift_vector_to_origin(fragments%mass, fragments%vc) + call fragments%set_coordinate_system() + + ! Try to put as much of the residual angular momentum into the spin of the fragments before the target body + call collider_local%get_energy_and_momentum(nbody_system, param, phase="after") + L_residual(:) = (collider_local%L_total(:,2) - collider_local%L_total(:,1)) + do i = istart, fragments%nbody fragments%L_spin(:,i) = fragments%L_spin(:,i) - L_residual(:) * fragments%mass(i) / fragments%mtot fragments%rot(:,i) = fragments%L_spin(:,i) / (fragments%mass(i) * fragments%radius(i)**2 * fragments%Ip(:,i)) end do - else - fragments%L_spin(:,1) = fragments%L_spin(:,1) - L_residual(:) - fragments%rot(:,1) = fragments%L_spin(:,1) / (fragments%mass(1) * fragments%radius(1)**2 * fragments%Ip(:,1)) + call collider_local%get_energy_and_momentum(nbody_system, param, phase="after") + L_residual(:) = (collider_local%L_total(:,2) - collider_local%L_total(:,1)) end if + + ! Put any remaining residual angular momentum into the spin of the target body + fragments%L_spin(:,1) = fragments%L_spin(:,1) - L_residual(:) + fragments%rot(:,1) = fragments%L_spin(:,1) / (fragments%mass(1) * fragments%radius(1)**2 * fragments%Ip(:,1)) fragments%rotmag(:) = .mag.fragments%rot(:,:) call collider_local%get_energy_and_momentum(nbody_system, param, phase="after") @@ -587,6 +615,7 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu dE = collider_local%te(2) - collider_local%te(1) E_residual = dE + impactors%Qloss + ! Check if we've converged on our constraints if ((abs(E_residual) < abs(E_residual_best)) .or. ((dE < 0.0_DP) .and. (dE_best >= 0.0_DP))) then ! This is our best case so far. Save it for posterity E_residual_best = E_residual dE_best = dE @@ -601,6 +630,8 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu end if if ((dE_best < 0.0_DP) .and. (dE_metric <= SUCCESS_METRIC * try)) exit outer ! As the tries increase, we relax the success metric. What was once a failure might become a success + ! If we have an offset of energy from our target, we'll remove it from both orbit and spin, keeping the proportion of each source roughly constant + ! Remove a constant amount of velocity from the bodies so we don't shift the center of mass and screw up the momentum f_spin = (fragments%ke_spin_tot )/ (fragments%ke_spin_tot + fragments%ke_orbit_tot) f_orbit = 1.0_DP - f_spin @@ -614,7 +645,7 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu ke_remove = min(f_spin * E_residual, 0.9_DP*fragments%ke_spin_tot) ke_rot_remove(:) = ke_remove * (fragments%ke_spin(:) / fragments%ke_spin_tot) where(ke_rot_remove(:) > fragments%ke_spin(:)) ke_rot_remove(:) = fragments%ke_spin(:) - do concurrent(i = 1:fragments%nbody, fragments%ke_spin(i) > 10*sqrt(tiny(1.0_DP))) + do concurrent(i = istart:fragments%nbody, fragments%ke_spin(i) > 10*sqrt(tiny(1.0_DP))) fscale = sqrt((fragments%ke_spin(i) - ke_rot_remove(i))/fragments%ke_spin(i)) fragments%rotmag(i) = fscale * fragments%rotmag(i) fragments%rot(:,i) = fscale * fragments%rot(:,i) @@ -648,7 +679,7 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu call fraggle_generate_rot_vec(collider_local, nbody_system, param) ! Increase the spatial size factor to get a less dense cloud - collider_local%fail_scale = collider_local%fail_scale*1.01_DP + collider_local%fail_scale = collider_local%fail_scale*1.001_DP ! Bring the minimum and maximum velocities closer together delta_v = 0.125_DP * (vmax_guess - vmin_guess)