From 3bb3b044b6683dc852a3e6253e50da80b4ff3101 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Mon, 26 Dec 2022 08:25:24 -0500 Subject: [PATCH] Lots of improvements to the simple collision model, which should improve Fraggle's ability to successfully find solutions --- examples/Fragmentation/Fragmentation_Movie.py | 40 +++-- src/collision/collision_generate.f90 | 154 ++++++++++-------- src/collision/collision_module.f90 | 52 +++--- src/collision/collision_resolve.f90 | 18 +- src/collision/collision_util.f90 | 46 +++++- src/fraggle/fraggle_generate.f90 | 2 +- src/fraggle/fraggle_util.f90 | 4 +- src/swiftest/swiftest_util.f90 | 22 +-- 8 files changed, 201 insertions(+), 137 deletions(-) diff --git a/examples/Fragmentation/Fragmentation_Movie.py b/examples/Fragmentation/Fragmentation_Movie.py index 381bca121..317ff7de7 100644 --- a/examples/Fragmentation/Fragmentation_Movie.py +++ b/examples/Fragmentation/Fragmentation_Movie.py @@ -36,30 +36,42 @@ # ---------------------------------------------------------------------------------------------------------------------- # Define the names and initial conditions of the various fragmentation simulation types # ---------------------------------------------------------------------------------------------------------------------- -available_movie_styles = ["disruption_headon", "supercatastrophic_off_axis", "hitandrun"] -movie_title_list = ["Head-on Disruption", "Off-axis Supercatastrophic", "Hit and Run"] +available_movie_styles = ["disruption_headon", "disruption_off_axis", "supercatastrophic_headon", "supercatastrophic_off_axis","hitandrun"] +movie_title_list = ["Head-on Disruption", "Off-axis Disruption", "Head-on Supercatastrophic", "Off-axis Supercatastrophic", "Hit and Run"] movie_titles = dict(zip(available_movie_styles, movie_title_list)) # These initial conditions were generated by trial and error names = ["Target","Projectile"] pos_vectors = {"disruption_headon" : [np.array([1.0, -5.0e-05, 0.0]), np.array([1.0, 5.0e-05 ,0.0])], + "disruption_off_axis" : [np.array([1.0, -5.0e-05, 0.0]), + np.array([1.0, 5.0e-05 ,0.0])], + "supercatastrophic_headon": [np.array([1.0, -5.0e-05, 0.0]), + np.array([1.0, 5.0e-05, 0.0])], "supercatastrophic_off_axis": [np.array([1.0, -5.0e-05, 0.0]), np.array([1.0, 5.0e-05, 0.0])], "hitandrun" : [np.array([1.0, -4.2e-05, 0.0]), np.array([1.0, 4.2e-05, 0.0])] } -vel_vectors = {"disruption_headon" : [np.array([-2.562596e-04, 6.280005, 0.0]), - np.array([-2.562596e-04, -6.280005, 0.0])], - "supercatastrophic_off_axis": [np.array([0.0, 6.28, 0.0]), - np.array([0.5, -6.28, 0.0])], - "hitandrun" : [np.array([0.0, 6.28, 0.0]), - np.array([-1.45, -6.28, 0.00])] +vel_vectors = {"disruption_headon" : [np.array([ 0.00, 6.280005, 0.0]), + np.array([ 0.00, -6.280005, 0.0])], + "disruption_off_axis" : [np.array([ 0.00, 6.280005, 0.0]), + np.array([ 0.50, -6.280005, 0.0])], + "supercatastrophic_headon": [np.array([ 0.00, 6.28, 0.0]), + np.array([ 0.00, -6.28, 0.0])], + "supercatastrophic_off_axis": [np.array([ 0.00, 6.28, 0.0]), + np.array([ 0.50, -6.28, 0.0])], + "hitandrun" : [np.array([ 0.00, 6.28, 0.0]), + np.array([-1.45, -6.28, 0.0])] } rot_vectors = {"disruption_headon" : [np.array([0.0, 0.0, 0.0]), np.array([0.0, 0.0, 0.0])], + "disruption_off_axis": [np.array([0.0, 0.0, 6.0e4]), + np.array([0.0, 0.0, 1.0e5])], + "supercatastrophic_headon": [np.array([0.0, 0.0, 0.0]), + np.array([0.0, 0.0, 0.0])], "supercatastrophic_off_axis": [np.array([0.0, 0.0, -6.0e4]), np.array([0.0, 0.0, 1.0e5])], "hitandrun" : [np.array([0.0, 0.0, 6.0e4]), @@ -67,6 +79,8 @@ } body_Gmass = {"disruption_headon" : [1e-7, 1e-10], + "disruption_off_axis" : [1e-7, 1e-10], + "supercatastrophic_headon": [1e-7, 1e-8], "supercatastrophic_off_axis": [1e-7, 1e-8], "hitandrun" : [1e-7, 7e-10] } @@ -185,12 +199,14 @@ def data_stream(self, frame=0): print("Select a fragmentation movie to generate.") print("1. Head-on disruption") - print("2. Off-axis supercatastrophic") - print("3. Hit and run") - print("4. All of the above") + print("2. Off-axis disruption") + print("3. Head-on supercatastrophic") + print("4. Off-axis supercatastrophic") + print("5. Hit and run") + print("6. All of the above") user_selection = int(input("? ")) - if user_selection > 0 and user_selection < 4: + if user_selection > 0 and user_selection < 6: movie_styles = [available_movie_styles[user_selection-1]] else: print("Generating all movie styles") diff --git a/src/collision/collision_generate.f90 b/src/collision/collision_generate.f90 index 0f6da8d79..16e0c39b7 100644 --- a/src/collision/collision_generate.f90 +++ b/src/collision/collision_generate.f90 @@ -32,6 +32,7 @@ module subroutine collision_generate_basic(self, nbody_system, param, t) return end subroutine collision_generate_basic + module subroutine collision_generate_bounce(self, nbody_system, param, t) !! author: David A. Minton !! @@ -200,8 +201,7 @@ module subroutine collision_generate_simple(self, nbody_system, param, t) class(base_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA additions real(DP), intent(in) :: t !! Time of collision ! Internals - integer(I4B) :: i, ibiggest, nfrag, status - logical :: lfailure + integer(I4B) :: i, ibiggest, nfrag character(len=STRMAX) :: message real(DP) :: dpe @@ -370,16 +370,14 @@ module subroutine collision_generate_disrupt(self, nbody_system, param, t, lfail implicit none ! Arguments class(collision_simple_disruption), intent(inout) :: self !! Simple fragment system 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 !! The time of the collision - logical, optional, intent(out) :: lfailure + 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 !! The time of the collision + logical, optional, intent(out) :: lfailure ! Internals - real(DP) :: r_max_start call self%get_energy_and_momentum(nbody_system, param, lbefore=.true.) - r_max_start = 1.5_DP * .mag.(self%impactors%rb(:,2) - self%impactors%rb(:,1)) - call collision_generate_simple_pos_vec(self, r_max_start) + call collision_generate_simple_pos_vec(self) call self%set_coordinate_system() call collision_generate_simple_vel_vec(self) call collision_generate_simple_rot_vec(self) @@ -388,7 +386,7 @@ module subroutine collision_generate_disrupt(self, nbody_system, param, t, lfail end subroutine collision_generate_disrupt - module subroutine collision_generate_simple_pos_vec(collider, r_max_start) + module subroutine collision_generate_simple_pos_vec(collider) !! Author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton !! !! Initializes the position vectors of the fragments around the center of mass based on the collision style. @@ -398,46 +396,51 @@ module subroutine collision_generate_simple_pos_vec(collider, r_max_start) implicit none ! Arguments class(collision_simple_disruption), intent(inout) :: collider !! Fraggle collision system object - real(DP), intent(in) :: r_max_start !! The maximum radial distance of fragments for disruptive collisions ! Internals - real(DP) :: dis, rad, r_max, fdistort - logical, dimension(:), allocatable :: loverlap - integer(I4B) :: i, j - logical :: lnoncat, lhitandrun + real(DP) :: dis + real(DP), dimension(NDIM,2) :: fragment_cloud_center + real(DP), dimension(2) :: fragment_cloud_radius + logical, dimension(collider%fragments%nbody) :: loverlap + integer(I4B) :: i, j, loop + logical :: lcat, lhitandrun + integer(I4B), parameter :: MAXLOOP = 10000 + real(DP) :: rdistance + real(DP), parameter :: fail_scale = 1.1_DP ! Scale factor to apply to cloud radius and distance if cloud generation fails - associate(fragments => collider%fragments, impactors => collider%impactors, nfrag => collider%fragments%nbody) - allocate(loverlap(nfrag)) - lnoncat = (impactors%regime /= COLLRESOLVE_REGIME_SUPERCATASTROPHIC) ! For non-catastrophic impacts, make the fragments act like ejecta and point away from the impact point - lhitandrun = (impactors%regime == COLLRESOLVE_REGIME_HIT_AND_RUN) ! Disruptive hit and runs have their own fragment distribution + associate(fragments => collider%fragments, impactors => collider%impactors, nfrag => collider%fragments%nbody) + lcat = (impactors%regime == COLLRESOLVE_REGIME_SUPERCATASTROPHIC) + lhitandrun = (impactors%regime == COLLRESOLVE_REGIME_HIT_AND_RUN) - ! Place the fragments into a region that is big enough that we should usually not have overlapping bodies - ! An overlapping bodies will collide in the next time step, so it's not a major problem if they do (it just slows the run down) - r_max = r_max_start - rad = sum(impactors%radius(:)) + ! We will treat the first two fragments of the list as special cases. + ! Place the first two bodies at the centers of the two fragment clouds, but be sure they are sufficiently far apart to avoid overlap + rdistance = .mag. (impactors%rc(:,2) - impactors%rc(:,1)) - sum(fragments%radius(1:2)) + rdistance = min(0.5_DP*rdistance, 1e-6_DP*impactors%radius(2)) - ! This is a factor that will "distort" the shape of the fragment cloud in the direction of the impact velocity - fdistort = .mag. (impactors%y_unit(:) .cross. impactors%v_unit(:)) + fragment_cloud_radius(:) = impactors%radius(:) - ! We will treat the first two fragments of the list as special cases. They get initialized at the original positions of the impactor bodies - fragments%rc(:, 1) = impactors%rb(:, 1) - impactors%rbcom(:) - fragments%rc(:, 2) = impactors%rb(:, 2) - impactors%rbcom(:) - call random_number(fragments%rc(:,3:nfrag)) loverlap(:) = .true. - do while (any(loverlap(3:nfrag))) - if (lhitandrun) then ! For a hit-and-run with disruption, the fragment cloud size is based on the radius of the disrupted body - r_max = 2 * impactors%radius(2) - else ! For disruptions, the the fragment cloud size is based on the mutual collision system - r_max = r_max + 0.1_DP * rad - end if - do i = 3, nfrag - if (loverlap(i)) then + do loop = 1, MAXLOOP + if (.not.any(loverlap(:))) exit + fragment_cloud_center(:,1) = impactors%rc(:,1) + rdistance * impactors%bounce_unit(:) + fragment_cloud_center(:,2) = impactors%rc(:,2) - rdistance * impactors%bounce_unit(:) + do concurrent(i = 1:nfrag, loverlap(i)) + if (i < 3) then + fragments%rc(:,i) = fragment_cloud_center(:,i) + else + ! Make a random cloud call random_number(fragments%rc(:,i)) - fragments%rc(:,i) = 2 * (fragments%rc(:, i) - 0.5_DP) - fragments%rc(:, i) = fragments%rc(:,i) + fdistort * impactors%v_unit(:) - fragments%rc(:, i) = r_max * fragments%rc(:, i) - fragments%rc(:, i) = fragments%rc(:, i) + (impactors%rbimp(:) - impactors%rbcom(:)) ! Shift the center of the fragment cloud to the impact point rather than the CoM - if (lnoncat .and. dot_product(fragments%rc(:,i), impactors%y_unit(:)) < 0.0_DP) fragments%rc(:, i) = -fragments%rc(:, i) ! Make sure the fragment cloud points away from the impact point + + ! Make the fragment cloud symmertic about 0 + fragments%rc(:,i) = 2 *(fragments%rc(:,i) - 0.5_DP) + + j = fragments%origin_body(i) + + ! Scale the cloud size + fragments%rc(:,i) = fragment_cloud_radius(j) * fragments%rc(:,i) + + ! Shift to the cloud center coordinates + fragments%rc(:,i) = fragments%rc(:,i) + fragment_cloud_center(:,j) end if end do @@ -449,6 +452,8 @@ module subroutine collision_generate_simple_pos_vec(collider, r_max_start) loverlap(i) = loverlap(i) .or. (dis <= (fragments%radius(i) + fragments%radius(j))) end do end do + rdistance = rdistance * fail_scale + fragment_cloud_radius(:) = fragment_cloud_radius(:) * fail_scale end do call collision_util_shift_vector_to_origin(fragments%mass, fragments%rc) call collider%set_coordinate_system() @@ -476,7 +481,7 @@ module subroutine collision_generate_simple_rot_vec(collider) ! Arguments class(collision_basic), intent(inout) :: collider !! Fraggle collision system object ! Internals - real(DP), dimension(NDIM) :: Ltotal, Lresidual + real(DP), dimension(NDIM) :: Lresidual integer(I4B) :: i associate(fragments => collider%fragments, impactors => collider%impactors, nfrag => collider%fragments%nbody) @@ -517,35 +522,52 @@ module subroutine collision_generate_simple_vel_vec(collider) ! Arguments class(collision_simple_disruption), intent(inout) :: collider !! Fraggle collision system object ! Internals - integer(I4B) :: i - logical :: lhr, lnoncat - real(DP), dimension(NDIM) :: vimp_unit, vcom, rnorm - real(DP), dimension(NDIM,collider%fragments%nbody) :: vnoise - real(DP), parameter :: VNOISE_MAG = 0.20_DP - real(DP) :: vmag + integer(I4B) :: i,j + logical :: lhitandrun, lnoncat + real(DP), dimension(NDIM) :: vimp_unit, rimp, vrot + real(DP), dimension(2) :: vimp + real(DP) :: vmag, vdisp + integer(I4B), dimension(collider%fragments%nbody) :: vsign + real(DP), dimension(collider%fragments%nbody) :: vscale associate(fragments => collider%fragments, impactors => collider%impactors, nfrag => collider%fragments%nbody) - lhr = (impactors%regime == COLLRESOLVE_REGIME_HIT_AND_RUN) + lhitandrun = (impactors%regime == COLLRESOLVE_REGIME_HIT_AND_RUN) lnoncat = (impactors%regime /= COLLRESOLVE_REGIME_SUPERCATASTROPHIC) ! For non-catastrophic impacts, make the fragments act like ejecta and point away from the impact point ! "Bounce" the first two bodies - do i = 1,2 - vcom(:) = impactors%vb(:,i) - impactors%vbcom(:) - rnorm(:) = impactors%y_unit(:) - ! Do the reflection - if (.not. lhr) vcom(:) = vcom(:) - 2 * dot_product(vcom(:),rnorm(:)) * rnorm(:) - fragments%vc(:,i) = vcom(:) - end do - - ! Compute the escape velocity and velocity dispursion - call random_number(vnoise) - if (lhr) then - vmag = 2 * sqrt(2 * impactors%Gmass(2) / impactors%radius(2)) + where(fragments%origin_body(:) == 1) + vsign(:) = -1 + elsewhere + vsign(:) = 1 + end where + + ! Compute the velocity dispersion based on the escape velocity + if (lhitandrun) then + vdisp = 2 * sqrt(2 * impactors%Gmass(2) / impactors%radius(2)) else - vmag = .mag. (impactors%vb(:,2) - impactors%vb(:,1)) + vdisp = 2 * sqrt(2 * sum(impactors%Gmass(:)) / sum(impactors%radius(:))) + !vmag = abs(dot_product(impactors%vb(:,2) - impactors%vb(:,1), impactors%y_unit(:))) end if - do i = 3, nfrag - vimp_unit(:) = .unit. (fragments%rc(:,i) - (impactors%rbimp(:) - impactors%rbcom(:))) - fragments%vc(:,i) = vmag * (vimp_unit(:) + vnoise(:,i)) + fragments%vc(:,2) + + vimp(:) = .mag.impactors%vc(:,:) + + ! Scale the magnitude of the velocity by the distance from the impact point and add a bit of shear + do concurrent(i = 1:nfrag) + rimp(:) = fragments%rc(:,i) - impactors%rbimp(:) + vscale(i) = .mag. rimp(:) / (.mag. (impactors%rb(:,2) - impactors%rb(:,1))) + end do + vscale(:) = vscale(:)/maxval(vscale(:)) + + fragments%vc(:,1) = .mag.impactors%vc(:,1) * impactors%bounce_unit(:) + do concurrent(i = 2:nfrag) + rimp(:) = fragments%rc(:,i) - impactors%rbimp(:) + vimp_unit(:) = .unit. rimp(:) + + j = fragments%origin_body(i) + vrot(:) = impactors%rot(:,j) .cross. (fragments%rc(:,i) - impactors%rb(:,j) + impactors%rbcom(:)) + + vmag = .mag.impactors%vc(:,j) * vscale(i) + + fragments%vc(:,i) = vmag * 0.5_DP * (impactors%bounce_unit(:) + vimp_unit(:)) * vsign(i) + vrot(:) end do do concurrent(i = 1:nfrag) fragments%vb(:,i) = fragments%vc(:,i) + impactors%vbcom(:) diff --git a/src/collision/collision_module.f90 b/src/collision/collision_module.f90 index ecc7aa3d5..c0df619c5 100644 --- a/src/collision/collision_module.f90 +++ b/src/collision/collision_module.f90 @@ -52,8 +52,10 @@ module collision type, extends(base_object) :: collision_impactors integer(I4B) :: ncoll !! Number of bodies involved in the collision integer(I4B), dimension(:), allocatable :: id !! Index of bodies involved in the collision - real(DP), dimension(NDIM,2) :: rb !! Two-body equivalent position vectors of the collider bodies prior to collision - real(DP), dimension(NDIM,2) :: vb !! Two-body equivalent velocity vectors of the collider bodies prior to collision + real(DP), dimension(NDIM,2) :: rb !! Two-body equivalent position vectors of the collider bodies prior to collision in system barycentric coordinates + real(DP), dimension(NDIM,2) :: vb !! Two-body equivalent velocity vectors of the collider bodies prior to collision in system barycentric coordinate + real(DP), dimension(NDIM,2) :: rc !! Two-body equivalent position vectors of the collider bodies prior to collision in collision center of mass coordinates + real(DP), dimension(NDIM,2) :: vc !! Two-body equivalent velocity vectors of the collider bodies prior to collision in collision center of mass coordinates real(DP), dimension(NDIM,2) :: rot !! Two-body equivalent principal axes moments of inertia the collider bodies prior to collision real(DP), dimension(NDIM,2) :: Lspin !! Two-body equivalent spin angular momentum vectors of the collider bodies prior to collision real(DP), dimension(NDIM,2) :: Lorbit !! Two-body equivalent orbital angular momentum vectors of the collider bodies prior to collision @@ -74,7 +76,7 @@ module collision real(DP), dimension(NDIM) :: rbcom !! Center of mass position vector of the collider nbody_system in nbody_system barycentric coordinates real(DP), dimension(NDIM) :: vbcom !! Velocity vector of the center of mass of the collider nbody_system in nbody_system barycentric coordinates real(DP), dimension(NDIM) :: rbimp !! Impact point position vector of the collider nbody_system in nbody_system barycentric coordinates - real(DP), dimension(NDIM) :: vbimp !! The impact point velocity vector is the component of the velocity in the distance vector direction + real(DP), dimension(NDIM) :: bounce_unit !! The impact point velocity vector is the component of the velocity in the distance vector direction contains procedure :: consolidate => collision_resolve_consolidate_impactors !! Consolidates a multi-body collision into an equivalent 2-body collision @@ -87,26 +89,27 @@ module collision !> Class definition for the variables that describe a collection of fragments in barycentric coordinates type, extends(base_multibody) :: collision_fragments - real(DP) :: mtot !! Total mass of fragments - class(base_particle_info), dimension(:), allocatable :: info !! Particle metadata information - integer(I4B), dimension(nbody) :: status !! An integrator-specific status indicator - real(DP), dimension(NDIM,nbody) :: rh !! Heliocentric position - real(DP), dimension(NDIM,nbody) :: vh !! Heliocentric velocity - real(DP), dimension(NDIM,nbody) :: rb !! Barycentric position - real(DP), dimension(NDIM,nbody) :: vb !! Barycentric velocity - real(DP), dimension(NDIM,nbody) :: rot !! rotation vectors of fragments - real(DP), dimension(NDIM,nbody) :: Ip !! Principal axes moment of inertia for fragments - real(DP), dimension(nbody) :: mass !! masses of fragments - real(DP), dimension(nbody) :: radius !! Radii of fragments - real(DP), dimension(nbody) :: density !! Radii of fragments - real(DP), dimension(NDIM,nbody) :: rc !! Position vectors in the collision coordinate frame - real(DP), dimension(NDIM,nbody) :: vc !! Velocity vectors in the collision coordinate frame - real(DP), dimension(nbody) :: rmag !! Array of radial distance magnitudes of individual fragments in the collisional coordinate frame - real(DP), dimension(nbody) :: vmag !! Array of radial distance magnitudes of individual fragments in the collisional coordinate frame - real(DP), dimension(nbody) :: rotmag !! Array of rotation magnitudes of individual fragments - real(DP), dimension(NDIM,nbody) :: v_r_unit !! Array of radial direction unit vectors of individual fragments in the collisional coordinate frame - real(DP), dimension(NDIM,nbody) :: v_t_unit !! Array of tangential direction unit vectors of individual fragments in the collisional coordinate frame - real(DP), dimension(NDIM,nbody) :: v_n_unit !! Array of normal direction unit vectors of individual fragments in the collisional coordinate frame + real(DP) :: mtot !! Total mass of fragments + class(base_particle_info), dimension(:), allocatable :: info !! Particle metadata information + integer(I4B), dimension(nbody) :: status !! An integrator-specific status indicator + real(DP), dimension(NDIM,nbody) :: rh !! Heliocentric position + real(DP), dimension(NDIM,nbody) :: vh !! Heliocentric velocity + real(DP), dimension(NDIM,nbody) :: rb !! Barycentric position + real(DP), dimension(NDIM,nbody) :: vb !! Barycentric velocity + real(DP), dimension(NDIM,nbody) :: rot !! rotation vectors of fragments + real(DP), dimension(NDIM,nbody) :: Ip !! Principal axes moment of inertia for fragments + real(DP), dimension(nbody) :: mass !! masses of fragments + real(DP), dimension(nbody) :: radius !! Radii of fragments + real(DP), dimension(nbody) :: density !! Radii of fragments + real(DP), dimension(NDIM,nbody) :: rc !! Position vectors in the collision coordinate frame + real(DP), dimension(NDIM,nbody) :: vc !! Velocity vectors in the collision coordinate frame + real(DP), dimension(nbody) :: rmag !! Array of radial distance magnitudes of individual fragments in the collisional coordinate frame + real(DP), dimension(nbody) :: vmag !! Array of radial distance magnitudes of individual fragments in the collisional coordinate frame + real(DP), dimension(nbody) :: rotmag !! Array of rotation magnitudes of individual fragments + real(DP), dimension(NDIM,nbody) :: v_r_unit !! Array of radial direction unit vectors of individual fragments in the collisional coordinate frame + real(DP), dimension(NDIM,nbody) :: v_t_unit !! Array of tangential direction unit vectors of individual fragments in the collisional coordinate frame + real(DP), dimension(NDIM,nbody) :: v_n_unit !! Array of normal direction unit vectors of individual fragments in the collisional coordinate frame + integer(I1B), dimension(nbody) :: origin_body !! Array of indices indicating which impactor body (1 or 2) the fragment originates from contains procedure :: reset => collision_util_reset_fragments !! Deallocates all allocatable arrays and sets everything else to 0 final :: collision_final_fragments !! Finalizer deallocates all allocatables @@ -253,10 +256,9 @@ module subroutine collision_generate_simple(self, nbody_system, param, t) real(DP), intent(in) :: t !! The time of the collision end subroutine collision_generate_simple - module subroutine collision_generate_simple_pos_vec(collider, r_max_start) + module subroutine collision_generate_simple_pos_vec(collider) implicit none class(collision_simple_disruption), intent(inout) :: collider !! Collision system object - real(DP), intent(in) :: r_max_start !! The maximum radial distance of fragments for disruptive collisions end subroutine collision_generate_simple_pos_vec module subroutine collision_generate_simple_rot_vec(collider) diff --git a/src/collision/collision_resolve.f90 b/src/collision/collision_resolve.f90 index f29e6abce..9e7e1e468 100644 --- a/src/collision/collision_resolve.f90 +++ b/src/collision/collision_resolve.f90 @@ -179,10 +179,10 @@ module subroutine collision_resolve_extract_plpl(self, nbody_system, param) class is (swiftest_nbody_system) select type (pl => nbody_system%pl) class is (swiftest_pl) - associate(plpl_encounter => self, idx1 => self%index1, idx2 => self%index2, plparent => pl%kin%parent) - nplplenc = plpl_encounter%nenc + associate(idx1 => self%index1, idx2 => self%index2, plparent => pl%kin%parent) + nplplenc = self%nenc allocate(lplpl_collision(nplplenc)) - lplpl_collision(:) = plpl_encounter%status(1:nplplenc) == COLLIDED + lplpl_collision(:) = self%status(1:nplplenc) == COLLIDED if (.not.any(lplpl_collision)) return ! Collisions have been detected in this step. So we need to determine which of them are between unique bodies. @@ -221,7 +221,7 @@ module subroutine collision_resolve_extract_plpl(self, nbody_system, param) ! Create a mask that contains only the pl-pl encounters that did not result in a collision, and then discard them lplpl_collision(:) = .false. lplpl_collision(collision_idx(:)) = .true. - call plpl_encounter%spill(nbody_system%plpl_collision, lplpl_collision, ldestructive=.true.) ! Extract any encounters that are not collisions from the list. + call self%spill(nbody_system%plpl_collision, lplpl_collision, ldestructive=.true.) ! Extract any encounters that are not collisions from the list. end associate end select end select @@ -483,9 +483,9 @@ module subroutine collision_resolve_plpl(self, nbody_system, param, t, dt, irec) !! implicit none ! Arguments - class(collision_list_plpl), intent(inout) :: self !! Swiftest pl-pl encounter list - class(base_nbody_system), intent(inout) :: nbody_system !! Swiftest nbody system object - class(base_parameters), intent(inout) :: param !! Current run configuration parameters with Swiftest additions + class(collision_list_plpl), intent(inout) :: self !! Swiftest pl-pl encounter list + class(base_nbody_system), intent(inout) :: nbody_system !! Swiftest nbody system object + class(base_parameters), intent(inout) :: param !! Current run configuration parameters with Swiftest additions real(DP), intent(in) :: t !! Current simulation time real(DP), intent(in) :: dt !! Current simulation step size integer(I4B), intent(in) :: irec !! Current recursion level @@ -504,7 +504,7 @@ module subroutine collision_resolve_plpl(self, nbody_system, param, t, dt, irec) class is (swiftest_pl) select type(param) class is (swiftest_parameters) - associate(plpl_encounter => self, plpl_collision => nbody_system%plpl_collision, & + associate(plpl_collision => nbody_system%plpl_collision, & collision_history => nbody_system%collision_history, pl => nbody_system%pl, cb => nbody_system%cb, & collider => nbody_system%collider, fragments => nbody_system%collider%fragments, impactors => nbody_system%collider%impactors) if (plpl_collision%nenc == 0) return ! No collisions to resolve @@ -562,7 +562,7 @@ module subroutine collision_resolve_plpl(self, nbody_system, param, t, dt, irec) call nbody_system%pl_adds%setup(0, param) ! Check whether or not any of the particles that were just added are themselves in a collision state. This will generate a new plpl_collision - call plpl_encounter%collision_check(nbody_system, param, t, dt, irec, lplpl_collision) + call self%collision_check(nbody_system, param, t, dt, irec, lplpl_collision) if (.not.lplpl_collision) exit if (loop == MAXCASCADE) then diff --git a/src/collision/collision_util.f90 b/src/collision/collision_util.f90 index d1cd1bf99..76b48ea9e 100644 --- a/src/collision/collision_util.f90 +++ b/src/collision/collision_util.f90 @@ -455,12 +455,23 @@ module subroutine collision_util_set_coordinate_impactors(self) mtot = sum(impactors%mass(:)) impactors%rbcom(:) = (impactors%mass(1) * impactors%rb(:,1) + impactors%mass(2) * impactors%rb(:,2)) / mtot impactors%vbcom(:) = (impactors%mass(1) * impactors%vb(:,1) + impactors%mass(2) * impactors%vb(:,2)) / mtot + + ! The center of mass coordinate position and velocities + impactors%rc(:,1) = impactors%rb(:,1) - impactors%rbcom(:) + impactors%rc(:,2) = impactors%rb(:,2) - impactors%rbcom(:) + impactors%vc(:,1) = impactors%vb(:,1) - impactors%vbcom(:) + impactors%vc(:,2) = impactors%vb(:,2) - impactors%vbcom(:) ! Find the point of impact between the two bodies - impactors%rbimp(:) = impactors%rb(:,1) + impactors%radius(1) * impactors%y_unit(:) + impactors%rbimp(:) = impactors%rb(:,1) + impactors%radius(1) * impactors%y_unit(:) - impactors%rbcom(:) + + ! Set the velocity direction as the "bounce" direction" for disruptions, and body 2's direction for hit and runs + if (impactors%regime == COLLRESOLVE_REGIME_HIT_AND_RUN) then + impactors%bounce_unit(:) = .unit. impactors%vc(:,2) + else + impactors%bounce_unit(:) = .unit. (impactors%vc(:,2) - 2 * dot_product(impactors%vc(:,2),impactors%y_unit(:)) * impactors%y_unit(:)) + end if - ! The "bounce" unit vector is the projection of the velocity vector into the distance vector - impactors%vbimp(:) = dot_product(delta_v(:),impactors%y_unit(:)) * impactors%y_unit(:) end associate return @@ -478,10 +489,10 @@ module subroutine collision_util_set_mass_dist(self, param) class(collision_simple_disruption), intent(inout) :: self !! Fraggle collision system object class(base_parameters), intent(in) :: param !! Current Swiftest run configuration parameters ! Internals - integer(I4B) :: i, jproj, jtarg, nfrag, istart + 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 + real(DP) :: mfrag, mremaining, min_mfrag, mtot, mcumul real(DP), parameter :: BETA = 2.85_DP 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) @@ -489,6 +500,7 @@ module subroutine collision_util_set_mass_dist(self, param) integer(I4B), parameter :: iMlr = 1 integer(I4B), parameter :: iMslr = 2 integer(I4B), parameter :: iMrem = 3 + logical :: flipper associate(impactors => self%impactors) ! Get mass weighted mean of Ip and density @@ -596,7 +608,31 @@ module subroutine collision_util_set_mass_dist(self, param) fragments%Ip(:, i) = Ip_avg(:) end do + ! For catastrophic impacts, we will assign each of the n>2 fragments to one of the two original bodies so that the fragment cloud occupies + ! roughly the same space as both original bodies. For all other disruption cases, we use body 2 as the center of the cloud. + fragments%origin_body(1) = 1 + fragments%origin_body(2) = 2 + if (impactors%regime == COLLRESOLVE_REGIME_SUPERCATASTROPHIC) then + mcumul = fragments%mass(1) + flipper = .true. + j = 2 + do i = 1, nfrag + if (flipper .and. (mcumul < impactors%mass(1))) then + flipper = .false. + j = 1 + else + j = 2 + flipper = .true. + end if + fragments%origin_body(i) = j + end do + else + fragments%origin_body(3:nfrag) = 2 + end if + end select + + end associate return diff --git a/src/fraggle/fraggle_generate.f90 b/src/fraggle/fraggle_generate.f90 index 0d8331970..21a297c08 100644 --- a/src/fraggle/fraggle_generate.f90 +++ b/src/fraggle/fraggle_generate.f90 @@ -94,7 +94,7 @@ module subroutine fraggle_generate_disrupt(self, nbody_system, param, t, lfailur lfailure_local = .false. call ieee_set_flag(ieee_all, .false.) ! Set all fpe flags to quiet - call collision_generate_simple_pos_vec(self, r_max_start) + call collision_generate_simple_pos_vec(self) call self%set_coordinate_system() ! Initial velocity guess will be the barycentric velocity of the colliding nbody_system so that the budgets are based on the much smaller collisional-frame velocities diff --git a/src/fraggle/fraggle_util.f90 b/src/fraggle/fraggle_util.f90 index 2e2e55961..43823642c 100644 --- a/src/fraggle/fraggle_util.f90 +++ b/src/fraggle/fraggle_util.f90 @@ -211,7 +211,7 @@ module subroutine fraggle_util_set_natural_scale_factors(self) impactors%rbcom(:) = impactors%rbcom(:) / collider%dscale impactors%vbcom(:) = impactors%vbcom(:) / collider%vscale impactors%rbimp(:) = impactors%rbimp(:) / collider%dscale - impactors%vbimp(:) = impactors%vbimp(:) / collider%vscale + impactors%bounce_unit(:) = impactors%bounce_unit(:) / collider%vscale impactors%rb(:,:) = impactors%rb(:,:) / collider%dscale impactors%vb(:,:) = impactors%vb(:,:) / collider%vscale impactors%mass(:) = impactors%mass(:) / collider%mscale @@ -254,7 +254,7 @@ module subroutine fraggle_util_set_original_scale_factors(self) impactors%rbcom(:) = impactors%rbcom(:) * collider%dscale impactors%vbcom(:) = impactors%vbcom(:) * collider%vscale impactors%rbimp(:) = impactors%rbimp(:) * collider%dscale - impactors%vbimp(:) = impactors%vbimp(:) * collider%vscale + impactors%bounce_unit(:) = impactors%bounce_unit(:) * collider%vscale impactors%mass = impactors%mass * collider%mscale impactors%radius = impactors%radius * collider%dscale diff --git a/src/swiftest/swiftest_util.f90 b/src/swiftest/swiftest_util.f90 index 2d207b1f7..b17dce303 100644 --- a/src/swiftest/swiftest_util.f90 +++ b/src/swiftest/swiftest_util.f90 @@ -1661,8 +1661,6 @@ module subroutine swiftest_util_rearray_pl(self, nbody_system, param) logical, dimension(:), allocatable :: lmask, ldump_mask class(encounter_list), allocatable :: plplenc_old logical :: lencounter - integer(I4B), dimension(:), allocatable :: levelg_orig_pl, levelm_orig_pl, levelg_orig_tp, levelm_orig_tp - integer(I4B), dimension(:), allocatable :: nplenc_orig_pl, nplenc_orig_tp, ntpenc_orig_pl associate(pl => self, tp => nbody_system%tp, pl_adds => nbody_system%pl_adds) @@ -1742,23 +1740,9 @@ module subroutine swiftest_util_rearray_pl(self, nbody_system, param) class is (symba_pl) select type(tp) class is (symba_tp) - ! allocate(levelg_orig_pl, source=pl%levelg) - ! allocate(levelm_orig_pl, source=pl%levelm) - ! allocate(nplenc_orig_tp, source=tp%nplenc) - ! call move_alloc(levelg_orig_pl, pl%levelg) - ! call move_alloc(levelm_orig_pl, pl%levelm) - ! call move_alloc(nplenc_orig_pl, pl%nplenc) lencounter = pl%encounter_check(param, nbody_system, param%dt, nbody_system%irec) if (tp%nbody > 0) then - ! allocate(levelg_orig_tp, source=tp%levelg) - ! allocate(levelm_orig_tp, source=tp%levelm) - ! allocate(nplenc_orig_tp, source=tp%nplenc) - ! allocate(ntpenc_orig_pl, source=pl%ntpenc) lencounter = tp%encounter_check(param, nbody_system, param%dt, nbody_system%irec) - ! call move_alloc(levelg_orig_tp, tp%levelg) - ! call move_alloc(levelm_orig_tp, tp%levelm) - ! call move_alloc(nplenc_orig_tp, tp%nplenc) - ! call move_alloc(ntpenc_orig_pl, pl%ntpenc) end if end select end select @@ -2523,7 +2507,11 @@ module subroutine swiftest_util_set_rhill(self,cb) if (self%nbody == 0) return call self%xv2el(cb) - self%rhill(1:self%nbody) = self%a(1:self%nbody) * (self%Gmass(1:self%nbody) / cb%Gmass / 3)**THIRD + where(self%a(1:self%nbody) > 0.0_DP) + self%rhill(1:self%nbody) = self%a(1:self%nbody) * (self%Gmass(1:self%nbody) / cb%Gmass / 3)**THIRD + elsewhere + self%rhill(1:self%nbody) = (.mag.self%rh(:,1:self%nbody)) * (self%Gmass(1:self%nbody) / cb%Gmass / 3)**THIRD + end where return end subroutine swiftest_util_set_rhill