From 52ecbefd61460c837178dddd2d5754f6424da7b8 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Tue, 13 Dec 2022 08:18:05 -0500 Subject: [PATCH] Restructuring so that the saving collisions is not a user option (collision info should always be saved when collisions are turned on) --- examples/Basic_Simulation/output_reader.py | 2 +- examples/Fragmentation/Fragmentation_Movie.py | 2 +- python/swiftest/swiftest/io.py | 5 +- python/swiftest/swiftest/simulation_class.py | 50 +++---------------- src/io/io.f90 | 2 +- src/modules/symba_classes.f90 | 10 ++-- src/symba/symba_collision.f90 | 19 ++++--- src/symba/symba_io.f90 | 12 +---- src/symba/symba_step.f90 | 6 +-- 9 files changed, 32 insertions(+), 76 deletions(-) diff --git a/examples/Basic_Simulation/output_reader.py b/examples/Basic_Simulation/output_reader.py index a41103ccd..fc332af0c 100644 --- a/examples/Basic_Simulation/output_reader.py +++ b/examples/Basic_Simulation/output_reader.py @@ -28,7 +28,7 @@ import matplotlib.pyplot as plt # Read in the simulation output and store it as an Xarray dataset. -sim = swiftest.Simulation(read_old_output_file=True) +sim = swiftest.Simulation(read_old_output=True) # Plot of the data and save the output plot. colors = ['white' if x == 'Massive Body' else 'black' for x in sim.data['particle_type']] diff --git a/examples/Fragmentation/Fragmentation_Movie.py b/examples/Fragmentation/Fragmentation_Movie.py index 090a4f4a3..a068205fb 100644 --- a/examples/Fragmentation/Fragmentation_Movie.py +++ b/examples/Fragmentation/Fragmentation_Movie.py @@ -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, collision_save="TRAJECTORY", encounter_save="TRAJECTORY", gmtiny=gmtiny, minimum_fragment_gmass=minimum_fragment_gmass, verbose=False) + sim.set_parameter(fragmentation=True, encounter_save="TRAJECTORY", 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/io.py b/python/swiftest/swiftest/io.py index be70cae50..029672ec7 100644 --- a/python/swiftest/swiftest/io.py +++ b/python/swiftest/swiftest/io.py @@ -32,8 +32,7 @@ "ENCOUNTER_CHECK", "TSTART", "DUMP_CADENCE", - "ENCOUNTER_SAVE", - "COLLISION_SAVE") + "ENCOUNTER_SAVE") @@ -55,7 +54,7 @@ float_param = ["T0", "TSTART", "TSTOP", "DT", "CHK_RMIN", "CHK_RMAX", "CHK_EJECT", "CHK_QMIN", "DU2M", "MU2KG", "TU2S", "MIN_GMFRAG", "GMTINY"] -upper_str_param = ["OUT_TYPE","OUT_FORM","OUT_STAT","IN_TYPE","IN_FORM","ENCOUNTER_SAVE","COLLISION_SAVE", "CHK_QMIN_COORD"] +upper_str_param = ["OUT_TYPE","OUT_FORM","OUT_STAT","IN_TYPE","IN_FORM","ENCOUNTER_SAVE", "CHK_QMIN_COORD"] lower_str_param = ["NC_IN", "PL_IN", "TP_IN", "CB_IN", "CHK_QMIN_RANGE"] param_keys = ['! VERSION'] + int_param + float_param + upper_str_param + lower_str_param+ bool_param diff --git a/python/swiftest/swiftest/simulation_class.py b/python/swiftest/swiftest/simulation_class.py index 1bedfeba4..a239b9293 100644 --- a/python/swiftest/swiftest/simulation_class.py +++ b/python/swiftest/swiftest/simulation_class.py @@ -43,7 +43,7 @@ class Simulation: This is a class that defines the basic Swift/Swifter/Swiftest simulation object """ - def __init__(self,read_param: bool = False, read_old_output_file: bool = False, simdir: os.PathLike | str = "simdata", **kwargs: Any): + def __init__(self,read_param: bool = False, read_old_output: bool = False, simdir: os.PathLike | str = "simdata", **kwargs: Any): """ Parameters @@ -65,7 +65,7 @@ def __init__(self,read_param: bool = False, read_old_output_file: bool = False, inside the current working directory, which can be changed by passing `param_file` as an argument. - The argument has an equivalent parameter or set of parameters in the parameter input file. 3. Default values (see below) - read_old_output_file : bool, default False + read_old_output : bool, default False If true, read in a pre-existing binary input file given by the argument `output_file_name` if it exists. Parameter input file equivalent: None simdir : PathLike, default `"simdir"` @@ -227,12 +227,6 @@ def __init__(self,read_param: bool = False, read_old_output_file: bool = False, If set to True, this turns on the Fraggle fragment generation code and `rotation` must also be True. This argument only applies to Swiftest-SyMBA simulations. It will be ignored otherwise. Parameter input file equivalent: `FRAGMENTATION` - collision_save : {"NONE","TRAJECTORY","CLOSEST"}, default "NONE" - Indicate if and how fragmentation data should be saved. If set to "TRAJECTORY" the full close encounter - trajectories associated with each collision are saved to file. If set to "CLOSEST" only the trajectories - at a the time the collision occurs are saved. If set to "NONE" no trajectory information is saved (collision - details are still logged fraggle.log). - *WARNING*: Enabling this feature could lead to very large files. minimum_fragment_gmass : float, optional If fragmentation is turned on, this sets the mimimum G*mass of a collisional fragment that can be generated. *Note.* Only set one of minimum_fragment_gmass or minimum_fragment_mass @@ -329,7 +323,7 @@ def __init__(self,read_param: bool = False, read_old_output_file: bool = False, msg += "\nDelete the file or change the location of param_file" raise NotADirectoryError(msg) else: - if read_old_output_file or read_param: + if read_old_output or read_param: raise NotADirectoryError(f"Cannot find directory {self.simdir.resolve()} ") else: self.simdir.mkdir(parents=True, exist_ok=False) @@ -354,8 +348,8 @@ def __init__(self,read_param: bool = False, read_old_output_file: bool = False, # If the user asks to read in an old parameter file or output file, override any default parameters with values from the file # If the file doesn't exist, flag it for now so we know to create it param_file_found = False - if read_param or read_old_output_file: - if self.read_param(read_init_cond = not read_old_output_file): + if read_param or read_old_output: + if self.read_param(read_init_cond = not read_old_output): # We will add the parameter file to the kwarg list. This will keep the set_parameter method from # overriding everything with defaults when there are no arguments passed to Simulation() kwargs['param_file'] = self.param_file @@ -375,7 +369,7 @@ def __init__(self,read_param: bool = False, read_old_output_file: bool = False, self.write_param() # Read in an old simulation file if requested - if read_old_output_file: + if read_old_output: binpath = os.path.join(self.simdir, self.param['BIN_OUT']) if os.path.exists(binpath): self.read_output_file() @@ -762,7 +756,7 @@ def set_parameter(self, verbose: bool = True, **kwargs): "init_cond_file_type": "NETCDF_DOUBLE", "init_cond_file_name": None, "init_cond_format": "EL", - "read_old_output_file": False, + "read_old_output": False, "output_file_type": "NETCDF_DOUBLE", "output_file_name": None, "output_format": "XVEL", @@ -794,8 +788,7 @@ def set_parameter(self, verbose: bool = True, **kwargs): "encounter_check_loops": "TRIANGULAR", "ephemeris_date": "MBCL", "restart": False, - "encounter_save" : "NONE", - "collision_save" : "NONE" + "encounter_save" : "NONE" } param_file = kwargs.pop("param_file",None) @@ -1032,7 +1025,6 @@ def set_feature(self, 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, - collision_save: Literal["NONE", "TRAJECTORY", "CLOSEST"] | None = None, verbose: bool | None = None, **kwargs: Any ): @@ -1054,12 +1046,6 @@ def set_feature(self, fragmentation : bool, optional If set to True, this turns on the Fraggle fragment generation code and `rotation` must also be True. This argument only applies to Swiftest-SyMBA simulations. It will be ignored otherwise. - collision_save : {"NONE","TRAJECTORY","CLOSEST"}, default "NONE" - Indicate if and how fragmentation data should be saved. If set to "TRAJECTORY" the full close encounter - trajectories associated with each collision are saved to file. If set to "CLOSEST" only the trajectories - at a the time the collision occurs are saved. If set to "NONE" no trajectory information is saved (collision - details are still logged fraggle.log). - *WARNING*: Enabling this feature could lead to very large files. minimum_fragment_gmass : float, optional If fragmentation is turned on, this sets the mimimum G*mass of a collisional fragment that can be generated. *Note.* Only set one of minimum_fragment_gmass or minimum_fragment_mass @@ -1226,20 +1212,6 @@ def set_feature(self, self.param["ENCOUNTER_SAVE"] = encounter_save update_list.append("encounter_save") - - if collision_save is not None: - collision_save = collision_save.upper() - valid_vals = ["NONE", "TRAJECTORY", "CLOSEST"] - if collision_save not in valid_vals: - msg = f"{collision_save} is not a valid option for collision_save." - msg += f"\nMust be one of {valid_vals}" - warnings.warn(msg,stacklevel=2) - if "COLLISION_SAVE" not in self.param: - self.param["COLLISION_SAVE"] = valid_vals[0] - else: - self.param["COLLISION_SAVE"] = collision_save - update_list.append("collision_save") - self.param["TIDES"] = False feature_dict = self.get_feature(update_list, verbose) @@ -1272,7 +1244,6 @@ def get_feature(self, arg_list: str | List[str] | None = None, verbose: bool | N valid_var = {"close_encounter_check": "CHK_CLOSE", "fragmentation": "FRAGMENTATION", "encounter_save": "ENCOUNTER_SAVE", - "collision_save": "COLLISION_SAVE", "minimum_fragment_gmass": "MIN_GMFRAG", "rotation": "ROTATION", "general_relativity": "GR", @@ -2740,11 +2711,6 @@ def read_output_file(self,read_init_cond : bool = True): else: read_encounters = False - if "COLLISION_SAVE" in self.param: - read_collisions = self.param["COLLISION_SAVE"] != "NONE" - else: - read_collisions = False - param_tmp = self.param.copy() param_tmp['BIN_OUT'] = os.path.join(self.simdir, self.param['BIN_OUT']) if self.codename == "Swiftest": diff --git a/src/io/io.f90 b/src/io/io.f90 index c58eb4dbb..f159e6ac7 100644 --- a/src/io/io.f90 +++ b/src/io/io.f90 @@ -682,7 +682,7 @@ module subroutine io_param_reader(self, unit, iotype, v_list, iostat, iomsg) param%lrestart = .true. end if ! Ignore SyMBA-specific, not-yet-implemented, or obsolete input parameters - case ("NPLMAX", "NTPMAX", "GMTINY", "MIN_GMFRAG", "FRAGMENTATION", "SEED", "YARKOVSKY", "YORP", "ENCOUNTER_SAVE", "COLLISION_SAVE") + case ("NPLMAX", "NTPMAX", "GMTINY", "MIN_GMFRAG", "FRAGMENTATION", "SEED", "YARKOVSKY", "YORP", "ENCOUNTER_SAVE") case default write(*,*) "Ignoring unknown parameter -> ",param_name end select diff --git a/src/modules/symba_classes.f90 b/src/modules/symba_classes.f90 index 47d9c9997..ef3297cde 100644 --- a/src/modules/symba_classes.f90 +++ b/src/modules/symba_classes.f90 @@ -31,8 +31,7 @@ module symba_classes integer(I4B), dimension(:), allocatable :: seed !! Random seeds logical :: lfragmentation = .false. !! Do fragmentation modeling instead of simple merger. character(STRMAX) :: encounter_save = "NONE" !! Indicate if and how encounter data should be saved - character(STRMAX) :: collision_save = "NONE" !! Indicate if and how fragmentation data should be saved - logical :: lencounter_save = .false. !! Turns on encounter saving + logical :: lencounter_save type(encounter_storage(nframes=:)), allocatable :: encounter_history !! Stores encounter history for later retrieval and saving to file type(collision_storage(nframes=:)), allocatable :: collision_history !! Stores encounter history for later retrieval and saving to file contains @@ -175,7 +174,7 @@ module symba_classes !> SyMBA class for tracking pl-pl close encounters in a step type, extends(symba_encounter) :: symba_plplenc contains - procedure :: extract_collisions => symba_collision_encounter_extract_collisions !! Processes the pl-pl encounter list remove only those encounters that led to a collision + 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 @@ -207,7 +206,7 @@ module symba_classes interface - module subroutine symba_collision_check_encounter(self, system, param, t, dt, irec, lany_collision, lany_closest) + module subroutine symba_collision_check_encounter(self, system, param, t, dt, irec, lany_collision) use swiftest_classes, only : swiftest_parameters implicit none class(symba_encounter), intent(inout) :: self !! SyMBA pl-tp encounter list object @@ -217,10 +216,9 @@ module subroutine symba_collision_check_encounter(self, system, param, t, dt, ir real(DP), intent(in) :: dt !! step size integer(I4B), intent(in) :: irec !! Current recursion level logical, intent(out) :: lany_collision !! Returns true if any pair of encounters resulted in a collision - logical, intent(out) :: lany_closest !! Returns true if any pair of encounters reached their closest approach without colliding end subroutine symba_collision_check_encounter - module subroutine symba_collision_encounter_extract_collisions(self, system, param) + module subroutine symba_collision_extract_collisions_from_encounters(self, system, param) implicit none class(symba_plplenc), intent(inout) :: self !! SyMBA pl-pl encounter list class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object diff --git a/src/symba/symba_collision.f90 b/src/symba/symba_collision.f90 index 8534237c2..96ab905a0 100644 --- a/src/symba/symba_collision.f90 +++ b/src/symba/symba_collision.f90 @@ -262,7 +262,7 @@ subroutine symba_collision_collider_message(pl, collidx, collider_message) end subroutine symba_collision_collider_message - module subroutine symba_collision_check_encounter(self, system, param, t, dt, irec, lany_collision, lany_closest) + module subroutine symba_collision_check_encounter(self, system, param, t, dt, irec, lany_collision) !! author: David A. Minton !! !! Check for merger between massive bodies and test particles in SyMBA @@ -278,7 +278,7 @@ module subroutine symba_collision_check_encounter(self, system, param, t, dt, ir real(DP), intent(in) :: t !! current time real(DP), intent(in) :: dt !! step size integer(I4B), intent(in) :: irec !! Current recursion level - logical, intent(out) :: lany_collision, lany_closest !! Returns true if cany pair of encounters resulted in a collision + logical, intent(out) :: lany_collision !! Returns true if cany pair of encounters resulted in a collision ! Internals logical, dimension(:), allocatable :: lcollision, lclosest, lmask real(DP), dimension(NDIM) :: xr, vr @@ -289,7 +289,6 @@ module subroutine symba_collision_check_encounter(self, system, param, t, dt, ir class(symba_encounter), allocatable :: tmp lany_collision = .false. - lany_closest = .false. if (self%nenc == 0) return select type(self) @@ -339,7 +338,6 @@ module subroutine symba_collision_check_encounter(self, system, param, t, dt, ir end if lany_collision = any(lcollision(:)) - lany_closest = any(lclosest(:)) if (lany_collision) then call pl%rh2rb(system%cb) ! Update the central body barycenteric position vector to get us out of DH and into bary @@ -391,6 +389,11 @@ module subroutine symba_collision_check_encounter(self, system, param, t, dt, ir end select end if + ! Take snapshots of pairs of bodies at close approach (but not collision) if requested + if (any(lclosest(:))) then + + end if + end select end select @@ -580,7 +583,7 @@ function symba_collision_consolidate_colliders(pl, cb, param, idx_parent, collid end function symba_collision_consolidate_colliders - module subroutine symba_collision_encounter_extract_collisions(self, system, param) + module subroutine symba_collision_extract_collisions_from_encounters(self, system, param) !! author: David A. Minton !! !! Processes the pl-pl encounter list remove only those encounters that led to a collision @@ -646,7 +649,7 @@ module subroutine symba_collision_encounter_extract_collisions(self, system, par end select return - end subroutine symba_collision_encounter_extract_collisions + end subroutine symba_collision_extract_collisions_from_encounters module subroutine symba_collision_make_colliders_pl(self, idx) @@ -984,7 +987,7 @@ module subroutine symba_resolve_collision_plplenc(self, system, param, t, dt, ir integer(I4B), intent(in) :: irec !! Current recursion level ! Internals real(DP) :: Eorbit_before, Eorbit_after - logical :: lplpl_collision, lplpl_closest + logical :: lplpl_collision character(len=STRMAX) :: timestr class(symba_parameters), allocatable :: tmp_param @@ -1037,7 +1040,7 @@ module subroutine symba_resolve_collision_plplenc(self, system, param, t, dt, ir deallocate(tmp_param) ! Check whether or not any of the particles that were just added are themselves in a collision state. This will generate a new plplcollision_list - call plplenc_list%collision_check(system, param, t, dt, irec, lplpl_collision, lplpl_closest) + call plplenc_list%collision_check(system, param, t, dt, irec, lplpl_collision) if (.not.lplpl_collision) exit end do diff --git a/src/symba/symba_io.f90 b/src/symba/symba_io.f90 index 4f19bfd30..27a455062 100644 --- a/src/symba/symba_io.f90 +++ b/src/symba/symba_io.f90 @@ -68,9 +68,6 @@ module subroutine symba_io_param_reader(self, unit, iotype, v_list, iostat, ioms case ("ENCOUNTER_SAVE") call io_toupper(param_value) read(param_value, *) param%encounter_save - case ("COLLISION_SAVE") - call io_toupper(param_value) - read(param_value, *) param%collision_save case("SEED") read(param_value, *) nseeds_from_file ! Because the number of seeds can vary between compilers/systems, we need to make sure we can handle cases in which the input file has a different @@ -128,14 +125,7 @@ module subroutine symba_io_param_reader(self, unit, iotype, v_list, iostat, ioms return end if - if ((param%collision_save /= "NONE") .and. (param%collision_save /= "TRAJECTORY") .and. (param%collision_save /= "CLOSEST")) then - write(iomsg,*) 'Invalid collision_save parameter: ',trim(adjustl(param%out_type)) - write(iomsg,*) 'Valid options are NONE, TRAJECTORY, or CLOSEST' - iostat = -1 - return - end if - param%lencounter_save = (param%encounter_save == "TRAJECTORY") .or. (param%encounter_save == "CLOSEST") .or. & - (param%collision_save == "TRAJECTORY") .or. (param%collision_save == "CLOSEST") + param%lencounter_save = (param%encounter_save == "TRAJECTORY") .or. (param%encounter_save == "CLOSEST") ! Call the base method (which also prints the contents to screen) call io_param_reader(param, unit, iotype, v_list, iostat, iomsg) diff --git a/src/symba/symba_step.f90 b/src/symba/symba_step.f90 index ba61844e5..eb37c718b 100644 --- a/src/symba/symba_step.f90 +++ b/src/symba/symba_step.f90 @@ -180,7 +180,7 @@ recursive module subroutine symba_step_recur_system(self, param, t, ireci) ! Internals integer(I4B) :: j, irecp, nloops real(DP) :: dtl, dth - logical :: lencounter, lplpl_closest, lpltp_closest + logical :: lencounter select type(param) class is (symba_parameters) @@ -241,8 +241,8 @@ recursive module subroutine symba_step_recur_system(self, param, t, ireci) end if if (param%lclose) then - call plplenc_list%collision_check(system, param, t+dtl, dtl, ireci, lplpl_collision, lplpl_closest) - call pltpenc_list%collision_check(system, param, t+dtl, dtl, ireci, lpltp_collision, lpltp_closest) + call plplenc_list%collision_check(system, param, t+dtl, dtl, ireci, lplpl_collision) + call pltpenc_list%collision_check(system, param, t+dtl, dtl, ireci, lpltp_collision) if (lplpl_collision) call plplenc_list%resolve_collision(system, param, t+dtl, dtl, ireci) if (lpltp_collision) call pltpenc_list%resolve_collision(system, param, t+dtl, dtl, ireci)