Fixed a bunch of stuff related to taking collision snapshots, includi…

…ng consolidation of collision NetCDF into a single file.
MintonGroup · Jan 5, 2023 · 785d8d8 · 785d8d8
1 parent b5a46fa
commit 785d8d8
Show file tree

Hide file tree

Showing 13 changed files with 442 additions and 404 deletions.
diff --git a/src/base/base_module.f90 b/src/base/base_module.f90
@@ -85,15 +85,15 @@ module base
       logical :: ltides         = .false. !! Include tidal dissipation 
 
       ! Initial values to pass to the energy report subroutine (usually only used in the case of a restart, otherwise these will be updated with initial conditions values)
-      real(DP)                  :: Eorbit_orig  = 0.0_DP  !! Initial orbital energy
+      real(DP)                  :: E_orbit_orig  = 0.0_DP  !! Initial orbital energy
       real(DP)                  :: GMtot_orig   = 0.0_DP  !! Initial system mass
-      real(DP), dimension(NDIM) :: Ltot_orig    = 0.0_DP  !! Initial total angular momentum vector
-      real(DP), dimension(NDIM) :: Lorbit_orig  = 0.0_DP  !! Initial orbital angular momentum
-      real(DP), dimension(NDIM) :: Lspin_orig   = 0.0_DP  !! Initial spin angular momentum vector
-      real(DP), dimension(NDIM) :: Lescape      = 0.0_DP  !! Angular momentum of bodies that escaped the system (used for bookeeping)
+      real(DP), dimension(NDIM) :: L_total_orig    = 0.0_DP  !! Initial total angular momentum vector
+      real(DP), dimension(NDIM) :: L_orbit_orig  = 0.0_DP  !! Initial orbital angular momentum
+      real(DP), dimension(NDIM) :: L_spin_orig   = 0.0_DP  !! Initial spin angular momentum vector
+      real(DP), dimension(NDIM) :: L_escape      = 0.0_DP  !! Angular momentum of bodies that escaped the system (used for bookeeping)
       real(DP)                  :: GMescape     = 0.0_DP  !! Mass of bodies that escaped the system (used for bookeeping)
-      real(DP)                  :: Ecollisions  = 0.0_DP  !! Energy lost from system due to collisions
-      real(DP)                  :: Euntracked   = 0.0_DP  !! Energy gained from system due to escaped bodies
+      real(DP)                  :: E_collisions  = 0.0_DP  !! Energy lost from system due to collisions
+      real(DP)                  :: E_untracked   = 0.0_DP  !! Energy gained from system due to escaped bodies
       logical                   :: lfirstenergy = .true.  !! This is the first time computing energe
       logical                   :: lfirstkick   = .true.  !! Initiate the first kick in a symplectic step
       logical                   :: lrestart     = .false. !! Indicates whether or not this is a restarted run

diff --git a/src/collision/collision_generate.f90 b/src/collision/collision_generate.f90
@@ -55,6 +55,13 @@ module subroutine collision_generate_bounce(self, nbody_system, param, t)
          associate(impactors => nbody_system%collider%impactors, fragments => nbody_system%collider%fragments)
             select case (impactors%regime) 
             case (COLLRESOLVE_REGIME_DISRUPTION, COLLRESOLVE_REGIME_SUPERCATASTROPHIC)
+
+               ! Manually save the before/after snapshots because this case doesn't use the mergeaddsub procedure
+               select type(before => self%before)
+               class is (swiftest_nbody_system)
+                  allocate(before%pl, source=pl) 
+               end select
+
                nfrag = size(impactors%id(:))
                do i = 1, nfrag
                   j = impactors%id(i)
@@ -68,13 +75,12 @@ module subroutine collision_generate_bounce(self, nbody_system, param, t)
                   pl%ldiscard(j) = .false.
                   pl%lcollision(j) = .false.
                end do
-               select type(before => self%before)
-               class is (swiftest_nbody_system)
+
                select type(after => self%after)
                class is (swiftest_nbody_system)
-                  allocate(after%pl, source=before%pl) ! Be sure to save the pl so that snapshots still work   
-               end select
+                  allocate(after%pl, source=pl) 
                end select
+
             case (COLLRESOLVE_REGIME_HIT_AND_RUN)
                call self%hitandrun(nbody_system, param, t)
             case (COLLRESOLVE_REGIME_MERGE, COLLRESOLVE_REGIME_GRAZE_AND_MERGE)
@@ -107,7 +113,6 @@ module subroutine collision_generate_hitandrun(self, nbody_system, param, t)
       ! Internals
       character(len=STRMAX) :: message
 
-
       select type(nbody_system)
       class is (swiftest_nbody_system)
       select type(pl => nbody_system%pl)
@@ -117,18 +122,21 @@ module subroutine collision_generate_hitandrun(self, nbody_system, param, t)
             call collision_io_collider_message(nbody_system%pl, impactors%id, message)
             call swiftest_io_log_one_message(COLLISION_LOG_OUT, trim(adjustl(message)))
 
+            ! Manually save the before/after snapshots because this case doesn't use the mergeaddsub procedure
+            select type(before => self%before)
+            class is (swiftest_nbody_system)
+               allocate(before%pl, source=pl) 
+            end select
+
             status = HIT_AND_RUN_PURE
             pl%status(impactors%id(:)) = ACTIVE
             pl%ldiscard(impactors%id(:)) = .false.
             pl%lcollision(impactors%id(:)) = .false.
-            ! Be sure to save the pl so that snapshots still work 
-            select type(before => self%before)
-            class is (swiftest_nbody_system)
+
             select type(after => self%after)
             class is (swiftest_nbody_system)
-               allocate(after%pl, source=before%pl) 
+               allocate(after%pl, source=pl) 
             end select
-            end select  
 
          end associate
       end select
@@ -154,7 +162,7 @@ module subroutine collision_generate_merge(self, nbody_system, param, t)
       real(DP),                 intent(in)    :: t            !! The time of the collision
       ! Internals
       integer(I4B)                              :: i, j, k, ibiggest
-      real(DP), dimension(NDIM)                 :: Lspin_new
+      real(DP), dimension(NDIM)                 :: L_spin_new
       real(DP)                                  :: volume, G
       character(len=STRMAX) :: message
 
@@ -193,12 +201,12 @@ module subroutine collision_generate_merge(self, nbody_system, param, t)
                   if (param%lrotation) then
                      do concurrent(i = 1:NDIM)
                         fragments%Ip(i,1) = sum(impactors%mass(:) * impactors%Ip(i,:)) 
-                        Lspin_new(i) = sum(impactors%Lorbit(i,:) + impactors%Lorbit(i,:))
+                        L_spin_new(i) = sum(impactors%L_orbit(i,:) + impactors%L_orbit(i,:))
                      end do
                      fragments%Ip(:,1) = fragments%Ip(:,1) / fragments%mass(1)
-                     fragments%rot(:,1) = Lspin_new(:) / (fragments%Ip(3,1) * fragments%mass(1) * fragments%radius(1)**2)
+                     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
-                     nbody_system%Lescape(:) = nbody_system%Lescape(:) + impactors%Lorbit(:,1) + impactors%Lorbit(:,2) 
+                     nbody_system%L_escape(:) = nbody_system%L_escape(:) + impactors%L_orbit(:,1) + impactors%L_orbit(:,2) 
                   end if
 
                   ! The fragment trajectory will be the barycentric trajectory