Added the structure necessary to snapshot collisions. Refactored frag…

…mentation_save to collision_save

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, fragmentation_save="TRAJECTORY", gmtiny=gmtiny, minimum_fragment_gmass=minimum_fragment_gmass, verbose=False)
+        sim.set_parameter(fragmentation=True, collision_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=0)
 
         print("Generating animation")

python/swiftest/swiftest/simulation_class.py

@@ -227,7 +227,7 @@ 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`
-        fragmentation_save : {"NONE","TRAJECTORY","CLOSEST"}, default "NONE"
+        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
@@ -794,7 +794,7 @@ def set_parameter(self, verbose: bool = True, **kwargs):
             "ephemeris_date": "MBCL",
             "restart": False,
             "encounter_save" : "NONE",
-            "fragmentation_save" : "NONE"
+            "collision_save" : "NONE"
         }
         param_file = kwargs.pop("param_file",None)
 
@@ -1031,7 +1031,7 @@ 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,
-                    fragmentation_save: Literal["NONE", "TRAJECTORY", "CLOSEST"] | None = None,
+                    collision_save: Literal["NONE", "TRAJECTORY", "CLOSEST"] | None = None,
                     verbose: bool | None = None,
                     **kwargs: Any
                     ):
@@ -1053,7 +1053,7 @@ 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.
-        fragmentation_save : {"NONE","TRAJECTORY","CLOSEST"}, default "NONE"
+        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
@@ -1226,18 +1226,18 @@ def set_feature(self,
                 update_list.append("encounter_save")
 
 
-        if fragmentation_save is not None:
-            fragmentation_save = fragmentation_save.upper()
+        if collision_save is not None:
+            collision_save = collision_save.upper()
             valid_vals = ["NONE", "TRAJECTORY", "CLOSEST"]
-            if fragmentation_save not in valid_vals:
-                msg = f"{fragmentation_save} is not a valid option for fragmentation_save."
+            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 "FRAGMENTATION_SAVE" not in self.param:
                     self.param["FRAGMENTATION_SAVE"] = valid_vals[0]
             else:
-                self.param["FRAGMENTATION_SAVE"] = fragmentation_save
-                update_list.append("fragmentation_save")
+                self.param["FRAGMENTATION_SAVE"] = collision_save
+                update_list.append("collision_save")
 
         self.param["TIDES"] = False
 
@@ -1271,7 +1271,7 @@ 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",
-                     "fragmentation_save": "FRAGMENTATION_SAVE",
+                     "collision_save": "FRAGMENTATION_SAVE",
                      "minimum_fragment_gmass": "MIN_GMFRAG",
                      "rotation": "ROTATION",
                      "general_relativity": "GR",

src/encounter/encounter_io.f90

@@ -28,7 +28,7 @@ module subroutine encounter_io_dump(self, param)
             param%ioutput = self%tslot(i)
 
             select type(snapshot => self%frame(i)%item)
-            class is (fraggle_encounter_snapshot)
+            class is (fraggle_collision_snapshot)
                call snapshot%write_frame(self%ncc,param)
                call snapshot%encounter_snapshot%write_frame(self%nce,param)
             class is (encounter_snapshot)

src/fraggle/fraggle_io.f90

@@ -154,7 +154,7 @@ module subroutine fraggle_io_write_frame(self, nc, param)
       use netcdf
       implicit none
       ! Arguments
-      class(fraggle_encounter_snapshot), intent(in)    :: self   !! Swiftest encounter structure
+      class(fraggle_collision_snapshot), intent(in)    :: self   !! Swiftest encounter structure
       class(encounter_io_parameters),    intent(inout) :: nc    !! Parameters used to identify a particular encounter io NetCDF dataset
       class(swiftest_parameters),        intent(inout) :: param  !! Current run configuration parameters
       ! Internals

src/fraggle/fraggle_util.f90

@@ -174,7 +174,7 @@ module subroutine fraggle_util_final_snapshot(self)
       !! Finalizer will deallocate all allocatables
       implicit none
       ! Arguments
-      type(fraggle_encounter_snapshot),  intent(inout) :: self !! Fraggle encountar storage object
+      type(fraggle_collision_snapshot),  intent(inout) :: self !! Fraggle encountar storage object
 
       call encounter_util_final_snapshot(self%encounter_snapshot)
 

src/modules/fraggle_classes.f90

@@ -130,14 +130,14 @@ module fraggle_classes
       procedure :: initialize => fraggle_io_initialize_output !! Initialize a set of parameters used to identify a NetCDF output object
    end type fraggle_io_parameters
 
-   type, extends(encounter_snapshot)  :: fraggle_encounter_snapshot
+   type, extends(encounter_snapshot)  :: fraggle_collision_snapshot
       logical                               :: lcollision !! Indicates that this snapshot contains at least one collision
       class(fraggle_colliders), allocatable :: colliders  !! Colliders object at this snapshot
       class(fraggle_fragments), allocatable :: fragments  !! Fragments object at this snapshot
    contains
       procedure :: write_frame => fraggle_io_write_frame !! Writes a frame of encounter data to file 
       final     ::                fraggle_util_final_snapshot
-   end type fraggle_encounter_snapshot
+   end type fraggle_collision_snapshot
 
    interface
       module subroutine fraggle_generate_fragments(self, colliders, system, param, lfailure)
@@ -158,7 +158,7 @@ end subroutine fraggle_io_initialize_output
 
       module subroutine fraggle_io_write_frame(self, nc, param)
          implicit none
-         class(fraggle_encounter_snapshot), intent(in)    :: self   !! Swiftest encounter structure
+         class(fraggle_collision_snapshot), intent(in)    :: self   !! Swiftest encounter structure
          class(encounter_io_parameters),    intent(inout) :: nc    !! Parameters used to identify a particular encounter io NetCDF dataset
          class(swiftest_parameters),        intent(inout) :: param  !! Current run configuration parameters
       end subroutine fraggle_io_write_frame
@@ -293,7 +293,7 @@ end subroutine fraggle_util_final_fragments
 
       module subroutine fraggle_util_final_snapshot(self)
          implicit none
-         type(fraggle_encounter_snapshot),  intent(inout) :: self !! Fraggle encountar storage object
+         type(fraggle_collision_snapshot),  intent(inout) :: self !! Fraggle encountar storage object
       end subroutine fraggle_util_final_snapshot
 
       module subroutine fraggle_util_get_energy_momentum(self, colliders, system, param, lbefore)

src/modules/symba_classes.f90

@@ -31,7 +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)                       :: fragmentation_save = "NONE"  !! Indicate if and how fragmentation 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
    contains
       procedure :: reader => symba_io_param_reader
@@ -200,7 +200,8 @@ module symba_classes
       procedure :: set_recur_levels => symba_step_set_recur_levels_system !! Sets recursion levels of bodies and encounter lists to the current system level
       procedure :: recursive_step   => symba_step_recur_system            !! Step interacting planets and active test particles ahead in democratic heliocentric coordinates at the current recursion level, if applicable, and descend to the next deeper level if necessary
       procedure :: reset            => symba_step_reset_system            !! Resets pl, tp,and encounter structures at the start of a new step 
-      procedure :: snapshot         => symba_util_take_encounter_snapshot !! Take a minimal snapshot of the system through an encounter
+      procedure :: encounter_snap   => symba_util_take_encounter_snapshot !! Take a minimal snapshot of the system through an encounter
+      procedure :: collision_snap   => symba_util_take_collision_snapshot !! Take a minimal snapshot of the system before and after a collision
       procedure :: dump_encounter   => symba_io_dump_encounter            !! Saves the encounter and/or fragmentation data to file(s)   
       final     ::                     symba_util_final_system            !! Finalizes the SyMBA nbody system object - deallocates all allocatables
    end type symba_nbody_system
@@ -371,6 +372,14 @@ module subroutine symba_util_set_renc(self, scale)
          integer(I4B),    intent(in)    :: scale !! Current recursion depth
       end subroutine symba_util_set_renc
 
+      module subroutine symba_util_take_collision_snapshot(self, param, t, stage)
+         implicit none
+         class(symba_nbody_system),  intent(inout) :: self  !! SyMBA nbody system object
+         class(swiftest_parameters), intent(in)    :: param !! Current run configuration parameters 
+         real(DP),                   intent(in)    :: t     !! current time
+         character(*),               intent(in)    :: stage !! Either before or afte
+      end subroutine symba_util_take_collision_snapshot
+
       module subroutine symba_util_take_encounter_snapshot(self, param, t)
          use swiftest_classes, only : swiftest_parameters
          implicit none

src/symba/symba_collision.f90

@@ -885,7 +885,7 @@ module subroutine symba_resolve_collision_fragmentations(self, system, param)
       logical                                     :: lgoodcollision
       integer(I4B)                                :: i
 
-      associate(plplcollision_list => self, ncollisions => self%nenc, idx1 => self%index1, idx2 => self%index2, colliders => system%colliders, fragments => system%fragments)
+      associate(plplcollision_list => self, ncollisions => self%nenc, idx1 => self%index1, idx2 => self%index2, colliders => system%colliders, fragments => system%fragments, t => system%t)
          select type(pl => system%pl)
          class is (symba_pl)
             select type (cb => system%cb)
@@ -897,7 +897,8 @@ module subroutine symba_resolve_collision_fragmentations(self, system, param)
                   if ((.not. lgoodcollision) .or. any(pl%status(idx_parent(:)) /= COLLISION)) cycle
 
                   call colliders%regime(fragments, system, param)
-
+
+                  if (param%lencounter_save) call system%collision_snap(param, t, "before") 
                   select case (fragments%regime)
                   case (COLLRESOLVE_REGIME_DISRUPTION, COLLRESOLVE_REGIME_SUPERCATASTROPHIC)
                      plplcollision_list%status(i) = symba_collision_casedisruption(system, param)
@@ -909,6 +910,7 @@ module subroutine symba_resolve_collision_fragmentations(self, system, param)
                      write(*,*) "Error in symba_collision, unrecognized collision regime"
                      call util_exit(FAILURE)
                   end select
+                  if (param%lencounter_save) call system%collision_snap(param, t, "after") 
                end do
             end select
          end select

src/symba/symba_io.f90

@@ -101,7 +101,7 @@ module subroutine symba_io_param_reader(self, unit, iotype, v_list, iostat, ioms
                   read(param_value, *) param%encounter_save
                case ("FRAGMENTATION_SAVE")
                   call io_toupper(param_value)
-                  read(param_value, *) param%fragmentation_save
+                  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
@@ -159,14 +159,14 @@ module subroutine symba_io_param_reader(self, unit, iotype, v_list, iostat, ioms
             return
          end if
 
-         if ((param%fragmentation_save /= "NONE") .and. (param%fragmentation_save /= "TRAJECTORY") .and. (param%fragmentation_save /= "CLOSEST")) then
-            write(iomsg,*) 'Invalid fragmentation_save parameter: ',trim(adjustl(param%out_type))
+         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%fragmentation_save == "TRAJECTORY") .or. (param%fragmentation_save == "CLOSEST") 
+                                 (param%collision_save == "TRAJECTORY") .or. (param%collision_save == "CLOSEST") 
 
          ! Call the base method (which also prints the contents to screen)
          call io_param_reader(param, unit, iotype, v_list, iostat, iomsg) 

src/symba/symba_step.f90

@@ -37,9 +37,9 @@ module subroutine symba_step_system(self, param, t, dt)
                call self%reset(param)
                lencounter = pl%encounter_check(param, self, dt, 0) .or. tp%encounter_check(param, self, dt, 0)
                if (lencounter) then
-                  if (param%lencounter_save) call self%snapshot(param, t) 
+                  if (param%lencounter_save) call self%encounter_snap(param, t) 
                   call self%interp(param, t, dt)
-                  if (param%lencounter_save) call self%snapshot(param, t+dt) 
+                  if (param%lencounter_save) call self%encounter_snap(param, t+dt) 
                else
                   self%irec = -1
                   call helio_step_system(self, param, t, dt)
@@ -244,7 +244,7 @@ recursive module subroutine symba_step_recur_system(self, param, t, ireci)
                         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)
                      end if
-                     if (param%lencounter_save) call system%snapshot(param, t+dtl)
+                     if (param%lencounter_save) call system%encounter_snap(param, t+dtl)
 
                      call self%set_recur_levels(ireci)
 

src/symba/symba_util.f90

@@ -1287,6 +1287,30 @@ module subroutine symba_util_spill_tp(self, discards, lspill_list, ldestructive)
    end subroutine symba_util_spill_tp
 
 
+   module subroutine symba_util_take_collision_snapshot(self, param, t, stage)
+      !! author: David A. Minton
+      !!
+      !! Takes a minimal snapshot of the state of the system during an encounter so that the trajectories
+      !! can be played back through the encounter
+      implicit none
+      ! Internals
+      class(symba_nbody_system),  intent(inout) :: self  !! SyMBA nbody system object
+      class(swiftest_parameters), intent(in)    :: param !! Current run configuration parameters 
+      real(DP),                   intent(in)    :: t     !! current time
+      character(*),               intent(in)    :: stage !! Either before or after
+      ! Arguments
+      class(fraggle_collision_snapshot), allocatable :: snapshot
+
+      select case(stage)
+      case("before")
+
+      case("after")
+
+      end select
+
+      return
+   end subroutine symba_util_take_collision_snapshot
+
    module subroutine symba_util_take_encounter_snapshot(self, param, t)
       !! author: David A. Minton
       !!
@@ -1304,7 +1328,7 @@ module subroutine symba_util_take_encounter_snapshot(self, param, t)
       associate(npl => self%pl%nbody,  ntp => self%tp%nbody)
 
          if (self%plplenc_list%lcollision) then
-            allocate(fraggle_encounter_snapshot :: snapshot)
+            allocate(fraggle_collision_snapshot :: snapshot)
          else
             allocate(encounter_snapshot :: snapshot)
          end if
@@ -1388,7 +1412,7 @@ module subroutine symba_util_take_encounter_snapshot(self, param, t)
                end select
 
                select type(snapshot)
-               class is (fraggle_encounter_snapshot)
+               class is (fraggle_collision_snapshot)
                   allocate(snapshot%colliders, source=self%colliders)
                   allocate(snapshot%fragments, source=self%fragments)
                end select