diff --git a/examples/Fragmentation/Fragmentation_Movie.py b/examples/Fragmentation/Fragmentation_Movie.py index 8f0d7d22b..d17191912 100644 --- a/examples/Fragmentation/Fragmentation_Movie.py +++ b/examples/Fragmentation/Fragmentation_Movie.py @@ -86,7 +86,7 @@ body_Gmass = {"disruption_headon" : [1e-7, 1e-10], "disruption_off_axis" : [1e-7, 1e-10], - "supercatastrophic_headon": [1e-7, 1e-8], + "supercatastrophic_headon" : [1e-7, 1e-8], "supercatastrophic_off_axis": [1e-7, 1e-8], "hitandrun_disrupt" : [1e-7, 7e-10], "hitandrun_pure" : [1e-7, 7e-10] @@ -236,7 +236,7 @@ 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="disruption", encounter_save="both", gmtiny=gmtiny, minimum_fragment_gmass=minimum_fragment_gmass, verbose=False) + sim.set_parameter(collision_model="fraggle", 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") diff --git a/python/swiftest/swiftest/simulation_class.py b/python/swiftest/swiftest/simulation_class.py index 530d8b034..e66848183 100644 --- a/python/swiftest/swiftest/simulation_class.py +++ b/python/swiftest/swiftest/simulation_class.py @@ -224,7 +224,7 @@ def __init__(self,read_param: bool = False, read_old_output: bool = False, simdi general_relativity : bool, default True Include the post-Newtonian correction in acceleration calculations. Parameter input file equivalent: `GR` - collision_model: {"MERGE","BOUNCE","SIMPLE","FRAGGLE"}, default "MERGE" + collision_model: {"MERGE","BOUNCE","FRAGGLE"}, default "MERGE" This is used to set the collision/fragmentation model. [TODO: DESCRIBE THESE] This argument only applies to Swiftest-SyMBA simulations. It will be ignored otherwise. Parameter input file equivalent: `COLLISION_MODEL` @@ -1015,7 +1015,7 @@ def get_integrator(self,arg_list: str | List[str] | None = None, verbose: bool | def set_feature(self, close_encounter_check: bool | None = None, general_relativity: bool | None = None, - collision_model: Literal["MERGE","BOUNCE","SIMPLE","FRAGGLE"] | None = None, + collision_model: Literal["MERGE","BOUNCE","FRAGGLE"] | None = None, minimum_fragment_gmass: float | None = None, minimum_fragment_mass: float | None = None, rotation: bool | None = None, @@ -1048,7 +1048,7 @@ def set_feature(self, *WARNING*: Enabling this feature could lead to very large files. general_relativity : bool, optional Include the post-Newtonian correction in acceleration calculations. - collision_model: {"MERGE","BOUNCE","SIMPLE","FRAGGLE"}, default "MERGE" + collision_model: {"MERGE","BOUNCE","FRAGGLE"}, default "MERGE" This is used to set the collision/fragmentation model. [TODO: DESCRIBE THESE] This argument only applies to Swiftest-SyMBA simulations. It will be ignored otherwise. Parameter input file equivalent: `COLLISION_MODEL` @@ -1129,7 +1129,7 @@ def set_feature(self, self.param["GR"] = general_relativity update_list.append("general_relativity") - fragmentation_models = ["FRAGGLE", "SIMPLE"] + fragmentation_models = ["FRAGGLE"] if collision_model is not None: collision_model = collision_model.upper() fragmentation = collision_model in fragmentation_models diff --git a/src/collision/collision_generate.f90 b/src/collision/collision_generate.f90 index 007ac8d0c..53900f353 100644 --- a/src/collision/collision_generate.f90 +++ b/src/collision/collision_generate.f90 @@ -105,10 +105,7 @@ module subroutine collision_generate_hitandrun(self, nbody_system, param, t) ! Result integer(I4B) :: status !! Status flag assigned to this outcome ! Internals - integer(I4B) :: i, ibiggest, nfrag, jtarg, jproj - logical :: lpure character(len=STRMAX) :: message - real(DP) :: dpe select type(nbody_system) @@ -120,69 +117,19 @@ module subroutine collision_generate_hitandrun(self, nbody_system, param, t) call collision_io_collider_message(nbody_system%pl, impactors%id, message) call swiftest_io_log_one_message(COLLISION_LOG_OUT, trim(adjustl(message))) - if (impactors%mass(1) > impactors%mass(2)) then - jtarg = 1 - jproj = 2 - else - jtarg = 2 - jproj = 1 - end if - - ! The simple disruption model (and its extended types allow for non-pure hit and run. - !For the basic merge model, all hit and runs are pure - select type(self) - class is (collision_disruption) - if (impactors%mass_dist(2) > 0.9_DP * impactors%mass(jproj)) then ! Pure hit and run, so we'll just keep the two bodies untouched - call swiftest_io_log_one_message(COLLISION_LOG_OUT, "Pure hit and run. No new fragments generated.") - nfrag = 0 - lpure = .true. - else ! Imperfect hit and run, so we'll keep the largest body and destroy the other - lpure = .false. - call self%set_mass_dist(param) - - ! Generate the position and velocity distributions of the fragments - call self%get_energy_and_momentum(nbody_system, param, lbefore=.true.) - call self%disrupt(nbody_system, param, t, lpure) - call self%get_energy_and_momentum(nbody_system, param, lbefore=.false.) - - dpe = self%pe(2) - self%pe(1) - nbody_system%Ecollisions = nbody_system%Ecollisions - dpe - nbody_system%Euntracked = nbody_system%Euntracked + dpe - - if (lpure) then - call swiftest_io_log_one_message(COLLISION_LOG_OUT, "Should have been a pure hit and run instead") - nfrag = 0 - else - nfrag = self%fragments%nbody - write(message, *) nfrag - call swiftest_io_log_one_message(COLLISION_LOG_OUT, "Generating " // trim(adjustl(message)) // " fragments") - end if - end if - class default - lpure = .true. - end select + status = HIT_AND_RUN_PURE + pl%status(impactors%id(:)) = ACTIVE + pl%ldiscard(impactors%id(:)) = .false. + pl%lcollision(impactors%id(:)) = .false. + ! Be sure to save the pl so that snapshots still work + select type(before => self%before) + class is (swiftest_nbody_system) + select type(after => self%after) + class is (swiftest_nbody_system) + allocate(after%pl, source=before%pl) + end select + end select - if (lpure) then ! Reset these bodies back to being active so that nothing further is done to them - status = HIT_AND_RUN_PURE - pl%status(impactors%id(:)) = ACTIVE - pl%ldiscard(impactors%id(:)) = .false. - pl%lcollision(impactors%id(:)) = .false. - ! Be sure to save the pl so that snapshots still work - select type(before => self%before) - class is (swiftest_nbody_system) - select type(after => self%after) - class is (swiftest_nbody_system) - allocate(after%pl, source=before%pl) - end select - end select - else - ibiggest = impactors%id(maxloc(pl%Gmass(impactors%id(:)), dim=1)) - self%fragments%id(1) = pl%id(ibiggest) - self%fragments%id(2:nfrag) = [(i, i = param%maxid + 1, param%maxid + nfrag - 1)] - param%maxid = self%fragments%id(nfrag) - status = HIT_AND_RUN_DISRUPT - call collision_resolve_mergeaddsub(nbody_system, param, t, status) - end if end associate end select end select @@ -191,81 +138,6 @@ module subroutine collision_generate_hitandrun(self, nbody_system, param, t) end subroutine collision_generate_hitandrun - module subroutine collision_generate_disruption(self, nbody_system, param, t) - !! author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton - !! - !! Create the fragments resulting from a non-catastrophic disruption collision - !! - implicit none - ! Arguments - class(collision_disruption), intent(inout) :: self - class(base_nbody_system), intent(inout) :: nbody_system !! Swiftest nbody system object - 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 - character(len=STRMAX) :: message - real(DP) :: dpe - - select type(nbody_system) - class is (swiftest_nbody_system) - select type(pl => nbody_system%pl) - class is (swiftest_pl) - associate(impactors => self%impactors, status => self%status) - - select case (impactors%regime) - case (COLLRESOLVE_REGIME_HIT_AND_RUN) - call self%hitandrun(nbody_system, param, t) - return - case (COLLRESOLVE_REGIME_MERGE, COLLRESOLVE_REGIME_GRAZE_AND_MERGE) - call self%merge(nbody_system, param, t) ! Use the default collision model, which is merge - return - case(COLLRESOLVE_REGIME_DISRUPTION) - message = "Disruption between" - case(COLLRESOLVE_REGIME_SUPERCATASTROPHIC) - message = "Supercatastrophic disruption between" - case default - write(*,*) "Error in swiftest_collision, unrecognized collision regime" - call util_exit(FAILURE) - end select - call self%set_mass_dist(param) - call self%get_energy_and_momentum(nbody_system, param, lbefore=.true.) - call self%set_budgets() - call impactors%set_coordinate_system() - call self%disrupt(nbody_system, param, t) - call self%get_energy_and_momentum(nbody_system, param, lbefore=.false.) - - dpe = self%pe(2) - self%pe(1) - nbody_system%Ecollisions = nbody_system%Ecollisions - dpe - nbody_system%Euntracked = nbody_system%Euntracked + dpe - - associate (fragments => self%fragments) - ! Populate the list of new bodies - nfrag = fragments%nbody - write(message, *) nfrag - call swiftest_io_log_one_message(COLLISION_LOG_OUT, "Generating " // trim(adjustl(message)) // " fragments") - select case(impactors%regime) - case(COLLRESOLVE_REGIME_DISRUPTION) - status = DISRUPTED - ibiggest = impactors%id(maxloc(pl%Gmass(impactors%id(:)), dim=1)) - fragments%id(1) = pl%id(ibiggest) - fragments%id(2:nfrag) = [(i, i = param%maxid + 1, param%maxid + nfrag - 1)] - param%maxid = fragments%id(nfrag) - case(COLLRESOLVE_REGIME_SUPERCATASTROPHIC) - status = SUPERCATASTROPHIC - fragments%id(1:nfrag) = [(i, i = param%maxid + 1, param%maxid + nfrag)] - param%maxid = fragments%id(nfrag) - end select - - call collision_resolve_mergeaddsub(nbody_system, param, t, status) - end associate - end associate - end select - end select - return - end subroutine collision_generate_disruption - - module subroutine collision_generate_merge(self, nbody_system, param, t) !! author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton !! @@ -370,273 +242,4 @@ module subroutine collision_generate_merge(self, nbody_system, param, t) end subroutine collision_generate_merge - module subroutine collision_generate_disrupt(self, nbody_system, param, t, lfailure) - !! 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. - implicit none - ! Arguments - class(collision_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 - ! Internals - - call collision_generate_disruption_pos_vec(self) - call collision_generate_disruption_rot_vec(self) - call collision_generate_disruption_vel_vec(self) - - return - end subroutine collision_generate_disrupt - - - module subroutine collision_generate_disruption_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. - !! 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_disruption), intent(inout) :: collider !! Fraggle collision system object - ! Internals - 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) - lcat = (impactors%regime == COLLRESOLVE_REGIME_SUPERCATASTROPHIC) - lhitandrun = (impactors%regime == COLLRESOLVE_REGIME_HIT_AND_RUN) - - ! 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)) - - fragment_cloud_radius(:) = impactors%radius(:) - - loverlap(:) = .true. - 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)) - - ! 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 - - ! 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 - 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() - - 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_disruption_pos_vec - - - module subroutine collision_generate_disruption_rot_vec(collider) - !! Author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton - !! - !! Calculates the spins of a collection of fragments such that they conserve angular momentum - implicit none - ! Arguments - class(collision_basic), intent(inout) :: collider !! Fraggle collision system object - ! Internals - real(DP), dimension(NDIM) :: Lbefore, Lafter, Lspin, rotdir - real(DP) :: v_init, v_final, mass_init, mass_final, rotmag - integer(I4B) :: i - - associate(fragments => collider%fragments, impactors => collider%impactors, nfrag => collider%fragments%nbody) - - ! Torque the first body based on the change in angular momentum betwen the pre- and post-impact system assuming a simple bounce - mass_init = impactors%mass(2) - mass_final = sum(fragments%mass(2:nfrag)) - v_init = .mag.(impactors%vb(:,2) - impactors%vb(:,1)) - v_final = sqrt(mass_init / mass_final * v_init**2 - 2 * impactors%Qloss / mass_final) - - Lbefore(:) = mass_init * (impactors%rb(:,2) - impactors%rb(:,1)) .cross. (impactors%vb(:,2) - impactors%vb(:,1)) - - Lafter(:) = mass_final * (impactors%rb(:,2) - impactors%rb(:,1)) .cross. (v_final * impactors%bounce_unit(:)) - Lspin(:) = impactors%Lspin(:,1) + (Lbefore(:) - Lafter(:)) - fragments%rot(:,1) = Lspin(:) / (fragments%mass(1) * fragments%radius(1)**2 * fragments%Ip(3,1)) - - ! Add in some random spin noise. The magnitude will be scaled by the before-after amount and the direction will be random - do concurrent(i = 2:nfrag) - call random_number(rotdir) - rotdir = rotdir - 0.5_DP - rotdir = .unit. rotdir - call random_number(rotmag) - rotmag = rotmag * .mag. (Lbefore(:) - Lafter(:)) / ((nfrag - 1) * fragments%mass(i) * fragments%radius(i)**2 * fragments%Ip(3,i)) - fragments%rot(:,i) = rotmag * rotdir - end do - - end associate - - return - end subroutine collision_generate_disruption_rot_vec - - - module subroutine collision_generate_disruption_vel_vec(collider) - !! Author: David A. Minton - !! - !! Generates an initial velocity distribution. For disruptions, the velocity magnitude is set to be - !! 2x the escape velocity of the colliding pair. For hit and runs the velocity magnitude is set to be - !! 2x the escape velocity of the smallest of the two bodies. - implicit none - ! Arguments - class(collision_disruption), intent(inout) :: collider !! Fraggle collision system object - ! Internals - integer(I4B) :: i,j, loop, istart - logical :: lhitandrun, lnoncat - real(DP), dimension(NDIM) :: vimp_unit, rimp, vrot, Lresidual - real(DP) :: vmag, vesc, rotmag, ke_residual, ke_per_dof - integer(I4B), dimension(collider%fragments%nbody) :: vsign - real(DP), dimension(collider%fragments%nbody) :: vscale, mass_vscale, ke_avail - integer(I4B), parameter :: MAXLOOP = 100 - - associate(fragments => collider%fragments, impactors => collider%impactors, nfrag => collider%fragments%nbody) - 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 - - ! The fragments will be divided into two "clouds" based on identified origin body. - ! These clouds will collectively travel like two impactors bouncing off of each other. - where(fragments%origin_body(:) == 1) - vsign(:) = -1 - elsewhere - vsign(:) = 1 - end where - - ! The minimum fragment velocity will be set by the escape velocity - if (lhitandrun) then - vesc = sqrt(2 * impactors%Gmass(2) / impactors%radius(2)) - else - vesc = sqrt(2 * sum(impactors%Gmass(:)) / sum(impactors%radius(:))) - end if - - ! Scale the magnitude of the velocity by the distance from the impact point - ! This will reduce the chances of fragments colliding with each other immediately, and is more physically correct - do concurrent(i = 1:nfrag) - rimp(:) = fragments%rc(:,i) - impactors%rbimp(:) - vscale(i) = .mag. rimp(:) / (.mag. (impactors%rb(:,2) - impactors%rb(:,1))) - end do - vscale(:) = vscale(:)/minval(vscale(:)) - - ! Give the fragment velocities a random value that is scaled with fragment mass - call random_number(mass_vscale) - mass_vscale(:) = (mass_vscale(:) + 1.0_DP) / 2 - mass_vscale(:) = mass_vscale(:) * (fragments%mtot / fragments%mass(:))**(0.125_DP) ! The 1/8 power is arbitrary. It just gives the velocity a small mass dependence - mass_vscale(:) = mass_vscale(:) / minval(mass_vscale(:)) - - ! Set the velocities of all fragments using all of the scale factors determined above - do concurrent(i = 1:nfrag) - j = fragments%origin_body(i) - vrot(:) = impactors%rot(:,j) .cross. (fragments%rc(:,i) - impactors%rc(:,j)) - if (lhitandrun) then - if (i == 1) then - fragments%vc(:,1) = impactors%vc(:,1) - else - vmag = .mag.impactors%vc(:,2) / (maxval(mass_vscale(:) * maxval(vscale(:)))) - fragments%vc(:,i) = vmag * mass_vscale(i) * vscale(i) * impactors%bounce_unit(:) * vsign(i) + vrot(:) - end if - else - ! Add more velocity dispersion to disruptions vs hit and runs. - vmag = vesc * vscale(i) * mass_vscale(i) - rimp(:) = fragments%rc(:,i) - impactors%rbimp(:) - vimp_unit(:) = .unit. rimp(:) - fragments%vc(:,i) = vmag * (impactors%bounce_unit(:) + vimp_unit(:)) * vsign(i) + vrot(:) - end if - end do - - if (lhitandrun) then - istart = 2 - else - istart = 1 - end if - do loop = 1, MAXLOOP - call collider%set_coordinate_system() - call fragments%get_kinetic_energy() - ke_residual = fragments%ke_budget - (fragments%ke_orbit + fragments%ke_spin) - ! Make sure we don't take away too much orbital kinetic energy, otherwise the fragment can't escape - ke_avail(:) = fragments%ke_orbit_frag(:) - impactors%Gmass(1)*impactors%mass(2)/fragments%vmag(:) - ke_per_dof = -ke_residual/(nfrag - istart + 1) - do concurrent(i = istart:nfrag, ke_avail(i) > ke_per_dof) - fragments%vmag(i) = sqrt(2 * (fragments%ke_orbit_frag(i) - ke_per_dof)/fragments%mass(i)) - fragments%vc(:,i) = fragments%vmag(i) * fragments%v_unit(:,i) - end do - call fragments%get_kinetic_energy() - ke_residual = fragments%ke_budget - (fragments%ke_orbit + fragments%ke_spin) - - ! Check for any residual angular momentum, and if there is any, put it into spin of the biggest body - call collider%set_coordinate_system() - call fragments%get_angular_momentum() - Lresidual(:) = fragments%L_budget(:) - (fragments%Lorbit(:) + fragments%Lspin(:)) - rotmag = .mag. fragments%rot(:,1) - fragments%rot(:,1) = fragments%rot(:,1) + Lresidual(:) / (fragments%mass(1) * fragments%radius(1)**2 * fragments%Ip(:,1)) - rotmag = .mag. fragments%rot(:,1) - - if (ke_residual >= 0.0_DP) exit - - end do - - do concurrent(i = 1:nfrag) - fragments%vb(:,i) = fragments%vc(:,i) + impactors%vbcom(:) - end do - - impactors%vbcom(:) = 0.0_DP - do concurrent(i = 1:nfrag) - impactors%vbcom(:) = impactors%vbcom(:) + fragments%mass(i) * fragments%vb(:,i) - end do - impactors%vbcom(:) = impactors%vbcom(:) / fragments%mtot - - end associate - return - end subroutine collision_generate_disruption_vel_vec - - end submodule s_collision_generate diff --git a/src/collision/collision_module.f90 b/src/collision/collision_module.f90 index c9d81322f..5fe3893cb 100644 --- a/src/collision/collision_module.f90 +++ b/src/collision/collision_module.f90 @@ -123,7 +123,6 @@ module collision procedure :: reset => collision_util_reset_fragments !! Deallocates all allocatable arrays and sets everything else to 0 procedure :: get_angular_momentum => collision_util_get_angular_momentum !! Calcualtes the current angular momentum of the fragments procedure :: get_kinetic_energy => collision_util_get_kinetic_energy !! Calcualtes the current kinetic energy of the fragments - procedure :: set_spins => collision_util_set_spins !! Calcualtes the spins of all fragments from the angular momentum budget and residual final :: collision_final_fragments !! Finalizer deallocates all allocatables end type collision_fragments @@ -168,13 +167,7 @@ module collision final :: collision_final_bounce !! Finalizer will deallocate all allocatables end type collision_bounce - type, extends(collision_basic) :: collision_disruption - contains - procedure :: generate => collision_generate_disruption !! A simple disruption models that does not constrain energy loss in collisions - procedure :: disrupt => collision_generate_disrupt !! Disrupt the colliders into the fragments - procedure :: set_mass_dist => collision_util_set_mass_dist !! Sets the distribution of mass among the fragments depending on the regime type - final :: collision_final_simple !! Finalizer will deallocate all allocatables - end type collision_disruption + !! NetCDF dimension and variable names for the enounter save object @@ -236,14 +229,6 @@ module subroutine collision_generate_bounce(self, nbody_system, param, t) real(DP), intent(in) :: t !! The time of the collision end subroutine collision_generate_bounce - module subroutine collision_generate_disrupt(self, nbody_system, param, t, lfailure) - implicit none - class(collision_disruption), intent(inout) :: self - class(base_nbody_system), intent(inout) :: nbody_system !! Swiftest nbody system object - class(base_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA additions - real(DP), intent(in) :: t !! Time of collision - logical, optional, intent(out) :: lfailure !! Disruption failed - end subroutine collision_generate_disrupt module subroutine collision_generate_hitandrun(self, nbody_system, param, t) implicit none @@ -261,28 +246,6 @@ module subroutine collision_generate_merge(self, nbody_system, param, t) real(DP), intent(in) :: t !! The time of the collision end subroutine collision_generate_merge - module subroutine collision_generate_disruption(self, nbody_system, param, t) - implicit none - class(collision_disruption), intent(inout) :: self !! Simple fragment nbody_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 - end subroutine collision_generate_disruption - - module subroutine collision_generate_disruption_pos_vec(collider) - implicit none - class(collision_disruption), intent(inout) :: collider !! Collision system object - end subroutine collision_generate_disruption_pos_vec - - module subroutine collision_generate_disruption_rot_vec(collider) - implicit none - class(collision_basic), intent(inout) :: collider !! Collision system object - end subroutine collision_generate_disruption_rot_vec - - module subroutine collision_generate_disruption_vel_vec(collider) - implicit none - class(collision_disruption), intent(inout) :: collider !! Collision system object - end subroutine collision_generate_disruption_vel_vec module subroutine collision_io_collider_message(pl, collidx, collider_message) implicit none @@ -446,17 +409,6 @@ module subroutine collision_util_set_coordinate_impactors(self) class(collision_impactors), intent(inout) :: self !! collisional system end subroutine collision_util_set_coordinate_impactors - module subroutine collision_util_set_mass_dist(self, param) - implicit none - class(collision_disruption), intent(inout) :: self !! Simple disruption collision object - class(base_parameters), intent(in) :: param !! Current Swiftest run configuration parameters - end subroutine collision_util_set_mass_dist - - module subroutine collision_util_set_spins(self) - implicit none - class(collision_fragments(*)), intent(inout) :: self !! Collision fragment system object - end subroutine collision_util_set_spins - module subroutine collision_util_setup_collider(self, nbody_system) implicit none class(collision_basic), intent(inout) :: self !! Encounter collision system object @@ -534,7 +486,6 @@ subroutine collision_final_fragments(self) return end subroutine collision_final_fragments - subroutine collision_final_impactors(self) !! author: David A. Minton !! @@ -548,7 +499,6 @@ subroutine collision_final_impactors(self) return end subroutine collision_final_impactors - subroutine collision_final_netcdf_parameters(self) !! author: David A. Minton !! @@ -562,7 +512,6 @@ subroutine collision_final_netcdf_parameters(self) return end subroutine collision_final_netcdf_parameters - subroutine collision_final_plpl(self) !! author: David A. Minton !! @@ -589,7 +538,6 @@ subroutine collision_final_pltp(self) return end subroutine collision_final_pltp - subroutine collision_final_snapshot(self) !! author: David A. Minton !! @@ -603,7 +551,6 @@ subroutine collision_final_snapshot(self) return end subroutine collision_final_snapshot - subroutine collision_final_storage(self) !! author: David A. Minton !! @@ -621,7 +568,6 @@ subroutine collision_final_storage(self) return end subroutine collision_final_storage - subroutine collision_final_bounce(self) !! author: David A. Minton !! @@ -637,21 +583,5 @@ subroutine collision_final_bounce(self) return end subroutine collision_final_bounce - - subroutine collision_final_simple(self) - !! author: David A. Minton - !! - !! Finalizer will deallocate all allocatables - implicit none - ! Arguments - type(collision_disruption), intent(inout) :: self !! Collision system object - - call self%reset() - if (allocated(self%impactors)) deallocate(self%impactors) - if (allocated(self%fragments)) deallocate(self%fragments) - - return - end subroutine collision_final_simple - end module collision diff --git a/src/collision/collision_util.f90 b/src/collision/collision_util.f90 index 0d76af301..ef582ee10 100644 --- a/src/collision/collision_util.f90 +++ b/src/collision/collision_util.f90 @@ -527,194 +527,6 @@ module subroutine collision_util_set_coordinate_impactors(self) end subroutine collision_util_set_coordinate_impactors - module subroutine collision_util_set_mass_dist(self, param) - !! author: David A. Minton - !! - !! Sets the mass of fragments based on the mass distribution returned by the regime calculation. - !! This subroutine must be run after the the setup routine has been run on the fragments - !! - implicit none - ! Arguments - class(collision_disruption), intent(inout) :: self !! Fraggle collision system object - class(base_parameters), intent(in) :: param !! Current Swiftest run configuration parameters - ! Internals - 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, 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) - integer(I4B), parameter :: NFRAG_SIZE_MULTIPLIER = 3 !! Log-space scale factor that scales the number of fragments by the collisional system mass - 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 - volume(1:2) = 4._DP / 3._DP * PI * impactors%radius(1:2)**3 - mtot = sum(impactors%mass(:)) - Ip_avg(:) = (impactors%mass(1) * impactors%Ip(:,1) + impactors%mass(2) * impactors%Ip(:,2)) / mtot - - if (impactors%mass(1) > impactors%mass(2)) then - jtarg = 1 - jproj = 2 - else - jtarg = 2 - jproj = 1 - end if - - select case(impactors%regime) - case(COLLRESOLVE_REGIME_DISRUPTION, COLLRESOLVE_REGIME_SUPERCATASTROPHIC, COLLRESOLVE_REGIME_HIT_AND_RUN) - ! The first two bins of the mass_dist are the largest and second-largest fragments that came out of collision_regime. - ! The remainder from the third bin will be distributed among nfrag-2 bodies. The following code will determine nfrag based on - ! the limits bracketed above and the model size distribution of fragments. - ! Check to see if our size distribution would give us a smaller number of fragments than the maximum number - - select type(param) - class is (swiftest_parameters) - min_mfrag = (param%min_GMfrag / param%GU) - ! The number of fragments we generate is bracked by the minimum required by fraggle_generate (7) and the - ! maximum set by the NFRAG_SIZE_MULTIPLIER which limits the total number of fragments to prevent the nbody - ! code from getting an overwhelmingly large number of fragments - nfrag = ceiling(NFRAG_SIZE_MULTIPLIER * log(mtot / min_mfrag)) - nfrag = max(min(nfrag, NFRAGMAX), NFRAGMIN) - class default - min_mfrag = 0.0_DP - nfrag = NFRAGMAX - end select - - i = iMrem - mremaining = impactors%mass_dist(iMrem) - do while (i <= nfrag) - mfrag = (1 + i - iMslr)**(-3._DP / BETA) * impactors%mass_dist(iMslr) - if (mremaining - mfrag < 0.0_DP) exit - mremaining = mremaining - mfrag - i = i + 1 - end do - if (i < nfrag) nfrag = max(i, NFRAGMIN) ! The sfd would actually give us fewer fragments than our maximum - call self%setup_fragments(nfrag) - - case (COLLRESOLVE_REGIME_MERGE, COLLRESOLVE_REGIME_GRAZE_AND_MERGE) - - call self%setup_fragments(1) - select type(fragments => self%fragments) - class is (collision_fragments(*)) - fragments%mass(1) = impactors%mass_dist(1) - fragments%radius(1) = impactors%radius(jtarg) - fragments%density(1) = impactors%mass_dist(1) / volume(jtarg) - if (param%lrotation) fragments%Ip(:, 1) = impactors%Ip(:,1) - end select - return - case default - write(*,*) "collision_util_set_mass_dist_fragments error: Unrecognized regime code",impactors%regime - end select - - select type(fragments => self%fragments) - class is (collision_fragments(*)) - fragments%mtot = mtot - - ! Make the first two bins the same as the Mlr and Mslr values that came from collision_regime - fragments%mass(1) = impactors%mass_dist(iMlr) - fragments%mass(2) = impactors%mass_dist(iMslr) - - ! Distribute the remaining mass the 3:nfrag bodies following the model SFD given by slope BETA - mremaining = impactors%mass_dist(iMrem) - do i = iMrem, nfrag - mfrag = (1 + i - iMslr)**(-3._DP / BETA) * impactors%mass_dist(iMslr) - fragments%mass(i) = mfrag - mremaining = mremaining - mfrag - end do - - ! If there is any residual mass (either positive or negative) we will distribute remaining mass proportionally among the the fragments - if (mremaining < 0.0_DP) then ! If the remainder is negative, this means that that the number of fragments required by the SFD is smaller than our lower limit set by fraggle_generate. - istart = iMrem ! We will reduce the mass of the 3:nfrag bodies to prevent the second-largest fragment from going smaller - else ! If the remainder is postiive, this means that the number of fragments required by the SFD is larger than our upper limit set by computational expediency. - istart = iMslr ! We will increase the mass of the 2:nfrag bodies to compensate, which ensures that the second largest fragment remains the second largest - end if - mfrag = 1._DP + mremaining / sum(fragments%mass(istart:nfrag)) - fragments%mass(istart:nfrag) = fragments%mass(istart:nfrag) * mfrag - - ! There may still be some small residual due to round-off error. If so, simply add it to the last bin of the mass distribution. - mremaining = fragments%mtot - sum(fragments%mass(1:nfrag)) - fragments%mass(nfrag) = fragments%mass(nfrag) + mremaining - - ! Compute physical properties of the new fragments - select case(impactors%regime) - case(COLLRESOLVE_REGIME_HIT_AND_RUN) ! The hit and run case always preserves the largest body intact, so there is no need to recompute the physical properties of the first fragment - fragments%radius(1) = impactors%radius(jtarg) - fragments%density(1) = impactors%mass_dist(iMlr) / volume(jtarg) - fragments%Ip(:, 1) = impactors%Ip(:,1) - istart = 2 - case default - istart = 1 - end select - - fragments%density(istart:nfrag) = fragments%mtot / sum(volume(:)) - fragments%radius(istart:nfrag) = (3 * fragments%mass(istart:nfrag) / (4 * PI * fragments%density(istart:nfrag)))**(1.0_DP / 3.0_DP) - do i = istart, nfrag - 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 - end subroutine collision_util_set_mass_dist - - - module subroutine collision_util_set_spins(self) - !! Author: David A. Minton - !! - !! Distributes any residual angular momentum into the spins of the n>1 fragments - implicit none - ! Arguments - class(collision_fragments(*)), intent(inout) :: self !! Collision fragment system object - ! Internals - integer(I4B) :: i - real(DP), dimension(NDIM) :: Lresidual - - call self%get_angular_momentum() - Lresidual(:) = self%L_budget(:) - (self%Lorbit(:) + self%Lspin(:)) - ! Distribute residual angular momentum amongst the fragments - if (.mag.(Lresidual(:)) > tiny(1.0_DP)) then - do i = 1,self%nbody - self%rot(:,i) = self%rot(:,i) + Lresidual(:) / (self%nbody * self%mass(i) * self%radius(i)**2 * self%Ip(3,i)) - end do - end if - - call self%get_angular_momentum() - Lresidual(:) = self%L_budget(:) - (self%Lorbit(:) + self%Lspin(:)) - - return - end subroutine collision_util_set_spins - - module subroutine collision_util_setup_collider(self, nbody_system) !! author: David A. Minton !! diff --git a/src/fraggle/fraggle_generate.f90 b/src/fraggle/fraggle_generate.f90 index f8361be19..08c184ab0 100644 --- a/src/fraggle/fraggle_generate.f90 +++ b/src/fraggle/fraggle_generate.f90 @@ -11,12 +11,86 @@ use swiftest use symba - integer(I4B), parameter :: NFRAG_MIN = 7 !! The minimum allowable number of fragments (set to 6 because that's how many unknowns are needed in the tangential velocity calculation) real(DP), parameter :: FRAGGLE_LTOL = 1e-4_DP !10 * epsilon(1.0_DP) real(DP), parameter :: FRAGGLE_ETOL = 1e-1_DP contains + module subroutine fraggle_generate(self, nbody_system, param, t) + !! author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton + !! + !! Create the fragments resulting from a non-catastrophic disruption collision + !! + implicit none + ! Arguments + class(collision_fraggle), intent(inout) :: self + class(base_nbody_system), intent(inout) :: nbody_system !! Swiftest nbody system object + 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 + character(len=STRMAX) :: message + real(DP) :: dpe + + select type(nbody_system) + class is (swiftest_nbody_system) + select type(pl => nbody_system%pl) + class is (swiftest_pl) + associate(impactors => self%impactors, status => self%status) + + select case (impactors%regime) + case (COLLRESOLVE_REGIME_HIT_AND_RUN) + call self%hitandrun(nbody_system, param, t) + return + case (COLLRESOLVE_REGIME_MERGE, COLLRESOLVE_REGIME_GRAZE_AND_MERGE) + call self%merge(nbody_system, param, t) ! Use the default collision model, which is merge + return + case(COLLRESOLVE_REGIME_DISRUPTION) + message = "Disruption between" + case(COLLRESOLVE_REGIME_SUPERCATASTROPHIC) + message = "Supercatastrophic disruption between" + case default + write(*,*) "Error in swiftest_collision, unrecognized collision regime" + call util_exit(FAILURE) + end select + call self%set_mass_dist(param) + call self%get_energy_and_momentum(nbody_system, param, lbefore=.true.) + call self%set_budgets() + call impactors%set_coordinate_system() + call self%disrupt(nbody_system, param, t) + call self%get_energy_and_momentum(nbody_system, param, lbefore=.false.) + + dpe = self%pe(2) - self%pe(1) + nbody_system%Ecollisions = nbody_system%Ecollisions - dpe + nbody_system%Euntracked = nbody_system%Euntracked + dpe + + associate (fragments => self%fragments) + ! Populate the list of new bodies + nfrag = fragments%nbody + write(message, *) nfrag + call swiftest_io_log_one_message(COLLISION_LOG_OUT, "Generating " // trim(adjustl(message)) // " fragments") + select case(impactors%regime) + case(COLLRESOLVE_REGIME_DISRUPTION) + status = DISRUPTED + ibiggest = impactors%id(maxloc(pl%Gmass(impactors%id(:)), dim=1)) + fragments%id(1) = pl%id(ibiggest) + fragments%id(2:nfrag) = [(i, i = param%maxid + 1, param%maxid + nfrag - 1)] + param%maxid = fragments%id(nfrag) + case(COLLRESOLVE_REGIME_SUPERCATASTROPHIC) + status = SUPERCATASTROPHIC + fragments%id(1:nfrag) = [(i, i = param%maxid + 1, param%maxid + nfrag)] + param%maxid = fragments%id(nfrag) + end select + + call collision_resolve_mergeaddsub(nbody_system, param, t, status) + end associate + end associate + end select + end select + return + end subroutine fraggle_generate + + module subroutine fraggle_generate_disrupt(self, nbody_system, param, t, lfailure) !! Author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton !! @@ -54,12 +128,6 @@ module subroutine fraggle_generate_disrupt(self, nbody_system, param, t, lfailur write(message,*) nfrag call swiftest_io_log_one_message(COLLISION_LOG_OUT, "Fraggle generating " // trim(adjustl(message)) // " fragments.") - if (nfrag < NFRAG_MIN) then - write(message,*) "Fraggle needs at least ",NFRAG_MIN," fragments, but only ",nfrag," were given." - call swiftest_io_log_one_message(COLLISION_LOG_OUT, message) - lfailure_local = .true. - return - end if if (param%lflatten_interactions) then lk_plpl = allocated(pl%k_plpl) @@ -88,14 +156,10 @@ module subroutine fraggle_generate_disrupt(self, nbody_system, param, t, lfailur ! Use the disruption collision model to generate initial conditions ! Compute the "before" energy/momentum and the budgets call self%get_energy_and_momentum(nbody_system, param, lbefore=.true.) - call self%collision_disruption%disrupt(nbody_system, param, t) + call fraggle_generate_pos_vec(self) + call fraggle_generate_rot_vec(self) + call fraggle_generate_vel_vec(self) call self%get_energy_and_momentum(nbody_system, param, lbefore=.false.) - call self%set_budgets() - - ! Minimize difference between energy/momentum and budgets - ! call fraggle_generate_minimize(self, lfailure_local) - ! call fraggle_generate_tan_vel(self, lfailure_local) - call fraggle_generate_rad_vel(self, lfailure_local) call self%get_energy_and_momentum(nbody_system, param, lbefore=.false.) dEtot = self%Etot(2) - self%Etot(1) @@ -151,433 +215,334 @@ module subroutine fraggle_generate_disrupt(self, nbody_system, param, t, lfailur end subroutine fraggle_generate_disrupt - subroutine fraggle_generate_minimize(collider, lfailure) - !! Author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton - !! - !! Minimizes the differences between the energy and momentum and the budget + module subroutine fraggle_generate_hitandrun(self, nbody_system, param, t) + !! author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton !! - !! A failure will trigger a restructuring of the fragments so we will try a new random distribution + !! Create the fragments resulting from a non-catastrophic hit-and-run collision + !! implicit none ! Arguments - class(collision_fraggle), intent(inout) :: collider !! Fraggle collision system object - logical, intent(out) :: lfailure !! Logical flag indicating whether this step fails or succeeds + class(collision_fraggle), intent(inout) :: self !! Collision system object + class(base_nbody_system), intent(inout) :: nbody_system !! Swiftest nbody system object + class(base_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA additions + real(DP), intent(in) :: t !! Time of collision + ! Result + integer(I4B) :: status !! Status flag assigned to this outcome ! Internals - real(DP), parameter :: TOL_INIT = 1e-5_DP - integer(I4B), parameter :: MAXLOOP = 100 - real(DP), dimension(collider%fragments%nbody) :: input_v - real(DP), dimension(:), allocatable :: output_v - type(lambda_obj) :: Efunc - real(DP) :: tol, fval - integer(I4B) :: loop,i, nelem - logical :: lfirst_Efunc - - associate(impactors => collider%impactors, nfrag => collider%fragments%nbody) - select type(fragments => collider%fragments) - class is (fraggle_fragments(*)) - - nelem = nfrag - lfailure = .false. - ! Find the local kinetic energy minimum for the nbody_system that conserves linear and angular momentum - Efunc = lambda_obj(E_objective_function) - tol = TOL_INIT - - fragments%r_unit(:,:) = .unit. fragments%vc(:,:) - fragments%vmag(:) = .mag. fragments%vc(:,1:nfrag) - input_v(:) = fragments%vmag(1:nfrag) - lfirst_Efunc = .true. - fval = E_objective_function(input_v) - - call minimize_bfgs(Efunc, nelem, input_v, tol, MAXLOOP, lfailure, output_v) - fval = E_objective_function(output_v) - fragments%vmag(1:nfrag) = output_v(1:nfrag) - do concurrent(i=1:nfrag) - fragments%vc(:,i) = abs(fragments%vmag(i)) * fragments%r_unit(:,i) - end do + integer(I4B) :: i, ibiggest, nfrag, jtarg, jproj + logical :: lpure + character(len=STRMAX) :: message + real(DP) :: dpe - ! Set spins in order to conserve angular momentum - call fragments%set_spins() - call collision_util_shift_vector_to_origin(fragments%mass, fragments%vc) - end select - end associate - return + select type(nbody_system) + class is (swiftest_nbody_system) + select type(pl => nbody_system%pl) + class is (swiftest_pl) + associate(impactors => self%impactors) + message = "Hit and run between" + call collision_io_collider_message(nbody_system%pl, impactors%id, message) + call swiftest_io_log_one_message(COLLISION_LOG_OUT, trim(adjustl(message))) + + if (impactors%mass(1) > impactors%mass(2)) then + jtarg = 1 + jproj = 2 + else + jtarg = 2 + jproj = 1 + end if - contains - - function E_objective_function(val_input) result(fval) - !! Author: David A. Minton - !! - !! Objective function for evaluating how close our fragment trajectories are to the energy budget - implicit none - ! Arguments - real(DP), dimension(:), intent(in) :: val_input !! Flattened velocity and rotation arrays - ! Result - real(DP) :: fval !! The objective function result, which is the sum of the squares of the difference between the calculated fragment kinetic energy and the components of angular and linear momentum - !! Minimizing this brings us closer to our objective - ! Internals - integer(I4B) :: i - type(fraggle_fragments(:)), allocatable :: tmp_frag - real(DP) :: deltaE - real(DP), save :: deltaE_scale = 1.0_DP - - associate(impactors => collider%impactors, nfrag => collider%fragments%nbody) - select type(fragments => collider%fragments) - class is (fraggle_fragments(*)) - allocate(tmp_frag, source=fragments) - tmp_frag%vmag(1:nfrag) = val_input(1:nfrag) - do concurrent(i=1:nfrag) - tmp_frag%vc(:,i) = abs(tmp_frag%vmag(i)) * tmp_frag%r_unit(:,i) - end do - - call collision_util_shift_vector_to_origin(tmp_frag%mass, tmp_frag%vc) - ! Set spins in order to conserve angular momentum - call collision_util_set_spins(tmp_frag) - - ! Get the current kinetic energy of the system - call tmp_frag%get_kinetic_energy() - deltaE = (tmp_frag%ke_budget - (tmp_frag%ke_orbit + tmp_frag%ke_spin)) / (tmp_frag%ke_budget) - - ! The result will be scaled so such that the initial deltaE is 1.0. This helps improve the success of the - ! minimizer, by keeping values from getting to large - if (lfirst_Efunc) then - deltaE_scale = deltaE - deltaE = 1.0_DP - lfirst_Efunc = .false. - else - deltaE = deltaE/deltaE_scale - end if - fval = deltaE**2 - deallocate(tmp_frag) - - end select - end associate - - return - end function E_objective_function + ! The frag disruption model (and its extended types allow for non-pure hit and run. + if (impactors%mass_dist(2) > 0.9_DP * impactors%mass(jproj)) then ! Pure hit and run, so we'll just keep the two bodies untouched + call swiftest_io_log_one_message(COLLISION_LOG_OUT, "Pure hit and run. No new fragments generated.") + nfrag = 0 + call self%collision_basic%hitandrun(nbody_system, param, t) + lpure = .true. + return + else ! Imperfect hit and run, so we'll keep the largest body and destroy the other + lpure = .false. + call self%set_mass_dist(param) + + ! Generate the position and velocity distributions of the fragments + call self%get_energy_and_momentum(nbody_system, param, lbefore=.true.) + call self%disrupt(nbody_system, param, t, lpure) + call self%get_energy_and_momentum(nbody_system, param, lbefore=.false.) + + dpe = self%pe(2) - self%pe(1) + nbody_system%Ecollisions = nbody_system%Ecollisions - dpe + nbody_system%Euntracked = nbody_system%Euntracked + dpe + + if (lpure) then + call swiftest_io_log_one_message(COLLISION_LOG_OUT, "Should have been a pure hit and run instead") + nfrag = 0 + else + nfrag = self%fragments%nbody + write(message, *) nfrag + call swiftest_io_log_one_message(COLLISION_LOG_OUT, "Generating " // trim(adjustl(message)) // " fragments") + end if + end if + + ibiggest = impactors%id(maxloc(pl%Gmass(impactors%id(:)), dim=1)) + self%fragments%id(1) = pl%id(ibiggest) + self%fragments%id(2:nfrag) = [(i, i = param%maxid + 1, param%maxid + nfrag - 1)] + param%maxid = self%fragments%id(nfrag) + status = HIT_AND_RUN_DISRUPT + call collision_resolve_mergeaddsub(nbody_system, param, t, status) + end associate + end select + end select + + return + end subroutine fraggle_generate_hitandrun - end subroutine fraggle_generate_minimize - subroutine fraggle_generate_tan_vel(collider, lfailure) + module subroutine fraggle_generate_pos_vec(collider) !! Author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton !! - !! Adjusts the tangential velocities and spins of a collection of fragments such that they conserve angular momentum without blowing the fragment kinetic energy budget. - !! This procedure works in several stages, with a goal to solve the angular and linear momentum constraints on the fragments, while still leaving a positive balance of - !! our fragment kinetic energy (fragments%ke_budget) that we can put into the radial velocity distribution. - !! - !! The first thing we'll try to do is solve for the tangential velocities of the first 6 fragments, using angular and linear momentum as constraints and an initial - !! tangential velocity distribution for the remaining bodies (if there are any) that distributes their angular momentum equally between them. - !! If that doesn't work and we blow our kinetic energy budget, we will attempt to find a tangential velocity distribution that minimizes the kinetic energy while - !! conserving momentum. - !! - !! A failure will trigger a restructuring of the fragments so we will try new values of the radial position distribution. + !! 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_fraggle), intent(inout) :: collider !! Fraggle fragment system object - logical, intent(out) :: lfailure !! Logical flag indicating whether this step fails or succeeds + class(collision_fraggle), intent(inout) :: collider !! Fraggle collision system object ! Internals - integer(I4B) :: i - real(DP), parameter :: TOL_MIN = 1e-1_DP ! This doesn't have to be very accurate, as we really just want a tangential velocity distribution with less kinetic energy than our initial guess. - real(DP), parameter :: TOL_INIT = 1e-6_DP - real(DP), parameter :: VNOISE_MAG = 1e-3_DP !! Magnitude of the noise to apply to initial conditions to help minimizer find a solution in case of failure - integer(I4B), parameter :: MAXLOOP = 10 - real(DP) :: tol, fval - real(DP), dimension(:), allocatable :: v_t_initial - real(DP), dimension(collider%fragments%nbody) :: kefrag, vnoise - type(lambda_obj_err) :: objective_function - real(DP), dimension(NDIM) :: Li, L_remainder, L_frag_tot - character(len=STRMAX) :: message - real(DP), dimension(:), allocatable :: v_t_output - logical :: lfirst_func - - associate(impactors => collider%impactors, nfrag => collider%fragments%nbody) - select type(fragments => collider%fragments) - class is (fraggle_fragments(*)) - lfailure = .false. - - ! Solve for the tangential component of the velocities that both conserves linear momentum and uses the remaining angular momentum not used in spin. - ! This will be done using a linear solver that solves for the tangential velocities of the first 6 fragments, constrained by the linear and angular momentum vectors, - ! which is embedded in a non-linear minimizer that will adjust the tangential velocities of the remaining i>6 fragments to minimize kinetic energy for a given momentum solution - ! The initial conditions fed to the minimizer for the fragments will be the remaining angular momentum distributed between the fragments. - tol = TOL_INIT - lfirst_func = .true. - do i = 1, nfrag - fragments%v_t_mag(i) = dot_product(fragments%vc(:,i), fragments%t_unit(:,i)) - fragments%v_r_mag(i) = dot_product(fragments%vc(:,i), fragments%r_unit(:,i)) + 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) + lcat = (impactors%regime == COLLRESOLVE_REGIME_SUPERCATASTROPHIC) + lhitandrun = (impactors%regime == COLLRESOLVE_REGIME_HIT_AND_RUN) + + ! 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)) + + fragment_cloud_radius(:) = impactors%radius(:) + + loverlap(:) = .true. + 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)) + + ! 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 - allocate(v_t_initial, source=fragments%v_t_mag) - do while(tol < TOL_MIN) - - ! ! Find the local kinetic energy minimum for the system that conserves linear and angular momentum - objective_function = lambda_obj(tangential_objective_function, lfailure) - fval = tangential_objective_function(v_t_output(:), lfailure) - - call minimize_bfgs(objective_function, nfrag-6, v_t_initial(7:nfrag), tol, MAXLOOP, lfailure, v_t_output) - fval = tangential_objective_function(v_t_initial(:), lfailure) - fragments%v_t_mag(7:nfrag) = v_t_output(:) - ! Now that the KE-minimized values of the i>6 fragments are found, calculate the momentum-conserving solution for tangential velociteis - v_t_initial(7:nfrag) = fragments%v_t_mag(7:nfrag) - - fragments%v_t_mag(1:nfrag) = solve_fragment_tan_vel(v_t_mag_input=v_t_initial(7:nfrag), lfailure=lfailure) - - ! Perform one final shift of the radial velocity vectors to align with the center of mass of the collisional system (the origin) - fragments%vb(:,1:nfrag) = fraggle_util_vmag_to_vb(fragments%v_r_mag(1:nfrag), fragments%r_unit(:,1:nfrag), fragments%v_t_mag(1:nfrag), & - fragments%t_unit(:,1:nfrag), fragments%mass(1:nfrag), impactors%vbcom(:)) - do concurrent (i = 1:nfrag) - fragments%vc(:,i) = fragments%vb(:,i) - impactors%vbcom(:) - end do - ! Now do a kinetic energy budget check to make sure we are still within the budget. - kefrag = 0.0_DP - do concurrent(i = 1:nfrag) - kefrag(i) = fragments%mass(i) * dot_product(fragments%vc(:,i), fragments%vc(:,i)) + ! 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 - fragments%ke_orbit = 0.5_DP * sum(kefrag(:)) - - ! If we are over the energy budget, flag this as a failure so we can try again - call fragments%get_angular_momentum() - lfailure = ((fragments%ke_budget - fragments%ke_spin - fragments%ke_orbit) < 0.0_DP) - if (.not.lfailure) exit - tol = tol * 2_DP ! Keep increasing the tolerance until we converge on a solution - ! Reduce fragment spins to try to get a better solution - fragments%rot(:,2:nfrag) = fragments%rot(:,2:nfrag) * 0.9_DP end do - if (lfailure) then - call swiftest_io_log_one_message(COLLISION_LOG_OUT, " ") - call swiftest_io_log_one_message(COLLISION_LOG_OUT, "Tangential velocity failure diagnostics") - call fragments%get_angular_momentum() - L_frag_tot = fragments%Lspin(:) + fragments%Lorbit(:) - write(message, *) .mag.(L_frag_tot(:)) / (.mag.fragments%L_budget(:)) - call swiftest_io_log_one_message(COLLISION_LOG_OUT, "|L_remainder| : " // trim(adjustl(message))) - write(message, *) fragments%ke_budget - call swiftest_io_log_one_message(COLLISION_LOG_OUT, "ke_budget : " // trim(adjustl(message))) - write(message, *) fragments%ke_spin - call swiftest_io_log_one_message(COLLISION_LOG_OUT, "ke_spin : " // trim(adjustl(message))) - write(message, *) fragments%ke_orbit - call swiftest_io_log_one_message(COLLISION_LOG_OUT, "ke_tangential : " // trim(adjustl(message))) - write(message, *) fragments%ke_budget - fragments%ke_spin - fragments%ke_orbit - call swiftest_io_log_one_message(COLLISION_LOG_OUT, "ke_radial : " // trim(adjustl(message))) - end if - end select + 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() + + 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 - contains - function solve_fragment_tan_vel(lfailure, v_t_mag_input) result(v_t_mag_output) - !! Author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton - !! - !! Adjusts the positions, velocities, and spins of a collection of fragments such that they conserve angular momentum - implicit none - ! Arguments - logical, intent(out) :: lfailure !! Error flag - real(DP), dimension(:), optional, intent(in) :: v_t_mag_input !! Unknown tangential velocities for fragments 7:nfrag - ! Internals - integer(I4B) :: i - ! Result - real(DP), dimension(:), allocatable :: v_t_mag_output - - real(DP), dimension(2 * NDIM, 2 * NDIM) :: A ! LHS of linear equation used to solve for momentum constraint in Gauss elimination code - real(DP), dimension(2 * NDIM) :: b ! RHS of linear equation used to solve for momentum constraint in Gauss elimination code - real(DP), dimension(NDIM) :: L_lin_others, L_orb_others, L, vtmp - - associate(impactors => collider%impactors, nfrag => collider%fragments%nbody) - select type(fragments => collider%fragments) - class is (fraggle_fragments(*)) - lfailure = .false. - ! We have 6 constraint equations (2 vector constraints in 3 dimensions each) - ! The first 3 are that the linear momentum of the fragments is zero with respect to the collisional barycenter - ! The second 3 are that the sum of the angular momentum of the fragments is conserved from the pre-impact state - L_lin_others(:) = 0.0_DP - L_orb_others(:) = 0.0_DP - do i = 1, nfrag - if (i <= 2 * NDIM) then ! The tangential velocities of the first set of bodies will be the unknowns we will solve for to satisfy the constraints - A(1:3, i) = fragments%mass(i) * fragments%t_unit(:, i) - A(4:6, i) = fragments%mass(i) * fragments%rmag(i) * (fragments%r_unit(:, i) .cross. fragments%t_unit(:, i)) - else if (present(v_t_mag_input)) then - vtmp(:) = v_t_mag_input(i - 6) * fragments%t_unit(:, i) - L_lin_others(:) = L_lin_others(:) + fragments%mass(i) * vtmp(:) - L(:) = fragments%mass(i) * (fragments%rc(:, i) .cross. vtmp(:)) - L_orb_others(:) = L_orb_others(:) + L(:) - end if - end do - b(1:3) = -L_lin_others(:) - b(4:6) = fragments%L_budget(:) - fragments%Lspin(:) - L_orb_others(:) - allocate(v_t_mag_output(nfrag)) - v_t_mag_output(1:6) = solve_linear_system(A, b, 6, lfailure) - if (present(v_t_mag_input)) v_t_mag_output(7:nfrag) = v_t_mag_input(:) - end select - end associate - return - end function solve_fragment_tan_vel - - - function tangential_objective_function(v_t_mag_input, lfailure) result(fval) - !! Author: David A. Minton - !! - !! Objective function for evaluating how close our fragment velocities get to minimizing KE error from our required value - implicit none - ! Arguments - real(DP), dimension(:), intent(in) :: v_t_mag_input !! Unknown tangential component of velocity vector set previously by angular momentum constraint - logical, intent(out) :: lfailure !! Error flag - ! Result - real(DP) :: fval - ! Internals - integer(I4B) :: i - real(DP), dimension(NDIM,collider%fragments%nbody) :: vc, vb - real(DP), dimension(collider%fragments%nbody) :: v_t_new, kearr - real(DP) :: keo - real(DP), save :: fval_scale = 1.0_DP - - associate(impactors => collider%impactors, nfrag => collider%fragments%nbody) - select type(fragments => collider%fragments) - class is (fraggle_fragments(*)) - lfailure = .false. - - v_t_new(:) = solve_fragment_tan_vel(v_t_mag_input=v_t_mag_input(:), lfailure=lfailure) - vb(:,1:nfrag) = fraggle_util_vmag_to_vb(fragments%v_r_mag(1:nfrag), fragments%r_unit(:,1:nfrag), v_t_new(1:nfrag), & - fragments%t_unit(:,1:nfrag), fragments%mass(1:nfrag), impactors%vbcom(:)) - do concurrent (i = 1:nfrag) - vc(:,i) = vb(:,i) - impactors%vbcom(:) - end do - kearr = 0.0_DP - do concurrent(i = 1:nfrag) - kearr(i) = fragments%mass(i) * dot_product(vc(:,i), vc(:,i)) - end do - keo = 0.5_DP * sum(kearr(:)) - fval = keo - if (lfirst_func) then - fval_scale = keo - lfirst_func = .false. - end if - fval = keo / fval_scale - lfailure = .false. - end select - end associate - return - end function tangential_objective_function + module subroutine fraggle_generate_rot_vec(collider) + !! Author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton + !! + !! Calculates the spins of a collection of fragments such that they conserve angular momentum + implicit none + ! Arguments + class(collision_fraggle), intent(inout) :: collider !! Fraggle collision system object + ! Internals + real(DP), dimension(NDIM) :: Lbefore, Lafter, Lspin, rotdir + real(DP) :: v_init, v_final, mass_init, mass_final, rotmag + integer(I4B) :: i + + associate(fragments => collider%fragments, impactors => collider%impactors, nfrag => collider%fragments%nbody) + + ! Torque the first body based on the change in angular momentum betwen the pre- and post-impact system assuming a simple bounce + mass_init = impactors%mass(2) + mass_final = sum(fragments%mass(2:nfrag)) + v_init = .mag.(impactors%vb(:,2) - impactors%vb(:,1)) + v_final = sqrt(mass_init / mass_final * v_init**2 - 2 * impactors%Qloss / mass_final) + + Lbefore(:) = mass_init * (impactors%rb(:,2) - impactors%rb(:,1)) .cross. (impactors%vb(:,2) - impactors%vb(:,1)) + + Lafter(:) = mass_final * (impactors%rb(:,2) - impactors%rb(:,1)) .cross. (v_final * impactors%bounce_unit(:)) + Lspin(:) = impactors%Lspin(:,1) + (Lbefore(:) - Lafter(:)) + fragments%rot(:,1) = Lspin(:) / (fragments%mass(1) * fragments%radius(1)**2 * fragments%Ip(3,1)) + + ! Add in some random spin noise. The magnitude will be scaled by the before-after amount and the direction will be random + do concurrent(i = 2:nfrag) + call random_number(rotdir) + rotdir = rotdir - 0.5_DP + rotdir = .unit. rotdir + call random_number(rotmag) + rotmag = rotmag * .mag. (Lbefore(:) - Lafter(:)) / ((nfrag - 1) * fragments%mass(i) * fragments%radius(i)**2 * fragments%Ip(3,i)) + fragments%rot(:,i) = rotmag * rotdir + end do - end subroutine fraggle_generate_tan_vel + end associate + return + end subroutine fraggle_generate_rot_vec - subroutine fraggle_generate_rad_vel(collider, lfailure) - !! Author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton + + module subroutine fraggle_generate_vel_vec(collider) + !! Author: David A. Minton !! - !! - !! Adjust the fragment velocities to set the fragment orbital kinetic energy. This will minimize the difference between the fragment kinetic energy and the energy budget + !! Generates an initial velocity distribution. For disruptions, the velocity magnitude is set to be + !! 2x the escape velocity of the colliding pair. For hit and runs the velocity magnitude is set to be + !! 2x the escape velocity of the smallest of the two bodies. implicit none ! Arguments class(collision_fraggle), intent(inout) :: collider !! Fraggle collision system object - logical, intent(out) :: lfailure !! Logical flag indicating whether this step fails or succeeds! ! Internals - real(DP), parameter :: TOL_MIN = FRAGGLE_ETOL ! This needs to be more accurate than the tangential step, as we are trying to minimize the total residual energy - real(DP), parameter :: TOL_INIT = 1e-14_DP - real(DP), parameter :: VNOISE_MAG = 1e-10_DP !! Magnitude of the noise to apply to initial conditions to help minimizer find a solution in case of failure - integer(I4B), parameter :: MAXLOOP = 100 - real(DP) :: ke_radial, tol, fval - integer(I4B) :: i - real(DP), dimension(:), allocatable :: v_r_initial - real(DP), dimension(collider%fragments%nbody) :: vnoise - type(lambda_obj) :: objective_function - character(len=STRMAX) :: message - real(DP), dimension(:), allocatable :: v_r_output - - associate(impactors => collider%impactors, nfrag => collider%fragments%nbody) - select type(fragments => collider%fragments) - class is (fraggle_fragments(*)) - ! Set the "target" ke for the radial component - ke_radial = fragments%ke_budget - fragments%ke_spin - fragments%ke_orbit - - do i = 1, nfrag - fragments%v_t_mag(i) = dot_product(fragments%vc(:,i), fragments%t_unit(:,i)) - fragments%v_r_mag(i) = dot_product(fragments%vc(:,i), fragments%r_unit(:,i)) - end do - - allocate(v_r_initial, source=fragments%v_r_mag) - - ! Initialize the lambda function using a structure constructor that calls the init method - ! Minimize the ke objective function using the BFGS optimizer - objective_function = lambda_obj(radial_objective_function) - tol = TOL_INIT - do while(tol < TOL_MIN) - fval = radial_objective_function(v_r_initial) - call minimize_bfgs(objective_function, nfrag, v_r_initial, tol, MAXLOOP, lfailure, v_r_output) - fragments%v_r_mag = v_r_output - if (.not.lfailure) exit - tol = tol * 2 ! Keep increasing the tolerance until we converge on a solution - v_r_initial(:) = fragments%v_r_mag(:) - call random_number(vnoise(1:nfrag)) ! Adding a bit of noise to the initial conditions helps it find a solution more often - vnoise(:) = 1.0_DP + VNOISE_MAG * (2 * vnoise(:) - 1._DP) - v_r_initial(:) = v_r_initial(:) * vnoise(:) - end do - - ! Shift the radial velocity vectors to align with the center of mass of the collisional system (the origin) - fragments%ke_orbit = 0.0_DP - fragments%vb(:,1:nfrag) = fraggle_util_vmag_to_vb(fragments%v_r_mag(1:nfrag), fragments%r_unit(:,1:nfrag), & - fragments%v_t_mag(1:nfrag), fragments%t_unit(:,1:nfrag), fragments%mass(1:nfrag), impactors%vbcom(:)) - do i = 1, nfrag - fragments%vc(:, i) = fragments%vb(:, i) - impactors%vbcom(:) - fragments%ke_orbit = fragments%ke_orbit + fragments%mass(i) * dot_product(fragments%vc(:, i), fragments%vc(:, i)) - end do - fragments%ke_orbit = 0.5_DP * fragments%ke_orbit - - lfailure = abs((fragments%ke_budget - (fragments%ke_orbit + fragments%ke_spin)) / fragments%ke_budget) > FRAGGLE_ETOL - if (lfailure) then - call swiftest_io_log_one_message(COLLISION_LOG_OUT, " ") - call swiftest_io_log_one_message(COLLISION_LOG_OUT, "Radial velocity failure diagnostics") - write(message, *) fragments%ke_budget - call swiftest_io_log_one_message(COLLISION_LOG_OUT, "ke_budget : " // trim(adjustl(message))) - write(message, *) fragments%ke_spin - call swiftest_io_log_one_message(COLLISION_LOG_OUT, "ke_spin : " // trim(adjustl(message))) - write(message, *) fragments%ke_orbit - call swiftest_io_log_one_message(COLLISION_LOG_OUT, "ke_orbit : " // trim(adjustl(message))) - write(message, *) fragments%ke_budget - (fragments%ke_orbit + fragments%ke_spin) - call swiftest_io_log_one_message(COLLISION_LOG_OUT, "ke_remainder : " // trim(adjustl(message))) + integer(I4B) :: i,j, loop, istart + logical :: lhitandrun, lnoncat + real(DP), dimension(NDIM) :: vimp_unit, rimp, vrot, Lresidual + real(DP) :: vmag, vesc, rotmag, ke_residual, ke_per_dof + integer(I4B), dimension(collider%fragments%nbody) :: vsign + real(DP), dimension(collider%fragments%nbody) :: vscale, mass_vscale, ke_avail + integer(I4B), parameter :: MAXLOOP = 100 + + associate(fragments => collider%fragments, impactors => collider%impactors, nfrag => collider%fragments%nbody) + 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 + + ! The fragments will be divided into two "clouds" based on identified origin body. + ! These clouds will collectively travel like two impactors bouncing off of each other. + where(fragments%origin_body(:) == 1) + vsign(:) = -1 + elsewhere + vsign(:) = 1 + end where + + ! The minimum fragment velocity will be set by the escape velocity + if (lhitandrun) then + vesc = sqrt(2 * impactors%Gmass(2) / impactors%radius(2)) + else + vesc = sqrt(2 * sum(impactors%Gmass(:)) / sum(impactors%radius(:))) + end if + + ! Scale the magnitude of the velocity by the distance from the impact point + ! This will reduce the chances of fragments colliding with each other immediately, and is more physically correct + do concurrent(i = 1:nfrag) + rimp(:) = fragments%rc(:,i) - impactors%rbimp(:) + vscale(i) = .mag. rimp(:) / (.mag. (impactors%rb(:,2) - impactors%rb(:,1))) + end do + vscale(:) = vscale(:)/minval(vscale(:)) + + ! Give the fragment velocities a random value that is scaled with fragment mass + call random_number(mass_vscale) + mass_vscale(:) = (mass_vscale(:) + 1.0_DP) / 2 + mass_vscale(:) = mass_vscale(:) * (fragments%mtot / fragments%mass(:))**(0.125_DP) ! The 1/8 power is arbitrary. It just gives the velocity a small mass dependence + mass_vscale(:) = mass_vscale(:) / minval(mass_vscale(:)) + + ! Set the velocities of all fragments using all of the scale factors determined above + do concurrent(i = 1:nfrag) + j = fragments%origin_body(i) + vrot(:) = impactors%rot(:,j) .cross. (fragments%rc(:,i) - impactors%rc(:,j)) + if (lhitandrun) then + if (i == 1) then + fragments%vc(:,1) = impactors%vc(:,1) + else + vmag = .mag.impactors%vc(:,2) / (maxval(mass_vscale(:) * maxval(vscale(:)))) + fragments%vc(:,i) = vmag * mass_vscale(i) * vscale(i) * impactors%bounce_unit(:) * vsign(i) + vrot(:) + end if + else + ! Add more velocity dispersion to disruptions vs hit and runs. + vmag = vesc * vscale(i) * mass_vscale(i) + rimp(:) = fragments%rc(:,i) - impactors%rbimp(:) + vimp_unit(:) = .unit. rimp(:) + fragments%vc(:,i) = vmag * (impactors%bounce_unit(:) + vimp_unit(:)) * vsign(i) + vrot(:) end if - end select + end do + + if (lhitandrun) then + istart = 2 + else + istart = 1 + end if + do loop = 1, MAXLOOP + call collider%set_coordinate_system() + call fragments%get_kinetic_energy() + ke_residual = fragments%ke_budget - (fragments%ke_orbit + fragments%ke_spin) + ! Make sure we don't take away too much orbital kinetic energy, otherwise the fragment can't escape + ke_avail(:) = fragments%ke_orbit_frag(:) - impactors%Gmass(1)*impactors%mass(2)/fragments%vmag(:) + ke_per_dof = -ke_residual/(nfrag - istart + 1) + do concurrent(i = istart:nfrag, ke_avail(i) > ke_per_dof) + fragments%vmag(i) = sqrt(2 * (fragments%ke_orbit_frag(i) - ke_per_dof)/fragments%mass(i)) + fragments%vc(:,i) = fragments%vmag(i) * fragments%v_unit(:,i) + end do + call fragments%get_kinetic_energy() + ke_residual = fragments%ke_budget - (fragments%ke_orbit + fragments%ke_spin) + + ! Check for any residual angular momentum, and if there is any, put it into spin of the biggest body + call collider%set_coordinate_system() + call fragments%get_angular_momentum() + Lresidual(:) = fragments%L_budget(:) - (fragments%Lorbit(:) + fragments%Lspin(:)) + rotmag = .mag. fragments%rot(:,1) + fragments%rot(:,1) = fragments%rot(:,1) + Lresidual(:) / (fragments%mass(1) * fragments%radius(1)**2 * fragments%Ip(:,1)) + rotmag = .mag. fragments%rot(:,1) + + if (ke_residual >= 0.0_DP) exit + + end do + + do concurrent(i = 1:nfrag) + fragments%vb(:,i) = fragments%vc(:,i) + impactors%vbcom(:) + end do + + impactors%vbcom(:) = 0.0_DP + do concurrent(i = 1:nfrag) + impactors%vbcom(:) = impactors%vbcom(:) + fragments%mass(i) * fragments%vb(:,i) + end do + impactors%vbcom(:) = impactors%vbcom(:) / fragments%mtot + end associate return + end subroutine fraggle_generate_vel_vec - contains - function radial_objective_function(v_r_mag_input) result(fval) - !! Author: David A. Minton - !! - !! Objective function for evaluating how close our fragment velocities get to minimizing KE error from our required value - implicit none - ! Arguments - real(DP), dimension(:), intent(in) :: v_r_mag_input !! Unknown radial component of fragment velocity vector - ! Result - real(DP) :: fval !! The objective function result, which is the square of the difference between the calculated fragment kinetic energy and our target - !! Minimizing this brings us closer to our objective - ! Internals - integer(I4B) :: i - real(DP), dimension(:,:), allocatable :: v_shift - real(DP), dimension(collider%fragments%nbody) :: kearr - real(DP) :: keo, ke_radial, rotmag2, vmag2 - - associate(impactors => collider%impactors, nfrag => collider%fragments%nbody) - select type(fragments => collider%fragments) - class is (fraggle_fragments(*)) - allocate(v_shift, mold=fragments%vb) - v_shift(:,:) = fraggle_util_vmag_to_vb(v_r_mag_input, fragments%r_unit, fragments%v_t_mag, fragments%t_unit, fragments%mass, impactors%vbcom) - do i = 1,fragments%nbody - v_shift(:,i) = v_shift(:,i) - impactors%vbcom(:) - rotmag2 = fragments%rot(1,i)**2 + fragments%rot(2,i)**2 + fragments%rot(3,i)**2 - vmag2 = v_shift(1,i)**2 + v_shift(2,i)**2 + v_shift(3,i)**2 - kearr(i) = fragments%mass(i) * (fragments%Ip(3, i) * fragments%radius(i)**2 * rotmag2 + vmag2) - end do - keo = 2 * fragments%ke_budget - sum(kearr(:)) - ke_radial = fragments%ke_budget - fragments%ke_orbit - fragments%ke_spin - ! The following ensures that fval = 0 is a local minimum, which is what the BFGS method is searching for - fval = (keo / (2 * ke_radial))**2 - end select - end associate - - return - end function radial_objective_function - - end subroutine fraggle_generate_rad_vel end submodule s_fraggle_generate diff --git a/src/fraggle/fraggle_module.f90 b/src/fraggle/fraggle_module.f90 index 2c5ee9672..130b3af5c 100644 --- a/src/fraggle/fraggle_module.f90 +++ b/src/fraggle/fraggle_module.f90 @@ -17,9 +17,6 @@ module fraggle !> Class definition for the variables that describe a collection of fragments by Fraggle barycentric coordinates type, extends(collision_fragments) :: fraggle_fragments - real(DP), dimension(nbody) :: v_r_mag !! Array of radial direction velocity magnitudes of individual fragments - real(DP), dimension(nbody) :: v_t_mag !! Array of tangential direction velocity magnitudes of individual fragments - real(DP), dimension(nbody) :: v_n_mag !! Array of normal direction velocity magnitudes of individual fragments contains procedure :: reset => fraggle_util_reset_fragments !! Resets all position and velocity-dependent fragment quantities in order to do a fresh calculation (does not reset mass, radius, or other values that get set prior to the call to fraggle_generate) @@ -27,7 +24,7 @@ module fraggle end type fraggle_fragments - type, extends(collision_disruption) :: collision_fraggle + type, extends(collision_basic) :: collision_fraggle ! Scale factors used to scale dimensioned quantities to a more "natural" system where important quantities (like kinetic energy, momentum) are of order ~1 real(DP) :: dscale = 1.0_DP !! Distance dimension scale factor real(DP) :: mscale = 1.0_DP !! Mass scale factor @@ -37,16 +34,27 @@ module fraggle real(DP) :: Lscale = 1.0_DP !! Angular momentum scale factor (a convenience unit that is derived from dscale, tscale, and mscale) contains procedure :: disrupt => fraggle_generate_disrupt !! Generates a system of fragments in barycentric coordinates that conserves energy and momentum. + procedure :: generate => fraggle_generate !! A simple disruption models that does not constrain energy loss in collisions + procedure :: hitandrun => fraggle_generate_hitandrun + procedure :: set_mass_dist => fraggle_util_set_mass_dist !! Sets the distribution of mass among the fragments depending on the regime type procedure :: set_natural_scale => fraggle_util_set_natural_scale_factors !! Scales dimenional quantities to ~O(1) with respect to the collisional system. procedure :: set_original_scale => fraggle_util_set_original_scale_factors !! Restores dimenional quantities back to the original system units procedure :: setup_fragments => fraggle_util_setup_fragments_system !! Initializer for the fragments of the collision system. - procedure :: construct_temporary_system => fraggle_util_construct_temporary_system !! Constructs temporary n-body system in order to compute pre- and post-impact energy and momentum procedure :: reset => fraggle_util_reset_system !! Deallocates all allocatables final :: fraggle_final_system !! Finalizer will deallocate all allocatables end type collision_fraggle - interface + + + module subroutine fraggle_generate(self, nbody_system, param, t) + implicit none + class(collision_fraggle), intent(inout) :: self !! Fraggle 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 + end subroutine fraggle_generate + module subroutine fraggle_generate_disrupt(self, nbody_system, param, t, lfailure) implicit none class(collision_fraggle), intent(inout) :: self !! Fraggle system object the outputs will be the fragmentation @@ -56,21 +64,35 @@ module subroutine fraggle_generate_disrupt(self, nbody_system, param, t, lfailur logical, optional, intent(out) :: lfailure !! Answers the question: Should this have been a merger instead? end subroutine fraggle_generate_disrupt + module subroutine fraggle_generate_hitandrun(self, nbody_system, param, t) + implicit none + class(collision_fraggle), intent(inout) :: self !! Collision system object + class(base_nbody_system), intent(inout) :: nbody_system !! Swiftest nbody system object + class(base_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA additions + real(DP), intent(in) :: t !! Time of collision + end subroutine fraggle_generate_hitandrun + + module subroutine fraggle_generate_pos_vec(collider) + implicit none + class(collision_fraggle), intent(inout) :: collider !! Fraggle ollision system object + end subroutine fraggle_generate_pos_vec + + module subroutine fraggle_generate_rot_vec(collider) + implicit none + class(collision_fraggle), intent(inout) :: collider !! Collision system object + end subroutine fraggle_generate_rot_vec + + module subroutine fraggle_generate_vel_vec(collider) + implicit none + class(collision_fraggle), intent(inout) :: collider !! Collision system object + end subroutine fraggle_generate_vel_vec + module subroutine fraggle_util_setup_fragments_system(self, nfrag) implicit none class(collision_fraggle), intent(inout) :: self !! Encounter collision system object integer(I4B), intent(in) :: nfrag !! Number of fragments to create end subroutine fraggle_util_setup_fragments_system - module subroutine fraggle_util_construct_temporary_system(self, nbody_system, param, tmpsys, tmpparam) - implicit none - class(collision_fraggle), intent(inout) :: self !! Fraggle collision system object - class(base_nbody_system), intent(in) :: nbody_system !! Original swiftest nbody system object - class(base_parameters), intent(in) :: param !! Current swiftest run configuration parameters - class(base_nbody_system), allocatable, intent(out) :: tmpsys !! Output temporary swiftest nbody system object - class(base_parameters), allocatable, intent(out) :: tmpparam !! Output temporary configuration run parameters - end subroutine fraggle_util_construct_temporary_system - module subroutine fraggle_util_reset_fragments(self) implicit none class(fraggle_fragments(*)), intent(inout) :: self @@ -81,6 +103,12 @@ module subroutine fraggle_util_reset_system(self) class(collision_fraggle), intent(inout) :: self !! Collision system object end subroutine fraggle_util_reset_system + module subroutine fraggle_util_set_mass_dist(self, param) + implicit none + class(collision_fraggle), intent(inout) :: self !! Fraggle collision object + class(base_parameters), intent(in) :: param !! Current Swiftest run configuration parameters + end subroutine fraggle_util_set_mass_dist + module subroutine fraggle_util_set_natural_scale_factors(self) implicit none class(collision_fraggle), intent(inout) :: self !! Fraggle collision system object @@ -91,15 +119,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 function fraggle_util_vmag_to_vb(v_r_mag, r_unit, v_t_mag, t_unit, m_frag, vcom) result(vb) - implicit none - real(DP), dimension(:), intent(in) :: v_r_mag !! Unknown radial component of fragment velocity vector - real(DP), dimension(:), intent(in) :: v_t_mag !! Tangential component of velocity vector set previously by angular momentum constraint - real(DP), dimension(:,:), intent(in) :: r_unit, t_unit !! Radial and tangential unit vectors for each fragment - real(DP), dimension(:), intent(in) :: m_frag !! Fragment masses - real(DP), dimension(:), intent(in) :: vcom !! Barycentric velocity of collisional system center of mass - real(DP), dimension(:,:), allocatable :: vb - end function fraggle_util_vmag_to_vb end interface contains diff --git a/src/fraggle/fraggle_util.f90 b/src/fraggle/fraggle_util.f90 index cff4fd832..e5cdaa20d 100644 --- a/src/fraggle/fraggle_util.f90 +++ b/src/fraggle/fraggle_util.f90 @@ -11,32 +11,6 @@ use swiftest contains - module subroutine fraggle_util_construct_temporary_system(self, nbody_system, param, tmpsys, tmpparam) - !! Author: David A. Minton - !! - !! Constructs a temporary internal system consisting of active bodies and additional fragments. This internal temporary system is used to calculate system energy with and without fragments - implicit none - ! Arguments - class(collision_fraggle), intent(inout) :: self !! Fraggle collision system object - class(base_nbody_system), intent(in) :: nbody_system !! Original swiftest nbody system object - class(base_parameters), intent(in) :: param !! Current swiftest run configuration parameters - class(base_nbody_system), allocatable, intent(out) :: tmpsys !! Output temporary swiftest nbody system object - class(base_parameters), allocatable, intent(out) :: tmpparam !! Output temporary configuration run parameters - - call collision_util_construct_temporary_system(self, nbody_system, param, tmpsys, tmpparam) - - select type(tmpsys) - class is (swiftest_nbody_system) - select type(tmpparam) - class is (swiftest_parameters) - call tmpsys%rescale(tmpparam, self%mscale, self%dscale, self%tscale) - end select - end select - - return - end subroutine fraggle_util_construct_temporary_system - - module subroutine fraggle_util_reset_fragments(self) !! author: David A. Minton !! @@ -58,9 +32,6 @@ module subroutine fraggle_util_reset_fragments(self) self%rmag(:) = 0.0_DP self%rotmag(:) = 0.0_DP - self%v_r_mag(:) = 0.0_DP - self%v_t_mag(:) = 0.0_DP - self%v_n_mag(:) = 0.0_DP return end subroutine fraggle_util_reset_fragments @@ -87,6 +58,167 @@ module subroutine fraggle_util_reset_system(self) end subroutine fraggle_util_reset_system + module subroutine fraggle_util_set_mass_dist(self, param) + !! author: David A. Minton + !! + !! Sets the mass of fragments based on the mass distribution returned by the regime calculation. + !! This subroutine must be run after the the setup routine has been run on the fragments + !! + implicit none + ! Arguments + class(collision_fraggle), intent(inout) :: self !! Fraggle collision system object + class(base_parameters), intent(in) :: param !! Current Swiftest run configuration parameters + ! Internals + 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, 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) + integer(I4B), parameter :: NFRAG_SIZE_MULTIPLIER = 3 !! Log-space scale factor that scales the number of fragments by the collisional system mass + 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 + volume(1:2) = 4._DP / 3._DP * PI * impactors%radius(1:2)**3 + mtot = sum(impactors%mass(:)) + Ip_avg(:) = (impactors%mass(1) * impactors%Ip(:,1) + impactors%mass(2) * impactors%Ip(:,2)) / mtot + + if (impactors%mass(1) > impactors%mass(2)) then + jtarg = 1 + jproj = 2 + else + jtarg = 2 + jproj = 1 + end if + + select case(impactors%regime) + case(COLLRESOLVE_REGIME_DISRUPTION, COLLRESOLVE_REGIME_SUPERCATASTROPHIC, COLLRESOLVE_REGIME_HIT_AND_RUN) + ! The first two bins of the mass_dist are the largest and second-largest fragments that came out of collision_regime. + ! The remainder from the third bin will be distributed among nfrag-2 bodies. The following code will determine nfrag based on + ! the limits bracketed above and the model size distribution of fragments. + ! Check to see if our size distribution would give us a smaller number of fragments than the maximum number + + select type(param) + class is (swiftest_parameters) + min_mfrag = (param%min_GMfrag / param%GU) + ! The number of fragments we generate is bracked by the minimum required by fraggle_generate (7) and the + ! maximum set by the NFRAG_SIZE_MULTIPLIER which limits the total number of fragments to prevent the nbody + ! code from getting an overwhelmingly large number of fragments + nfrag = ceiling(NFRAG_SIZE_MULTIPLIER * log(mtot / min_mfrag)) + nfrag = max(min(nfrag, NFRAGMAX), NFRAGMIN) + class default + min_mfrag = 0.0_DP + nfrag = NFRAGMAX + end select + + i = iMrem + mremaining = impactors%mass_dist(iMrem) + do while (i <= nfrag) + mfrag = (1 + i - iMslr)**(-3._DP / BETA) * impactors%mass_dist(iMslr) + if (mremaining - mfrag < 0.0_DP) exit + mremaining = mremaining - mfrag + i = i + 1 + end do + if (i < nfrag) nfrag = max(i, NFRAGMIN) ! The sfd would actually give us fewer fragments than our maximum + call self%setup_fragments(nfrag) + + case (COLLRESOLVE_REGIME_MERGE, COLLRESOLVE_REGIME_GRAZE_AND_MERGE) + + call self%setup_fragments(1) + select type(fragments => self%fragments) + class is (collision_fragments(*)) + fragments%mass(1) = impactors%mass_dist(1) + fragments%radius(1) = impactors%radius(jtarg) + fragments%density(1) = impactors%mass_dist(1) / volume(jtarg) + if (param%lrotation) fragments%Ip(:, 1) = impactors%Ip(:,1) + end select + return + case default + write(*,*) "collision_util_set_mass_dist_fragments error: Unrecognized regime code",impactors%regime + end select + + select type(fragments => self%fragments) + class is (collision_fragments(*)) + fragments%mtot = mtot + + ! Make the first two bins the same as the Mlr and Mslr values that came from collision_regime + fragments%mass(1) = impactors%mass_dist(iMlr) + fragments%mass(2) = impactors%mass_dist(iMslr) + + ! Distribute the remaining mass the 3:nfrag bodies following the model SFD given by slope BETA + mremaining = impactors%mass_dist(iMrem) + do i = iMrem, nfrag + mfrag = (1 + i - iMslr)**(-3._DP / BETA) * impactors%mass_dist(iMslr) + fragments%mass(i) = mfrag + mremaining = mremaining - mfrag + end do + + ! If there is any residual mass (either positive or negative) we will distribute remaining mass proportionally among the the fragments + if (mremaining < 0.0_DP) then ! If the remainder is negative, this means that that the number of fragments required by the SFD is smaller than our lower limit set by fraggle_generate. + istart = iMrem ! We will reduce the mass of the 3:nfrag bodies to prevent the second-largest fragment from going smaller + else ! If the remainder is postiive, this means that the number of fragments required by the SFD is larger than our upper limit set by computational expediency. + istart = iMslr ! We will increase the mass of the 2:nfrag bodies to compensate, which ensures that the second largest fragment remains the second largest + end if + mfrag = 1._DP + mremaining / sum(fragments%mass(istart:nfrag)) + fragments%mass(istart:nfrag) = fragments%mass(istart:nfrag) * mfrag + + ! There may still be some small residual due to round-off error. If so, simply add it to the last bin of the mass distribution. + mremaining = fragments%mtot - sum(fragments%mass(1:nfrag)) + fragments%mass(nfrag) = fragments%mass(nfrag) + mremaining + + ! Compute physical properties of the new fragments + select case(impactors%regime) + case(COLLRESOLVE_REGIME_HIT_AND_RUN) ! The hit and run case always preserves the largest body intact, so there is no need to recompute the physical properties of the first fragment + fragments%radius(1) = impactors%radius(jtarg) + fragments%density(1) = impactors%mass_dist(iMlr) / volume(jtarg) + fragments%Ip(:, 1) = impactors%Ip(:,1) + istart = 2 + case default + istart = 1 + end select + + fragments%density(istart:nfrag) = fragments%mtot / sum(volume(:)) + fragments%radius(istart:nfrag) = (3 * fragments%mass(istart:nfrag) / (4 * PI * fragments%density(istart:nfrag)))**(1.0_DP / 3.0_DP) + do i = istart, nfrag + 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 + end subroutine fraggle_util_set_mass_dist + + module subroutine fraggle_util_set_natural_scale_factors(self) !! author: David A. Minton !! @@ -226,37 +358,4 @@ module subroutine fraggle_util_setup_fragments_system(self, nfrag) end subroutine fraggle_util_setup_fragments_system - module function fraggle_util_vmag_to_vb(v_r_mag, r_unit, v_t_mag, t_unit, m_frag, vcom) result(vb) - !! Author: David A. Minton - !! - !! Converts radial and tangential velocity magnitudes into barycentric velocity - implicit none - ! Arguments - real(DP), dimension(:), intent(in) :: v_r_mag !! Unknown radial component of fragment velocity vector - real(DP), dimension(:), intent(in) :: v_t_mag !! Tangential component of velocity vector set previously by angular momentum constraint - real(DP), dimension(:,:), intent(in) :: r_unit, t_unit !! Radial and tangential unit vectors for each fragment - real(DP), dimension(:), intent(in) :: m_frag !! Fragment masses - real(DP), dimension(:), intent(in) :: vcom !! Barycentric velocity of collisional system center of mass - ! Result - real(DP), dimension(:,:), allocatable :: vb - ! Internals - integer(I4B) :: i, nfrag - - allocate(vb, mold=r_unit) - ! Make sure the velocity magnitude stays positive - nfrag = size(m_frag) - do i = 1, nfrag - vb(:,i) = abs(v_r_mag(i)) * 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 collision_util_shift_vector_to_origin(m_frag, vb) - - do i = 1, nfrag - vb(:, i) = vb(:, i) + v_t_mag(i) * t_unit(:, i) + vcom(:) - end do - - return - end function fraggle_util_vmag_to_vb - - end submodule s_fraggle_util diff --git a/src/swiftest/swiftest_io.f90 b/src/swiftest/swiftest_io.f90 index 2ba53b5e5..5f1a90276 100644 --- a/src/swiftest/swiftest_io.f90 +++ b/src/swiftest/swiftest_io.f90 @@ -2086,15 +2086,14 @@ module subroutine swiftest_io_param_reader(self, unit, iotype, v_list, iostat, i if ((param%collision_model /= "MERGE") .and. & (param%collision_model /= "BOUNCE") .and. & - (param%collision_model /= "DISRUPTION") .and. & (param%collision_model /= "FRAGGLE")) then write(iomsg,*) 'Invalid collision_model parameter: ',trim(adjustl(param%out_type)) - write(iomsg,*) 'Valid options are MERGE, BOUNCE, DISRUPTION, or FRAGGLE' + write(iomsg,*) 'Valid options are MERGE, BOUNCE, or FRAGGLE' iostat = -1 return end if - if (param%collision_model == "FRAGGLE" .or. param%collision_model == "DISRUPTION") then + if (param%collision_model == "FRAGGLE" ) then if (seed_set) then call random_seed(put = param%seed) else diff --git a/src/swiftest/swiftest_util.f90 b/src/swiftest/swiftest_util.f90 index 49c2d0362..e3e0fccee 100644 --- a/src/swiftest/swiftest_util.f90 +++ b/src/swiftest/swiftest_util.f90 @@ -2642,8 +2642,6 @@ module subroutine swiftest_util_setup_construct_system(nbody_system, param) allocate(collision_basic :: nbody_system%collider) case("BOUNCE") allocate(collision_bounce :: nbody_system%collider) - case("DISRUPTION") - allocate(collision_disruption :: nbody_system%collider) case("FRAGGLE") allocate(collision_fraggle :: nbody_system%collider) end select