Fixed bugs in Fraggle causing arithmetic exception when initial spins…

… were 0
MintonGroup · Nov 29, 2022 · 8b10f65 · 8b10f65
1 parent 3175f87
commit 8b10f65
Show file tree

Hide file tree

Showing 10 changed files with 42 additions and 92 deletions.
diff --git a/examples/Fragmentation/Fragmentation_Movie.py b/examples/Fragmentation/Fragmentation_Movie.py
@@ -43,6 +43,7 @@
 movie_title_list = ["Head-on Disrutption", "Off-axis Supercatastrophic", "Hit and Run"]
 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, -1.807993e-05, 0.0]),
                                               np.array([1.0,  1.807993e-05 ,0.0])],
                "supercatastrophic_off_axis": [np.array([1.0, -4.2e-05,      0.0]),
@@ -67,7 +68,7 @@
                                               np.array([0.0, 0.0, 1.0e5])]
                }
 
-body_Gmass = {"disruption_headon"         : [1e-7, 7e-10],
+body_Gmass = {"disruption_headon"        : [1e-7, 1e-10],
              "supercatastrophic_off_axis": [1e-7, 1e-8],
              "hitandrun"                 : [1e-7, 7e-10]
                }
@@ -131,8 +132,10 @@ def animate(i,ds,movie_title):
     sim.add_body(Gmass=body_Gmass[style], radius=body_radius[style], xh=pos_vectors[style], vh=vel_vectors[style], rot=rot_vectors[style])
 
     # Set fragmentation parameters
-    sim.set_parameter(fragmentation = True, gmtiny=1e-11, minimum_fragment_gmass=1e-11)
-    sim.run(dt=1e-8, tstop=1e-5)
+    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, gmtiny=gmtiny, minimum_fragment_gmass=minimum_fragment_gmass)
+    sim.run(dt=1e-8, tstop=1.e-5)
 
     # Calculate the number of frames in the dataset.
     nframes = int(sim.data['time'].size)
@@ -144,5 +147,5 @@ def animate(i,ds,movie_title):
     # Set up the figure and the animation.
     fig, ax = plt.subplots(figsize=(4,4))
     # Generate the movie.
-    ani = animation.FuncAnimation(fig, animate, interval=1, frames=nframes, repeat=False)
+    ani = animation.FuncAnimation(fig, animate, fargs=(sim.data, movie_titles[style]), interval=1, frames=nframes, repeat=False)
     ani.save(movie_filename, fps=60, dpi=300, extra_args=['-vcodec', 'libx264'])
diff --git a/examples/Fragmentation/cb.in b/examples/Fragmentation/cb.in
diff --git a/examples/Fragmentation/disruption_headon.in b/examples/Fragmentation/disruption_headon.in
diff --git a/examples/Fragmentation/hitandrun.in b/examples/Fragmentation/hitandrun.in
diff --git a/examples/Fragmentation/supercatastrophic_off_axis.in b/examples/Fragmentation/supercatastrophic_off_axis.in
diff --git a/examples/Fragmentation/tp.in b/examples/Fragmentation/tp.in
diff --git a/python/swiftest/swiftest/simulation_class.py b/python/swiftest/swiftest/simulation_class.py
@@ -31,6 +31,7 @@
     Literal,
     Dict,
     List,
+    Tuple,
     Any
 )
 
@@ -387,7 +388,7 @@ def _type_scrub(output_data):
         if "nplm" in self.data:
             post_message += f" nplm: {self.data['nplm'].values[0]}"
         if self.param['ENERGY']:
