From da12dc98adb755bb6baaacac8cadcfacf74a9ba4 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Sat, 17 Dec 2022 10:29:09 -0500 Subject: [PATCH] More cleanup. It compiles again! --- src/CMakeLists.txt | 2 +- .../collision_regime.f90} | 25 ++++---- src/collision/collision_setup.f90 | 12 ++-- src/collision/collision_util.f90 | 19 +++--- src/fraggle/fraggle_set.f90 | 4 +- src/fraggle/fraggle_util.f90 | 6 +- src/modules/collision_classes.f90 | 49 ++++++++++----- src/modules/symba_classes.f90 | 5 +- src/setup/setup.f90 | 11 ++++ src/symba/symba_collision.f90 | 63 +++++++++---------- 10 files changed, 116 insertions(+), 80 deletions(-) rename src/{fraggle/fraggle_regime.f90 => collision/collision_regime.f90} (95%) diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 90ace7141..95ee7e9a0 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -27,6 +27,7 @@ SET(FAST_MATH_FILES ${SRC}/modules/whm_classes.f90 ${SRC}/modules/swiftest.f90 ${SRC}/collision/collision_io.f90 + ${SRC}/collision/collision_regime.f90 ${SRC}/collision/collision_setup.f90 ${SRC}/collision/collision_util.f90 ${SRC}/discard/discard.f90 @@ -37,7 +38,6 @@ SET(FAST_MATH_FILES ${SRC}/encounter/encounter_io.f90 ${SRC}/fraggle/fraggle_generate.f90 ${SRC}/fraggle/fraggle_io.f90 - ${SRC}/fraggle/fraggle_regime.f90 ${SRC}/fraggle/fraggle_set.f90 ${SRC}/fraggle/fraggle_setup.f90 ${SRC}/fraggle/fraggle_util.f90 diff --git a/src/fraggle/fraggle_regime.f90 b/src/collision/collision_regime.f90 similarity index 95% rename from src/fraggle/fraggle_regime.f90 rename to src/collision/collision_regime.f90 index 09015c66c..89fdeb269 100644 --- a/src/fraggle/fraggle_regime.f90 +++ b/src/collision/collision_regime.f90 @@ -7,20 +7,19 @@ !! You should have received a copy of the GNU General Public License along with Swiftest. !! If not, see: https://www.gnu.org/licenses. -submodule(fraggle_classes) s_encounter_regime +submodule(collision_classes) s_collision_regime use swiftest contains - module subroutine encounter_regime_impactors(self, fragments, system, param) + module subroutine collision_regime_impactors(self, system, param) !! Author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton !! !! Determine which fragmentation regime the set of impactors will be. This subroutine is a wrapper for the non-polymorphic raggle_regime_collresolve subroutine. !! It converts to SI units prior to calling implicit none ! Arguments - class(collision_impactors), intent(inout) :: self !! Fraggle impactors object - class(fraggle_fragments), intent(inout) :: fragments !! Fraggle fragment system object + class(collision_impactors), intent(inout) :: self !! Collision system impactors object class(swiftest_nbody_system), intent(in) :: system !! Swiftest nbody system object class(swiftest_parameters), intent(in) :: param !! Current Swiftest run configuration parameters ! Internals @@ -58,7 +57,7 @@ module subroutine encounter_regime_impactors(self, fragments, system, param) dentot = sum(mass_si(:) * density_si(:)) / mtot !! Use the positions and velocities of the parents from indside the step (at collision) to calculate the collisional regime - call encounter_regime_collresolve(Mcb_si, mass_si(jtarg), mass_si(jproj), radius_si(jtarg), radius_si(jproj), & + call collision_regime_collresolve(Mcb_si, mass_si(jtarg), mass_si(jproj), radius_si(jtarg), radius_si(jproj), & x1_si(:), x2_si(:), v1_si(:), v2_si(:), density_si(jtarg), density_si(jproj), & min_mfrag_si, impactors%regime, mlr, mslr, impactors%Qloss) @@ -67,9 +66,9 @@ module subroutine encounter_regime_impactors(self, fragments, system, param) impactors%mass_dist(3) = min(max(mtot - mlr - mslr, 0.0_DP), mtot) ! Find the center of mass of the collisional system - fragments%mtot = sum(impactors%mass(:)) - impactors%rbcom(:) = (impactors%mass(1) * impactors%rb(:,1) + impactors%mass(2) * impactors%rb(:,2)) / fragments%mtot - impactors%vbcom(:) = (impactors%mass(1) * impactors%vb(:,1) + impactors%mass(2) * impactors%vb(:,2)) / fragments%mtot + mtot = sum(impactors%mass(:)) + impactors%rbcom(:) = (impactors%mass(1) * impactors%rb(:,1) + impactors%mass(2) * impactors%rb(:,2)) / mtot + impactors%vbcom(:) = (impactors%mass(1) * impactors%vb(:,1) + impactors%mass(2) * impactors%vb(:,2)) / mtot ! Find the point of impact between the two bodies runit(:) = impactors%rb(:,2) - impactors%rb(:,1) @@ -80,14 +79,14 @@ module subroutine encounter_regime_impactors(self, fragments, system, param) impactors%mass_dist(:) = (impactors%mass_dist(:) / param%MU2KG) impactors%Qloss = impactors%Qloss * (param%TU2S / param%DU2M)**2 / param%MU2KG - call fraggle_io_log_regime(impactors, fragments) + !call fraggle_io_log_regime(impactors, fragments) end associate return - end subroutine encounter_regime_impactors + end subroutine collision_regime_impactors - subroutine encounter_regime_collresolve(Mcb, m1, m2, rad1, rad2, rh1, rh2, vb1, vb2, den1, den2, min_mfrag, & + subroutine collision_regime_collresolve(Mcb, m1, m2, rad1, rad2, rh1, rh2, vb1, vb2, den1, den2, min_mfrag, & regime, Mlr, Mslr, Qloss) !! Author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton !! @@ -379,6 +378,6 @@ function calc_c_star(Rc1) result(c_star) return end function calc_c_star - end subroutine encounter_regime_collresolve + end subroutine collision_regime_collresolve -end submodule s_encounter_regime \ No newline at end of file +end submodule s_collision_regime \ No newline at end of file diff --git a/src/collision/collision_setup.f90 b/src/collision/collision_setup.f90 index 287d74109..d98883ca4 100644 --- a/src/collision/collision_setup.f90 +++ b/src/collision/collision_setup.f90 @@ -11,21 +11,25 @@ use swiftest contains - module subroutine collision_setup_system(self, system, param) + module subroutine collision_setup_system(self, param) !! author: David A. Minton !! !! Initializer for the encounter collision system. Allocates the collider and fragments classes and the before/after snapshots implicit none ! Arguments - class(collision_system), intent(inout) :: self !! Encounter collision system object - class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object - class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters + class(collision_system), intent(inout) :: self !! Encounter collision system object + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters ! Internals + ! TODO: Check parameter file for fragmentation model in SyMBA + allocate(collision_impactors :: self%impactors) + allocate(fraggle_fragments :: self%fragments) + return end subroutine collision_setup_system + module subroutine collision_setup_fragments(self, n, param) !! author: David A. Minton !! diff --git a/src/collision/collision_util.f90 b/src/collision/collision_util.f90 index 0c435cfc6..819f9df76 100644 --- a/src/collision/collision_util.f90 +++ b/src/collision/collision_util.f90 @@ -81,8 +81,12 @@ module subroutine collision_util_final_system(self) implicit none ! Arguments type(collision_system), intent(inout) :: self !! Collision system object + ! Internals + type(swiftest_parameters) :: tmp_param - call self%reset() + call self%reset(tmp_param) + if (allocated(self%impactors)) deallocate(self%impactors) + if (allocated(self%fragments)) deallocate(self%fragments) return end subroutine collision_util_final_system @@ -298,17 +302,15 @@ module subroutine collision_util_reset_impactors(self) return end subroutine collision_util_reset_impactors - - module subroutine collision_util_reset_system(self) + module subroutine collision_util_reset_system(self, param) !! author: David A. Minton !! !! Resets the collider system and deallocates all allocatables implicit none ! Arguments - class(collision_system), intent(inout) :: self + class(collision_system), intent(inout) :: self !! Collision system object + class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters - if (allocated(self%impactors)) deallocate(self%impactors) - if (allocated(self%fragments)) deallocate(self%fragments) if (allocated(self%before)) deallocate(self%before) if (allocated(self%after)) deallocate(self%after) @@ -320,12 +322,13 @@ module subroutine collision_util_reset_system(self) self%pe(:) = 0.0_DP self%Etot(:) = 0.0_DP + call self%impactors%reset() + call self%fragments%setup(0, param) + return end subroutine collision_util_reset_system - - module subroutine collision_util_set_coordinate_system(self) !! author: David A. Minton !! diff --git a/src/fraggle/fraggle_set.f90 b/src/fraggle/fraggle_set.f90 index 20770a1d0..7c405db90 100644 --- a/src/fraggle/fraggle_set.f90 +++ b/src/fraggle/fraggle_set.f90 @@ -77,7 +77,7 @@ module subroutine fraggle_set_mass_dist(self, param) 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 encounter_regime. + ! 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 @@ -117,7 +117,7 @@ module subroutine fraggle_set_mass_dist(self, param) write(*,*) "fraggle_set_mass_dist_fragments error: Unrecognized regime code",impactors%regime end select - ! Make the first two bins the same as the Mlr and Mslr values that came from encounter_regime + ! 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) diff --git a/src/fraggle/fraggle_util.f90 b/src/fraggle/fraggle_util.f90 index d550b7166..c1e2fccb5 100644 --- a/src/fraggle/fraggle_util.f90 +++ b/src/fraggle/fraggle_util.f90 @@ -185,8 +185,12 @@ module subroutine fraggle_util_final_system(self) implicit none ! Arguments type(fraggle_system), intent(inout) :: self !! Collision impactors storage object + ! Internals + type(swiftest_parameters) :: tmp_param - call self%reset() + call self%reset(tmp_param) + if (allocated(self%impactors)) deallocate(self%impactors) + if (allocated(self%fragments)) deallocate(self%fragments) return end subroutine fraggle_util_final_system diff --git a/src/modules/collision_classes.f90 b/src/modules/collision_classes.f90 index 262920269..5c31ab592 100644 --- a/src/modules/collision_classes.f90 +++ b/src/modules/collision_classes.f90 @@ -55,6 +55,7 @@ module collision_classes real(DP), dimension(NDIM) :: rbimp !! Impact point position vector of the collider system in system barycentric coordinates contains + procedure :: get_regime => collision_regime_impactors !! Determine which fragmentation regime the set of impactors will be procedure :: set_coordinate_system => collision_set_coordinate_impactors !! Defines the collisional coordinate system, including the unit vectors of both the system and individual fragments. procedure :: setup => collision_setup_impactors !! Allocates arrays for n fragments in a fragment system. Passing n = 0 deallocates all arrays. procedure :: reset => collision_util_reset_impactors !! Resets the collider object variables to 0 and deallocates the index and mass distributions @@ -102,7 +103,8 @@ module collision_classes real(DP), dimension(2) :: pe !! Before/after potential energy real(DP), dimension(2) :: Etot !! Before/after total system energy contains - procedure :: regime => collision_regime_system !! Determine which fragmentation regime the set of impactors will be + procedure :: generate_fragments => abstract_generate_fragments !! Generates a system of fragments + procedure :: set_mass_dist => abstract_set_mass_dist !! Sets the distribution of mass among the fragments depending on the regime type procedure :: setup => collision_setup_system !! Initializer for the encounter collision system. Allocates the collider and fragments classes and the before/after snapshots procedure :: get_energy_and_momentum => collision_util_get_energy_momentum !! Calculates total 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 @@ -110,6 +112,24 @@ module collision_classes final :: collision_util_final_system !! Finalizer will deallocate all allocatables end type collision_system + abstract interface + subroutine abstract_generate_fragments(self, system, param, lfailure) + import collision_system, swiftest_nbody_system, swiftest_parameters + implicit none + class(collision_system), intent(inout) :: self !! Fraggle fragment system object + class(swiftest_nbody_system), intent(inout) :: 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 abstract_generate_fragments + + subroutine abstract_set_mass_dist(self, param) + import collision_system, swiftest_parameters + implicit none + class(collision_system), intent(inout) :: self !! Collision system object + class(swiftest_parameters), intent(in) :: param !! Current Swiftest run configuration parameters + end subroutine abstract_set_mass_dist + end interface + !! NetCDF dimension and variable names for the enounter save object type, extends(encounter_io_parameters) :: collision_io_parameters integer(I4B) :: stage_dimid !! ID for the stage dimension @@ -201,12 +221,12 @@ module subroutine collision_util_placeholder_step(self, system, param, t, dt) real(DP), intent(in) :: dt !! Stepsiz end subroutine collision_util_placeholder_step - module subroutine collision_regime_system(self, system, param) + module subroutine collision_regime_impactors(self, system, param) implicit none - class(collision_system), intent(inout) :: self !! Collision system object - class(swiftest_nbody_system), intent(in) :: system !! Swiftest nbody system object - class(swiftest_parameters), intent(in) :: param !! Current Swiftest run configuration parameters - end subroutine collision_regime_system + class(collision_impactors), intent(inout) :: self !! Collision system impactors object + class(swiftest_nbody_system), intent(in) :: system !! Swiftest nbody system object + class(swiftest_parameters), intent(in) :: param !! Current Swiftest run configuration parameters + end subroutine collision_regime_impactors module subroutine collision_set_coordinate_impactors(self) implicit none @@ -223,7 +243,6 @@ module subroutine collision_util_set_coordinate_system(self) class(collision_system), intent(inout) :: self !! Collisional system end subroutine collision_util_set_coordinate_system - module subroutine collision_setup_fragments(self, n, param) implicit none class(collision_fragments), intent(inout) :: self !! Fragment system object @@ -238,12 +257,11 @@ module subroutine collision_setup_impactors(self, system, param) class(swiftest_parameters), intent(in) :: param !! Current swiftest run configuration parameters end subroutine collision_setup_impactors - module subroutine collision_setup_system(self, system, param) - use swiftest_classes, only : swiftest_nbody_system, swiftest_parameters + module subroutine collision_setup_system(self, param) + use swiftest_classes, only : swiftest_parameters implicit none - class(collision_system), intent(inout) :: self !! Encounter collision system object - class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object - class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters + class(collision_system), intent(inout) :: self !! Collision system object + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters end subroutine collision_setup_system module subroutine collision_util_dealloc_fragments(self) @@ -293,12 +311,13 @@ end subroutine collision_util_index_map module subroutine collision_util_reset_impactors(self) implicit none - class(collision_impactors), intent(inout) :: self + class(collision_impactors), intent(inout) :: self !! Collision system object end subroutine collision_util_reset_impactors - module subroutine collision_util_reset_system(self) + module subroutine collision_util_reset_system(self, param) implicit none - class(collision_system), intent(inout) :: self + class(collision_system), intent(inout) :: self !! Collision system object + class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters end subroutine collision_util_reset_system module subroutine collision_util_snapshot(self, param, system, t, arg) diff --git a/src/modules/symba_classes.f90 b/src/modules/symba_classes.f90 index 98143ed08..86988fc9e 100644 --- a/src/modules/symba_classes.f90 +++ b/src/modules/symba_classes.f90 @@ -17,7 +17,7 @@ module symba_classes use helio_classes, only : helio_cb, helio_pl, helio_tp, helio_nbody_system use fraggle_classes, only : collision_impactors, fraggle_fragments use encounter_classes, only : encounter_list, encounter_storage - use collision_classes, only : collision_storage + use collision_classes, only : collision_storage, collision_system implicit none public @@ -190,8 +190,7 @@ module symba_classes class(symba_plplenc), allocatable :: plplenc_list !! List of massive body-massive body encounters in a single step class(symba_plplenc), allocatable :: plplcollision_list !! List of massive body-massive body collisions in a single step integer(I4B) :: irec !! System recursion level - class(collision_impactors), allocatable :: impactors !! Fraggle impactors object - class(fraggle_fragments), allocatable :: fragments !! Fraggle fragmentation system object + class(collision_system), allocatable :: collision_system !! Collision system object contains procedure :: write_discard => symba_io_write_discard !! Write out information about discarded and merged planets and test particles in SyMBA procedure :: initialize => symba_setup_initialize_system !! Performs SyMBA-specific initilization steps diff --git a/src/setup/setup.f90 b/src/setup/setup.f90 index 353f60fb2..9b27416ae 100644 --- a/src/setup/setup.f90 +++ b/src/setup/setup.f90 @@ -75,6 +75,15 @@ module subroutine setup_construct_system(system, param) select type(param) class is (symba_parameters) + + if (param%lfragmentation) then + allocate(fraggle_system :: system%collision_system) + call system%collision_system%setup(param) + + allocate(symba_nbody_system :: system%collision_system%before) + allocate(symba_nbody_system :: system%collision_system%after) + end if + if (param%lenc_save_trajectory .or. param%lenc_save_closest) then allocate(encounter_storage :: param%encounter_history) associate (encounter_history => param%encounter_history) @@ -98,6 +107,8 @@ module subroutine setup_construct_system(system, param) end associate end select + + end select case (RINGMOONS) write(*,*) 'RINGMOONS-SyMBA integrator not yet enabled' diff --git a/src/symba/symba_collision.f90 b/src/symba/symba_collision.f90 index 663e40b55..3af2c7e99 100644 --- a/src/symba/symba_collision.f90 +++ b/src/symba/symba_collision.f90 @@ -30,8 +30,7 @@ module function symba_collision_casedisruption(system, param, t) result(status) character(len=STRMAX) :: message real(DP) :: dpe - associate(impactors => system%impactors, fragments => system%fragments) - + associate(collision_system => system%collision_system, impactors => system%collision_system%impactors, fragments => system%collision_system%fragments) select case(impactors%regime) case(COLLRESOLVE_REGIME_DISRUPTION) message = "Disruption between" @@ -42,12 +41,12 @@ module function symba_collision_casedisruption(system, param, t) result(status) call io_log_one_message(FRAGGLE_LOG_OUT, message) ! Collisional fragments will be uniformly distributed around the pre-impact barycenter - call fragments%set_mass_dist(impactors, param) + call collision_system%set_mass_dist(param) ! Generate the position and velocity distributions of the fragments - call fragments%generate_fragments(impactors, system, param, lfailure) + call collision_system%generate_fragments(system, param, lfailure) - dpe = fragments%pe_after - fragments%pe_before + dpe = collision_system%pe(2) - collision_system%pe(1) system%Ecollisions = system%Ecollisions - dpe system%Euntracked = system%Euntracked + dpe @@ -105,7 +104,7 @@ module function symba_collision_casehitandrun(system, param, t) result(status) character(len=STRMAX) :: message real(DP) :: dpe - associate(impactors => system%impactors, fragments => system%fragments) + associate(collision_system => system%collision_system, impactors => system%collision_system%impactors, fragments => system%collision_system%fragments) message = "Hit and run between" call symba_collision_collider_message(system%pl, impactors%idx, message) call io_log_one_message(FRAGGLE_LOG_OUT, trim(adjustl(message))) @@ -124,12 +123,12 @@ module function symba_collision_casehitandrun(system, param, t) result(status) lpure = .true. else ! Imperfect hit and run, so we'll keep the largest body and destroy the other lpure = .false. - call fragments%set_mass_dist(impactors, param) + call collision_system%set_mass_dist(param) ! Generate the position and velocity distributions of the fragments - call fragments%generate_fragments(impactors, system, param, lpure) + call collision_system%generate_fragments(system, param, lpure) - dpe = fragments%pe_after - fragments%pe_before + dpe = collision_system%pe(2) - collision_system%pe(1) system%Ecollisions = system%Ecollisions - dpe system%Euntracked = system%Euntracked + dpe @@ -150,7 +149,7 @@ module function symba_collision_casehitandrun(system, param, t) result(status) pl%ldiscard(impactors%idx(:)) = .false. pl%lcollision(impactors%idx(:)) = .false. end select - allocate(system%fragments%pl, source=system%impactors%pl) ! Be sure to save the pl so that snapshots still work + allocate(collision_system%after%pl, source=collision_system%before%pl) ! Be sure to save the pl so that snapshots still work else ibiggest = impactors%idx(maxloc(system%pl%Gmass(impactors%idx(:)), dim=1)) fragments%id(1) = system%pl%id(ibiggest) @@ -187,7 +186,7 @@ module function symba_collision_casemerge(system, param, t) result(status) real(DP) :: dpe character(len=STRMAX) :: message - associate(impactors => system%impactors, fragments => system%fragments) + associate(collision_system => system%collision_system, impactors => system%collision_system%impactors, fragments => system%collision_system%fragments) message = "Merging" call symba_collision_collider_message(system%pl, impactors%idx, message) call io_log_one_message(FRAGGLE_LOG_OUT, message) @@ -195,10 +194,10 @@ module function symba_collision_casemerge(system, param, t) result(status) select type(pl => system%pl) class is (symba_pl) - call fragments%set_mass_dist(impactors, param) + call collision_system%set_mass_dist(param) ! Calculate the initial energy of the system without the collisional family - call fragments%get_energy_and_momentum(impactors, system, param, lbefore=.true.) + call collision_system%get_energy_and_momentum(system, param, lbefore=.true.) ibiggest = impactors%idx(maxloc(pl%Gmass(impactors%idx(:)), dim=1)) fragments%id(1) = pl%id(ibiggest) @@ -217,8 +216,8 @@ module function symba_collision_casemerge(system, param, t) result(status) ! Keep track of the component of potential energy due to the pre-impact impactors%idx for book-keeping ! Get the energy of the system after the collision - call fragments%get_energy_and_momentum(impactors, system, param, lbefore=.false.) - dpe = fragments%pe_after - fragments%pe_before + call collision_system%get_energy_and_momentum(system, param, lbefore=.false.) + dpe = collision_system%pe(2) - collision_system%pe(1) system%Ecollisions = system%Ecollisions - dpe system%Euntracked = system%Euntracked + dpe @@ -757,7 +756,8 @@ subroutine symba_collision_mergeaddsub(system, param, t, status) class is (symba_pl) select type(pl_discards => system%pl_discards) class is (symba_merger) - associate(info => pl%info, pl_adds => system%pl_adds, cb => system%cb, npl => pl%nbody, impactors => system%impactors, fragments => system%fragments) + associate(info => pl%info, pl_adds => system%pl_adds, cb => system%cb, npl => pl%nbody, & + collision_system => system%collision_system, impactors => system%collision_system%impactors, fragments => system%collision_system%fragments) ! Add the impactors%idx bodies to the subtraction list nimpactors = impactors%ncoll nfrag = fragments%nbody @@ -864,7 +864,7 @@ subroutine symba_collision_mergeaddsub(system, param, t, status) end where ! Log the properties of the new bodies - allocate(system%fragments%pl, source=plnew) + allocate(collision_system%after%pl, source=plnew) ! Append the new merged body to the list nstart = pl_adds%nbody + 1 @@ -920,35 +920,32 @@ subroutine symba_resolve_collision(plplcollision_list , system, param, t) logical :: lgoodcollision integer(I4B) :: i - associate(ncollisions => plplcollision_list%nenc, idx1 => plplcollision_list%index1, idx2 => plplcollision_list%index2, collision_history => param%collision_history) + associate(ncollisions => plplcollision_list%nenc, idx1 => plplcollision_list%index1, idx2 => plplcollision_list%index2, collision_history => param%collision_history, & + collision_system => system%collision_system, impactors => system%collision_system%impactors, fragments => system%collision_system%fragments) select type(pl => system%pl) class is (symba_pl) select type (cb => system%cb) class is (symba_cb) do i = 1, ncollisions - allocate(collision_impactors :: system%impactors) - allocate(fraggle_fragments :: system%fragments) idx_parent(1) = pl%kin(idx1(i))%parent idx_parent(2) = pl%kin(idx2(i))%parent - lgoodcollision = symba_collision_consolidate_impactors(pl, cb, param, idx_parent, system%impactors) + lgoodcollision = symba_collision_consolidate_impactors(pl, cb, param, idx_parent, impactors) if ((.not. lgoodcollision) .or. any(pl%status(idx_parent(:)) /= COLLISION)) cycle if (param%lfragmentation) then - call system%impactors%regime(system%fragments, system, param) + call impactors%get_regime(system, param) else - associate(fragments => system%fragments, impactors => system%impactors) - impactors%regime = COLLRESOLVE_REGIME_MERGE - fragments%mtot = sum(impactors%mass(:)) - impactors%mass_dist(1) = fragments%mtot - impactors%mass_dist(2) = 0.0_DP - impactors%mass_dist(3) = 0.0_DP - impactors%rbcom(:) = (impactors%mass(1) * impactors%rb(:,1) + impactors%mass(2) * impactors%rb(:,2)) / fragments%mtot - impactors%vbcom(:) = (impactors%mass(1) * impactors%vb(:,1) + impactors%mass(2) * impactors%vb(:,2)) / fragments%mtot - end associate + impactors%regime = COLLRESOLVE_REGIME_MERGE + fragments%mtot = sum(impactors%mass(:)) + impactors%mass_dist(1) = fragments%mtot + impactors%mass_dist(2) = 0.0_DP + impactors%mass_dist(3) = 0.0_DP + impactors%rbcom(:) = (impactors%mass(1) * impactors%rb(:,1) + impactors%mass(2) * impactors%rb(:,2)) / fragments%mtot + impactors%vbcom(:) = (impactors%mass(1) * impactors%vb(:,1) + impactors%mass(2) * impactors%vb(:,2)) / fragments%mtot end if call collision_history%take_snapshot(param,system, t, "before") - select case (system%impactors%regime) + select case (impactors%regime) case (COLLRESOLVE_REGIME_DISRUPTION, COLLRESOLVE_REGIME_SUPERCATASTROPHIC) plplcollision_list%status(i) = symba_collision_casedisruption(system, param, t) case (COLLRESOLVE_REGIME_HIT_AND_RUN) @@ -960,7 +957,7 @@ subroutine symba_resolve_collision(plplcollision_list , system, param, t) call util_exit(FAILURE) end select call collision_history%take_snapshot(param,system, t, "after") - deallocate(system%impactors,system%fragments) + call impactors%reset() end do end select end select