Skip to content
This repository was archived by the owner on Aug 28, 2024. It is now read-only.

Commit

Permalink
Merge branch 'debug'
Browse files Browse the repository at this point in the history
  • Loading branch information
@daminton
daminton committed Jan 13, 2023
2 parents 3a6d1f0 + 7530265 commit 4c51dd3
Show file tree
Hide file tree
Showing 6 changed files with 103 additions and 79 deletions.
4 changes: 4 additions & 0 deletions cmake/Modules/SetFortranFlags.cmake
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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}"
Expand Down
30 changes: 15 additions & 15 deletions examples/Fragmentation/Fragmentation_Movie.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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]),
Expand All @@ -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
Expand Down
9 changes: 6 additions & 3 deletions src/collision/collision_generate.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand All @@ -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
Expand Down Expand Up @@ -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
!!
Expand Down
57 changes: 27 additions & 30 deletions src/collision/collision_util.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down
52 changes: 36 additions & 16 deletions src/fraggle/fraggle_generate.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand All @@ -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
Expand Down Expand Up @@ -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.
Expand Down Expand Up @@ -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))
Expand All @@ -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)

Expand Down Expand Up @@ -454,17 +468,17 @@ 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
integer(I4B), dimension(collider%fragments%nbody) :: vsign
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
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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
Expand Down
30 changes: 15 additions & 15 deletions src/swiftest/swiftest_util.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down

0 comments on commit 4c51dd3

Please sign in to comment.