Adjusted Fraggle so that if there is a failure, it does a bounce inst…

…ead of a merge

cmake/Modules/SetFortranFlags.cmake

@@ -121,6 +121,10 @@ SET_COMPILE_FLAG(CMAKE_Fortran_FLAGS_DEBUG "${CMAKE_Fortran_FLAGS_DEBUG}"
                          "/check"  # Intel Windows
                          "-fcheck=all" # GNU 
                 )
+SET_COMPILE_FLAG(CMAKE_Fortran_FLAGS_DEBUG "${CMAKE_Fortran_FLAGS_DEBUG}"
+                 Fortran "-fstack-check" # GNU 
+                )
+
 
 # Initializes matrices/arrays with NaN values
 SET_COMPILE_FLAG(CMAKE_Fortran_FLAGS_DEBUG "${CMAKE_Fortran_FLAGS_DEBUG}"

src/collision/collision_generate.f90

@@ -102,7 +102,7 @@ module subroutine collision_generate_bounce(self, nbody_system, param, t)
       return
    end subroutine collision_generate_bounce
 
-
+   
    module subroutine collision_generate_hitandrun(self, nbody_system, param, t) 
       !! author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton
       !!

src/collision/collision_util.f90

@@ -486,39 +486,36 @@ module subroutine collision_util_setup_fragments(self, n)
 
       if (n < 0) return
 
-      call self%dealloc()
-
       self%nbody = n
       if (n == 0) return
 
-      allocate(swiftest_particle_info :: self%info(n))
-
-      allocate(self%id(n))
-      allocate(self%status(n))
-      allocate(self%rh(NDIM, n))
-      allocate(self%vh(NDIM, n))
-      allocate(self%rb(NDIM, n))
-      allocate(self%vb(NDIM, n))
-      allocate(self%rc(NDIM, n))
-      allocate(self%vc(NDIM, n))
-      allocate(self%r_unit(NDIM, n))
-      allocate(self%v_unit(NDIM, n))
-      allocate(self%t_unit(NDIM, n))
-      allocate(self%n_unit(NDIM, n))
-      allocate(self%rot(NDIM, n))
-      allocate(self%Ip(NDIM, n))
-      allocate(self%Gmass(n))
-      allocate(self%mass(n))
-      allocate(self%radius(n))
-      allocate(self%density(n))
-      allocate(self%rmag(n))
-      allocate(self%vmag(n))
-      allocate(self%rotmag(n))
-      allocate(self%origin_body(n))
-      allocate(self%L_orbit(NDIM, n))
-      allocate(self%L_spin(NDIM, n))
-      allocate(self%ke_orbit(n))
-      allocate(self%ke_spin(n))
+      if(allocated(self%info)) deallocate(self%info); allocate(swiftest_particle_info :: self%info(n))
+      if (allocated(self%id)) deallocate(self%id); allocate(self%id(n))
+      if (allocated(self%status)) deallocate(self%status); allocate(self%status(n))
+      if (allocated(self%rh)) deallocate(self%rh); allocate(self%rh(NDIM, n))
+      if (allocated(self%vh)) deallocate(self%vh); allocate(self%vh(NDIM, n))
+      if (allocated(self%rb)) deallocate(self%rb); allocate(self%rb(NDIM, n))
+      if (allocated(self%vb)) deallocate(self%vb); allocate(self%vb(NDIM, n))
+      if (allocated(self%rc)) deallocate(self%rc); allocate(self%rc(NDIM, n))
+      if (allocated(self%vc)) deallocate(self%vc); allocate(self%vc(NDIM, n))
+      if (allocated(self%r_unit)) deallocate(self%r_unit); allocate(self%r_unit(NDIM, n))
+      if (allocated(self%v_unit)) deallocate(self%v_unit); allocate(self%v_unit(NDIM, n))
+      if (allocated(self%t_unit)) deallocate(self%t_unit); allocate(self%t_unit(NDIM, n))
+      if (allocated(self%n_unit)) deallocate(self%n_unit); allocate(self%n_unit(NDIM, n))
+      if (allocated(self%rot)) deallocate(self%rot); allocate(self%rot(NDIM, n))
+      if (allocated(self%Ip)) deallocate(self%Ip); allocate(self%Ip(NDIM, n))
+      if (allocated(self%Gmass)) deallocate(self%Gmass); allocate(self%Gmass(n))
+      if (allocated(self%mass)) deallocate(self%mass); allocate(self%mass(n))
+      if (allocated(self%radius)) deallocate(self%radius); allocate(self%radius(n))
+      if (allocated(self%density)) deallocate(self%density); allocate(self%density(n))
+      if (allocated(self%rmag)) deallocate(self%rmag); allocate(self%rmag(n))
+      if (allocated(self%vmag)) deallocate(self%vmag); allocate(self%vmag(n))
+      if (allocated(self%rotmag)) deallocate(self%rotmag); allocate(self%rotmag(n))
+      if (allocated(self%origin_body)) deallocate(self%origin_body); allocate(self%origin_body(n))
+      if (allocated(self%L_orbit)) deallocate(self%L_orbit); allocate(self%L_orbit(NDIM, n))
+      if (allocated(self%L_spin)) deallocate(self%L_spin); allocate(self%L_spin(NDIM, n))
+      if (allocated(self%ke_orbit)) deallocate(self%ke_orbit); allocate(self%ke_orbit(n))
+      if (allocated(self%ke_spin)) deallocate(self%ke_spin); allocate(self%ke_spin(n))
 
       self%id(:) = 0
       select type(info => self%info)

src/fraggle/fraggle_generate.f90

@@ -25,7 +25,8 @@ module subroutine fraggle_generate(self, nbody_system, param, t)
       class(base_parameters),   intent(inout) :: param        !! Current run configuration parameters 
       real(DP),                 intent(in)    :: t            !! Time of collision
       ! Internals
-      integer(I4B)          :: i, ibiggest, nfrag
+      integer(I4B)          :: i, j, ibiggest, nfrag, nimp
+      real(DP), dimension(NDIM) :: vcom, rnorm
       character(len=STRMAX) :: message 
       logical               :: lfailure
 
@@ -52,15 +53,29 @@ module subroutine fraggle_generate(self, nbody_system, param, t)
             call self%set_mass_dist(param) 
             call self%disrupt(nbody_system, param, t, lfailure)
             if (lfailure) then
-               call swiftest_io_log_one_message(COLLISION_LOG_OUT, "Fraggle failed to find an energy-losing solution. Treating this as a merger.") 
-               call self%merge(nbody_system, param, t) 
-               return
+               call swiftest_io_log_one_message(COLLISION_LOG_OUT, "Fraggle failed to find an energy-losing solution. Treating this as a bounce.") 
+
+               nimp = size(impactors%id(:))
+               call self%fragments%setup(nimp)
+               do i = 1, nimp
+                  j = impactors%id(i)
+                  vcom(:) = pl%vb(:,j) - impactors%vbcom(:)
+                  rnorm(:) = .unit. (impactors%rb(:,2) - impactors%rb(:,1))
+                  ! Do the reflection
+                  vcom(:) = vcom(:) - 2 * dot_product(vcom(:),rnorm(:)) * rnorm(:)
+                  self%fragments%vb(:,i) = impactors%vbcom(:) + vcom(:)
+                  self%fragments%rb(:,i) = pl%rb(:,j)
+                  self%fragments%mass(i) = pl%mass(j)
+                  self%fragments%Gmass(i) = pl%Gmass(j)
+                  self%fragments%radius(i) = pl%radius(j)
+                  self%fragments%rot(:,i) = pl%rot(:,j)
+                  self%fragments%Ip(:,i) = pl%Ip(:,j)
+               end do
             end if
 
             associate (fragments => self%fragments)
                ! Populate the list of new bodies
                nfrag = fragments%nbody
-               write(message, *) nfrag
                select case(impactors%regime)
                case(COLLRESOLVE_REGIME_DISRUPTION)
                   status = DISRUPTED
@@ -454,7 +469,7 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu
       class(swiftest_parameters),   intent(inout) :: param        !! Current run configuration parameters 
       logical,                      intent(out)   :: lfailure     !! Did the velocity computation fail?
       ! Internals
-      integer(I4B) :: i, j, loop, try, istart, nfrag, nlast
+      integer(I4B) :: i, j, loop, try, istart, nfrag, nlast, nsteps
       logical :: lhitandrun, lsupercat
       real(DP), dimension(NDIM) :: vimp_unit, rimp, vrot, L_residual
       real(DP) :: vmag, vesc, dE, E_residual, ke_min, ke_avail, ke_remove, dE_best, E_residual_best, fscale, f_spin, f_orbit, dE_metric
@@ -504,7 +519,7 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu
             lfailure = .false.
             dE_metric = huge(1.0_DP)
 
-            outer: do try = 1, maxtry
+            outer: do try = 1, nfrag
                ! Scale the magnitude of the velocity by the distance from the impact point
                ! This will reduce the chances of fragments colliding with each other immediately, and is more physically correct  
                do concurrent(i = 1:nfrag)
@@ -543,6 +558,7 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu
                ke_min = 0.5_DP * fragments%mtot * vesc**2
 
                do loop = 1, MAXLOOP
+                  nsteps = loop * try
                   call collider_local%get_energy_and_momentum(nbody_system, param, phase="after")
                   ke_avail = max(fragments%ke_orbit_tot - ke_min, 0.0_DP)
                   ! Check for any residual angular momentum, and if there is any, put it into spin of the largest body
@@ -601,8 +617,9 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu
                   call fragments%set_coordinate_system()
 
                end do
+               if (dE_best < 0.0_DP) exit outer
                ! We didn't converge. Reset the fragment positions and velocities and try a new configuration with some slightly different parameters
-               if (fragments%nbody == 2) exit
+               if (fragments%nbody == 2) exit outer
                ! Reduce the number of fragments by one
                nlast = fragments%nbody
                fragments%Ip(:,1) = fragments%mass(1) * fragments%Ip(:,1) + fragments%mass(nlast) * fragments%Ip(:,nlast)
@@ -634,7 +651,7 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu
             end do outer
             lfailure = dE_best > 0.0_DP
 
-            write(message, *) try*loop
+            write(message, *) nsteps
             if (lfailure) then
                call swiftest_io_log_one_message(COLLISION_LOG_OUT, "Fraggle velocity calculation failed to converge after " // trim(adjustl(message)) // " steps. This collision would add energy.")
             else