Started implementing the CLOSEST vs. TRAJECTORY user options

MintonGroup · Dec 13, 2022 · c0adfd8 · c0adfd8
1 parent 92f1d8e
commit c0adfd8
Show file tree

Hide file tree

Showing 5 changed files with 18 additions and 15 deletions.
diff --git a/src/modules/symba_classes.f90 b/src/modules/symba_classes.f90
@@ -26,12 +26,13 @@ module symba_classes
    real(DP),     private, parameter :: RSHELL  = 0.48075_DP
 
    type, extends(swiftest_parameters) :: symba_parameters
-      real(DP)                                :: GMTINY             = -1.0_DP !! Smallest G*mass that is fully gravitating
-      real(DP)                                :: min_GMfrag         = -1.0_DP !! Smallest G*mass that can be produced in a fragmentation event
-      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
-      logical                                 :: lencounter_save   
+      real(DP)                                :: GMTINY               = -1.0_DP !! Smallest G*mass that is fully gravitating
+      real(DP)                                :: min_GMfrag           = -1.0_DP !! Smallest G*mass that can be produced in a fragmentation event
+      integer(I4B), dimension(:), allocatable :: seed                           !! Random seeds for fragmentation modeling
+      logical                                 :: lfragmentation       = .false. !! Do fragmentation modeling instead of simple merger.
+      character(STRMAX)                       :: encounter_save       = "NONE"  !! Indicate if and how encounter data should be saved
+      logical                                 :: lenc_trajectory_save = .false. !! Indicates that when encounters are saved, the full trajectory through recursion steps are saved
+      logical                                 :: lenc_closest_save    = .false. !! Indicates that when encounters are saved, the closest approach distance between pairs of bodies is saved
       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

diff --git a/src/setup/setup.f90 b/src/setup/setup.f90
@@ -75,7 +75,7 @@ module subroutine setup_construct_system(system, param)
 
             select type(param)
             class is (symba_parameters)
-               if (param%lencounter_save) then
+               if (param%lenc_trajectory_save .or. param%lenc_closest_save) then
                   allocate(encounter_storage :: param%encounter_history)
                   associate (encounter_history => param%encounter_history)
                      allocate(encounter_io_parameters :: encounter_history%nc)
@@ -85,7 +85,9 @@ module subroutine setup_construct_system(system, param)
                         nc%file_number = param%iloop / param%dump_cadence
                      end select
                   end associate
-
+               end if
+
+               if (param%lclose) then
                   allocate(collision_storage :: param%collision_history)
                   associate (collision_history => param%collision_history)
                      allocate(fraggle_io_parameters :: collision_history%nc)
@@ -98,7 +100,6 @@ module subroutine setup_construct_system(system, param)
                end if
             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
@@ -905,7 +905,7 @@ module subroutine symba_resolve_collision_fragmentations(self, system, param)
 
                   call system%colliders%regime(system%fragments, system, param)
 
-                  if (param%lencounter_save) call collision_history%take_snapshot(param,system, t, "before") 
+                  if (param%lenc_trajectory_save) call collision_history%take_snapshot(param,system, t, "before") 
                   select case (system%fragments%regime)
                   case (COLLRESOLVE_REGIME_DISRUPTION, COLLRESOLVE_REGIME_SUPERCATASTROPHIC)
                      plplcollision_list%status(i) = symba_collision_casedisruption(system, param)
@@ -917,7 +917,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 collision_history%take_snapshot(param,system, t, "after") 
+                  if (param%lenc_trajectory_save) call collision_history%take_snapshot(param,system, t, "after") 
                   deallocate(system%colliders,system%fragments)
                end do
             end select

diff --git a/src/symba/symba_io.f90 b/src/symba/symba_io.f90
@@ -125,7 +125,8 @@ module subroutine symba_io_param_reader(self, unit, iotype, v_list, iostat, ioms
             return
          end if
 
-         param%lencounter_save = (param%encounter_save == "TRAJECTORY") .or. (param%encounter_save == "CLOSEST")
+         param%lenc_trajectory_save = (param%encounter_save == "TRAJECTORY") 
+         param%lenc_closest_save = (param%encounter_save == "CLOSEST") .or. param%lenc_trajectory_save ! Closest approaches are always saved when trajectories are saved
 
          ! 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
@@ -38,9 +38,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 encounter_history%take_snapshot(param, self, t) 
+                     if (param%lenc_trajectory_save) call encounter_history%take_snapshot(param, self, t) 
                      call self%interp(param, t, dt)
-                     if (param%lencounter_save) call encounter_history%take_snapshot(param, self, t+dt) 
+                     if (param%lenc_trajectory_save) call encounter_history%take_snapshot(param, self, t+dt) 
                   else
                      self%irec = -1
                      call helio_step_system(self, param, t, dt)
@@ -247,7 +247,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 encounter_history%take_snapshot(param, self, t+dtl) 
+                     if (param%lenc_trajectory_save) call encounter_history%take_snapshot(param, self, t+dtl) 
 
                      call self%set_recur_levels(ireci)