diff --git a/src/collision/collision_generate.f90 b/src/collision/collision_generate.f90 index aa2e79761..333e5d059 100644 --- a/src/collision/collision_generate.f90 +++ b/src/collision/collision_generate.f90 @@ -12,6 +12,26 @@ use swiftest contains + module subroutine collision_generate_basic(self, nbody_system, param, t) + !! author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton + !! + !! Merge massive bodies no matter the regime + !! + !! Adapted from David E. Kaufmann's Swifter routines symba_merge_pl.f90 and symba_discard_merge_pl.f90 + !! + !! Adapted from Hal Levison's Swift routines symba5_merge.f and discard_mass_merge.f + implicit none + ! Arguments + class(collision_basic), intent(inout) :: self !! Merge 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 + + call self%merge(nbody_system, param, t) + + return + end subroutine collision_generate_basic + module subroutine collision_generate_bounce(self, nbody_system, param, t) !! author: David A. Minton !! @@ -55,9 +75,9 @@ module subroutine collision_generate_bounce(self, nbody_system, param, t) end select end select case (COLLRESOLVE_REGIME_HIT_AND_RUN) - call collision_generate_hitandrun(self, nbody_system, param, t) + call self%hitandrun(nbody_system, param, t) case (COLLRESOLVE_REGIME_MERGE, COLLRESOLVE_REGIME_GRAZE_AND_MERGE) - call self%collision_merge%generate(nbody_system, param, t) ! Use the default collision model, which is merge + call self%merge(nbody_system, param, t) ! Use the default collision model, which is merge case default write(*,*) "Error in swiftest_collision, unrecognized collision regime" call util_exit(FAILURE) @@ -70,54 +90,178 @@ module subroutine collision_generate_bounce(self, nbody_system, param, t) end subroutine collision_generate_bounce - module subroutine collision_generate_hitandrun(collider, nbody_system, param, t) + module subroutine collision_generate_hitandrun(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 hit-and-run collision !! implicit none ! Arguments - class(collision_merge), intent(inout) :: collider !! Fraggle 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 + class(collision_basic), 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 integer(I4B) :: i, ibiggest, nfrag, jtarg, jproj logical :: lpure character(len=STRMAX) :: message real(DP) :: dpe + select type(nbody_system) class is (swiftest_nbody_system) - select type (pl => nbody_system%pl) + select type(pl => nbody_system%pl) class is (swiftest_pl) - - associate(impactors => collider%impactors, fragments => collider%fragments, pl => nbody_system%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))) - collider%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 => collider%before) - class is (swiftest_nbody_system) - select type(after => collider%after) - class is (swiftest_nbody_system) - allocate(after%pl, source=before%pl) - end select - end select + 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_simple_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%disrupt(nbody_system, param, t, lpure) + + 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 + + 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 + return + end subroutine collision_generate_hitandrun + + + module subroutine collision_generate_simple(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_simple_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, status + logical :: lfailure + 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, fragments => self%fragments, 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 collision_io_collider_message(nbody_system%pl, impactors%id, message) + call swiftest_io_log_one_message(COLLISION_LOG_OUT, message) + + ! Collisional fragments will be uniformly distributed around the pre-impact barycenter + call self%set_mass_dist(param) + + ! Generate the position and velocity distributions of the fragments + call self%disrupt(nbody_system, param, t) + + dpe = self%pe(2) - self%pe(1) + nbody_system%Ecollisions = nbody_system%Ecollisions - dpe + nbody_system%Euntracked = nbody_system%Euntracked + dpe + + ! 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 select + end select return - end subroutine collision_generate_hitandrun + end subroutine collision_generate_simple module subroutine collision_generate_merge(self, nbody_system, param, t) @@ -130,7 +274,7 @@ module subroutine collision_generate_merge(self, nbody_system, param, t) !! Adapted from Hal Levison's Swift routines symba5_merge.f and discard_mass_merge.f implicit none ! Arguments - class(collision_merge), intent(inout) :: self !! Merge fragment system object + class(collision_basic), intent(inout) :: self !! Merge 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 @@ -220,7 +364,7 @@ module subroutine collision_generate_merge(self, nbody_system, param, t) end subroutine collision_generate_merge - module subroutine collision_generate_simple_disruption(self, nbody_system, param, t) + 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 @@ -231,18 +375,16 @@ module subroutine collision_generate_simple_disruption(self, nbody_system, param class(base_nbody_system), intent(inout) :: nbody_system !! Swiftest nbody system object class(base_parameters), intent(inout) :: param !! Current run configuration parameters real(DP), intent(in) :: t !! The time of the collision + logical, optional, intent(out) :: lfailure ! Internals real(DP) :: r_max_start - associate(impactors => self%impactors, fragments => self%fragments, nfrag => self%fragments%nbody) - r_max_start = 1.1_DP * .mag.(impactors%rb(:,2) - impactors%rb(:,1)) - call collision_generate_simple_pos_vec(self, r_max_start) - call self%set_coordinate_system() - call collision_generate_simple_vel_vec(self) - end associate - + r_max_start = 1.1_DP * .mag.(self%impactors%rb(:,2) - self%impactors%rb(:,1)) + call collision_generate_simple_pos_vec(self, r_max_start) + call self%set_coordinate_system() + call collision_generate_simple_vel_vec(self) return - end subroutine collision_generate_simple_disruption + end subroutine collision_generate_disrupt module subroutine collision_generate_simple_pos_vec(collider, r_max_start) @@ -371,6 +513,4 @@ module subroutine collision_generate_simple_vel_vec(collider) end subroutine collision_generate_simple_vel_vec - - end submodule s_collision_generate \ No newline at end of file diff --git a/src/collision/collision_module.f90 b/src/collision/collision_module.f90 index 046341a93..51be502a4 100644 --- a/src/collision/collision_module.f90 +++ b/src/collision/collision_module.f90 @@ -113,11 +113,11 @@ module collision end type collision_fragments - type :: collision_merge + type :: collision_basic !! This class defines a collisional system that stores impactors and fragments. This is written so that various collision models (i.e. Fraggle) could potentially be used !! to resolve collision by defining extended types of encounters_impactors and/or encounetr_fragments !! - !! The generate method for this class is the merge model. This allows any extended type to have access to the merge procedure by selecting the collision_merge parent class + !! The generate method for this class is the merge model. This allows any extended type to have access to the merge procedure by selecting the collision_basic parent class class(collision_fragments(:)), allocatable :: fragments !! Object containing information on the pre-collision system class(collision_impactors), allocatable :: impactors !! Object containing information on the post-collision system class(base_nbody_system), allocatable :: before !! A snapshot of the subset of the nbody_system involved in the collision @@ -141,18 +141,21 @@ module collision procedure :: get_energy_and_momentum => collision_util_get_energy_momentum !! Calculates total nbody_system energy in either the pre-collision outcome state (lbefore = .true.) or the post-collision outcome state (lbefore = .false.) procedure :: reset => collision_util_reset_system !! Deallocates all allocatables procedure :: set_coordinate_system => collision_util_set_coordinate_collider !! Sets the coordinate system of the collisional system - procedure :: generate => collision_generate_merge !! Merges the impactors to make a single final body - end type collision_merge + procedure :: generate => collision_generate_basic !! Merges the impactors to make a single final body + procedure :: hitandrun => collision_generate_hitandrun !! Merges the impactors to make a single final body + procedure :: merge => collision_generate_merge !! Merges the impactors to make a single final body + end type collision_basic - type, extends(collision_merge) :: collision_bounce + type, extends(collision_basic) :: collision_bounce contains procedure :: generate => collision_generate_bounce !! If a collision would result in a disruption, "bounce" the bodies instead. final :: collision_final_bounce !! Finalizer will deallocate all allocatables end type collision_bounce - type, extends(collision_merge) :: collision_simple_disruption + type, extends(collision_basic) :: collision_simple_disruption contains - procedure :: generate => collision_generate_simple_disruption !! If a collision would result in a disruption [TODO: SOMETHING LIKE CHAMBERS 2012] + procedure :: generate => collision_generate_simple !! 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_simple_disruption @@ -182,7 +185,7 @@ module collision type, extends(encounter_snapshot) :: collision_snapshot logical :: lcollision !! Indicates that this snapshot contains at least one collision - class(collision_merge), allocatable :: collider !! Collider object at this snapshot + class(collision_basic), allocatable :: collider !! Collider object at this snapshot contains procedure :: write_frame => collision_io_netcdf_write_frame_snapshot !! Writes a frame of encounter data to file procedure :: get_idvals => collision_util_get_idvalues_snapshot !! Gets an array of all id values saved in this snapshot @@ -201,10 +204,34 @@ module collision interface + module subroutine collision_generate_basic(self, nbody_system, param, t) + implicit none + class(collision_basic), intent(inout) :: self !! Merge 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_basic + + module subroutine collision_generate_bounce(self, nbody_system, param, t) + implicit none + class(collision_bounce), intent(inout) :: self !! Bounce 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_bounce + + module subroutine collision_generate_disrupt(self, nbody_system, param, t, lfailure) + implicit none + class(collision_simple_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(collider, nbody_system, param, t) + module subroutine collision_generate_hitandrun(self, nbody_system, param, t) implicit none - class(collision_merge), intent(inout) :: collider !! Merge (or extended) collision system object + class(collision_basic), 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 @@ -212,27 +239,19 @@ end subroutine collision_generate_hitandrun module subroutine collision_generate_merge(self, nbody_system, param, t) implicit none - class(collision_merge), intent(inout) :: self !! Merge fragment nbody_system object + class(collision_basic), intent(inout) :: self !! Merge 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_merge - module subroutine collision_generate_bounce(self, nbody_system, param, t) - implicit none - class(collision_bounce), intent(inout) :: self !! Bounce 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_bounce - - module subroutine collision_generate_simple_disruption(self, nbody_system, param, t) + module subroutine collision_generate_simple(self, nbody_system, param, t) implicit none class(collision_simple_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_simple_disruption + end subroutine collision_generate_simple module subroutine collision_generate_simple_pos_vec(collider, r_max_start) implicit none @@ -364,14 +383,14 @@ end subroutine collision_resolve_pltp module subroutine collision_util_add_fragments_to_collider(self, nbody_system, param) implicit none - class(collision_merge), intent(in) :: self !! Collision system object + class(collision_basic), intent(in) :: self !! Collision system object class(base_nbody_system), intent(inout) :: nbody_system !! Swiftest nbody system object class(base_parameters), intent(in) :: param !! Current swiftest run configuration parameters end subroutine collision_util_add_fragments_to_collider module subroutine collision_util_construct_temporary_system(self, nbody_system, param, tmpsys, tmpparam) implicit none - class(collision_merge), intent(inout) :: self !! Collision system object + class(collision_basic), intent(inout) :: self !! 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 @@ -391,7 +410,7 @@ end subroutine collision_util_set_mass_dist module subroutine collision_util_set_coordinate_collider(self) implicit none - class(collision_merge), intent(inout) :: self !! collisional system + class(collision_basic), intent(inout) :: self !! collisional system end subroutine collision_util_set_coordinate_collider module subroutine collision_util_set_coordinate_impactors(self) @@ -401,18 +420,18 @@ end subroutine collision_util_set_coordinate_impactors module subroutine collision_util_setup_collider(self, nbody_system) implicit none - class(collision_merge), intent(inout) :: self !! Encounter collision system object + class(collision_basic), intent(inout) :: self !! Encounter collision system object class(base_nbody_system), intent(in) :: nbody_system !! Current nbody system. Used as a mold for the before/after snapshots end subroutine collision_util_setup_collider module subroutine collision_util_setup_impactors_collider(self) implicit none - class(collision_merge), intent(inout) :: self !! Encounter collision system object + class(collision_basic), intent(inout) :: self !! Encounter collision system object end subroutine collision_util_setup_impactors_collider module subroutine collision_util_setup_fragments_collider(self, nfrag) implicit none - class(collision_merge), intent(inout) :: self !! Encounter collision system object + class(collision_basic), intent(inout) :: self !! Encounter collision system object integer(I4B), intent(in) :: nfrag !! Number of fragments to create end subroutine collision_util_setup_fragments_collider @@ -431,7 +450,7 @@ end subroutine collision_util_get_idvalues_snapshot module subroutine collision_util_get_energy_momentum(self, nbody_system, param, lbefore) use base, only : base_nbody_system, base_parameters implicit none - class(collision_merge), intent(inout) :: self !! Encounter collision system object + class(collision_basic), intent(inout) :: self !! Encounter collision system object class(base_nbody_system), intent(inout) :: nbody_system !! Swiftest nbody system object class(base_parameters), intent(inout) :: param !! Current swiftest run configuration parameters logical, intent(in) :: lbefore !! Flag indicating that this the "before" state of the nbody_system, with impactors included and fragments excluded or vice versa @@ -449,7 +468,7 @@ end subroutine collision_util_reset_impactors module subroutine collision_util_reset_system(self) implicit none - class(collision_merge), intent(inout) :: self !! Collision system object + class(collision_basic), intent(inout) :: self !! Collision system object end subroutine collision_util_reset_system module subroutine collision_util_snapshot(self, param, nbody_system, t, arg) diff --git a/src/collision/collision_resolve.f90 b/src/collision/collision_resolve.f90 index af4e4583c..9e91b1ecd 100644 --- a/src/collision/collision_resolve.f90 +++ b/src/collision/collision_resolve.f90 @@ -327,7 +327,7 @@ module subroutine collision_resolve_mergeaddsub(nbody_system, param, t, status) select type(param) class is (swiftest_parameters) associate(pl => nbody_system%pl, pl_discards => nbody_system%pl_discards, info => nbody_system%pl%info, pl_adds => nbody_system%pl_adds, cb => nbody_system%cb, npl => pl%nbody, & - collision_merge => nbody_system%collider, impactors => nbody_system%collider%impactors,fragments => nbody_system%collider%fragments) + collision_basic => nbody_system%collider, impactors => nbody_system%collider%impactors,fragments => nbody_system%collider%fragments) ! Add the impactors%id bodies to the subtraction list nimpactors = impactors%ncoll @@ -437,7 +437,7 @@ module subroutine collision_resolve_mergeaddsub(nbody_system, param, t, status) end where ! Log the properties of the new bodies - select type(after => collision_merge%after) + select type(after => collision_basic%after) class is (swiftest_nbody_system) allocate(after%pl, source=plnew) end select diff --git a/src/collision/collision_util.f90 b/src/collision/collision_util.f90 index da68b42d6..d1cd1bf99 100644 --- a/src/collision/collision_util.f90 +++ b/src/collision/collision_util.f90 @@ -17,7 +17,7 @@ module subroutine collision_util_add_fragments_to_collider(self, nbody_system, p !! Adds fragments to the temporary system pl object implicit none ! Arguments - class(collision_merge), intent(in) :: self !! Collision system system object + class(collision_basic), intent(in) :: self !! Collision system system object class(base_nbody_system), intent(inout) :: nbody_system !! Swiftest nbody system object class(base_parameters), intent(in) :: param !! Current swiftest run configuration parameters ! Internals @@ -67,7 +67,7 @@ module subroutine collision_util_construct_temporary_system(self, nbody_system, !! 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_merge), intent(inout) :: self !! Fraggle collision system object + class(collision_basic), 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 @@ -184,7 +184,7 @@ module subroutine collision_util_get_energy_momentum(self, nbody_system, param, !! This will temporarily expand the massive body object in a temporary system object called tmpsys to feed it into symba_energy implicit none ! Arguments - class(collision_merge), intent(inout) :: self !! Encounter collision system object + class(collision_basic), intent(inout) :: self !! Encounter collision system object class(base_nbody_system), intent(inout) :: nbody_system !! Swiftest nbody system object class(base_parameters), intent(inout) :: param !! Current swiftest run configuration parameters logical, intent(in) :: lbefore !! Flag indicating that this the "before" state of the nbody_system, with impactors included and fragments excluded or vice versa @@ -353,7 +353,7 @@ module subroutine collision_util_reset_system(self) !! Resets the collider nbody_system and deallocates all allocatables implicit none ! Arguments - class(collision_merge), intent(inout) :: self !! Collision system object + class(collision_basic), intent(inout) :: self !! Collision system object select type(before => self%before) class is (swiftest_nbody_system) @@ -387,7 +387,7 @@ module subroutine collision_util_set_coordinate_collider(self) !! Defines the collisional coordinate nbody_system, including the unit vectors of both the nbody_system and individual fragments. implicit none ! Arguments - class(collision_merge), intent(inout) :: self !! Collisional nbody_system + class(collision_basic), intent(inout) :: self !! Collisional nbody_system ! Internals integer(I4B) :: i real(DP), dimension(NDIM, self%fragments%nbody) :: L_sigma @@ -610,7 +610,7 @@ module subroutine collision_util_setup_collider(self, nbody_system) !! but not fragments. Those are setup later when the number of fragments is known. implicit none ! Arguments - class(collision_merge), intent(inout) :: self !! Encounter collision system object + class(collision_basic), intent(inout) :: self !! Encounter collision system object class(base_nbody_system), intent(in) :: nbody_system !! Current nbody system. Used as a mold for the before/after snapshots call self%setup_impactors() @@ -630,7 +630,7 @@ module subroutine collision_util_setup_impactors_collider(self) !! Initializer for the impactors for the encounter collision system. Deallocates old impactors before creating new ones implicit none ! Arguments - class(collision_merge), intent(inout) :: self !! Encounter collision system object + class(collision_basic), intent(inout) :: self !! Encounter collision system object if (allocated(self%impactors)) deallocate(self%impactors) allocate(collision_impactors :: self%impactors) @@ -645,7 +645,7 @@ module subroutine collision_util_setup_fragments_collider(self, nfrag) !! Initializer for the fragments of the collision system. implicit none ! Arguments - class(collision_merge), intent(inout) :: self !! Encounter collision system object + class(collision_basic), intent(inout) :: self !! Encounter collision system object integer(I4B), intent(in) :: nfrag !! Number of fragments to create if (allocated(self%fragments)) deallocate(self%fragments) diff --git a/src/fraggle/fraggle_generate.f90 b/src/fraggle/fraggle_generate.f90 index 4424b2412..3c48b71bf 100644 --- a/src/fraggle/fraggle_generate.f90 +++ b/src/fraggle/fraggle_generate.f90 @@ -18,207 +18,19 @@ contains - module subroutine fraggle_generate_disruption(collider, 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) :: collider - class(swiftest_nbody_system), intent(inout) :: nbody_system !! Swiftest nbody system object - class(swiftest_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 - logical :: lfailure - character(len=STRMAX) :: message - real(DP) :: dpe - - associate(impactors => collider%impactors, fragments => collider%fragments, pl => nbody_system%pl, status => collider%status) - select case(impactors%regime) - case(COLLRESOLVE_REGIME_DISRUPTION) - message = "Disruption between" - case(COLLRESOLVE_REGIME_SUPERCATASTROPHIC) - message = "Supercatastrophic disruption between" - end select - call collision_io_collider_message(nbody_system%pl, impactors%id, message) - call swiftest_io_log_one_message(COLLISION_LOG_OUT, message) - - ! Collisional fragments will be uniformly distributed around the pre-impact barycenter - call collider%set_mass_dist(param) - - ! Generate the position and velocity distributions of the fragments - call fraggle_generate_fragments(collider, nbody_system, param, lfailure) - - dpe = collider%pe(2) - collider%pe(1) - nbody_system%Ecollisions = nbody_system%Ecollisions - dpe - nbody_system%Euntracked = nbody_system%Euntracked + dpe - - if (lfailure) then - call swiftest_io_log_one_message(COLLISION_LOG_OUT, "No fragment solution found, so treat as a pure hit-and-run") - status = ACTIVE - nfrag = 0 - pl%status(impactors%id(:)) = status - pl%ldiscard(impactors%id(:)) = .false. - pl%lcollision(impactors%id(:)) = .false. - select type(before => collider%before) - class is (swiftest_nbody_system) - select type(after => collider%after) - class is (swiftest_nbody_system) - allocate(after%pl, source=before%pl) ! Be sure to save the pl so that snapshots still work - end select - end select - else - ! 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 if - end associate - - return - end subroutine fraggle_generate_disruption - - - module subroutine fraggle_generate_hitandrun(collider, nbody_system, param, t) - !! author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton - !! - !! 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 - class(swiftest_nbody_system), intent(inout) :: nbody_system !! Swiftest nbody system object - class(swiftest_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 - integer(I4B) :: i, ibiggest, nfrag, jtarg, jproj - logical :: lpure - character(len=STRMAX) :: message - real(DP) :: dpe - - select type(before => collider%before) - class is (swiftest_nbody_system) - select type(after => collider%after) - class is (swiftest_nbody_system) - associate(impactors => collider%impactors,pl => nbody_system%pl) - 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 - - 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 collider%set_mass_dist(param) - - ! Generate the position and velocity distributions of the fragments - call fraggle_generate_fragments(collider, nbody_system, param, lpure) - - dpe = collider%pe(2) - collider%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 = collider%fragments%nbody - write(message, *) nfrag - call swiftest_io_log_one_message(COLLISION_LOG_OUT, "Generating " // trim(adjustl(message)) // " fragments") - end if - end if - 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. - allocate(after%pl, source=before%pl) ! Be sure to save the pl so that snapshots still work - else - ibiggest = impactors%id(maxloc(pl%Gmass(impactors%id(:)), dim=1)) - collider%fragments%id(1) = pl%id(ibiggest) - collider%fragments%id(2:nfrag) = [(i, i = param%maxid + 1, param%maxid + nfrag - 1)] - param%maxid = collider%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 - - - return - end subroutine fraggle_generate_hitandrun - - - module subroutine fraggle_generate_system(self, nbody_system, param, t) - implicit none - class(collision_fraggle), intent(inout) :: self !! Fraggle 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 - - select type(nbody_system) - class is (swiftest_nbody_system) - select type(param) - class is (swiftest_parameters) - select case (self%impactors%regime) - case (COLLRESOLVE_REGIME_DISRUPTION, COLLRESOLVE_REGIME_SUPERCATASTROPHIC) - call fraggle_generate_disruption(self, nbody_system, param, t) - case (COLLRESOLVE_REGIME_HIT_AND_RUN) - call fraggle_generate_hitandrun(self, nbody_system, param, t) - case (COLLRESOLVE_REGIME_MERGE, COLLRESOLVE_REGIME_GRAZE_AND_MERGE) - call self%collision_merge%generate(nbody_system, param, t) - case default - write(*,*) "Error in swiftest_collision, unrecognized collision regime" - call util_exit(FAILURE) - end select - end select - end select - - return - end subroutine fraggle_generate_system - - - module subroutine fraggle_generate_fragments(collider, nbody_system, param, lfailure) + module subroutine fraggle_generate_disrupt(self, nbody_system, param, t, lfailure) !! Author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton !! !! Generates a nbody_system of fragments in barycentric coordinates that conserves energy and momentum. use, intrinsic :: ieee_exceptions implicit none ! Arguments - class(collision_fraggle), intent(inout) :: collider !! Fraggle nbody_system object the outputs will be the fragmentation - class(swiftest_nbody_system), intent(inout) :: nbody_system !! Swiftest nbody system object - class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters - logical, intent(out) :: lfailure !! Answers the question: Should this have been a merger instead? - ! Internals + class(collision_fraggle), intent(inout) :: self !! Fraggle system object the outputs will be the fragmentation + 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 !! Time of collision + logical, optional, intent(out) :: lfailure !! Answers the question: Should this have been a merger instead? + ! Internals integer(I4B) :: i integer(I4B) :: try real(DP) :: r_max_start, f_spin, dEtot, dLmag @@ -238,9 +50,9 @@ module subroutine fraggle_generate_fragments(collider, nbody_system, param, lfai class is (swiftest_nbody_system) select type(param) class is (swiftest_parameters) - select type(fragments => collider%fragments) + select type(fragments => self%fragments) class is (fraggle_fragments(*)) - associate(impactors => collider%impactors, nfrag => fragments%nbody, pl => nbody_system%pl) + associate(impactors => self%impactors, nfrag => fragments%nbody, pl => nbody_system%pl) write(message,*) nfrag call swiftest_io_log_one_message(COLLISION_LOG_OUT, "Fraggle generating " // trim(adjustl(message)) // " fragments.") @@ -258,12 +70,12 @@ module subroutine fraggle_generate_fragments(collider, nbody_system, param, lfai lk_plpl = .false. end if - call collider%set_natural_scale() + call self%set_natural_scale() call fragments%reset() ! Calculate the initial energy of the nbody_system without the collisional family - call collider%get_energy_and_momentum(nbody_system, param, lbefore=.true.) + call self%get_energy_and_momentum(nbody_system, param, lbefore=.true.) ! Start out the fragments close to the initial separation distance. This will be increased if there is any overlap or we fail to find a solution r_max_start = 1.1_DP * .mag.(impactors%rb(:,2) - impactors%rb(:,1)) @@ -282,40 +94,40 @@ module subroutine fraggle_generate_fragments(collider, nbody_system, param, lfai lfailure = .false. call ieee_set_flag(ieee_all, .false.) ! Set all fpe flags to quiet - call collision_generate_simple_pos_vec(collider, r_max_start) - call collider%set_coordinate_system() + call collision_generate_simple_pos_vec(self, r_max_start) + call self%set_coordinate_system() ! Initial velocity guess will be the barycentric velocity of the colliding nbody_system so that the budgets are based on the much smaller collisional-frame velocities do concurrent (i = 1:nfrag) fragments%vb(:, i) = impactors%vbcom(:) end do - call collider%get_energy_and_momentum(nbody_system, param, lbefore=.false.) - call collider%set_budgets() + call self%get_energy_and_momentum(nbody_system, param, lbefore=.false.) + call self%set_budgets() - call collision_generate_simple_vel_vec(collider) + call collision_generate_simple_vel_vec(self) - call fraggle_generate_tan_vel(collider, lfailure) + call fraggle_generate_tan_vel(self, lfailure) if (lfailure) then call swiftest_io_log_one_message(COLLISION_LOG_OUT, "Fraggle failed to find tangential velocities") cycle end if - call fraggle_generate_rad_vel(collider, lfailure) + call fraggle_generate_rad_vel(self, lfailure) if (lfailure) then call swiftest_io_log_one_message(COLLISION_LOG_OUT, "Fraggle failed to find radial velocities") cycle end if - call fraggle_generate_spins(collider, f_spin, lfailure) + call fraggle_generate_spins(self, f_spin, lfailure) if (lfailure) then call swiftest_io_log_one_message(COLLISION_LOG_OUT, "Fraggle failed to find spins") cycle end if - call collider%get_energy_and_momentum(nbody_system, param, lbefore=.false.) - dEtot = collider%Etot(2) - collider%Etot(1) - dLmag = .mag. (collider%Ltot(:,2) - collider%Ltot(:,1)) + call self%get_energy_and_momentum(nbody_system, param, lbefore=.false.) + dEtot = self%Etot(2) - self%Etot(1) + dLmag = .mag. (self%Ltot(:,2) - self%Ltot(:,1)) exit lfailure = ((abs(dEtot + impactors%Qloss) > FRAGGLE_ETOL) .or. (dEtot > 0.0_DP)) @@ -326,9 +138,9 @@ module subroutine fraggle_generate_fragments(collider, nbody_system, param, lfai cycle end if - lfailure = ((abs(dLmag) / (.mag.collider%Ltot(:,1))) > FRAGGLE_LTOL) + lfailure = ((abs(dLmag) / (.mag.self%Ltot(:,1))) > FRAGGLE_LTOL) if (lfailure) then - write(message,*) dLmag / (.mag.collider%Ltot(:,1)) + write(message,*) dLmag / (.mag.self%Ltot(:,1)) call swiftest_io_log_one_message(COLLISION_LOG_OUT, "Fraggle failed due to high angular momentum error: " // & trim(adjustl(message))) cycle @@ -351,7 +163,7 @@ module subroutine fraggle_generate_fragments(collider, nbody_system, param, lfai trim(adjustl(message)) // " tries") end if - call collider%set_original_scale() + call self%set_original_scale() ! Restore the big array if (lk_plpl) call pl%flatten(param) @@ -362,7 +174,8 @@ module subroutine fraggle_generate_fragments(collider, nbody_system, param, lfai call ieee_set_halting_mode(IEEE_ALL,fpe_halting_modes) ! Save the current halting modes so we can turn them off temporarily return - end subroutine fraggle_generate_fragments + end subroutine fraggle_generate_disrupt + subroutine fraggle_generate_spins(collider, f_spin, lfailure) !! Author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton @@ -373,8 +186,8 @@ subroutine fraggle_generate_spins(collider, f_spin, lfailure) implicit none ! Arguments class(collision_fraggle), intent(inout) :: collider !! Fraggle collision system object - real(DP), intent(in) :: f_spin !! Fraction of energy or momentum that goes into spin (whichever gives the lowest kinetic energy) - logical, intent(out) :: lfailure !! Logical flag indicating whether this step fails or succeeds! + real(DP), intent(in) :: f_spin !! Fraction of energy or momentum that goes into spin (whichever gives the lowest kinetic energy) + logical, intent(out) :: lfailure !! Logical flag indicating whether this step fails or succeeds! ! Internals real(DP), dimension(NDIM) :: L_remainder, rot_L, rot_ke, L real(DP), dimension(NDIM,collider%fragments%nbody) :: frot_rand ! The random rotation factor applied to fragments diff --git a/src/fraggle/fraggle_module.f90 b/src/fraggle/fraggle_module.f90 index 7e364dc6f..fe0eaf8a3 100644 --- a/src/fraggle/fraggle_module.f90 +++ b/src/fraggle/fraggle_module.f90 @@ -43,11 +43,11 @@ module fraggle real(DP) :: Escale = 1.0_DP !! Energy scale factor (a convenience unit that is derived from dscale, tscale, and mscale) real(DP) :: Lscale = 1.0_DP !! Angular momentum scale factor (a convenience unit that is derived from dscale, tscale, and mscale) contains - procedure :: generate => fraggle_generate_system !! Generates a system of fragments in barycentric coordinates that conserves energy and momentum. - procedure :: set_budgets => fraggle_util_set_budgets !! Sets the energy and momentum budgets of the fragments based on the collider value - 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 :: disrupt => fraggle_generate_disrupt !! Generates a system of fragments in barycentric coordinates that conserves energy and momentum. + procedure :: set_budgets => fraggle_util_set_budgets !! Sets the energy and momentum budgets of the fragments based on the collider value + 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 @@ -55,39 +55,14 @@ module fraggle interface - - - module subroutine fraggle_generate_disruption(collider, nbody_system, param, t) - implicit none - class(collision_fraggle), intent(inout) :: collider !! Fraggle collision system object - class(swiftest_nbody_system), intent(inout) :: nbody_system !! Swiftest nbody system object - class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA additions - real(DP), intent(in) :: t !! Time of collision - end subroutine fraggle_generate_disruption - - module subroutine fraggle_generate_fragments(collider, nbody_system, param, lfailure) - implicit none - class(collision_fraggle), intent(inout) :: collider !! Fraggle system object the outputs will be the fragmentation - class(swiftest_nbody_system), intent(inout) :: nbody_system !! Swiftest nbody system object - class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters - logical, intent(out) :: lfailure !! Answers the question: Should this have been a merger instead? - end subroutine fraggle_generate_fragments - - module subroutine fraggle_generate_hitandrun(collider, nbody_system, param, t) - implicit none - class(collision_fraggle), intent(inout) :: collider !! Fraggle collision system object - class(swiftest_nbody_system), intent(inout) :: nbody_system !! Swiftest nbody system object - class(swiftest_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_system(self, nbody_system, param, t) + module subroutine fraggle_generate_disrupt(self, nbody_system, param, t, lfailure) implicit none - class(collision_fraggle), intent(inout) :: self !! Fraggle fragment system object + class(collision_fraggle), intent(inout) :: self !! Fraggle system object the outputs will be the fragmentation 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 !! Time of collision - end subroutine fraggle_generate_system + real(DP), intent(in) :: t !! Time of collision + logical, optional, intent(out) :: lfailure !! Answers the question: Should this have been a merger instead? + end subroutine fraggle_generate_disrupt module subroutine fraggle_util_setup_fragments_system(self, nfrag) implicit none diff --git a/src/fraggle/fraggle_util.f90 b/src/fraggle/fraggle_util.f90 index 928e3655e..2e2e55961 100644 --- a/src/fraggle/fraggle_util.f90 +++ b/src/fraggle/fraggle_util.f90 @@ -197,36 +197,36 @@ module subroutine fraggle_util_set_natural_scale_factors(self) ! Internals integer(I4B) :: i - associate(collision_merge => self, fragments => self%fragments, impactors => self%impactors) + associate(collider => self, fragments => self%fragments, impactors => self%impactors) ! Set scale factors - collision_merge%Escale = 0.5_DP * ( impactors%mass(1) * dot_product(impactors%vb(:,1), impactors%vb(:,1)) & + collider%Escale = 0.5_DP * ( impactors%mass(1) * dot_product(impactors%vb(:,1), impactors%vb(:,1)) & + impactors%mass(2) * dot_product(impactors%vb(:,2), impactors%vb(:,2))) - collision_merge%dscale = sum(impactors%radius(:)) - collision_merge%mscale = fragments%mtot - collision_merge%vscale = sqrt(collision_merge%Escale / collision_merge%mscale) - collision_merge%tscale = collision_merge%dscale / collision_merge%vscale - collision_merge%Lscale = collision_merge%mscale * collision_merge%dscale * collision_merge%vscale + collider%dscale = sum(impactors%radius(:)) + collider%mscale = fragments%mtot + collider%vscale = sqrt(collider%Escale / collider%mscale) + collider%tscale = collider%dscale / collider%vscale + collider%Lscale = collider%mscale * collider%dscale * collider%vscale ! Scale all dimensioned quantities of impactors and fragments - impactors%rbcom(:) = impactors%rbcom(:) / collision_merge%dscale - impactors%vbcom(:) = impactors%vbcom(:) / collision_merge%vscale - impactors%rbimp(:) = impactors%rbimp(:) / collision_merge%dscale - impactors%vbimp(:) = impactors%vbimp(:) / collision_merge%vscale - impactors%rb(:,:) = impactors%rb(:,:) / collision_merge%dscale - impactors%vb(:,:) = impactors%vb(:,:) / collision_merge%vscale - impactors%mass(:) = impactors%mass(:) / collision_merge%mscale - impactors%radius(:) = impactors%radius(:) / collision_merge%dscale - impactors%Lspin(:,:) = impactors%Lspin(:,:) / collision_merge%Lscale - impactors%Lorbit(:,:) = impactors%Lorbit(:,:) / collision_merge%Lscale + impactors%rbcom(:) = impactors%rbcom(:) / collider%dscale + impactors%vbcom(:) = impactors%vbcom(:) / collider%vscale + impactors%rbimp(:) = impactors%rbimp(:) / collider%dscale + impactors%vbimp(:) = impactors%vbimp(:) / collider%vscale + impactors%rb(:,:) = impactors%rb(:,:) / collider%dscale + impactors%vb(:,:) = impactors%vb(:,:) / collider%vscale + impactors%mass(:) = impactors%mass(:) / collider%mscale + impactors%radius(:) = impactors%radius(:) / collider%dscale + impactors%Lspin(:,:) = impactors%Lspin(:,:) / collider%Lscale + impactors%Lorbit(:,:) = impactors%Lorbit(:,:) / collider%Lscale do i = 1, 2 impactors%rot(:,i) = impactors%Lspin(:,i) / (impactors%mass(i) * impactors%radius(i)**2 * impactors%Ip(3, i)) end do - fragments%mtot = fragments%mtot / collision_merge%mscale - fragments%mass = fragments%mass / collision_merge%mscale - fragments%radius = fragments%radius / collision_merge%dscale - impactors%Qloss = impactors%Qloss / collision_merge%Escale + fragments%mtot = fragments%mtot / collider%mscale + fragments%mass = fragments%mass / collider%mscale + fragments%radius = fragments%radius / collider%dscale + impactors%Qloss = impactors%Qloss / collider%Escale end associate return @@ -248,51 +248,51 @@ module subroutine fraggle_util_set_original_scale_factors(self) call ieee_get_halting_mode(IEEE_ALL,fpe_halting_modes) ! Save the current halting modes so we can turn them off temporarily call ieee_set_halting_mode(IEEE_ALL,.false.) - associate(collision_merge => self, fragments => self%fragments, impactors => self%impactors) + associate(collider => self, fragments => self%fragments, impactors => self%impactors) ! Restore scale factors - impactors%rbcom(:) = impactors%rbcom(:) * collision_merge%dscale - impactors%vbcom(:) = impactors%vbcom(:) * collision_merge%vscale - impactors%rbimp(:) = impactors%rbimp(:) * collision_merge%dscale - impactors%vbimp(:) = impactors%vbimp(:) * collision_merge%vscale + impactors%rbcom(:) = impactors%rbcom(:) * collider%dscale + impactors%vbcom(:) = impactors%vbcom(:) * collider%vscale + impactors%rbimp(:) = impactors%rbimp(:) * collider%dscale + impactors%vbimp(:) = impactors%vbimp(:) * collider%vscale - impactors%mass = impactors%mass * collision_merge%mscale - impactors%radius = impactors%radius * collision_merge%dscale - impactors%rb = impactors%rb * collision_merge%dscale - impactors%vb = impactors%vb * collision_merge%vscale - impactors%Lspin = impactors%Lspin * collision_merge%Lscale + impactors%mass = impactors%mass * collider%mscale + impactors%radius = impactors%radius * collider%dscale + impactors%rb = impactors%rb * collider%dscale + impactors%vb = impactors%vb * collider%vscale + impactors%Lspin = impactors%Lspin * collider%Lscale do i = 1, 2 impactors%rot(:,i) = impactors%Lspin(:,i) * (impactors%mass(i) * impactors%radius(i)**2 * impactors%Ip(3, i)) end do - fragments%mtot = fragments%mtot * collision_merge%mscale - fragments%mass = fragments%mass * collision_merge%mscale - fragments%radius = fragments%radius * collision_merge%dscale - fragments%rot = fragments%rot / collision_merge%tscale - fragments%rc = fragments%rc * collision_merge%dscale - fragments%vc = fragments%vc * collision_merge%vscale + fragments%mtot = fragments%mtot * collider%mscale + fragments%mass = fragments%mass * collider%mscale + fragments%radius = fragments%radius * collider%dscale + fragments%rot = fragments%rot / collider%tscale + fragments%rc = fragments%rc * collider%dscale + fragments%vc = fragments%vc * collider%vscale do i = 1, fragments%nbody fragments%rb(:, i) = fragments%rc(:, i) + impactors%rbcom(:) fragments%vb(:, i) = fragments%vc(:, i) + impactors%vbcom(:) end do - impactors%Qloss = impactors%Qloss * collision_merge%Escale + impactors%Qloss = impactors%Qloss * collider%Escale - collision_merge%Lorbit(:,:) = collision_merge%Lorbit(:,:) * collision_merge%Lscale - collision_merge%Lspin(:,:) = collision_merge%Lspin(:,:) * collision_merge%Lscale - collision_merge%Ltot(:,:) = collision_merge%Ltot(:,:) * collision_merge%Lscale - collision_merge%ke_orbit(:) = collision_merge%ke_orbit(:) * collision_merge%Escale - collision_merge%ke_spin(:) = collision_merge%ke_spin(:) * collision_merge%Escale - collision_merge%pe(:) = collision_merge%pe(:) * collision_merge%Escale - collision_merge%Etot(:) = collision_merge%Etot(:) * collision_merge%Escale + collider%Lorbit(:,:) = collider%Lorbit(:,:) * collider%Lscale + collider%Lspin(:,:) = collider%Lspin(:,:) * collider%Lscale + collider%Ltot(:,:) = collider%Ltot(:,:) * collider%Lscale + collider%ke_orbit(:) = collider%ke_orbit(:) * collider%Escale + collider%ke_spin(:) = collider%ke_spin(:) * collider%Escale + collider%pe(:) = collider%pe(:) * collider%Escale + collider%Etot(:) = collider%Etot(:) * collider%Escale - collision_merge%mscale = 1.0_DP - collision_merge%dscale = 1.0_DP - collision_merge%vscale = 1.0_DP - collision_merge%tscale = 1.0_DP - collision_merge%Lscale = 1.0_DP - collision_merge%Escale = 1.0_DP + collider%mscale = 1.0_DP + collider%dscale = 1.0_DP + collider%vscale = 1.0_DP + collider%tscale = 1.0_DP + collider%Lscale = 1.0_DP + collider%Escale = 1.0_DP end associate call ieee_set_halting_mode(IEEE_ALL,fpe_halting_modes) diff --git a/src/swiftest/swiftest_module.f90 b/src/swiftest/swiftest_module.f90 index 82e233bb3..d4f88412b 100644 --- a/src/swiftest/swiftest_module.f90 +++ b/src/swiftest/swiftest_module.f90 @@ -322,7 +322,7 @@ module swiftest class(encounter_list), allocatable :: plpl_encounter !! List of massive body-massive body encounters in a single step class(collision_list_plpl), allocatable :: plpl_collision !! List of massive body-massive body collisions in a single step class(collision_list_plpl), allocatable :: pltp_collision !! List of massive body-massive body collisions in a single step - class(collision_merge), allocatable :: collider !! Collision system object + class(collision_basic), allocatable :: collider !! Collision system object class(encounter_storage(nframes=:)), allocatable :: encounter_history !! Stores encounter history for later retrieval and saving to file class(collision_storage(nframes=:)), allocatable :: collision_history !! Stores encounter history for later retrieval and saving to file diff --git a/src/swiftest/swiftest_util.f90 b/src/swiftest/swiftest_util.f90 index 32a137c50..eae7556e0 100644 --- a/src/swiftest/swiftest_util.f90 +++ b/src/swiftest/swiftest_util.f90 @@ -2653,7 +2653,7 @@ module subroutine swiftest_util_setup_construct_system(nbody_system, param) select case(param%collision_model) case("MERGE") - allocate(collision_merge :: nbody_system%collider) + allocate(collision_basic :: nbody_system%collider) case("BOUNCE") allocate(collision_bounce :: nbody_system%collider) case("SIMPLE")