diff --git a/examples/Fragmentation/Fragmentation_Movie.py b/examples/Fragmentation/Fragmentation_Movie.py index c77c69da6..89e4e9cfa 100755 --- a/examples/Fragmentation/Fragmentation_Movie.py +++ b/examples/Fragmentation/Fragmentation_Movie.py @@ -36,8 +36,8 @@ # ---------------------------------------------------------------------------------------------------------------------- # Define the names and initial conditions of the various fragmentation simulation types # ---------------------------------------------------------------------------------------------------------------------- -available_movie_styles = ["disruption_headon", "disruption_off_axis", "supercatastrophic_headon", "supercatastrophic_off_axis","hitandrun_disrupt", "hitandrun_pure", "merge"] -movie_title_list = ["Head-on Disruption", "Off-axis Disruption", "Head-on Supercatastrophic", "Off-axis Supercatastrophic", "Hit and Run w/ Runner Disruption", "Pure Hit and Run", "Merge"] +available_movie_styles = ["disruption_headon", "disruption_off_axis", "supercatastrophic_headon", "supercatastrophic_off_axis","hitandrun_disrupt", "hitandrun_pure", "merge", "merge_spinner"] +movie_title_list = ["Head-on Disruption", "Off-axis Disruption", "Head-on Supercatastrophic", "Off-axis Supercatastrophic", "Hit and Run w/ Runner Disruption", "Pure Hit and Run", "Merge", "Merge crossing the spin barrier"] movie_titles = dict(zip(available_movie_styles, movie_title_list)) num_movie_frames = 1000 @@ -56,6 +56,8 @@ "hitandrun_pure" : [np.array([1.0, -4.2e-05, 0.0]), np.array([1.0, 4.2e-05, 0.0])], "merge" : [np.array([1.0, -5.0e-05, 0.0]), + np.array([1.0, 5.0e-05 ,0.0])], + "merge_spinner" : [np.array([1.0, -5.0e-05, 0.0]), np.array([1.0, 5.0e-05 ,0.0])] } @@ -72,6 +74,8 @@ "hitandrun_pure" : [np.array([ 0.00, 6.28, 0.0]), np.array([-1.52, -6.28, 0.0])], "merge" : [np.array([ 0.04, 6.28, 0.0]), + np.array([ 0.05, 6.18, 0.0])], + "merge_spinner" : [np.array([ 0.04, 6.28, 0.0]), np.array([ 0.05, 6.18, 0.0])] } @@ -88,7 +92,9 @@ "hitandrun_pure" : [np.array([0.0, 0.0, 6.0e3]), np.array([0.0, 0.0, 1.0e4])], "merge" : [np.array([0.0, 0.0, 0.0]), - np.array([0.0, 0.0, 0.0])] + np.array([0.0, 0.0, 0.0])], + "merge_spinner" : [np.array([0.0, 0.0, -1.2e6]), + np.array([0.0, 0.0, 0.0])], } body_Gmass = {"disruption_headon" : [1e-7, 1e-9], @@ -97,7 +103,8 @@ "supercatastrophic_off_axis": [1e-7, 1e-8], "hitandrun_disrupt" : [1e-7, 7e-10], "hitandrun_pure" : [1e-7, 7e-10], - "merge" : [1e-7, 1e-8] + "merge" : [1e-7, 1e-8], + "merge_spinner" : [1e-7, 1e-8] } tstop = {"disruption_headon" : 2.0e-3, @@ -107,6 +114,7 @@ "hitandrun_disrupt" : 2.0e-4, "hitandrun_pure" : 2.0e-4, "merge" : 5.0e-3, + "merge_spinner" : 5.0e-3, } density = 3000 * swiftest.AU2M**3 / swiftest.MSun @@ -312,10 +320,11 @@ def vec_props(self, c): print("5. Hit and run with disruption of the runner") print("6. Pure hit and run") print("7. Merge") - print("8. All of the above") + print("8. Merge crossing the spin barrier") + print("9. All of the above") user_selection = int(input("? ")) - if user_selection > 0 and user_selection < 8: + if user_selection > 0 and user_selection < 9: movie_styles = [available_movie_styles[user_selection-1]] else: print("Generating all movie styles") diff --git a/src/collision/collision_module.f90 b/src/collision/collision_module.f90 index 207de2772..21b606202 100644 --- a/src/collision/collision_module.f90 +++ b/src/collision/collision_module.f90 @@ -33,6 +33,7 @@ module collision character(len=NAMELEN),parameter :: REGIME_NAME_GRAZE_AND_MERGE = "Graze and Merge" character(len=NAMELEN),parameter :: REGIME_NAME_HIT_AND_RUN = "Hit and Run" character(len=NAMELEN),dimension(5), parameter :: REGIME_NAMES = [REGIME_NAME_MERGE, REGIME_NAME_DISRUPTION, REGIME_NAME_SUPERCATASTROPHIC, REGIME_NAME_GRAZE_AND_MERGE, REGIME_NAME_HIT_AND_RUN] + real(DP), parameter :: MAX_ROT_SI = 7.108e-4 !! Spin limit in rad/s of cohesionless body from Holsapple (2007) !> Swiftest class for tracking pl-pl close encounters in a step when collisions are possible type, extends(encounter_list) :: collision_list_plpl @@ -88,7 +89,6 @@ module collision contains procedure :: consolidate => collision_resolve_consolidate_impactors !! Consolidates a multi-body collision into an equivalent 2-body collision - procedure :: get_regime => collision_regime_impactors !! Determine which fragmentation regime the set of impactors will be procedure :: dealloc => collision_util_dealloc_impactors !! Resets the collider object variables to 0 and deallocates the index and mass distributions procedure :: set_coordinate_system => collision_util_set_coordinate_impactors !! Sets the coordinate system of the impactors final :: collision_final_impactors !! Finalizer will deallocate all allocatables @@ -152,8 +152,9 @@ module collision 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 + real(DP) :: max_rot !! Maximum rotation rate (in system or natural units, depending on ) - ! Scale factors used to scale dimensioned quantities to a more "natural" system where important quantities (like kinetic energy, momentum) are of order ~1 + ! Scale factors used to scale dimensioned quantities to a more "natural" system where escape velocity is 1 and body masses are of order 1 real(DP) :: dscale = 1.0_DP !! Distance dimension scale factor real(DP) :: mscale = 1.0_DP !! Mass scale factor real(DP) :: tscale = 1.0_DP !! Time scale factor @@ -178,6 +179,7 @@ module collision procedure :: add_fragments => collision_util_add_fragments_to_collider !! Add fragments to nbody_system procedure :: get_energy_and_momentum => collision_util_get_energy_and_momentum !! Calculates total nbody_system energy in either the pre-collision outcome state (lbefore = .true.) or the post-collision outcome state (lbefore = .false.) procedure :: dealloc => collision_util_dealloc_basic !! Deallocates all allocatables + procedure :: get_regime => collision_regime_collider !! Determine which fragmentation regime the set of impactors will be procedure :: setup => collision_util_setup_collider !! Initializer for the encounter collision system and the before/after snapshots procedure :: setup_impactors => collision_util_setup_impactors_collider !! Initializer for the impactors for the encounter collision system. Deallocates old impactors before creating new ones procedure :: setup_fragments => collision_util_setup_fragments_collider !! Initializer for the fragments of the collision system. @@ -303,16 +305,16 @@ module subroutine collision_io_netcdf_write_frame_snapshot(self, history, param) class(base_parameters), intent(inout) :: param !! Current run configuration parameters end subroutine collision_io_netcdf_write_frame_snapshot - module subroutine collision_regime_impactors(self, nbody_system, param) + module subroutine collision_regime_collider(self, nbody_system, param) implicit none - class(collision_impactors), intent(inout) :: self !! Collision system impactors object - class(base_nbody_system), intent(in) :: nbody_system !! Swiftest nbody system object - class(base_parameters), intent(in) :: param !! Current Swiftest run configuration parameters - end subroutine collision_regime_impactors + class(collision_basic), intent(inout) :: self !! Collision system object + class(base_nbody_system), intent(in) :: nbody_system !! Swiftest nbody system object + class(base_parameters), intent(in) :: param !! Current Swiftest run configuration parameters + end subroutine collision_regime_collider module subroutine collision_check_plpl(self, nbody_system, param, t, dt, irec, lany_collision) implicit none - class(collision_list_plpl), intent(inout) :: self !! encounter list object + class(collision_list_plpl), intent(inout) :: self !! encounter list object class(base_nbody_system), intent(inout) :: nbody_system !! Swiftest nbody system object class(base_parameters), intent(inout) :: param !! Current run configuration parameters real(DP), intent(in) :: t !! current time diff --git a/src/collision/collision_regime.f90 b/src/collision/collision_regime.f90 index 966e5ae76..14d1286e1 100644 --- a/src/collision/collision_regime.f90 +++ b/src/collision/collision_regime.f90 @@ -13,20 +13,20 @@ contains - module subroutine collision_regime_impactors(self, nbody_system, param) + module subroutine collision_regime_collider(self, nbody_system, param) !! Author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton !! !! Determine which fragmentation regime the set of impactors will be. This subroutine is a wrapper for the non-polymorphic raggle_regime_collresolve subroutine. !! It converts to SI units prior to calling implicit none ! Arguments - class(collision_impactors), intent(inout) :: self !! Collision system impactors object - class(base_nbody_system), intent(in) :: nbody_system !! Swiftest nbody system object - class(base_parameters), intent(in) :: param !! Current Swiftest run configuration parameters + class(collision_basic), intent(inout) :: self !! Collision system impactors object + class(base_nbody_system), intent(in) :: nbody_system !! Swiftest nbody system object + class(base_parameters), intent(in) :: param !! Current Swiftest run configuration parameters ! Internals real (DP) :: mtot - associate(impactors => self) + associate(impactors => self%impactors) select type (nbody_system) class is (swiftest_nbody_system) select type(param) @@ -39,18 +39,18 @@ module subroutine collision_regime_impactors(self, nbody_system, param) allocate(impactors%mass_dist(1)) impactors%mass_dist(1) = mtot case default - call collision_regime_LS12(impactors, nbody_system, param) - call collision_io_log_regime(self) + call collision_regime_LS12(self, nbody_system, param) + call collision_io_log_regime(self%impactors) end select end select end select end associate return - end subroutine collision_regime_impactors + end subroutine collision_regime_collider - subroutine collision_regime_LS12(impactors, nbody_system, param) + subroutine collision_regime_LS12(collider, nbody_system, param) !! Author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton !! !! Determine the collisional regime of two colliding bodies based on the model by Leinhard and Stewart (2012) @@ -58,60 +58,93 @@ subroutine collision_regime_LS12(impactors, nbody_system, param) !! This is a wrapper subroutine that converts quantities to SI units and calls the main LS12 subroutine implicit none ! Arguments - class(collision_impactors), intent(inout) :: impactors !! The impactors to determine the regime for + class(collision_basic), intent(inout) :: collider !! The impactors to determine the regime for class(swiftest_nbody_system), intent(in) :: nbody_system !! Swiftest n-body system object class(swiftest_parameters), intent(in) :: param !! The current parameters ! Internals - integer(I4B) :: jtarg, jproj + integer(I4B) :: i,jtarg, jproj real(DP), dimension(2) :: radius_si, mass_si, density_si real(DP) :: min_mfrag_si, Mcb_si real(DP), dimension(NDIM) :: x1_si, v1_si, x2_si, v2_si - real(DP) :: mlr, mslr, mtot, dentot + real(DP) :: mlr, mslr, mslr_hr, mtot, dentot, Qloss, Qmerge integer(I4B), parameter :: NMASS_DIST = 3 !! Number of mass bins returned by the regime calculation (largest fragment, second largest, and remainder) + real(DP), dimension(NDIM) :: Ip, rot, L_spin + real(DP) :: radius, volume + + associate(impactors => collider%impactors) + + ! Convert all quantities to SI units and determine which of the pair is the projectile vs. target before sending them to the regime determination subroutine + if (impactors%mass(1) > impactors%mass(2)) then + jtarg = 1 + jproj = 2 + else + jtarg = 2 + jproj = 1 + end if + mass_si(:) = impactors%mass([jtarg, jproj]) * param%MU2KG !! The two-body equivalent masses of the collider system + radius_si(:) = impactors%radius([jtarg, jproj]) * param%DU2M !! The two-body equivalent radii of the collider system + density_si(:) = mass_si(:) / (4.0_DP / 3._DP * PI * radius_si(:)**3) !! The two-body equivalent density of the collider system + x1_si(:) = impactors%rb(:,jtarg) * param%DU2M !! The first body of the two-body equivalent position vector the collider system + v1_si(:) = impactors%vb(:,jtarg) * param%DU2M / param%TU2S !! The first body of the two-body equivalent velocity vector the collider system + x2_si(:) = impactors%rb(:,jproj) * param%DU2M !! The second body of the two-body equivalent position vector the collider system + v2_si(:) = impactors%vb(:,jproj) * param%DU2M / param%TU2S !! The second body of the two-body equivalent velocity vector the collider system + Mcb_si = nbody_system%cb%mass * param%MU2KG !! The central body mass of the system + min_mfrag_si = (param%min_GMfrag / param%GU) * param%MU2KG !! The minimum fragment mass to generate. Collider systems that would otherwise generate less massive fragments than this value will be forced to merge instead + + mtot = sum(mass_si(:)) + dentot = sum(mass_si(:) * density_si(:)) / mtot + + !! Use the positions and velocities of the parents from indside the step (at collision) to calculate the collisional regime + call collision_regime_LS12_SI(Mcb_si, mass_si(jtarg), mass_si(jproj), radius_si(jtarg), radius_si(jproj), & + x1_si(:), x2_si(:), v1_si(:), v2_si(:), density_si(jtarg), density_si(jproj), & + min_mfrag_si, impactors%regime, mlr, mslr, mslr_hr, Qloss, Qmerge) + + ! Convert back from SI to system units + mlr = mlr / param%MU2KG + mslr = mslr / param%MU2kg + mslr_hr = mslr_hr / param%MU2kg + Qloss = Qloss * (param%TU2S / param%DU2M)**2 / param%MU2KG + Qmerge = Qmerge * (param%TU2S / param%DU2M)**2 / param%MU2KG + mtot = mtot / param%MU2kg + + ! If this is came back as a merger, check to make sure that the rotation of the merged body doesn't exceed the spin barrier + if (impactors%regime == COLLRESOLVE_REGIME_MERGE) then + volume = 4._DP / 3._DP * PI * sum(impactors%radius(:)**3) + radius = (3._DP * volume / (4._DP * PI))**(THIRD) + do concurrent(i = 1:NDIM) + Ip(i) = sum(impactors%mass(:) * impactors%Ip(i,:)) + L_spin(i) = sum(impactors%L_orbit(i,:) + impactors%L_spin(i,:)) + end do + Ip(:) = Ip(:) / mtot + rot(:) = L_spin(:) / (Ip(3) * mtot * radius**2) + if (.mag.rot(:) > collider%max_rot) then ! The merged body would spin too fast, so reclasify this as a hit and run + mlr = impactors%mass(jtarg) + mslr = mslr_hr + impactors%regime = COLLRESOLVE_REGIME_HIT_AND_RUN + impactors%Qloss = Qloss + else + mlr = mtot + mslr = 0.0_DP + impactors%Qloss = Qmerge + end if + else + impactors%Qloss = Qloss + end if - ! Convert all quantities to SI units and determine which of the pair is the projectile vs. target before sending them to the regime determination subroutine - if (impactors%mass(1) > impactors%mass(2)) then - jtarg = 1 - jproj = 2 - else - jtarg = 2 - jproj = 1 - end if - mass_si(:) = impactors%mass([jtarg, jproj]) * param%MU2KG !! The two-body equivalent masses of the collider system - radius_si(:) = impactors%radius([jtarg, jproj]) * param%DU2M !! The two-body equivalent radii of the collider system - density_si(:) = mass_si(:) / (4.0_DP / 3._DP * PI * radius_si(:)**3) !! The two-body equivalent density of the collider system - x1_si(:) = impactors%rb(:,jtarg) * param%DU2M !! The first body of the two-body equivalent position vector the collider system - v1_si(:) = impactors%vb(:,jtarg) * param%DU2M / param%TU2S !! The first body of the two-body equivalent velocity vector the collider system - x2_si(:) = impactors%rb(:,jproj) * param%DU2M !! The second body of the two-body equivalent position vector the collider system - v2_si(:) = impactors%vb(:,jproj) * param%DU2M / param%TU2S !! The second body of the two-body equivalent velocity vector the collider system - Mcb_si = nbody_system%cb%mass * param%MU2KG !! The central body mass of the system - min_mfrag_si = (param%min_GMfrag / param%GU) * param%MU2KG !! The minimum fragment mass to generate. Collider systems that would otherwise generate less massive fragments than this value will be forced to merge instead - - mtot = sum(mass_si(:)) - dentot = sum(mass_si(:) * density_si(:)) / mtot - - !! Use the positions and velocities of the parents from indside the step (at collision) to calculate the collisional regime - call collision_regime_LS12_SI(Mcb_si, mass_si(jtarg), mass_si(jproj), radius_si(jtarg), radius_si(jproj), & - x1_si(:), x2_si(:), v1_si(:), v2_si(:), density_si(jtarg), density_si(jproj), & - min_mfrag_si, impactors%regime, mlr, mslr, impactors%Qloss) - - if (allocated(impactors%mass_dist)) deallocate(impactors%mass_dist) - allocate(impactors%mass_dist(NMASS_DIST)) - impactors%mass_dist(1) = min(max(mlr, 0.0_DP), mtot) - impactors%mass_dist(2) = min(max(mslr, 0.0_DP), mtot) - impactors%mass_dist(3) = min(max(mtot - mlr - mslr, 0.0_DP), mtot) - + if (allocated(impactors%mass_dist)) deallocate(impactors%mass_dist) + allocate(impactors%mass_dist(NMASS_DIST)) + impactors%mass_dist(1) = min(max(mlr, 0.0_DP), mtot) + impactors%mass_dist(2) = min(max(mslr, 0.0_DP), mtot) + impactors%mass_dist(3) = min(max(mtot - mlr - mslr, 0.0_DP), mtot) - ! Convert quantities back to the system units and save them into the fragment system - impactors%mass_dist(:) = (impactors%mass_dist(:) / param%MU2KG) - impactors%Qloss = impactors%Qloss * (param%TU2S / param%DU2M)**2 / param%MU2KG + end associate return end subroutine collision_regime_LS12 subroutine collision_regime_LS12_SI(Mcb, m1, m2, rad1, rad2, rh1, rh2, vb1, vb2, den1, den2, min_mfrag, & - regime, Mlr, Mslr, Qloss) + regime, Mlr, Mslr, Mslr_hitandrun, Qloss, Qmerge) !! Author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton !! !! Determine the collisional regime of two colliding bodies. @@ -133,8 +166,9 @@ subroutine collision_regime_LS12_SI(Mcb, m1, m2, rad1, rad2, rh1, rh2, vb1, vb2, real(DP), intent(in) :: Mcb, m1, m2, rad1, rad2, den1, den2, min_mfrag real(DP), dimension(:), intent(in) :: rh1, rh2, vb1, vb2 integer(I4B), intent(out) :: regime - real(DP), intent(out) :: Mlr, Mslr - real(DP), intent(out) :: Qloss !! The residual energy after the collision + real(DP), intent(out) :: Mlr, Mslr, Mslr_hitandrun ! Largest and second-largest remnant defined for various regimes + real(DP), intent(out) :: Qloss !! The energy lost in the collision if it was a fragmentation event + real(DP), intent(out) :: Qmerge !! The energy lost in the collision if it was a perfect merger ! Constants integer(I4B), parameter :: N1 = 1 !number of objects with mass equal to the largest remnant from LS12 integer(I4B), parameter :: N2 = 2 !number of objects with mass larger than second largest remnant from LS12 @@ -149,7 +183,7 @@ subroutine collision_regime_LS12_SI(Mcb, m1, m2, rad1, rad2, rh1, rh2, vb1, vb2, real(DP), parameter :: CRUFU = 2.0_DP - 3 * MU_BAR ! central potential variable from Rufu and Aharonson (2019) real(DP), parameter :: SUPERCAT_QRATIO = 1.8_DP ! See Section 4.1 of LS12 ! Internals - real(DP) :: a1, alpha, aint, b, bcrit, c_star, egy, zeta, l, lint, mu, phi, theta, ke, pe, Qmerge + real(DP) :: a1, alpha, aint, b, bcrit, c_star, egy, zeta, l, lint, mu, phi, theta, ke, pe real(DP) :: Qr, Qrd_pstar, Qr_erosion, Qr_supercat real(DP) :: Vhr, Verosion, Vescp, Vhill, Vimp, Vsupercat real(DP) :: Mint, Mtot, Mtmp @@ -201,8 +235,6 @@ subroutine collision_regime_LS12_SI(Mcb, m1, m2, rad1, rad2, rh1, rh2, vb1, vb2, Qr = mu*(Vimp**2) / Mtot / 2.0_DP !calculate mass largest remnant Mlr - Mlr = max((1.0_DP - Qr / Qrd_pstar / 2.0_DP) * Mtot, min_mfrag) ! [kg] # LS12 eq (5) - !calculate Vsupercat Qr_supercat = SUPERCAT_QRATIO * Qrd_pstar ! See LS12 Section 4.1 Vsupercat = sqrt(2 * Qr_supercat * Mtot / mu) @@ -212,64 +244,52 @@ subroutine collision_regime_LS12_SI(Mcb, m1, m2, rad1, rad2, rh1, rh2, vb1, vb2, theta = 1.0_DP - b Vhr = Vescp * (C1 * zeta**2 * theta**(2.5_DP) + C2 * zeta**2 + C3 * theta**(2.5_DP) + C4) ! Kokubo & Genda (2010) eq. (3) bcrit = rad1 / (rad1 + rad2) - Qloss = 0.0_DP ! Specific binding energy U_binding = (3 * GC * Mtot) / (5 * Rp) ! LS12 eq. 27 ke = 0.5_DP * Vimp**2 pe = - GC * m1 * m2 / (Mtot * norm2(rh2 - rh1)) - Qmerge = ke + pe + U_binding ! The specific energy lost if this were a perfect merger if ((m1 < min_mfrag).or.(m2 < min_mfrag)) then regime = COLLRESOLVE_REGIME_MERGE !perfect merging regime - Mlr = Mtot - Mslr = 0.0_DP - Qloss = Qmerge call swiftest_io_log_one_message(COLLISION_LOG_OUT, & "Fragments would have mass below the minimum. Converting this collision into a merger.") else if( Vimp < Vescp) then regime = COLLRESOLVE_REGIME_MERGE !perfect merging regime - Mlr = Mtot - Mslr = 0.0_DP - Qloss = Qmerge else if (Vimp < Verosion) then if (b < bcrit) then regime = COLLRESOLVE_REGIME_MERGE !partial accretion regime" - Mlr = Mtot - Mslr = 0.0_DP - Qloss = Qmerge else if ((b > bcrit) .and. (Vimp < Vhr)) then regime = COLLRESOLVE_REGIME_MERGE ! graze and merge - Mlr = Mtot - Mslr = 0.0_DP - Qloss = Qmerge else - Mlr = m1 - Mslr = max(calc_Qrd_rev(m2, m1, Mint, den1, den2, Vimp, c_star), min_mfrag) regime = COLLRESOLVE_REGIME_HIT_AND_RUN !hit and run - Qloss = (c_star - 1.0_DP) * U_binding end if else if (Vimp > Verosion .and. Vimp < Vsupercat) then if (m2 < 0.001_DP * m1) then regime = COLLRESOLVE_REGIME_MERGE !cratering regime" - Mlr = Mtot - Mslr = 0.0_DP - Qloss = Qmerge else - Mslr = max(Mtot * (3.0_DP - BETA) * (1.0_DP - N1 * Mlr / Mtot) / (N2 * BETA), min_mfrag) ! LS12 eq (37) regime = COLLRESOLVE_REGIME_DISRUPTION !disruption - Qloss = (c_star - 1.0_DP) * U_binding end if else if (Vimp > Vsupercat) then - Mlr = max(Mtot * 0.1_DP * (Qr / (Qrd_pstar * SUPERCAT_QRATIO))**(ETA), min_mfrag) !LS12 eq (44) - Mslr = max(Mtot * (3.0_DP - BETA) * (1.0_DP - N1 * Mlr / Mtot) / (N2 * BETA), min_mfrag) !LS12 eq (37) regime = COLLRESOLVE_REGIME_SUPERCATASTROPHIC ! supercatastrophic - Qloss = (c_star - 1.0_DP) * U_binding else call swiftest_io_log_one_message(COLLISION_LOG_OUT,"Error no regime found in symba_regime") end if end if + if (regime == COLLRESOLVE_REGIME_SUPERCATASTROPHIC) then + Mlr = max(Mtot * 0.1_DP * (Qr / (Qrd_pstar * SUPERCAT_QRATIO))**(ETA), min_mfrag) !LS12 eq (44) + else + Mlr = max((1.0_DP - Qr / Qrd_pstar / 2.0_DP) * Mtot, min_mfrag) ! [kg] # LS12 eq (5) + end if + + Mslr_hitandrun = max(calc_Qrd_rev(m2, m1, Mint, den1, den2, Vimp, c_star), min_mfrag) + if (regime == COLLRESOLVE_REGIME_HIT_AND_RUN ) then + Mslr = Mslr_hitandrun + else + Mslr = max(Mtot * (3.0_DP - BETA) * (1.0_DP - N1 * Mlr / Mtot) / (N2 * BETA), min_mfrag) !LS12 eq (37) + end if + Mresidual = Mtot - Mlr - Mslr if (Mresidual < 0.0_DP) then ! prevents final masses from going negative Mlr = Mlr + Mresidual @@ -281,8 +301,8 @@ subroutine collision_regime_LS12_SI(Mcb, m1, m2, rad1, rad2, rh1, rh2, vb1, vb2, Mslr = Mtmp end if - Qloss = Qloss * Mtot ! Convert specific energy loss to total energy loss in the system - + Qloss = (c_star - 1.0_DP) * U_binding * Mtot ! Convert specific energy loss to total energy loss in the system + Qmerge = (ke + pe + U_binding) * Mtot ! The energy lost if this were a perfect merger return diff --git a/src/collision/collision_resolve.f90 b/src/collision/collision_resolve.f90 index 525b53969..5e0f2296a 100644 --- a/src/collision/collision_resolve.f90 +++ b/src/collision/collision_resolve.f90 @@ -566,7 +566,7 @@ module subroutine collision_resolve_plpl(self, nbody_system, param, t, dt, irec) call swiftest_io_log_one_message(COLLISION_LOG_OUT, "collision_id " // trim(adjustl(idstr))) ! Get the collision regime - call impactors%get_regime(nbody_system, param) + call collider%get_regime(nbody_system, param) call collision_history%take_snapshot(param,nbody_system, t, "before") diff --git a/src/collision/collision_util.f90 b/src/collision/collision_util.f90 index 33b4dddbf..a73238ef7 100644 --- a/src/collision/collision_util.f90 +++ b/src/collision/collision_util.f90 @@ -923,6 +923,7 @@ module subroutine collision_util_set_natural_scale_factors(self) impactors%Qloss = impactors%Qloss / collider%Escale collider%min_mfrag = collider%min_mfrag / collider%mscale + collider%max_rot = collider%max_rot * collider%tscale end associate return @@ -989,6 +990,7 @@ module subroutine collision_util_set_original_scale_factors(self) collider%be(:) = collider%be(:) * collider%Escale collider%te(:) = collider%te(:) * collider%Escale collider%min_mfrag = collider%min_mfrag * collider%mscale + collider%max_rot = collider%max_rot / collider%tscale 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 d0aa3bcb3..49f104503 100644 --- a/src/fraggle/fraggle_generate.f90 +++ b/src/fraggle/fraggle_generate.f90 @@ -191,6 +191,7 @@ module subroutine fraggle_generate_disrupt(self, nbody_system, param, t, lfailur call swiftest_io_log_one_message(COLLISION_LOG_OUT, message) call self%set_original_scale() + self%max_rot = MAX_ROT_SI * param%TU2S ! Re-compute the spin limit from scratch so it doesn't drift due to floating point errors every time we convert ! Restore the big array if (lk_plpl) call pl%flatten(param) @@ -449,6 +450,10 @@ module subroutine fraggle_generate_rot_vec(collider, nbody_system, param) rotmag = random_scale_factor * .mag. fragments%rot(:,i) fragments%rot(:,i) = fragments%rot(:,i) + rotmag * rotdir fragments%rotmag(i) = .mag.fragments%rot(:,i) + if (fragments%rotmag(i) > collider%max_rot) then + fragments%rotmag(i) = collider%max_rot + fragments%rot(:,i) = fragments%rotmag(i) * .unit.fragments%rot(:,i) + end if end do end associate @@ -471,17 +476,18 @@ 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, vdir, L_residual_unit, Li, Lrat, r_lever - real(DP) :: rmag, vimp, 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 + real(DP), dimension(NDIM) :: vimp_unit, rimp, vrot, L_residual, L_residual_unit, dL, drot, rot_new + real(DP) :: vimp, vmag, vesc, dE, E_residual, ke_min, ke_avail, ke_remove, dE_best, E_residual_best, fscale, dE_metric, dM, mfrag, drotmag integer(I4B), dimension(collider%fragments%nbody) :: vsign - real(DP), dimension(collider%fragments%nbody) :: vscale, ke_rot_remove, volume + real(DP), dimension(collider%fragments%nbody) :: vscale, 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) :: vmax_guess real(DP) :: delta_v, GC 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 :: MAX_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 :: ROT_MAX_FRAC = 0.05_DP !! Fraction of difference between current rotation and maximum to add when angular momentum budget gets too high real(DP), parameter :: SUCCESS_METRIC = 1.0e-2_DP class(collision_fraggle), allocatable :: collider_local character(len=STRMAX) :: message @@ -522,7 +528,7 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu end if vimp = .mag. (impactors%vc(:,2) - impactors%vc(:,1)) - vmax_guess = 2 * vimp + vmax_guess = 1.1_DP * vimp E_residual_best = huge(1.0_DP) lfailure = .false. @@ -533,7 +539,7 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu ! 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) - rimp(:) = fragments%rc(:,i) !- impactors%rbimp(:) + rimp(:) = fragments%rc(:,i) - impactors%rbimp(:) vscale(i) = .mag. rimp(:) / (.mag. (impactors%rb(:,2) - impactors%rb(:,1))) end do @@ -550,10 +556,10 @@ 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) + fragments%vc(:,i) = vsign(i) * impactors%bounce_unit(:) * vscale(i) end if else - vmag = vesc * vscale(i) + vmag = vscale(i) rimp(:) = fragments%rc(:,i) - impactors%rbimp(:) vimp_unit(:) = .unit. (rimp(:) + vsign(i) * impactors%bounce_unit(:)) fragments%vc(:,i) = vmag * vimp_unit(:) + vrot(:) @@ -569,13 +575,33 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu do loop = 1, MAXLOOP nsteps = loop * try - ! 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 + ! 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 + angmtm: do j = 1, MAXTRY + do i = istart, fragments%nbody + call collider_local%get_energy_and_momentum(nbody_system, param, phase="after") + L_residual(:) = (collider_local%L_total(:,2) - collider_local%L_total(:,1)) + if (.mag.L_residual(:) / .mag.collider_local%L_total(:,1) <= epsilon(1.0_DP)) exit angmtm + L_residual_unit(:) = .unit. L_residual(:) + mfrag = sum(fragments%mass(i:fragments%nbody)) + drot(:) = -L_residual(:) / (mfrag * 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. + drotmag = collider_local%max_rot - fragments%rotmag(i) + drot(:) = ROT_MAX_FRAC * drotmag * L_residual_unit(:) ! Put a fraction of the difference between the spin barrier and the current spin into the new rotation + fragments%rot(:,i) = fragments%rot(:,i) + drot(:) + dL(:) = -L_residual(:) * fragments%mass(i) / mfrag + drot(:) * fragments%Ip(3,i) * fragments%mass(i) * fragments%radius(i)**2 + call fraggle_generate_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) + fragments%rotmag(i) = .mag.fragments%rot(:,i) + 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 + end do + end do angmtm call collider_local%get_energy_and_momentum(nbody_system, param, phase="after") ke_avail = 0.0_DP @@ -583,40 +609,6 @@ 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 put it into spin and shear - L_residual(:) = (collider_local%L_total(:,2) - collider_local%L_total(:,1)) - 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() - call collider_local%get_energy_and_momentum(nbody_system, param, phase="after") - L_residual(:) = (collider_local%L_total(:,2) - collider_local%L_total(:,1)) - - ! 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") - L_residual(:) = (collider_local%L_total(:,2) - collider_local%L_total(:,1)) - end if - dE = collider_local%te(2) - collider_local%te(1) E_residual = dE + impactors%Qloss @@ -637,20 +629,9 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu ! Remove a constant amount of velocity from the bodies so we don't shift the center of mass and screw up the momentum ke_remove = min(E_residual, ke_avail) - f_orbit = ke_remove / E_residual fscale = sqrt((max(fragments%ke_orbit_tot - ke_remove, 0.0_DP))/fragments%ke_orbit_tot) fragments%vc(:,:) = fscale * fragments%vc(:,:) - ! f_spin = 1.0_DP - f_orbit - ! 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 = 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) - ! end do - ! Update the unit vectors and magnitudes for the fragments based on their new orbits and rotations call collision_util_shift_vector_to_origin(fragments%mass, fragments%vc) call fragments%set_coordinate_system() @@ -661,7 +642,7 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu 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) + dM = min(MAX_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 @@ -702,4 +683,33 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu end subroutine fraggle_generate_vel_vec + subroutine fraggle_generate_velocity_torque(dL, mass, r, v) + !! author: David A. Minton + !! + !! Applies a torque to a body's center of mass velocity given a change in angular momentum + implicit none + ! Arguments + real(DP), dimension(:), intent(in) :: dL !! Change in angular momentum to apply + real(DP), intent(in) :: mass !! Mass of body + real(DP), dimension(:), intent(in) :: r !! Position of body wrt system center of mass + real(DP), dimension(:), intent(inout) :: v !! Velocity of body wrt system center of mass + ! Internals + real(DP), dimension(NDIM) :: dL_unit, r_unit, r_lever, vapply + real(DP) :: rmag, r_lever_mag + + dL_unit(:) = .unit. dL + r_unit(:) = .unit.r(:) + rmag = .mag.r(:) + ! Project the position vector onto the plane defined by the angular momentum vector and the origin to get the "lever arm" distance + r_lever(:) = dL_unit(:) .cross. (r(:) .cross. dL_unit(:)) + r_lever_mag = .mag.r_lever(:) + if (r_lever_mag > epsilon(1.0_DP)) then + vapply(:) = (dL(:) .cross. r(:)) / (mass * rmag**2) + v(:) = v(:) + vapply(:) + end if + + return + end subroutine fraggle_generate_velocity_torque + + end submodule s_fraggle_generate diff --git a/src/symba/symba_util.f90 b/src/symba/symba_util.f90 index a8e28c477..a30c08bac 100644 --- a/src/symba/symba_util.f90 +++ b/src/symba/symba_util.f90 @@ -331,6 +331,8 @@ module subroutine symba_util_setup_initialize_system(self, param) allocate(collision_fraggle :: nbody_system%collider) end select call nbody_system%collider%setup(nbody_system) + + nbody_system%collider%max_rot = MAX_ROT_SI * param%TU2S select type(nc => collision_history%nc) class is (collision_netcdf_parameters) nbody_system%collider%maxid_collision = nc%max_idslot