diff --git a/examples/Fragmentation/Fragmentation_Movie.py b/examples/Fragmentation/Fragmentation_Movie.py index 374574e9c..381bca121 100644 --- a/examples/Fragmentation/Fragmentation_Movie.py +++ b/examples/Fragmentation/Fragmentation_Movie.py @@ -206,8 +206,8 @@ def data_stream(self, frame=0): # Set fragmentation parameters minimum_fragment_gmass = 0.2 * body_Gmass[style][1] # Make the minimum fragment mass a fraction of the smallest body gmtiny = 0.99 * body_Gmass[style][1] # Make GMTINY just smaller than the smallest original body. This will prevent runaway collisional cascades - sim.set_parameter(collision_model="fraggle", encounter_save="both", gmtiny=gmtiny, minimum_fragment_gmass=minimum_fragment_gmass, verbose=False) + sim.set_parameter(collision_model="simple", encounter_save="both", gmtiny=gmtiny, minimum_fragment_gmass=minimum_fragment_gmass, verbose=False) sim.run(dt=1e-3, tstop=1.0e-3, istep_out=1, dump_cadence=0) print("Generating animation") - anim = AnimatedScatter(sim,movie_filename,movie_titles[style],style,nskip=1) \ No newline at end of file + anim = AnimatedScatter(sim,movie_filename,movie_titles[style],style,nskip=1) diff --git a/src/collision/collision_generate.f90 b/src/collision/collision_generate.f90 index 2be2408bb..a3ff0cbca 100644 --- a/src/collision/collision_generate.f90 +++ b/src/collision/collision_generate.f90 @@ -123,7 +123,7 @@ end subroutine collision_generate_hitandrun module subroutine collision_generate_merge(self, nbody_system, param, t) !! author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton !! - !! Merge massive bodies in any collisionals ystem. + !! Merge massive bodies in any collisional system. !! !! Adapted from David E. Kaufmann's Swifter routines symba_merge_pl.f90 and symba_discard_merge_pl.f90 !! @@ -221,11 +221,146 @@ end subroutine collision_generate_merge module subroutine collision_generate_simple_disruption(self, nbody_system, param, t) + !! author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton + !! + !! Generates a simple fragment position and velocity distribution based on the collision + !! regime. It makes no attempt to constrain the energy of the collision 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 + ! Internals + real(DP) :: r_max_start + + associate(impactors => self%impactors, fragments => self%fragments, nfrag => self%fragments%nbody) + r_max_start = 1.1_DP * .mag.(impactors%rb(:,2) - impactors%rb(:,1)) + call collision_generate_simple_pos_vec(self, r_max_start) + call self%set_coordinate_system() + call collision_generate_simple_vel_vec(self) + end associate + + return end subroutine collision_generate_simple_disruption + + module subroutine collision_generate_simple_pos_vec(collider, r_max_start) + !! 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. + !! For hit and run with disruption, the fragments are generated in a random cloud around the smallest of the two colliders (body 2) + !! For disruptive collisions, the fragments are generated in a random cloud around the impact point. Bodies are checked for overlap and + !! regenerated if they overlap. + 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 + + 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 + + ! 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(:)) + + ! 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(:)) + + ! 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 + 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 + end if + end do + + ! Check for any overlapping bodies. + loverlap(:) = .false. + do j = 1, nfrag + do i = j + 1, nfrag + dis = .mag.(fragments%rc(:,j) - fragments%rc(:,i)) + loverlap(i) = loverlap(i) .or. (dis <= (fragments%radius(i) + fragments%radius(j))) + end do + end do + end do + call collision_util_shift_vector_to_origin(fragments%mass, fragments%rc) + call collider%set_coordinate_system() + + do concurrent(i = 1:nfrag) + fragments%rb(:,i) = fragments%rc(:,i) + impactors%rbcom(:) + end do + + impactors%rbcom(:) = 0.0_DP + do concurrent(i = 1:nfrag) + impactors%rbcom(:) = impactors%rbcom(:) + fragments%mass(i) * fragments%rb(:,i) + end do + impactors%rbcom(:) = impactors%rbcom(:) / fragments%mtot + end associate + + return + end subroutine collision_generate_simple_pos_vec + + + module subroutine collision_generate_simple_vel_vec(collider) + !! Author: David A. Minton + !! + !! Generates an initial "guess" for the velocitity vectors + implicit none + ! Arguments + class(collision_simple_disruption), intent(inout) :: collider !! Fraggle collision system object + ! Internals + integer(I4B) :: i, j + logical :: lcat + real(DP), dimension(NDIM) :: vimp_unit + real(DP), dimension(NDIM,collider%fragments%nbody) :: vnoise + real(DP), parameter :: VNOISE_MAG = 0.10_DP + real(DP) :: vmag + + associate(fragments => collider%fragments, impactors => collider%impactors, nfrag => collider%fragments%nbody) + lcat = (impactors%regime == COLLRESOLVE_REGIME_SUPERCATASTROPHIC) + + ! "Bounce" the first two bodies + do i = 1,2 + fragments%vc(:,i) = impactors%vb(:,i) - impactors%vbcom(:) - 2 * impactors%vbimp(:) + end do + + vmag = .mag.impactors%vbcom(:) + vimp_unit(:) = .unit. impactors%vbimp(:) + call random_number(vnoise) + vnoise = (2 * vnoise - 1.0_DP) * vmag + do i = 3, nfrag + vimp_unit(:) = .unit. (fragments%rc(:,i) + impactors%rbcom(:) - impactors%rbimp(:)) + fragments%vc(:,i) = vmag * vimp_unit(:) + vnoise(:,i) + end do + end associate + return + end subroutine collision_generate_simple_vel_vec + + + + end submodule s_collision_generate \ No newline at end of file diff --git a/src/collision/collision_module.f90 b/src/collision/collision_module.f90 index 570e1cccc..ef75b56e9 100644 --- a/src/collision/collision_module.f90 +++ b/src/collision/collision_module.f90 @@ -233,6 +233,17 @@ module subroutine collision_generate_simple_disruption(self, nbody_system, param real(DP), intent(in) :: t !! The time of the collision end subroutine collision_generate_simple_disruption + module subroutine collision_generate_simple_pos_vec(collider, r_max_start) + implicit none + 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 + end subroutine collision_generate_simple_pos_vec + + module subroutine collision_generate_simple_vel_vec(collider) + implicit none + class(collision_simple_disruption), intent(inout) :: collider !! Fraggle collision system object + end subroutine collision_generate_simple_vel_vec + module subroutine collision_io_collider_message(pl, collidx, collider_message) implicit none class(base_object), intent(in) :: pl !! Swiftest massive body object @@ -349,32 +360,6 @@ module subroutine collision_resolve_pltp(self, nbody_system, param, t, dt, irec) integer(I4B), intent(in) :: irec !! Current recursion level end subroutine collision_resolve_pltp - module subroutine collision_util_set_coordinate_collider(self) - implicit none - class(collision_merge), intent(inout) :: self !! collisional system - end subroutine collision_util_set_coordinate_collider - - module subroutine collision_util_set_coordinate_impactors(self) - implicit none - class(collision_impactors), intent(inout) :: self !! collisional system - end subroutine collision_util_set_coordinate_impactors - - module subroutine collision_util_setup_collider(self, nbody_system) - implicit none - class(collision_merge), intent(inout) :: self !! Encounter collision system object - class(base_nbody_system), intent(in) :: nbody_system !! Current nbody system. Used as a mold for the before/after snapshots - end subroutine collision_util_setup_collider - - module subroutine collision_util_setup_impactors_collider(self) - implicit none - class(collision_merge), intent(inout) :: self !! Encounter collision system object - end subroutine collision_util_setup_impactors_collider - - module subroutine collision_util_setup_fragments_collider(self, nfrag) - implicit none - class(collision_merge), intent(inout) :: self !! Encounter collision system object - integer(I4B), intent(in) :: nfrag !! Number of fragments to create - end subroutine collision_util_setup_fragments_collider module subroutine collision_util_add_fragments_to_collider(self, nbody_system, param) implicit none @@ -403,6 +388,39 @@ module subroutine collision_util_set_mass_dist(self, param) class(base_parameters), intent(in) :: param !! Current Swiftest run configuration parameters end subroutine collision_util_set_mass_dist + module subroutine collision_util_set_coordinate_collider(self) + implicit none + class(collision_merge), intent(inout) :: self !! collisional system + end subroutine collision_util_set_coordinate_collider + + module subroutine collision_util_set_coordinate_impactors(self) + implicit none + class(collision_impactors), intent(inout) :: self !! collisional system + end subroutine collision_util_set_coordinate_impactors + + module subroutine collision_util_setup_collider(self, nbody_system) + implicit none + class(collision_merge), intent(inout) :: self !! Encounter collision system object + class(base_nbody_system), intent(in) :: nbody_system !! Current nbody system. Used as a mold for the before/after snapshots + end subroutine collision_util_setup_collider + + module subroutine collision_util_setup_impactors_collider(self) + implicit none + class(collision_merge), intent(inout) :: self !! Encounter collision system object + end subroutine collision_util_setup_impactors_collider + + module subroutine collision_util_setup_fragments_collider(self, nfrag) + implicit none + class(collision_merge), intent(inout) :: self !! Encounter collision system object + integer(I4B), intent(in) :: nfrag !! Number of fragments to create + end subroutine collision_util_setup_fragments_collider + + module subroutine collision_util_shift_vector_to_origin(m_frag, vec_frag) + implicit none + real(DP), dimension(:), intent(in) :: m_frag !! Fragment masses + real(DP), dimension(:,:), intent(inout) :: vec_frag !! Fragment positions or velocities in the center of mass frame + end subroutine + module subroutine collision_util_get_idvalues_snapshot(self, idvals) implicit none class(collision_snapshot), intent(in) :: self !! Collision snapshot object diff --git a/src/collision/collision_util.f90 b/src/collision/collision_util.f90 index 06a39f4ed..da68b42d6 100644 --- a/src/collision/collision_util.f90 +++ b/src/collision/collision_util.f90 @@ -459,8 +459,8 @@ module subroutine collision_util_set_coordinate_impactors(self) ! Find the point of impact between the two bodies impactors%rbimp(:) = impactors%rb(:,1) + impactors%radius(1) * impactors%y_unit(:) - ! The "bounce" unit vector is the projection of - impactors%vbimp(:) = dot_product(delta_v(:),impactors%x_unit(:)) * impactors%x_unit(:) + ! 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 @@ -656,6 +656,36 @@ module subroutine collision_util_setup_fragments_collider(self, nfrag) end subroutine collision_util_setup_fragments_collider + module subroutine collision_util_shift_vector_to_origin(m_frag, vec_frag) + !! Author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton + !! + !! Adjusts the position or velocity of the fragments as needed to align them with the origin + implicit none + ! Arguments + real(DP), dimension(:), intent(in) :: m_frag !! Fragment masses + real(DP), dimension(:,:), intent(inout) :: vec_frag !! Fragment positions or velocities in the center of mass frame + + ! Internals + real(DP), dimension(NDIM) :: mvec_frag, COM_offset + integer(I4B) :: i, nfrag + real(DP) :: mtot + + mvec_frag(:) = 0.0_DP + mtot = sum(m_frag) + nfrag = size(m_frag) + + do i = 1, nfrag + mvec_frag = mvec_frag(:) + vec_frag(:,i) * m_frag(i) + end do + COM_offset(:) = -mvec_frag(:) / mtot + do i = 1, nfrag + vec_frag(:, i) = vec_frag(:, i) + COM_offset(:) + end do + + return + end subroutine collision_util_shift_vector_to_origin + + subroutine collision_util_save_snapshot(collision_history, snapshot) !! author: David A. Minton !! diff --git a/src/fraggle/fraggle_generate.f90 b/src/fraggle/fraggle_generate.f90 index b809e45f4..4424b2412 100644 --- a/src/fraggle/fraggle_generate.f90 +++ b/src/fraggle/fraggle_generate.f90 @@ -282,7 +282,7 @@ module subroutine fraggle_generate_fragments(collider, nbody_system, param, lfai lfailure = .false. call ieee_set_flag(ieee_all, .false.) ! Set all fpe flags to quiet - call fraggle_generate_pos_vec(collider, r_max_start) + call collision_generate_simple_pos_vec(collider, r_max_start) call collider%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 @@ -293,46 +293,46 @@ module subroutine fraggle_generate_fragments(collider, nbody_system, param, lfai call collider%get_energy_and_momentum(nbody_system, param, lbefore=.false.) call collider%set_budgets() - call fraggle_generate_vel_vec_guess(collider) - - ! call fraggle_generate_tan_vel(collider, lfailure) - ! if (lfailure) then - ! call swiftest_io_log_one_message(COLLISION_LOG_OUT, "Fraggle failed to find tangential velocities") - ! cycle - ! end if - - ! call fraggle_generate_rad_vel(collider, lfailure) - ! if (lfailure) then - ! call swiftest_io_log_one_message(COLLISION_LOG_OUT, "Fraggle failed to find radial velocities") - ! cycle - ! end if - - ! call fraggle_generate_spins(collider, f_spin, lfailure) - ! if (lfailure) then - ! call swiftest_io_log_one_message(COLLISION_LOG_OUT, "Fraggle failed to find spins") - ! cycle - ! end if - - ! call collider%get_energy_and_momentum(nbody_system, param, lbefore=.false.) - ! dEtot = collider%Etot(2) - collider%Etot(1) - ! dLmag = .mag. (collider%Ltot(:,2) - collider%Ltot(:,1)) - ! exit - - ! lfailure = ((abs(dEtot + impactors%Qloss) > FRAGGLE_ETOL) .or. (dEtot > 0.0_DP)) - ! if (lfailure) then - ! write(message, *) dEtot, abs(dEtot + impactors%Qloss) / FRAGGLE_ETOL - ! call swiftest_io_log_one_message(COLLISION_LOG_OUT, "Fraggle failed due to high energy error: " // & - ! trim(adjustl(message))) - ! cycle - ! end if - - ! lfailure = ((abs(dLmag) / (.mag.collider%Ltot(:,1))) > FRAGGLE_LTOL) - ! if (lfailure) then - ! write(message,*) dLmag / (.mag.collider%Ltot(:,1)) - ! call swiftest_io_log_one_message(COLLISION_LOG_OUT, "Fraggle failed due to high angular momentum error: " // & - ! trim(adjustl(message))) - ! cycle - ! end if + call collision_generate_simple_vel_vec(collider) + + call fraggle_generate_tan_vel(collider, lfailure) + if (lfailure) then + call swiftest_io_log_one_message(COLLISION_LOG_OUT, "Fraggle failed to find tangential velocities") + cycle + end if + + call fraggle_generate_rad_vel(collider, lfailure) + if (lfailure) then + call swiftest_io_log_one_message(COLLISION_LOG_OUT, "Fraggle failed to find radial velocities") + cycle + end if + + call fraggle_generate_spins(collider, f_spin, lfailure) + if (lfailure) then + call swiftest_io_log_one_message(COLLISION_LOG_OUT, "Fraggle failed to find spins") + cycle + end if + + call collider%get_energy_and_momentum(nbody_system, param, lbefore=.false.) + dEtot = collider%Etot(2) - collider%Etot(1) + dLmag = .mag. (collider%Ltot(:,2) - collider%Ltot(:,1)) + exit + + lfailure = ((abs(dEtot + impactors%Qloss) > FRAGGLE_ETOL) .or. (dEtot > 0.0_DP)) + if (lfailure) then + write(message, *) dEtot, abs(dEtot + impactors%Qloss) / FRAGGLE_ETOL + call swiftest_io_log_one_message(COLLISION_LOG_OUT, "Fraggle failed due to high energy error: " // & + trim(adjustl(message))) + cycle + end if + + lfailure = ((abs(dLmag) / (.mag.collider%Ltot(:,1))) > FRAGGLE_LTOL) + if (lfailure) then + write(message,*) dLmag / (.mag.collider%Ltot(:,1)) + call swiftest_io_log_one_message(COLLISION_LOG_OUT, "Fraggle failed due to high angular momentum error: " // & + trim(adjustl(message))) + cycle + end if ! Check if any of the usual floating point exceptions happened, and fail the try if so call ieee_get_flag(ieee_usual, fpe_flag) @@ -364,122 +364,6 @@ module subroutine fraggle_generate_fragments(collider, nbody_system, param, lfai return end subroutine fraggle_generate_fragments - - subroutine fraggle_generate_pos_vec(collider, r_max_start) - !! Author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton - !! - !! Initializes the orbits of the fragments around the center of mass. The fragments are initially placed on a plane defined by the - !! pre-impact angular momentum. They are distributed on an ellipse surrounding the center of mass. - !! The initial positions do not conserve energy or momentum, so these need to be adjusted later. - implicit none - ! Arguments - class(collision_fraggle), intent(inout) :: collider !! Fraggle collision system object - real(DP), intent(in) :: r_max_start !! Initial guess for the starting maximum radial distance of fragments - ! Internals - real(DP) :: dis, rad, r_max, fdistort - logical, dimension(:), allocatable :: loverlap - integer(I4B) :: i, j - logical :: lnoncat, lhitandrun - - 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 - - ! 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(:)) - - ! 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(:)) - - ! 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 - 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 - end if - end do - - ! Check for any overlapping bodies. - loverlap(:) = .false. - do j = 1, nfrag - do i = j + 1, nfrag - dis = .mag.(fragments%rc(:,j) - fragments%rc(:,i)) - loverlap(i) = loverlap(i) .or. (dis <= (fragments%radius(i) + fragments%radius(j))) - end do - end do - end do - call fraggle_util_shift_vector_to_origin(fragments%mass, fragments%rc) - call collider%set_coordinate_system() - - do concurrent(i = 1:nfrag) - fragments%rb(:,i) = fragments%rc(:,i) + impactors%rbcom(:) - end do - - impactors%rbcom(:) = 0.0_DP - do concurrent(i = 1:nfrag) - impactors%rbcom(:) = impactors%rbcom(:) + fragments%mass(i) * fragments%rb(:,i) - end do - impactors%rbcom(:) = impactors%rbcom(:) / fragments%mtot - end associate - - return - end subroutine fraggle_generate_pos_vec - - - subroutine fraggle_generate_vel_vec_guess(collider) - !! Author: David A. Minton - !! - !! Generates an initial "guess" for the velocitity vectors - implicit none - ! Arguments - class(collision_fraggle), intent(inout) :: collider !! Fraggle collision system object - ! Internals - integer(I4B) :: i, j - logical :: lcat - real(DP), dimension(NDIM) :: vimp_unit - real(DP), dimension(NDIM,collider%fragments%nbody) :: vnoise - real(DP), parameter :: VNOISE_MAG = 0.10_DP - real(DP) :: vmag - - associate(fragments => collider%fragments, impactors => collider%impactors, nfrag => collider%fragments%nbody) - lcat = (impactors%regime == COLLRESOLVE_REGIME_SUPERCATASTROPHIC) - - ! "Bounce" the first two bodies - do i = 1,2 - fragments%vc(:,i) = impactors%vb(:,i) - impactors%vbcom(:) - 2 * impactors%vbimp(:) - end do - - vmag = .mag.impactors%vbcom(:) - vimp_unit(:) = .unit. impactors%vbimp(:) - call random_number(vnoise) - vnoise = (2 * vnoise - 1.0_DP) * vmag - do i = 3, nfrag - vimp_unit(:) = .unit. (fragments%rc(:,i) + impactors%rbcom(:) - impactors%rbimp(:)) - fragments%vc(:,i) = vmag * vimp_unit(:) + vnoise(:,i) - end do - end associate - return - end subroutine fraggle_generate_vel_vec_guess - - subroutine fraggle_generate_spins(collider, f_spin, lfailure) !! Author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton !! diff --git a/src/fraggle/fraggle_module.f90 b/src/fraggle/fraggle_module.f90 index 72c0e5dc3..7e364dc6f 100644 --- a/src/fraggle/fraggle_module.f90 +++ b/src/fraggle/fraggle_module.f90 @@ -143,13 +143,6 @@ module subroutine fraggle_util_set_original_scale_factors(self) class(collision_fraggle), intent(inout) :: self !! Fraggle collision system object end subroutine fraggle_util_set_original_scale_factors - - module subroutine fraggle_util_shift_vector_to_origin(m_frag, vec_frag) - implicit none - real(DP), dimension(:), intent(in) :: m_frag !! Fragment masses - real(DP), dimension(:,:), intent(inout) :: vec_frag !! Fragment positions or velocities in the center of mass frame - end subroutine - module function fraggle_util_vmag_to_vb(v_r_mag, v_r_unit, v_t_mag, v_t_unit, m_frag, vcom) result(vb) implicit none real(DP), dimension(:), intent(in) :: v_r_mag !! Unknown radial component of fragment velocity vector diff --git a/src/fraggle/fraggle_util.f90 b/src/fraggle/fraggle_util.f90 index 22691f17a..928e3655e 100644 --- a/src/fraggle/fraggle_util.f90 +++ b/src/fraggle/fraggle_util.f90 @@ -317,36 +317,6 @@ module subroutine fraggle_util_setup_fragments_system(self, nfrag) end subroutine fraggle_util_setup_fragments_system - module subroutine fraggle_util_shift_vector_to_origin(m_frag, vec_frag) - !! Author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton - !! - !! Adjusts the position or velocity of the fragments as needed to align them with the origin - implicit none - ! Arguments - real(DP), dimension(:), intent(in) :: m_frag !! Fragment masses - real(DP), dimension(:,:), intent(inout) :: vec_frag !! Fragment positions or velocities in the center of mass frame - - ! Internals - real(DP), dimension(NDIM) :: mvec_frag, COM_offset - integer(I4B) :: i, nfrag - real(DP) :: mtot - - mvec_frag(:) = 0.0_DP - mtot = sum(m_frag) - nfrag = size(m_frag) - - do i = 1, nfrag - mvec_frag = mvec_frag(:) + vec_frag(:,i) * m_frag(i) - end do - COM_offset(:) = -mvec_frag(:) / mtot - do i = 1, nfrag - vec_frag(:, i) = vec_frag(:, i) + COM_offset(:) - end do - - return - end subroutine fraggle_util_shift_vector_to_origin - - module function fraggle_util_vmag_to_vb(v_r_mag, v_r_unit, v_t_mag, v_t_unit, m_frag, vcom) result(vb) !! Author: David A. Minton !! @@ -370,7 +340,7 @@ module function fraggle_util_vmag_to_vb(v_r_mag, v_r_unit, v_t_mag, v_t_unit, m_ vb(:,i) = abs(v_r_mag(i)) * v_r_unit(:, i) end do ! In order to keep satisfying the kinetic energy constraint, we must shift the origin of the radial component of the velocities to the center of mass - call fraggle_util_shift_vector_to_origin(m_frag, vb) + call collision_util_shift_vector_to_origin(m_frag, vb) do i = 1, nfrag vb(:, i) = vb(:, i) + v_t_mag(i) * v_t_unit(:, i) + vcom(:)