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/examples/Fragmentation/Fragmentation_Movie.py b/examples/Fragmentation/Fragmentation_Movie.py index 90199fd4b..7d774f983 100755 --- a/examples/Fragmentation/Fragmentation_Movie.py +++ b/examples/Fragmentation/Fragmentation_Movie.py @@ -59,29 +59,29 @@ } vel_vectors = {"disruption_headon" : [np.array([ 0.00, 6.280005, 0.0]), - np.array([ 0.00, -6.280005, 0.0])], + np.array([ 0.00, 4.28, 0.0])], "disruption_off_axis" : [np.array([ 0.00, 6.280005, 0.0]), - np.array([ 0.50, -6.280005, 0.0])], + np.array([ 0.05, 4.28, 0.0])], "supercatastrophic_headon": [np.array([ 0.00, 6.28, 0.0]), - np.array([ 0.00, -6.28, 0.0])], + np.array([ 0.00, 4.28, 0.0])], "supercatastrophic_off_axis": [np.array([ 0.00, 6.28, 0.0]), - np.array([ 0.50, -6.28, 0.0])], + np.array([ 0.05, 4.28, 0.0])], "hitandrun_disrupt" : [np.array([ 0.00, 6.28, 0.0]), np.array([-1.45, -6.28, 0.0])], "hitandrun_pure" : [np.array([ 0.00, 6.28, 0.0]), np.array([-1.52, -6.28, 0.0])], - "merge" : [np.array([ 0.00, 0.0, 0.0]), - np.array([ 0.01, -0.100005, 0.0])] + "merge" : [np.array([ 0.05, 6.28, 0.0]), + np.array([ 0.05, 6.18, 0.0])] } rot_vectors = {"disruption_headon" : [np.array([0.0, 0.0, 0.0]), np.array([0.0, 0.0, 0.0])], - "disruption_off_axis": [np.array([0.0, 0.0, -6.0e3]), - np.array([0.0, 0.0, 1.0e4])], + "disruption_off_axis": [np.array([0.0, 0.0, 0.0]), + np.array([0.0, 0.0, 0.0])], "supercatastrophic_headon": [np.array([0.0, 0.0, 0.0]), np.array([0.0, 0.0, 0.0])], - "supercatastrophic_off_axis": [np.array([0.0, 0.0, -6.0e3]), - np.array([0.0, 0.0, 1.0e4])], + "supercatastrophic_off_axis": [np.array([0.0, 0.0, 0.0]), + np.array([0.0, 0.0, 0.0])], "hitandrun_disrupt" : [np.array([0.0, 0.0, 6.0e3]), np.array([0.0, 0.0, 1.0e4])], "hitandrun_pure" : [np.array([0.0, 0.0, 6.0e3]), @@ -99,13 +99,13 @@ "merge" : [1e-7, 1e-8] } -tstop = {"disruption_headon" : 5.0e-4, - "disruption_off_axis" : 5.0e-4, - "supercatastrophic_headon" : 5.0e-4, - "supercatastrophic_off_axis": 5.0e-4, +tstop = {"disruption_headon" : 2.0e-3, + "disruption_off_axis" : 2.0e-3, + "supercatastrophic_headon" : 2.0e-3, + "supercatastrophic_off_axis": 2.0e-3, "hitandrun_disrupt" : 2.0e-4, "hitandrun_pure" : 2.0e-4, - "merge" : 2.0e-3, + "merge" : 5.0e-3, } density = 3000 * swiftest.AU2M**3 / swiftest.MSun diff --git a/src/collision/collision_generate.f90 b/src/collision/collision_generate.f90 index e5c3035b9..b3ef6303e 100644 --- a/src/collision/collision_generate.f90 +++ b/src/collision/collision_generate.f90 @@ -46,7 +46,7 @@ module subroutine collision_generate_bounce(self, nbody_system, param, t) real(DP), intent(in) :: t !! The time of the collision ! Internals integer(I4B) :: i,j,nimp - real(DP), dimension(NDIM) :: vcom, rnorm + real(DP), dimension(NDIM) :: rcom, vcom, rnorm logical, dimension(:), allocatable :: lmask select type(nbody_system) @@ -70,10 +70,13 @@ module subroutine collision_generate_bounce(self, nbody_system, param, t) nimp = size(impactors%id(:)) do i = 1, nimp j = impactors%id(i) + rcom(:) = pl%rb(:,j) - impactors%rbcom(:) vcom(:) = pl%vb(:,j) - impactors%vbcom(:) - rnorm(:) = .unit. (impactors%rb(:,2) - impactors%rb(:,1)) + rnorm(:) = .unit. rcom(:) ! Do the reflection vcom(:) = vcom(:) - 2 * dot_product(vcom(:),rnorm(:)) * rnorm(:) + ! Shift the positions so that the collision doesn't immediately occur again + pl%rb(:,j) = pl%rb(:,j) + 0.5_DP * pl%radius(j) * rnorm(:) pl%vb(:,j) = impactors%vbcom(:) + vcom(:) self%status = DISRUPTED pl%status(j) = ACTIVE @@ -102,7 +105,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..cf30bea7c 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) :: rcom, vcom, rnorm character(len=STRMAX) :: message logical :: lfailure @@ -52,15 +53,31 @@ 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(:) + rcom(:) = pl%rb(:,j) - impactors%rbcom(:) + rnorm(:) = .unit. rcom(:) + ! Do the reflection + vcom(:) = vcom(:) - 2 * dot_product(vcom(:),rnorm(:)) * rnorm(:) + self%fragments%vb(:,i) = impactors%vbcom(:) + vcom(:) + 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) + ! Ensure that the bounce doesn't happen again + self%fragments%rb(:,i) = pl%rb(:,j) + 0.5_DP * self%fragments%radius(i) * rnorm(:) + 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 @@ -287,7 +304,7 @@ module subroutine fraggle_generate_pos_vec(collider) else if (lsupercat) then rdistance = 0.5_DP * sum(impactors%radius(:)) else - rdistance = 10 * impactors%radius(2) + rdistance = 2 * impactors%radius(2) end if ! Give the fragment positions a random value that is scaled with fragment mass so that the more massive bodies tend to be closer to the impact point ! Later, velocities will be scaled such that the farther away a fragment is placed from the impact point, the higher will its velocity be. @@ -318,7 +335,6 @@ module subroutine fraggle_generate_pos_vec(collider) istart = 2 end if - do i = istart, nfrag if (loverlap(i)) then call random_number(phi(i)) @@ -327,8 +343,6 @@ module subroutine fraggle_generate_pos_vec(collider) end if end do - ! Make the fragment cloud symmertic about 0 - do concurrent(i = istart:nfrag, loverlap(i)) j = fragments%origin_body(i) @@ -454,7 +468,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 @@ -462,9 +476,9 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu real(DP), dimension(collider%fragments%nbody) :: vscale, ke_rot_remove ! 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) :: vmax_guess = 4.0_DP real(DP) :: delta_v, volume - integer(I4B), parameter :: MAXLOOP = 100 + integer(I4B), parameter :: MAXLOOP = 200 integer(I4B), parameter :: MAXTRY = 1000 real(DP), parameter :: SUCCESS_METRIC = 1.0e-2_DP class(collision_fraggle), allocatable :: collider_local @@ -504,7 +518,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,8 +557,13 @@ 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) + ke_avail = 0.0_DP + do i = 1, fragments%nbody + ke_avail = ke_avail + 0.5_DP * fragments%mass(i) * max(fragments%vmag(i) - vesc,0.0_DP)**2 + end do + ! Check for any residual angular momentum, and if there is any, put it into spin of the largest body L_residual(:) = collider_local%L_total(:,2) - collider_local%L_total(:,1) if (ke_avail < epsilon(1.0_DP)) then @@ -601,8 +620,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 +654,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 diff --git a/src/swiftest/swiftest_util.f90 b/src/swiftest/swiftest_util.f90 index d79a9b042..0c25abbd5 100644 --- a/src/swiftest/swiftest_util.f90 +++ b/src/swiftest/swiftest_util.f90 @@ -1865,21 +1865,21 @@ module subroutine swiftest_util_rearray_pl(self, nbody_system, param) call plplenc_old%copy(nbody_system%plpl_encounter) end if - ! Re-build the encounter list - ! Be sure to get the level info if this is a SyMBA nbody_system - select type(nbody_system) - class is (symba_nbody_system) - select type(pl) - class is (symba_pl) - select type(tp) - class is (symba_tp) - lencounter = pl%encounter_check(param, nbody_system, param%dt, nbody_system%irec) - if (tp%nbody > 0) then - lencounter = tp%encounter_check(param, nbody_system, param%dt, nbody_system%irec) - end if - end select - end select - end select + ! ! Re-build the encounter list + ! ! Be sure to get the level info if this is a SyMBA nbody_system + ! select type(nbody_system) + ! class is (symba_nbody_system) + ! select type(pl) + ! class is (symba_pl) + ! select type(tp) + ! class is (symba_tp) + ! lencounter = pl%encounter_check(param, nbody_system, param%dt, nbody_system%irec) + ! if (tp%nbody > 0) then + ! lencounter = tp%encounter_check(param, nbody_system, param%dt, nbody_system%irec) + ! end if + ! end select + ! end select + ! end select ! Re-index the encounter list as the index values may have changed if (nenc_old > 0) then