-            post_message += f" dL/L0: {0.0:.5e} dE/|E0|: {0.0:.5e}"
+            post_message += f" dL/L0: {0.0:.5e} dE/|E0|: {0.0:+.5e}"
         pbar = tqdm(total=noutput, desc=pre_message, postfix=post_message, bar_format='{l_bar}{bar}{postfix}')
         try:
             with subprocess.Popen(shlex.split(cmd),
@@ -410,7 +411,7 @@ def _type_scrub(output_data):
                         if "LTOTERR" in output_data:
                             post_message += f" dL/L0: {output_data['LTOTERR']:.5e}"
                         if "ETOTERR" in output_data:
-                            post_message += f" dE/|E0|: {output_data['ETOTERR']:.5e}"
+                            post_message += f" dE/|E0|: {output_data['ETOTERR']:+.5e}"
                         interval = output_data['ILOOP'] - iloop
                         if interval > 0:
                            pbar.update(interval)
@@ -2278,8 +2279,8 @@ def add_body(self,
                  v4: float | List[float] | npt.NDArray[np.float_] | None = None,
                  v5: float | List[float] | npt.NDArray[np.float_] | None = None,
                  v6: float | List[float] | npt.NDArray[np.float_] | None = None,
-                 xh: List[float] | npt.NDArray[np.float_] | None = None,
-                 vh: List[float] | npt.NDArray[np.float_] | None = None,
+                 xh: List[float] | List[npt.NDArray[np.float_]] | npt.NDArray[np.float_] | None = None,
+                 vh: List[float] | List[npt.NDArray[np.float_]] | npt.NDArray[np.float_] | None = None,
                  mass: float | List[float] | npt.NDArray[np.float_] | None=None,
                  Gmass: float | List[float] | npt.NDArray[np.float_] | None=None,
                  radius: float | List[float] | npt.NDArray[np.float_] | None=None,
@@ -2291,7 +2292,7 @@ def add_body(self,
                  rotx: float | List[float] | npt.NDArray[np.float_] | None=None,
                  roty: float | List[float] | npt.NDArray[np.float_] | None=None,
                  rotz: float | List[float] | npt.NDArray[np.float_] | None=None,
-                 rot: List[float] | npt.NDArray[np.float_] | None=None,
+                 rot: List[float] | List[npt.NDArray[np.float_]] | npt.NDArray[np.float_] | None=None,
                  J2: float | List[float] | npt.NDArray[np.float_] | None=None,
                  J4: float | List[float] | npt.NDArray[np.float_] | None=None):
         """

diff --git a/src/fraggle/fraggle_generate.f90 b/src/fraggle/fraggle_generate.f90
@@ -57,8 +57,12 @@ module subroutine fraggle_generate_fragments(self, colliders, system, param, lfa
          end if
          f_spin = F_SPIN_FIRST
 
-         lk_plpl = allocated(pl%k_plpl)
-         if (lk_plpl) deallocate(pl%k_plpl)
+         if (param%lflatten_interactions) then
+            lk_plpl = allocated(pl%k_plpl)
+            if (lk_plpl) deallocate(pl%k_plpl)
+         else 
+            lk_plpl = .false.
+         end if
 
          call frag%set_natural_scale(colliders)
 
@@ -250,19 +254,21 @@ subroutine fraggle_generate_spins(frag, f_spin, lfailure)
          frag%rot(:,:) = 0.0_DP
 
          frag%ke_spin = 0.0_DP
-         do i = 1, nfrag
-            ! Convert a fraction (f_spin) of either the remaining angular momentum or kinetic energy budget into spin, whichever gives the smaller rotation so as not to blow any budgets
-            rot_ke(:) = sqrt(2 * f_spin * frag%ke_budget / (nfrag * frag%mass(i) * frag%radius(i)**2 * frag%Ip(3, i))) &
-                        * L_remainder(:) / norm2(L_remainder(:))
-            rot_L(:) = f_spin * L_remainder(:) / (nfrag * frag%mass(i) * frag%radius(i)**2 * frag%Ip(3, i))
-            if (norm2(rot_ke) < norm2(rot_L)) then
-               frag%rot(:,i) = rot_ke(:)
-            else
-               frag%rot(:, i) = rot_L(:)
-            end if
-            frag%ke_spin = frag%ke_spin + frag%mass(i) * frag%Ip(3, i) * frag%radius(i)**2 &
-                                                       * dot_product(frag%rot(:, i), frag%rot(:, i))
-         end do
+         if (norm2(L_remainder(:)) > FRAGGLE_LTOL) then 
+            do i = 1, nfrag
+               ! Convert a fraction (f_spin) of either the remaining angular momentum or kinetic energy budget into spin, whichever gives the smaller rotation so as not to blow any budgets
+               rot_ke(:) = sqrt(2 * f_spin * frag%ke_budget / (nfrag * frag%mass(i) * frag%radius(i)**2 * frag%Ip(3, i))) &
+                           * L_remainder(:) / norm2(L_remainder(:))
+               rot_L(:) = f_spin * L_remainder(:) / (nfrag * frag%mass(i) * frag%radius(i)**2 * frag%Ip(3, i))
+               if (norm2(rot_ke) < norm2(rot_L)) then
+                  frag%rot(:,i) = rot_ke(:)
+               else
+                  frag%rot(:, i) = rot_L(:)
+               end if
+               frag%ke_spin = frag%ke_spin + frag%mass(i) * frag%Ip(3, i) * frag%radius(i)**2 &
+                                                         * dot_product(frag%rot(:, i), frag%rot(:, i))
+            end do
+         end if
          frag%ke_spin = 0.5_DP * frag%ke_spin
 
          lfailure = ((frag%ke_budget - frag%ke_spin - frag%ke_orbit) < 0.0_DP)

diff --git a/src/fraggle/fraggle_util.f90 b/src/fraggle/fraggle_util.f90
@@ -87,7 +87,6 @@ module subroutine fraggle_util_construct_temporary_system(frag, system, param, t
       !! Author: David A. Minton
       !!
       !! Constructs a temporary internal system consisting of active bodies and additional fragments. This internal temporary system is used to calculate system energy with and without fragments
-      !! and optionally including fragments.
       implicit none
       ! Arguments
       class(fraggle_fragments),                   intent(in)  :: frag     !! Fraggle fragment system object
@@ -147,7 +146,6 @@ module subroutine fraggle_util_get_energy_momentum(self, colliders, system, para
       class(swiftest_parameters),   intent(inout) :: param     !! Current swiftest run configuration parameters
       logical,                      intent(in)    :: lbefore   !! Flag indicating that this the "before" state of the system, with colliders included and fragments excluded or vice versa
       ! Internals
-      logical, dimension(:), allocatable :: lexclude
       class(swiftest_nbody_system), allocatable, save :: tmpsys
       class(swiftest_parameters), allocatable, save   :: tmpparam
       integer(I4B)  :: npl_before, npl_after
@@ -162,24 +160,24 @@ module subroutine fraggle_util_get_energy_momentum(self, colliders, system, para
          npl_before = pl%nbody
          npl_after = npl_before + nfrag
 
-         ! Build the exluded body logical mask
-         allocate(lexclude(npl_after))
          if (lbefore) then
-            lexclude(1:npl_before) = .false.
-            lexclude(npl_before+1:npl_after) = .true.
             call fraggle_util_construct_temporary_system(frag, system, param, tmpsys, tmpparam)
+            ! Build the exluded body logical mask for the *before* case: Only the original bodies are used to compute energy and momentum
+            tmpsys%pl%status(colliders%idx(1:colliders%ncoll)) = ACTIVE
+            tmpsys%pl%status(npl_before+1:npl_after) = INACTIVE
          else
-            lexclude(1:npl_after) = .false.
-            lexclude(colliders%idx(1:colliders%ncoll)) = .true.
             if (.not.allocated(tmpsys)) then
                write(*,*) "Error in fraggle_util_get_energy_momentum. " // &
                          " This must be called with lbefore=.true. at least once before calling it with lbefore=.false."
                call util_exit(FAILURE)
             end if
+            ! Build the exluded body logical mask for the *after* case: Only the new bodies are used to compute energy and momentum
             call fraggle_util_add_fragments_to_system(frag, colliders, tmpsys, tmpparam)
+            tmpsys%pl%status(colliders%idx(1:colliders%ncoll)) = INACTIVE
+            tmpsys%pl%status(npl_before+1:npl_after) = ACTIVE
          end if 
 
-         call tmpsys%pl%flatten(param)
+         if (param%lflatten_interactions) call tmpsys%pl%flatten(param)
 
          call tmpsys%get_energy_and_momentum(param) 
 

diff --git a/src/util/util_get_energy_momentum.f90 b/src/util/util_get_energy_momentum.f90
@@ -23,7 +23,6 @@ module subroutine util_get_energy_momentum_system(self, param)
       class(swiftest_parameters),   intent(in)    :: param    !! Current run configuration parameters
       ! Internals
       integer(I4B) :: i
-      integer(I8B) :: nplpl
       real(DP) :: kecb, kespincb
       real(DP), dimension(self%pl%nbody) :: kepl, kespinpl
       real(DP), dimension(self%pl%nbody) :: Lplorbitx, Lplorbity, Lplorbitz
@@ -32,7 +31,6 @@ module subroutine util_get_energy_momentum_system(self, param)
       real(DP) :: hx, hy, hz
 
       associate(system => self, pl => self%pl, npl => self%pl%nbody, cb => self%cb)
-         nplpl = pl%nplpl
          system%Lorbit(:) = 0.0_DP
          system%Lspin(:) = 0.0_DP
          system%Ltot(:) = 0.0_DP
@@ -89,7 +87,7 @@ module subroutine util_get_energy_momentum_system(self, param)
          end if
 
          if (param%lflatten_interactions) then
-            call util_get_energy_potential_flat(npl, nplpl, pl%k_plpl, pl%lmask, cb%Gmass, pl%Gmass, pl%mass, pl%xb, system%pe)
+            call util_get_energy_potential_flat(npl, pl%nplpl, pl%k_plpl, pl%lmask, cb%Gmass, pl%Gmass, pl%mass, pl%xb, system%pe)
          else
             call util_get_energy_potential_triangular(npl, pl%lmask, cb%Gmass, pl%Gmass, pl%mass, pl%xb, system%pe)
          end if