Added new runtime checks to the debug compile mode to catch calls to …

…unallocated variables
MintonGroup · Jan 13, 2023 · bf4fc81 · bf4fc81
1 parent ff3d7bb
commit bf4fc81
Show file tree

Hide file tree

Showing 3 changed files with 30 additions and 10 deletions.
diff --git a/cmake/Modules/SetFortranFlags.cmake b/cmake/Modules/SetFortranFlags.cmake
@@ -110,10 +110,10 @@ SET_COMPILE_FLAG(CMAKE_Fortran_FLAGS_DEBUG "${CMAKE_Fortran_FLAGS_DEBUG}"
                          "-ftrace=full" # GNU (g95)
                 )
 
-# Check array bounds
+# Check everything
 SET_COMPILE_FLAG(CMAKE_Fortran_FLAGS_DEBUG "${CMAKE_Fortran_FLAGS_DEBUG}"
-                 Fortran "-check bounds"  # Intel
-                         "/check:bounds"  # Intel Windows
+                 Fortran "-check"  # Intel
+                         "/check"  # Intel Windows
                          "-fcheck=bounds" # GNU (New style)
                          "-fbounds-check" # GNU (Old style)
                          "-Mbounds"       # Portland Group

diff --git a/src/fraggle/fraggle_generate.f90 b/src/fraggle/fraggle_generate.f90
@@ -447,7 +447,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
+      integer(I4B) :: i, j, loop, try, istart, nfrag, nlast
       logical :: lhitandrun, lsupercat
       real(DP), dimension(NDIM) :: vimp_unit, rimp, vrot, L_residual
       real(DP) :: vmag, vesc, dE, E_residual, ke_avail, ke_remove, dE_best, E_residual_best, fscale, f_spin, f_orbit, dE_metric
@@ -456,7 +456,7 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu
       ! For the initial "guess" of fragment velocities, this is the minimum and maximum velocity relative to escape velocity that the fragments will have
       real(DP)                :: vmin_guess = 1.5_DP 
       real(DP)                :: vmax_guess = 10.0_DP
-      real(DP)                :: delta_v
+      real(DP)                :: delta_v, volume
       integer(I4B), parameter :: MAXLOOP = 20
       integer(I4B), parameter :: MAXTRY  = 100
       real(DP), parameter :: SUCCESS_METRIC = 1.0e-3_DP
@@ -497,7 +497,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 - 1
                ! 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)
@@ -551,7 +551,7 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu
                      E_residual_best = E_residual
                      dE_best = dE
 
-                     do concurrent(i = 1:nfrag)
+                     do concurrent(i = 1:fragments%nbody)
                         fragments%vb(:,i) = fragments%vc(:,i) + impactors%vbcom(:)
                      end do
 
@@ -576,7 +576,7 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu
                   ke_remove = min(f_spin * E_residual, fragments%ke_spin_tot)
                   ke_rot_remove(:) = ke_remove * (fragments%ke_spin(:) / fragments%ke_spin_tot)
                   where(ke_rot_remove(:) > fragments%ke_spin(:)) ke_rot_remove(:) = fragments%ke_spin(:) 
-                  do concurrent(i = 1:nfrag, fragments%ke_spin(i) > 10*sqrt(tiny(1.0_DP)))
+                  do concurrent(i = 1:fragments%nbody, fragments%ke_spin(i) > 10*sqrt(tiny(1.0_DP)))
                      fscale = sqrt((fragments%ke_spin(i) - ke_rot_remove(i))/fragments%ke_spin(i))
                      fragments%rot(:,i) = fscale * fragments%rot(:,i)
                   end do
@@ -586,16 +586,34 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu
                   call fragments%set_coordinate_system()
 
                end do
-               ! We didn't converge. Try another configuration and see if we get a better result
+               ! We didn't converge. Reset the fragment positions and velocities and try a new configuration with some slightly different parameters
+               ! Reduce the number of fragments by one
+               nlast = fragments%nbody
+               fragments%Ip(:,1) = fragments%mass(1) * impactors%Ip(:,1) + fragments%mass(nlast) * impactors%Ip(:,nlast)
+               fragments%mass(1) = fragments%mass(1) + fragments%mass(nlast)
+               fragments%Ip(:,1) = fragments%Ip(:,1) / fragments%mass(1)
+               fragments%Gmass(1) = fragments%Gmass(1) + fragments%mass(nlast)
+               volume = 4.0_DP / 3.0_DP * PI * (fragments%radius(1) + fragments%radius(nlast))
+               fragments%density(1) = fragments%mass(1) / volume
+               fragments%radius(1) = (3._DP * volume / (4._DP * PI))**(THIRD)
+               fragments%Ip(:,nlast) = 0.0_DP
+               fragments%mass(nlast) = 0.0_DP
+               fragments%Gmass(nlast) = 0.0_DP
+               fragments%radius(nlast) = 0.0_DP
+               fragments%nbody = nlast - 1
+
                call fragments%reset()
                call fraggle_generate_pos_vec(collider_local)
                call fraggle_generate_rot_vec(collider_local, nbody_system, param)
+
+               ! Increase the spatial size factor to get a less dense cloud
                collider_local%fail_scale = collider_local%fail_scale*1.01_DP
 
-               ! Bring the minimum and maximum velocities closer together for the next round
+               ! Bring the minimum and maximum velocities closer together 
                delta_v = 0.125_DP * (vmax_guess - vmin_guess)
                vmin_guess = vmin_guess + delta_v
                vmax_guess = vmax_guess - delta_v
+
             end do outer
             lfailure = dE_best > 0.0_DP
 

diff --git a/src/swiftest/swiftest_util.f90 b/src/swiftest/swiftest_util.f90
@@ -1770,6 +1770,8 @@ module subroutine swiftest_util_peri_body(self, nbody_system, param)
       ! Internals
       integer(I4B) :: i
 
+      if (self%nbody == 0) return
+
       select type(self)
       class is (swiftest_pl)
          if (self%lfirst) self%isperi(:) = 0