From f8066b91f8689adce20a32fdd3f82b2448a915e4 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Mon, 13 Feb 2023 09:30:27 -0500 Subject: [PATCH 1/3] Put a limiter on the outer loop to keep Fraggle from churning away needlessly when nfrag is very large --- src/fraggle/fraggle_generate.f90 | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/src/fraggle/fraggle_generate.f90 b/src/fraggle/fraggle_generate.f90 index 7c0b922f7..f34ce0fe2 100644 --- a/src/fraggle/fraggle_generate.f90 +++ b/src/fraggle/fraggle_generate.f90 @@ -528,7 +528,8 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu real(DP), parameter :: hitandrun_vscale = 0.25_DP real(DP) :: vmin_guess real(DP) :: vmax_guess - integer(I4B), parameter :: MAXINNER = 20 + integer(I4B), parameter :: MAXLOOP = 20 + integer(I4B), parameter :: MAXTRY = 10 integer(I4B), parameter :: MAXANGMTM = 1000 class(collision_fraggle), allocatable :: collider_local character(len=STRMAX) :: message @@ -555,7 +556,7 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu dE_best = huge(1.0_DP) nsteps_best = 0 nsteps = 0 - outer: do try = 1, nfrag - istart - 1 + outer: do try = 1, min(nfrag - istart - 1, MAXTRY) associate(fragments => collider_local%fragments) if (allocated(vsign)) deallocate(vsign); allocate(vsign(fragments%nbody)) if (allocated(vscale)) deallocate(vscale); allocate(vscale(fragments%nbody)) @@ -611,7 +612,7 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu call fragments%set_coordinate_system() E_residual = huge(1.0_DP) - inner: do loop = 1, MAXINNER + inner: do loop = 1, MAXLOOP nsteps = nsteps + 1 mfrag = sum(fragments%mass(istart:fragments%nbody)) From 23f92deedd2b2f87a2d6ea561d53fa3c8f1a1604 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Mon, 13 Feb 2023 09:31:11 -0500 Subject: [PATCH 2/3] Improved disruptive hit and run fragment positioning and target rotation --- src/fraggle/fraggle_generate.f90 | 48 ++++++++++++++++++++++---------- 1 file changed, 33 insertions(+), 15 deletions(-) diff --git a/src/fraggle/fraggle_generate.f90 b/src/fraggle/fraggle_generate.f90 index f34ce0fe2..75ab1521b 100644 --- a/src/fraggle/fraggle_generate.f90 +++ b/src/fraggle/fraggle_generate.f90 @@ -324,8 +324,10 @@ module subroutine fraggle_generate_pos_vec(collider, nbody_system, param, lfailu ! 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 loverlap(:) = .true. + istart = 3 if (lhitandrun) then - rdistance = 1.0_DP + rdistance = impactors%radius(2) + istart = 2 else if (lsupercat) then rdistance = 0.5_DP * sum(impactors%radius(:)) else @@ -342,20 +344,15 @@ module subroutine fraggle_generate_pos_vec(collider, nbody_system, param, lfailu 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 + if (lhitandrun) then ! Keep the target unchanged and set the 2nd fragment cloud to be centered on the projectile 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(:) + fragment_cloud_center(:,1) = impactors%rc(:,1) + fragment_cloud_center(:,2) = impactors%rc(:,2) + fragments%rc(:,1) = fragment_cloud_center(:,1) 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(:) @@ -379,18 +376,20 @@ module subroutine fraggle_generate_pos_vec(collider, nbody_system, param, lfailu ! Make a random cloud phi(i) = TWOPI * phi(i) theta(i) = acos(2 * theta(i) - 1.0_DP) + ! Scale the cloud size fragments%rmag(i) = fragment_cloud_radius(j) * mass_rscale(i) * u(i)**(THIRD) + ! Position the fragment in a random point within the cloud 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)) ! Shift to the cloud center coordinates fragments%rc(:,i) = fragments%rc(:,i) + fragment_cloud_center(:,j) - if (lhitandrun) then ! Stretch out the hit and run cloud along the flight trajectory - fragments%rc(:,i) = fragments%rc(:,i) + cloud_size_scale_factor * rdistance * fragments%rmag(i) * impactors%bounce_unit(:) - end if + + ! Stretch out the hit and run cloud along the flight trajectory + if (lhitandrun) fragments%rc(:,i) = fragments%rc(:,i) + cloud_size_scale_factor * rdistance * u(i) * impactors%bounce_unit(:) ! Make sure that the fragments are positioned away from the impact point direction = dot_product(fragments%rc(:,i) - impactors%rcimp(:), fragment_cloud_center(:,j) - impactors%rcimp(:)) @@ -461,11 +460,14 @@ module subroutine fraggle_generate_rot_vec(collider, nbody_system, param) ! 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), parameter :: hitandrun_momentum_transfer = 0.01_DP ! Fraction of projectile momentum transfered to target in a hit and run real(DP) :: mass_fac real(DP), dimension(NDIM) :: drot, dL integer(I4B), parameter :: MAXLOOP = 10 + logical :: lhitandrun associate(fragments => collider%fragments, impactors => collider%impactors, nfrag => collider%fragments%nbody) + lhitandrun = (impactors%regime == COLLRESOLVE_REGIME_HIT_AND_RUN) ! 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 @@ -488,6 +490,22 @@ module subroutine fraggle_generate_rot_vec(collider, nbody_system, param) end do fragments%rot(:,1) = fragments%rot(:,1) + drot(:) end if + + if (lhitandrun) then + dL(:) = hitandrun_momentum_transfer * impactors%mass(2) * (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)) + 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)) do concurrent (i = 2:nfrag) mass_fac = fragments%mass(i) / impactors%mass(2) @@ -564,8 +582,8 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu if (lhitandrun) then vesc = sqrt(2 * sum(fragments%Gmass(istart:fragments%nbody)) / sum(fragments%radius(istart:fragments%nbody))) - vmin_guess = .mag.impactors%vc(:,2) - vimp * hitandrun_vscale - vmax_guess = .mag.impactors%vc(:,2) + vimp * hitandrun_vscale + vmin_guess = .mag.impactors%vc(:,2) * (1.0_DP - hitandrun_vscale) + vmax_guess = .mag.impactors%vc(:,2) * (1.0_DP + hitandrun_vscale) else vesc = sqrt(2 * sum(impactors%Gmass(:)) / sum(impactors%radius(:))) vmin_guess = 1.001_DP * vesc From ee1ccb7da0f3693a2974b3a704a07f0525ca5896 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Mon, 13 Feb 2023 09:31:43 -0500 Subject: [PATCH 3/3] Adjusted initial conditions of hit and run to showcase momentum transfer --- examples/Fragmentation/Fragmentation_Movie.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/examples/Fragmentation/Fragmentation_Movie.py b/examples/Fragmentation/Fragmentation_Movie.py index 33178436d..5afb9484d 100755 --- a/examples/Fragmentation/Fragmentation_Movie.py +++ b/examples/Fragmentation/Fragmentation_Movie.py @@ -84,10 +84,10 @@ np.array([0.0, 0.0, -5.0e5])], "supercatastrophic_off_axis": [np.array([0.0, 0.0, 1.0e5]), np.array([0.0, 0.0, -1.0e4])], - "hitandrun_disrupt" : [np.array([0.0, 0.0, -1.0e5]), + "hitandrun_disrupt" : [np.array([0.0, 0.0, 0.0]), + np.array([0.0, 0.0, 1.0e5])], + "hitandrun_pure" : [np.array([0.0, 0.0, 0.0]), np.array([0.0, 0.0, 1.0e5])], - "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])], "merge_spinner" : [np.array([0.0, 0.0, -1.2e6]),