diff --git a/examples/Fragmentation/Fragmentation_Movie.py b/examples/Fragmentation/Fragmentation_Movie.py index 88e275d1d..43822bf25 100644 --- a/examples/Fragmentation/Fragmentation_Movie.py +++ b/examples/Fragmentation/Fragmentation_Movie.py @@ -54,7 +54,7 @@ "supercatastrophic_off_axis": [np.array([0.0, 6.28, 0.0]), np.array([0.5, -6.28, 0.0])], "hitandrun" : [np.array([0.0, 6.28, 0.0]), - np.array([-0.1, -6.28, 0.0])] + np.array([-0.9, -6.28, 0.0])] } rot_vectors = {"disruption_headon" : [np.array([0.0, 0.0, 0.0]), @@ -67,7 +67,7 @@ body_Gmass = {"disruption_headon" : [1e-7, 1e-10], "supercatastrophic_off_axis": [1e-7, 1e-8], - "hitandrun" : [1e-7, 7e-10] + "hitandrun" : [1e-7, 1e-10] } density = 3000 * swiftest.AU2M**3 / swiftest.MSun @@ -139,7 +139,7 @@ def setup_plot(self): scale_frame = abs(rhy1) + abs(rhy2) ax = plt.Axes(fig, [0.1, 0.1, 0.8, 0.8]) - self.ax_pt_size = self.figsize[0] * 0.8 * 72 / (np.sqrt(2)*scale_frame) + self.ax_pt_size = self.figsize[0] * 0.7 * 72 / scale_frame ax.set_xlim(-scale_frame, scale_frame) ax.set_ylim(-scale_frame, scale_frame) ax.set_xticks([]) @@ -203,8 +203,8 @@ def data_stream(self, frame=0): # Set fragmentation parameters minimum_fragment_gmass = 0.2 * body_Gmass[style][1] # Make the minimum fragment mass a fraction of the smallest body gmtiny = 0.99 * body_Gmass[style][1] # Make GMTINY just smaller than the smallest original body. This will prevent runaway collisional cascades - sim.set_parameter(fragmentation=True, encounter_save="trajectory", gmtiny=gmtiny, minimum_fragment_gmass=minimum_fragment_gmass, verbose=False) - sim.run(dt=1e-3, tstop=1.0e-3, istep_out=1, dump_cadence=1) + sim.set_parameter(fragmentation=True, encounter_save="both", gmtiny=gmtiny, minimum_fragment_gmass=minimum_fragment_gmass, verbose=False) + sim.run(dt=1e-3, tstop=1.0e-3, istep_out=1, dump_cadence=0) print("Generating animation") anim = AnimatedScatter(sim,movie_filename,movie_titles[style],style,nskip=1) \ No newline at end of file diff --git a/python/swiftest/swiftest/simulation_class.py b/python/swiftest/swiftest/simulation_class.py index eea75c1ed..e517e7366 100644 --- a/python/swiftest/swiftest/simulation_class.py +++ b/python/swiftest/swiftest/simulation_class.py @@ -2713,7 +2713,8 @@ def read_output_file(self,read_init_cond : bool = True): param_tmp['BIN_OUT'] = os.path.join(self.simdir, self.param['BIN_OUT']) if self.codename == "Swiftest": self.data = io.swiftest2xr(param_tmp, verbose=self.verbose) - if self.verbose: print('Swiftest simulation data stored as xarray DataSet .data') + if self.verbose: + print('Swiftest simulation data stored as xarray DataSet .data') if read_init_cond: if self.verbose: print("Reading initial conditions file as .init_cond") @@ -2725,6 +2726,8 @@ def read_output_file(self,read_init_cond : bool = True): self.read_encounters() self.read_collisions() + if self.verbose: + print("Finished reading Swiftest dataset files.") elif self.codename == "Swifter": self.data = io.swifter2xr(param_tmp, verbose=self.verbose) diff --git a/src/fraggle/fraggle_io.f90 b/src/fraggle/fraggle_io.f90 index 6f75c7d01..61650e700 100644 --- a/src/fraggle/fraggle_io.f90 +++ b/src/fraggle/fraggle_io.f90 @@ -98,23 +98,24 @@ module subroutine fraggle_io_initialize_output(self, param) call check( nf90_def_var(nc%id, nc%rot_varname, nc%out_type,& [ nc%space_dimid, nc%name_dimid, nc%stage_dimid, nc%event_dimid], nc%rot_varid), "fraggle_io_initialize nf90_def_var rot_varid") + + if (param%lenergy) then - call check( nf90_def_var(nc%id, nc%ke_orb_varname, nc%out_type,& - [ nc%stage_dimid, nc%event_dimid], nc%KE_orb_varid), "fraggle_io_initialize_output nf90_def_var KE_orb_varid") - - call check( nf90_def_var(nc%id, nc%ke_spin_varname, nc%out_type,& - [ nc%stage_dimid, nc%event_dimid], nc%KE_spin_varid), "fraggle_io_initialize_output nf90_def_var KE_spin_varid" ) - - call check( nf90_def_var(nc%id, nc%pe_varname, nc%out_type,& - [ nc%stage_dimid, nc%event_dimid], nc%PE_varid), "fraggle_io_initialize_output nf90_def_var PE_varid" ) + call check( nf90_def_var(nc%id, nc%ke_orb_varname, nc%out_type,& + [ nc%stage_dimid, nc%event_dimid], nc%KE_orb_varid), "fraggle_io_initialize_output nf90_def_var KE_orb_varid") - call check( nf90_def_var(nc%id, nc%L_orb_varname, nc%out_type, & - [ nc%space_dimid, nc%stage_dimid, nc%event_dimid], nc%L_orb_varid), "fraggle_io_initialize_output nf90_def_var L_orb_varid" ) + call check( nf90_def_var(nc%id, nc%ke_spin_varname, nc%out_type,& + [ nc%stage_dimid, nc%event_dimid], nc%KE_spin_varid), "fraggle_io_initialize_output nf90_def_var KE_spin_varid" ) - call check( nf90_def_var(nc%id, nc%L_spin_varname, nc%out_type,& - [ nc%space_dimid, nc%stage_dimid, nc%event_dimid], nc%L_spin_varid), "fraggle_io_initialize_output nf90_def_var L_spin_varid" ) + call check( nf90_def_var(nc%id, nc%pe_varname, nc%out_type,& + [ nc%stage_dimid, nc%event_dimid], nc%PE_varid), "fraggle_io_initialize_output nf90_def_var PE_varid" ) + call check( nf90_def_var(nc%id, nc%L_orb_varname, nc%out_type, & + [ nc%space_dimid, nc%stage_dimid, nc%event_dimid], nc%L_orb_varid), "fraggle_io_initialize_output nf90_def_var L_orb_varid" ) + call check( nf90_def_var(nc%id, nc%L_spin_varname, nc%out_type,& + [ nc%space_dimid, nc%stage_dimid, nc%event_dimid], nc%L_spin_varid), "fraggle_io_initialize_output nf90_def_var L_spin_varid" ) + end if call check( nf90_inquire(nc%id, nVariables=nvar), "fraggle_io_initialize nf90_inquire nVariables" ) do varid = 1, nvar @@ -202,14 +203,18 @@ module subroutine fraggle_io_write_frame(self, nc, param) call check( nf90_put_var(nc%id, nc%rot_varid, pl%rot(:,i), start=[1, idslot, stage, eslot], count=[NDIM,1,1,1]), "fraggle_io_write_frame nf90_put_var rotx_varid" ) end do end do - call check( nf90_put_var(nc%id, nc%ke_orb_varid, fragments%ke_orbit_before, start=[ 1, eslot]), "fraggle_io_write_frame nf90_put_var ke_orb_varid before" ) - call check( nf90_put_var(nc%id, nc%ke_orb_varid, fragments%ke_orbit_after, start=[ 2, eslot]), "fraggle_io_write_frame nf90_put_var ke_orb_varid after" ) - call check( nf90_put_var(nc%id, nc%pe_varid, fragments%pe_before, start=[ 1, eslot]), "fraggle_io_write_frame nf90_put_var pe_varid before" ) - call check( nf90_put_var(nc%id, nc%pe_varid, fragments%pe_after, start=[ 2, eslot]), "fraggle_io_write_frame nf90_put_var pe_varid after" ) - call check( nf90_put_var(nc%id, nc%L_orb_varid, fragments%Lorbit_before(:), start=[1, 1, eslot], count=[NDIM, 1, 1]), "fraggle_io_write_frame nf90_put_var L_orb_varid before" ) - call check( nf90_put_var(nc%id, nc%L_orb_varid, fragments%Lorbit_after(:), start=[1, 2, eslot], count=[NDIM, 1, 1]), "fraggle_io_write_frame nf90_put_var L_orb_varid after" ) - call check( nf90_put_var(nc%id, nc%L_spin_varid, fragments%Lspin_before(:), start=[1, 1, eslot], count=[NDIM, 1, 1]), "fraggle_io_write_frame nf90_put_var L_spin_varid before" ) - call check( nf90_put_var(nc%id, nc%L_spin_varid, fragments%Lspin_after(:), start=[1, 2, eslot], count=[NDIM, 1, 1]), "fraggle_io_write_frame nf90_put_var L_spin_varid after" ) + if (param%lenergy) then + call check( nf90_put_var(nc%id, nc%ke_orb_varid, fragments%ke_orbit_before, start=[ 1, eslot]), "fraggle_io_write_frame nf90_put_var ke_orb_varid before" ) + call check( nf90_put_var(nc%id, nc%ke_orb_varid, fragments%ke_orbit_after, start=[ 2, eslot]), "fraggle_io_write_frame nf90_put_var ke_orb_varid after" ) + call check( nf90_put_var(nc%id, nc%ke_spin_varid, fragments%ke_spin_before, start=[ 1, eslot]), "fraggle_io_write_frame nf90_put_var ke_spin_varid before" ) + call check( nf90_put_var(nc%id, nc%ke_spin_varid, fragments%ke_spin_after, start=[ 2, eslot]), "fraggle_io_write_frame nf90_put_var ke_spin_varid after" ) + call check( nf90_put_var(nc%id, nc%pe_varid, fragments%pe_before, start=[ 1, eslot]), "fraggle_io_write_frame nf90_put_var pe_varid before" ) + call check( nf90_put_var(nc%id, nc%pe_varid, fragments%pe_after, start=[ 2, eslot]), "fraggle_io_write_frame nf90_put_var pe_varid after" ) + call check( nf90_put_var(nc%id, nc%L_orb_varid, fragments%Lorbit_before(:), start=[1, 1, eslot], count=[NDIM, 1, 1]), "fraggle_io_write_frame nf90_put_var L_orb_varid before" ) + call check( nf90_put_var(nc%id, nc%L_orb_varid, fragments%Lorbit_after(:), start=[1, 2, eslot], count=[NDIM, 1, 1]), "fraggle_io_write_frame nf90_put_var L_orb_varid after" ) + call check( nf90_put_var(nc%id, nc%L_spin_varid, fragments%Lspin_before(:), start=[1, 1, eslot], count=[NDIM, 1, 1]), "fraggle_io_write_frame nf90_put_var L_spin_varid before" ) + call check( nf90_put_var(nc%id, nc%L_spin_varid, fragments%Lspin_after(:), start=[1, 2, eslot], count=[NDIM, 1, 1]), "fraggle_io_write_frame nf90_put_var L_spin_varid after" ) + end if call check( nf90_set_fill(nc%id, old_mode, old_mode) ) end associate @@ -219,24 +224,6 @@ module subroutine fraggle_io_write_frame(self, nc, param) end subroutine fraggle_io_write_frame - module subroutine fraggle_io_log_pl(pl, param) - !! author: David A. Minton - !! - !! Writes a single message to the fraggle log file - implicit none - ! Arguments - class(swiftest_pl), intent(in) :: pl !! Swiftest massive body object (only the new bodies generated in a collision) - class(swiftest_parameters), intent(in) :: param !! Current swiftest run configuration parameters - ! Internals - integer(I4B) :: i - character(STRMAX) :: errmsg - - return - 667 continue - write(*,*) "Error writing Fraggle message to log file: " // trim(adjustl(errmsg)) - end subroutine fraggle_io_log_pl - - module subroutine fraggle_io_log_regime(colliders, frag) !! author: David A. Minton !! diff --git a/src/fraggle/fraggle_set.f90 b/src/fraggle/fraggle_set.f90 index 4a70130b6..a8e61130d 100644 --- a/src/fraggle/fraggle_set.f90 +++ b/src/fraggle/fraggle_set.f90 @@ -287,7 +287,6 @@ module subroutine fraggle_set_original_scale_factors(self, colliders) do i = 1, 2 colliders%rot(:,i) = colliders%L_spin(:,i) * (colliders%mass(i) * colliders%radius(i)**2 * colliders%Ip(3, i)) end do - frag%Qloss = frag%Qloss * frag%Escale frag%mtot = frag%mtot * frag%mscale frag%mass = frag%mass * frag%mscale @@ -300,6 +299,24 @@ module subroutine fraggle_set_original_scale_factors(self, colliders) frag%rb(:, i) = frag%x_coll(:, i) + frag%rbcom(:) frag%vb(:, i) = frag%v_coll(:, i) + frag%vbcom(:) end do + + frag%Qloss = frag%Qloss * frag%Escale + + frag%Lorbit_before(:) = frag%Lorbit_before * frag%Lscale + frag%Lspin_before(:) = frag%Lspin_before * frag%Lscale + frag%Ltot_before(:) = frag%Ltot_before * frag%Lscale + frag%ke_orbit_before = frag%ke_orbit_before * frag%Escale + frag%ke_spin_before = frag%ke_spin_before * frag%Escale + frag%pe_before = frag%pe_before * frag%Escale + frag%Etot_before = frag%Etot_before * frag%Escale + + frag%Lorbit_after(:) = frag%Lorbit_after * frag%Lscale + frag%Lspin_after(:) = frag%Lspin_after * frag%Lscale + frag%Ltot_after(:) = frag%Ltot_after * frag%Lscale + frag%ke_orbit_after = frag%ke_orbit_after * frag%Escale + frag%ke_spin_after = frag%ke_spin_after * frag%Escale + frag%pe_after = frag%pe_after * frag%Escale + frag%Etot_after = frag%Etot_after * frag%Escale frag%mscale = 1.0_DP frag%dscale = 1.0_DP diff --git a/src/modules/fraggle_classes.f90 b/src/modules/fraggle_classes.f90 index 8c75a3fc6..8f6cd4310 100644 --- a/src/modules/fraggle_classes.f90 +++ b/src/modules/fraggle_classes.f90 @@ -86,12 +86,12 @@ module fraggle_classes real(DP) :: Etot_before, Etot_after !! Before/after total system energy ! 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 !! Distance dimension scale factor - real(DP) :: mscale !! Mass scale factor - real(DP) :: tscale !! Time scale factor - real(DP) :: vscale !! Velocity scale factor (a convenience unit that is derived from dscale and tscale) - real(DP) :: Escale !! Energy scale factor (a convenience unit that is derived from dscale, tscale, and mscale) - real(DP) :: Lscale !! Angular momentum scale factor (a convenience unit that is derived from dscale, tscale, and mscale) + real(DP) :: dscale = 1.0_DP !! Distance dimension scale factor + real(DP) :: mscale = 1.0_DP !! Mass scale factor + real(DP) :: tscale = 1.0_DP !! Time scale factor + real(DP) :: vscale = 1.0_DP !! Velocity scale factor (a convenience unit that is derived from dscale and tscale) + 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_fragments => fraggle_generate_fragments !! Generates a system of fragments in barycentric coordinates that conserves energy and momentum. procedure :: accel => fraggle_placeholder_accel !! Placeholder subroutine to fulfill requirement for an accel method @@ -166,12 +166,6 @@ module subroutine fraggle_io_write_frame(self, nc, param) class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters end subroutine fraggle_io_write_frame - module subroutine fraggle_io_log_pl(pl, param) - implicit none - class(swiftest_pl), intent(in) :: pl !! Swiftest massive body object (only the new bodies generated in a collision) - class(swiftest_parameters), intent(in) :: param !! Current swiftest run configuration parameters - end subroutine fraggle_io_log_pl - module subroutine fraggle_io_log_regime(colliders, frag) implicit none class(fraggle_colliders), intent(in) :: colliders diff --git a/src/modules/symba_classes.f90 b/src/modules/symba_classes.f90 index 1a2c2ef8e..507b16745 100644 --- a/src/modules/symba_classes.f90 +++ b/src/modules/symba_classes.f90 @@ -176,8 +176,6 @@ module symba_classes type, extends(symba_encounter) :: symba_plplenc contains procedure :: extract_collisions => symba_collision_extract_collisions_from_encounters !! Processes the pl-pl encounter list remove only those encounters that led to a collision - procedure :: resolve_fragmentations => symba_resolve_collision_fragmentations !! Process list of collisions, determine the collisional regime, and then create fragments - procedure :: resolve_mergers => symba_resolve_collision_mergers !! Process list of collisions and merge colliding bodies together procedure :: resolve_collision => symba_resolve_collision_plplenc !! Process the pl-pl collision list, then modifiy the massive bodies based on the outcome of the c end type symba_plplenc @@ -232,22 +230,6 @@ module subroutine symba_collision_make_colliders_pl(self,idx) integer(I4B), dimension(2), intent(in) :: idx !! Array holding the indices of the two bodies involved in the collision end subroutine symba_collision_make_colliders_pl - module subroutine symba_resolve_collision_fragmentations(self, system, param, t) - implicit none - class(symba_plplenc), intent(inout) :: self !! SyMBA pl-pl encounter list - class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object - class(symba_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA additions - real(DP), intent(in) :: t !! Time of collision - end subroutine symba_resolve_collision_fragmentations - - module subroutine symba_resolve_collision_mergers(self, system, param, t) - implicit none - class(symba_plplenc), intent(inout) :: self !! SyMBA pl-pl encounter list - class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object - class(symba_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA additions - real(DP), intent(in) :: t !! Time of collision - end subroutine symba_resolve_collision_mergers - module subroutine symba_resolve_collision_plplenc(self, system, param, t, dt, irec) implicit none class(symba_plplenc), intent(inout) :: self !! SyMBA pl-pl encounter list diff --git a/src/symba/symba_collision.f90 b/src/symba/symba_collision.f90 index 0353cc7d5..c6ad292c4 100644 --- a/src/symba/symba_collision.f90 +++ b/src/symba/symba_collision.f90 @@ -185,6 +185,10 @@ module function symba_collision_casemerge(system, param, t) result(status) class is (symba_pl) call fragments%set_mass_dist(colliders, param) + + ! Calculate the initial energy of the system without the collisional family + call fragments%get_energy_and_momentum(colliders, system, param, lbefore=.true.) + ibiggest = colliders%idx(maxloc(pl%Gmass(colliders%idx(:)), dim=1)) fragments%id(1) = pl%id(ibiggest) fragments%rb(:,1) = fragments%rbcom(:) @@ -197,18 +201,16 @@ module function symba_collision_casemerge(system, param, t) result(status) ! Assume prinicpal axis rotation on 3rd Ip axis fragments%rot(:,1) = L_spin_new(:) / (fragments%Ip(3,1) * fragments%mass(1) * fragments%radius(1)**2) else ! If spin is not enabled, we will consider the lost pre-collision angular momentum as "escaped" and add it to our bookkeeping variable - param%Lescape(:) = param%Lescape(:) + colliders%L_orbit(:,1) + colliders%L_orbit(:,2) + system%Lescape(:) = system%Lescape(:) + colliders%L_orbit(:,1) + colliders%L_orbit(:,2) end if ! Keep track of the component of potential energy due to the pre-impact colliders%idx for book-keeping - pe = 0.0_DP - do j = 1, colliders%ncoll - do i = j + 1, colliders%ncoll - pe = pe - pl%Gmass(i) * pl%mass(j) / norm2(pl%rb(:, i) - pl%rb(:, j)) - end do - end do - system%Ecollisions = system%Ecollisions + pe - system%Euntracked = system%Euntracked - pe + ! Get the energy of the system after the collision + call fragments%get_energy_and_momentum(colliders, system, param, lbefore=.false.) + pe = fragments%pe_after - fragments%pe_before + system%Ecollisions = system%Ecollisions - pe + system%Euntracked = system%Euntracked + pe + ! Update any encounter lists that have the removed bodies in them so that they instead point to the new do k = 1, system%plplenc_list%nenc @@ -843,8 +845,14 @@ subroutine symba_collision_mergeaddsub(system, param, t, status) plnew%levelg(1:nfrag) = pl%levelg(ibiggest) plnew%levelm(1:nfrag) = pl%levelm(ibiggest) + plnew%lmtiny(1:nfrag) = plnew%Gmass(1:nfrag) < param%GMTINY + where(plnew%lmtiny(1:nfrag)) + plnew%info(1:nfrag)%particle_type = PL_TINY_TYPE_NAME + elsewhere + plnew%info(1:nfrag)%particle_type = PL_TYPE_NAME + end where + ! Log the properties of the new bodies - call fraggle_io_log_pl(plnew, param) allocate(system%fragments%pl, source=plnew) ! Append the new merged body to the list @@ -884,24 +892,24 @@ subroutine symba_collision_mergeaddsub(system, param, t, status) end subroutine symba_collision_mergeaddsub - module subroutine symba_resolve_collision_fragmentations(self, system, param, t) + subroutine symba_resolve_collision(plplcollision_list , system, param, t) !! author: David A. Minton !! !! Process list of collisions, determine the collisional regime, and then create fragments. !! implicit none ! Arguments - class(symba_plplenc), intent(inout) :: self !! SyMBA pl-pl encounter list - class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object - class(symba_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA additions - real(DP), intent(in) :: t !! Time of collision + class(symba_plplenc), intent(inout) :: plplcollision_list !! SyMBA pl-pl encounter list + class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object + class(symba_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA additions + real(DP), intent(in) :: t !! Time of collision ! Internals ! Internals integer(I4B), dimension(2) :: idx_parent !! Index of the two bodies considered the "parents" of the collision logical :: lgoodcollision integer(I4B) :: i - associate(plplcollision_list => self, ncollisions => self%nenc, idx1 => self%index1, idx2 => self%index2, collision_history => param%collision_history) + associate(ncollisions => plplcollision_list%nenc, idx1 => plplcollision_list%index1, idx2 => plplcollision_list%index2, collision_history => param%collision_history) select type(pl => system%pl) class is (symba_pl) select type (cb => system%cb) @@ -914,7 +922,19 @@ module subroutine symba_resolve_collision_fragmentations(self, system, param, t) lgoodcollision = symba_collision_consolidate_colliders(pl, cb, param, idx_parent, system%colliders) if ((.not. lgoodcollision) .or. any(pl%status(idx_parent(:)) /= COLLISION)) cycle - call system%colliders%regime(system%fragments, system, param) + if (param%lfragmentation) then + call system%colliders%regime(system%fragments, system, param) + else + associate(fragments => system%fragments, colliders => system%colliders) + fragments%regime = COLLRESOLVE_REGIME_MERGE + fragments%mtot = sum(colliders%mass(:)) + fragments%mass_dist(1) = fragments%mtot + fragments%mass_dist(2) = 0.0_DP + fragments%mass_dist(3) = 0.0_DP + fragments%rbcom(:) = (colliders%mass(1) * colliders%rb(:,1) + colliders%mass(2) * colliders%rb(:,2)) / fragments%mtot + fragments%vbcom(:) = (colliders%mass(1) * colliders%vb(:,1) + colliders%mass(2) * colliders%vb(:,2)) / fragments%mtot + end associate + end if if (param%lenc_save_trajectory) call collision_history%take_snapshot(param,system, t, "before") select case (system%fragments%regime) @@ -936,52 +956,7 @@ module subroutine symba_resolve_collision_fragmentations(self, system, param, t) end associate return - end subroutine symba_resolve_collision_fragmentations - - - module subroutine symba_resolve_collision_mergers(self, system, param, t) - !! author: David A. Minton - !! - !! Process list of collisions and merge colliding bodies together. - !! - implicit none - ! Arguments - class(symba_plplenc), intent(inout) :: self !! SyMBA pl-pl encounter list - class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object - class(symba_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA additions - real(DP), intent(in) :: t !! Time of collision - ! Internals - integer(I4B), dimension(2) :: idx_parent !! Index of the two bodies considered the "parents" of the collision - logical :: lgoodcollision - integer(I4B) :: i - - associate(plplcollision_list => self, ncollisions => self%nenc, idx1 => self%index1, idx2 => self%index2, fragments => system%fragments, colliders => system%colliders) - select type(pl => system%pl) - class is (symba_pl) - select type(cb => system%cb) - class is (symba_cb) - do i = 1, ncollisions - idx_parent(1) = pl%kin(idx1(i))%parent - idx_parent(2) = pl%kin(idx2(i))%parent - lgoodcollision = symba_collision_consolidate_colliders(pl, cb, param, idx_parent, colliders) - if (.not. lgoodcollision) cycle - if (any(pl%status(idx_parent(:)) /= COLLISION)) cycle ! One of these two bodies has already been resolved - - fragments%regime = COLLRESOLVE_REGIME_MERGE - fragments%mtot = sum(colliders%mass(:)) - fragments%mass_dist(1) = fragments%mtot - fragments%mass_dist(2) = 0.0_DP - fragments%mass_dist(3) = 0.0_DP - fragments%rbcom(:) = (colliders%mass(1) * colliders%rb(:,1) + colliders%mass(2) * colliders%rb(:,2)) / fragments%mtot - fragments%vbcom(:) = (colliders%mass(1) * colliders%vb(:,1) + colliders%mass(2) * colliders%vb(:,2)) / fragments%mtot - plplcollision_list%status(i) = symba_collision_casemerge(system, param, t) - end do - end select - end select - end associate - - return - end subroutine symba_resolve_collision_mergers + end subroutine symba_resolve_collision module subroutine symba_resolve_collision_plplenc(self, system, param, t, dt, irec) @@ -1029,11 +1004,8 @@ module subroutine symba_resolve_collision_plplenc(self, system, param, t, dt, ir call io_log_one_message(FRAGGLE_LOG_OUT, "***********************************************************" // & "***********************************************************") allocate(tmp_param, source=param) - if (param%lfragmentation) then - call plplcollision_list%resolve_fragmentations(system, param, t) - else - call plplcollision_list%resolve_mergers(system, param, t) - end if + + call symba_resolve_collision(plplcollision_list, system, param, t) ! Destroy the collision list now that the collisions are resolved call plplcollision_list%setup(0_I8B)