diff --git a/cmake/Modules/SetFortranFlags.cmake b/cmake/Modules/SetFortranFlags.cmake index cb5ab68ad..2903375bb 100644 --- a/cmake/Modules/SetFortranFlags.cmake +++ b/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}" diff --git a/src/collision/collision_generate.f90 b/src/collision/collision_generate.f90 index e5c3035b9..a7ef92363 100644 --- a/src/collision/collision_generate.f90 +++ b/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 !! diff --git a/src/collision/collision_util.f90 b/src/collision/collision_util.f90 index 4a089dfa6..772d6f988 100644 --- a/src/collision/collision_util.f90 +++ b/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) diff --git a/src/fraggle/fraggle_generate.f90 b/src/fraggle/fraggle_generate.f90 index aea492a63..fd8f79014 100644 --- a/src/fraggle/fraggle_generate.f90 +++ b/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