diff --git a/src/collision/collision_module.f90 b/src/collision/collision_module.f90 index 4c1f78971..207de2772 100644 --- a/src/collision/collision_module.f90 +++ b/src/collision/collision_module.f90 @@ -151,6 +151,7 @@ module collision integer(I4B) :: status !! Status flag to pass to the collision list once the collision has been resolved integer(I4B) :: collision_id !! ID number of this collision event integer(I4B) :: maxid_collision = 0 !! The current maximum collision id number + real(DP) :: min_mfrag !! Minimum fragment mass ! Scale factors used to scale dimensioned quantities to a more "natural" system where important quantities (like kinetic energy, momentum) are of order ~1 real(DP) :: dscale = 1.0_DP !! Distance dimension scale factor diff --git a/src/collision/collision_util.f90 b/src/collision/collision_util.f90 index 772d6f988..33b4dddbf 100644 --- a/src/collision/collision_util.f90 +++ b/src/collision/collision_util.f90 @@ -921,6 +921,8 @@ module subroutine collision_util_set_natural_scale_factors(self) fragments%Gmass(:) = fragments%Gmass(:) / (collider%dscale**3/collider%tscale**2) fragments%radius(:) = fragments%radius(:) / collider%dscale impactors%Qloss = impactors%Qloss / collider%Escale + + collider%min_mfrag = collider%min_mfrag / collider%mscale end associate return @@ -986,6 +988,7 @@ module subroutine collision_util_set_original_scale_factors(self) collider%pe(:) = collider%pe(:) * collider%Escale collider%be(:) = collider%be(:) * collider%Escale collider%te(:) = collider%te(:) * collider%Escale + collider%min_mfrag = collider%min_mfrag * collider%mscale collider%mscale = 1.0_DP collider%dscale = 1.0_DP diff --git a/src/fraggle/fraggle_generate.f90 b/src/fraggle/fraggle_generate.f90 index daae97a48..a55ed2741 100644 --- a/src/fraggle/fraggle_generate.f90 +++ b/src/fraggle/fraggle_generate.f90 @@ -468,18 +468,19 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters logical, intent(out) :: lfailure !! Did the velocity computation fail? ! Internals - integer(I4B) :: i, j, loop, try, istart, nfrag, nlast, nsteps + 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 + 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 integer(I4B), dimension(collider%fragments%nbody) :: vsign - real(DP), dimension(collider%fragments%nbody) :: vscale, ke_rot_remove + 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.1_DP - real(DP) :: vmax_guess = 10.0_DP - real(DP) :: delta_v, volume - integer(I4B), parameter :: MAXLOOP = 200 - integer(I4B), parameter :: MAXTRY = 1000 + 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 + 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 character(len=STRMAX) :: message @@ -489,9 +490,14 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu lhitandrun = (impactors%regime == COLLRESOLVE_REGIME_HIT_AND_RUN) lsupercat = (impactors%regime == COLLRESOLVE_REGIME_SUPERCATASTROPHIC) + GC = impactors%Gmass(1) / impactors%mass(1) + allocate(collider_local, source=collider) associate(fragments => collider_local%fragments) + volume(:) = 4.0_DP / 3.0_DP * PI * (fragments%radius(:))**3 + fragments%density(:) = fragments%mass(:) / volume(:) + ! The fragments will be divided into two "clouds" based on identified origin body. ! These clouds will collectively travel like two impactors bouncing off of each other. where(fragments%origin_body(:) == 1) @@ -517,8 +523,9 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu E_residual_best = huge(1.0_DP) lfailure = .false. dE_metric = huge(1.0_DP) + dE_best = huge(1.0_DP) - outer: do try = 1, nfrag + outer: do try = 1, MAXTRY ! Scale the magnitude of the velocity by the distance from the impact point ! This will reduce the chances of fragments colliding with each other immediately, and is more physically correct do concurrent(i = 1:nfrag) @@ -582,7 +589,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 - if ((abs(E_residual) < abs(E_residual_best)) .or. ((dE < 0.0_DP) .and. (E_residual_best >= 0.0_DP))) then ! This is our best case so far. Save it for posterity + 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 @@ -620,24 +627,23 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu call fragments%set_coordinate_system() end do - if (dE_best < 0.0_DP) exit outer ! We didn't converge. Reset the fragment positions and velocities and try a new configuration with some slightly different parameters - if (fragments%nbody == 2) exit outer - ! Reduce the number of fragments by one - nlast = fragments%nbody - fragments%Ip(:,1) = fragments%mass(1) * fragments%Ip(:,1) + fragments%mass(nlast) * fragments%Ip(:,nlast) - fragments%mass(1) = fragments%mass(1) + fragments%mass(nlast) - fragments%Ip(:,1) = fragments%Ip(:,1) / fragments%mass(1) - fragments%Gmass(1) = fragments%Gmass(1) + fragments%mass(nlast) - volume = 4.0_DP / 3.0_DP * PI * ((fragments%radius(1))**3 + (fragments%radius(nlast))**3) - fragments%density(1) = fragments%mass(1) / volume - fragments%radius(1) = (3._DP * volume / (4._DP * PI))**(THIRD) - fragments%Ip(:,nlast) = 0.0_DP - fragments%mass(nlast) = 0.0_DP - fragments%Gmass(nlast) = 0.0_DP - fragments%radius(nlast) = 0.0_DP - fragments%status(nlast) = INACTIVE - fragments%nbody = nlast - 1 + ! Reduce the fragment masses and add it to the largest remenant and try again + if (any(fragments%mass(2:nfrag) > collider%min_mfrag)) then + do i = 2, nfrag + if (fragments%mass(i) > collider%min_mfrag) then + dM = min(mass_reduction_ratio * fragments%mass(i), fragments%mass(i) - collider%min_mfrag) + fragments%mass(i) = fragments%mass(i) - dM + fragments%mass(1) = fragments%mass(1) + dM + end if + end do + else + exit outer + end if + fragments%Gmass(:) = GC * fragments%mass(:) + + volume(:) = fragments%mass(:) / fragments%density(:) + fragments%radius(:) = (3._DP * volume(:) / (4._DP * PI))**(THIRD) call fragments%reset() call fraggle_generate_pos_vec(collider_local) diff --git a/src/fraggle/fraggle_util.f90 b/src/fraggle/fraggle_util.f90 index 58ff84d4d..18ad14759 100644 --- a/src/fraggle/fraggle_util.f90 +++ b/src/fraggle/fraggle_util.f90 @@ -25,7 +25,7 @@ module subroutine fraggle_util_set_mass_dist(self, param) integer(I4B) :: i, j, jproj, jtarg, nfrag, istart real(DP), dimension(2) :: volume real(DP), dimension(NDIM) :: Ip_avg - real(DP) :: mfrag, mremaining, min_mfrag, mtot, mcumul, G + real(DP) :: mfrag, mremaining, mtot, mcumul, G real(DP), dimension(:), allocatable :: mass real(DP), parameter :: BETA = 2.85_DP integer(I4B), parameter :: NFRAGMAX = 100 !! Maximum number of fragments that can be generated @@ -37,7 +37,7 @@ module subroutine fraggle_util_set_mass_dist(self, param) integer(I4B), dimension(:), allocatable :: ind logical :: flipper - associate(impactors => self%impactors) + associate(impactors => self%impactors, min_mfrag => self%min_mfrag) ! Get mass weighted mean of Ip and density volume(1:2) = 4._DP / 3._DP * PI * impactors%radius(1:2)**3 mtot = sum(impactors%mass(:)) diff --git a/src/swiftest/swiftest_util.f90 b/src/swiftest/swiftest_util.f90 index ecd3708ed..c1efe888f 100644 --- a/src/swiftest/swiftest_util.f90 +++ b/src/swiftest/swiftest_util.f90 @@ -1869,21 +1869,21 @@ module subroutine swiftest_util_rearray_pl(self, nbody_system, param) call plplenc_old%copy(nbody_system%plpl_encounter) end if - ! ! Re-build the encounter list - ! ! Be sure to get the level info if this is a SyMBA nbody_system - ! select type(nbody_system) - ! class is (symba_nbody_system) - ! select type(pl) - ! class is (symba_pl) - ! select type(tp) - ! class is (symba_tp) - ! lencounter = pl%encounter_check(param, nbody_system, param%dt, nbody_system%irec) - ! if (tp%nbody > 0) then - ! lencounter = tp%encounter_check(param, nbody_system, param%dt, nbody_system%irec) - ! end if - ! end select - ! end select - ! end select + ! Re-build the encounter list + ! Be sure to get the level info if this is a SyMBA nbody_system + select type(nbody_system) + class is (symba_nbody_system) + select type(pl) + class is (symba_pl) + select type(tp) + class is (symba_tp) + lencounter = pl%encounter_check(param, nbody_system, param%dt, nbody_system%irec) + if (tp%nbody > 0) then + lencounter = tp%encounter_check(param, nbody_system, param%dt, nbody_system%irec) + end if + end select + end select + end select ! Re-index the encounter list as the index values may have changed if (nenc_old > 0) then