diff --git a/examples/Fragmentation/Fragmentation_Movie.py b/examples/Fragmentation/Fragmentation_Movie.py index a068205fb..84488ca14 100644 --- a/examples/Fragmentation/Fragmentation_Movie.py +++ b/examples/Fragmentation/Fragmentation_Movie.py @@ -41,7 +41,7 @@ movie_titles = dict(zip(available_movie_styles, movie_title_list)) # These initial conditions were generated by trial and error -pos_vectors = {"disruption_headon" : [np.array([1.0, -5.0e-05, 0.0]), +pos_vectors = {"disruption_headon" : [np.array([1.0005, -5.0e-05, 0.0]), np.array([1.0, 5.0e-05 ,0.0])], "supercatastrophic_off_axis": [np.array([1.0, -4.2e-05, 0.0]), np.array([1.0, 4.2e-05, 0.0])], @@ -202,7 +202,7 @@ def data_stream(self, frame=0): # Set fragmentation parameters minimum_fragment_gmass = 0.2 * body_Gmass[style][1] # Make the minimum fragment mass a fraction of the smallest body gmtiny = 0.99 * body_Gmass[style][1] # Make GMTINY just smaller than the smallest original body. This will prevent runaway collisional cascades - sim.set_parameter(fragmentation=True, encounter_save="TRAJECTORY", gmtiny=gmtiny, minimum_fragment_gmass=minimum_fragment_gmass, verbose=False) + sim.set_parameter(fragmentation=True, encounter_save="both", gmtiny=gmtiny, minimum_fragment_gmass=minimum_fragment_gmass, verbose=False) sim.run(dt=1e-4, tstop=1.0e-3, istep_out=1, dump_cadence=1) print("Generating animation") diff --git a/python/swiftest/swiftest/simulation_class.py b/python/swiftest/swiftest/simulation_class.py index f5517678f..9edf289b3 100644 --- a/python/swiftest/swiftest/simulation_class.py +++ b/python/swiftest/swiftest/simulation_class.py @@ -214,10 +214,12 @@ def __init__(self,read_param: bool = False, read_old_output: bool = False, simdi Check for close encounters between bodies. If set to True, then the radii of massive bodies must be included in initial conditions. Parameter input file equivalent: `CHK_CLOSE` - encounter_save : {"NONE","TRAJECTORY","CLOSEST"}, default "NONE" - Indicate if and how encounter data should be saved. If set to "TRAJECTORY" the full close encounter - trajectories are saved to file. If set to "CLOSEST" only the trajectories at the time of closest approach - are saved. If set to "NONE" no trajectory information is saved. + encounter_save : {"NONE","TRAJECTORY","CLOSEST", "BOTH"}, default "NONE" + Indicate if and how encounter data should be saved. If set to "TRAJECTORY", the position and velocity vectors + of all bodies undergoing close encounters are saved at each intermediate step to the encounter files. + If set to "CLOSEST", the position and velocities at the point of closest approach between pairs of bodies are + computed and stored to the encounter files. If set to "BOTH", then this stores the values that would be computed + in "TRAJECTORY" and "CLOSEST". If set to "NONE" no trajectory information is saved. *WARNING*: Enabling this feature could lead to very large files. general_relativity : bool, default True Include the post-Newtonian correction in acceleration calculations. @@ -1023,7 +1025,7 @@ def set_feature(self, tides: bool | None = None, interaction_loops: Literal["TRIANGULAR", "FLAT", "ADAPTIVE"] | None = None, encounter_check_loops: Literal["TRIANGULAR", "SORTSWEEP", "ADAPTIVE"] | None = None, - encounter_save: Literal["NONE", "TRAJECTORY", "CLOSEST"] | None = None, + encounter_save: Literal["NONE", "TRAJECTORY", "CLOSEST", "BOTH"] | None = None, verbose: bool | None = None, **kwargs: Any ): @@ -1035,10 +1037,12 @@ def set_feature(self, close_encounter_check : bool, optional Check for close encounters between bodies. If set to True, then the radii of massive bodies must be included in initial conditions. - encounter_save : {"NONE","TRAJECTORY","CLOSEST"}, default "NONE" - Indicate if and how encounter data should be saved. If set to "TRAJECTORY" the full close encounter - trajectories are saved to file. If set to "CLOSEST" only the trajectories at the time of closest approach - are saved. If set to "NONE" no trajectory information is saved. + encounter_save : {"NONE","TRAJECTORY","CLOSEST","BOTH"}, default "NONE" + Indicate if and how encounter data should be saved. If set to "TRAJECTORY", the position and velocity vectors + of all bodies undergoing close encounters are saved at each intermediate step to the encounter files. + If set to "CLOSEST", the position and velocities at the point of closest approach between pairs of bodies are + computed and stored to the encounter files. If set to "BOTH", then this stores the values that would be computed + in "TRAJECTORY" and "CLOSEST". If set to "NONE" no trajectory information is saved. *WARNING*: Enabling this feature could lead to very large files. general_relativity : bool, optional Include the post-Newtonian correction in acceleration calculations. @@ -1199,7 +1203,7 @@ def set_feature(self, update_list.append("encounter_check_loops") if encounter_save is not None: - valid_vals = ["NONE", "TRAJECTORY", "CLOSEST"] + valid_vals = ["NONE", "TRAJECTORY", "CLOSEST", "BOTH"] encounter_save = encounter_save.upper() if encounter_save not in valid_vals: msg = f"{encounter_save} is not a valid option for encounter_save." diff --git a/src/encounter/encounter_util.f90 b/src/encounter/encounter_util.f90 index 5d9e93a2e..33689dba2 100644 --- a/src/encounter/encounter_util.f90 +++ b/src/encounter/encounter_util.f90 @@ -570,9 +570,9 @@ module subroutine encounter_util_snapshot_encounter(self, param, system, t, arg) character(*), intent(in), optional :: arg !! Optional argument (needed for extended storage type used in collision snapshots) ! Arguments class(encounter_snapshot), allocatable :: snapshot - integer(I4B) :: i, j, k, npl_snap, ntp_snap + integer(I4B) :: i, j, k, npl_snap, ntp_snap, iflag real(DP), dimension(NDIM) :: rrel, vrel, rcom, vcom - real(DP) :: mu, a, q, capm, tperi + real(DP) :: Gmtot, a, q, capm, tperi real(DP), dimension(NDIM,2) :: rb,vb if (.not.present(t)) then @@ -621,6 +621,8 @@ module subroutine encounter_util_snapshot_encounter(self, param, system, t, arg) ntp_snap = count(tp%lmask(1:ntp)) end if + if (npl_snap + ntp_snap == 0) return ! Nothing to snapshot + pl_snap%nbody = npl_snap ! Take snapshot of the currently encountering massive bodies @@ -674,7 +676,7 @@ module subroutine encounter_util_snapshot_encounter(self, param, system, t, arg) pl_snap%id(:) = pl%id([i,j]) pl_snap%info(:) = pl%info([i,j]) pl_snap%Gmass(:) = pl%Gmass([i,j]) - mu = sum(pl_snap%Gmass(:)) + Gmtot = sum(pl_snap%Gmass(:)) if (param%lclose) pl_snap%radius(:) = pl%radius([i,j]) if (param%lrotation) then do i = 1, NDIM @@ -686,15 +688,36 @@ module subroutine encounter_util_snapshot_encounter(self, param, system, t, arg) ! Compute pericenter passage time to get the closest approach parameters rrel(:) = plplenc_list%r2(:,k) - plplenc_list%r1(:,k) vrel(:) = plplenc_list%v2(:,k) - plplenc_list%v1(:,k) - call orbel_xv2aqt(mu, rrel(1), rrel(2), rrel(3), vrel(1), vrel(2), vrel(3), a, q, capm, tperi) + call orbel_xv2aqt(Gmtot, rrel(1), rrel(2), rrel(3), vrel(1), vrel(2), vrel(3), a, q, capm, tperi) snapshot%t = t + tperi ! Computer the center mass of the pair - - ! do i = 1, NDIM - ! pl_snap%rh(i,:) = pl%rh(i,[i,j]) - ! pl_snap%vh(i,:) = pl%vb(i,1:npl), pl%lmask(1:npl)) - ! end do + rcom(:) = (plplenc_list%r1(:,k) * pl_snap%Gmass(1) + plplenc_list%r2(:,k) * pl_snap%Gmass(2)) / Gmtot + vcom(:) = (plplenc_list%v1(:,k) * pl_snap%Gmass(1) + plplenc_list%v2(:,k) * pl_snap%Gmass(2)) / Gmtot + rb(:,1) = plplenc_list%r1(:,k) - rcom(:) + rb(:,2) = plplenc_list%r2(:,k) - rcom(:) + vb(:,1) = plplenc_list%v1(:,k) - vcom(:) + vb(:,2) = plplenc_list%v2(:,k) - vcom(:) + + ! Drift the relative orbit to get the new relative position and velocity + call drift_one(Gmtot, rrel(1), rrel(2), rrel(3), vrel(1), vrel(2), vrel(3), tperi, iflag) + if (iflag /= 0) write(*,*) "Danby error in encounter_util_snapshot_encounter. Closest approach positions and vectors may not be accurate." + + ! Get the new position and velocity vectors + rb(:,1) = -(pl_snap%Gmass(2) / Gmtot) * rrel(:) + rb(:,2) = (pl_snap%Gmass(1)) / Gmtot * rrel(:) + + vb(:,1) = -(pl_snap%Gmass(2) / Gmtot) * vrel(:) + vb(:,2) = (pl_snap%Gmass(1)) / Gmtot * vrel(:) + + ! Move the CoM assuming constant velocity over the time it takes to reach periapsis + !rcom(:) = rcom(:) + vcom(:) * tperi + + ! Compute the heliocentric position and velocity vector at periapsis + pl_snap%rh(:,1) = rb(:,1) + rcom(:) + pl_snap%rh(:,2) = rb(:,2) + rcom(:) + pl_snap%vh(:,1) = vb(:,1) + vcom(:) + pl_snap%vh(:,2) = vb(:,2) + vcom(:) call pl_snap%sort("id", ascending=.true.) call encounter_util_save_encounter(param%encounter_history,snapshot,snapshot%t) diff --git a/src/symba/symba_io.f90 b/src/symba/symba_io.f90 index 916fb1a4c..29f1c1fbe 100644 --- a/src/symba/symba_io.f90 +++ b/src/symba/symba_io.f90 @@ -118,15 +118,18 @@ module subroutine symba_io_param_reader(self, unit, iotype, v_list, iostat, ioms ! All reporting of collision information in SyMBA (including mergers) is now recorded in the Fraggle logfile call io_log_start(param, FRAGGLE_LOG_OUT, "Fraggle logfile") - if ((param%encounter_save /= "NONE") .and. (param%encounter_save /= "TRAJECTORY") .and. (param%encounter_save /= "CLOSEST")) then + if ((param%encounter_save /= "NONE") .and. & + (param%encounter_save /= "TRAJECTORY") .and. & + (param%encounter_save /= "CLOSEST") .and. & + (param%encounter_save /= "BOTH")) then write(iomsg,*) 'Invalid encounter_save parameter: ',trim(adjustl(param%out_type)) - write(iomsg,*) 'Valid options are NONE, TRAJECTORY, or CLOSEST' + write(iomsg,*) 'Valid options are NONE, TRAJECTORY, CLOSEST, or BOTH' iostat = -1 return end if - param%lenc_save_trajectory = (param%encounter_save == "TRAJECTORY") - param%lenc_save_closest = (param%encounter_save == "CLOSEST") .or. param%lenc_save_trajectory ! Closest approaches are always saved when trajectories are saved + param%lenc_save_trajectory = (param%encounter_save == "TRAJECTORY") .or. (param%encounter_save == "BOTH") + param%lenc_save_closest = (param%encounter_save == "CLOSEST") .or. (param%encounter_save == "BOTH") ! Call the base method (which also prints the contents to screen) call io_param_reader(param, unit, iotype, v_list, iostat, iomsg)