From ba80e3ca4d8586d18f3b7518fbdca5e85fe24882 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Fri, 23 Dec 2022 12:07:58 -0500 Subject: [PATCH] Moved the fragment position and initial velocity distribution code into the "simple disruption" model. This will allow these simpler distributions to be tested before adding in the energy and momentum constraints that Fraggle uses. --- examples/Fragmentation/Fragmentation_Movie.py | 4 +- src/collision/collision_generate.f90 | 137 +++++++++++- src/collision/collision_module.f90 | 70 ++++--- src/collision/collision_util.f90 | 34 ++- src/fraggle/fraggle_generate.f90 | 198 ++++-------------- src/fraggle/fraggle_module.f90 | 7 - src/fraggle/fraggle_util.f90 | 32 +-- 7 files changed, 256 insertions(+), 226 deletions(-) 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(:)