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

Commit

Permalink
cleaned up some missing bits from the append utility
Browse files Browse the repository at this point in the history
  • Loading branch information
@daminton
daminton committed Dec 23, 2022
1 parent ca90311 commit 9a3723f
Show file tree
Hide file tree
Showing 6 changed files with 108 additions and 62 deletions.
5 changes: 4 additions & 1 deletion src/collision/collision_generate.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -149,12 +149,15 @@ module subroutine collision_generate_merge(self, nbody_system, param, t)

select type(pl => nbody_system%pl)
class is (swiftest_pl)
! Get coordinate system
call self%set_coordinate_system()

! Generate the merged body as a single fragment
call self%setup_fragments(1)

! Calculate the initial energy of the nbody_system without the collisional family
call self%get_energy_and_momentum(nbody_system, param, lbefore=.true.)

! The new body's metadata will be taken from the largest of the two impactor bodies, so we need
! its index in the main pl structure
ibiggest = impactors%id(maxloc(pl%Gmass(impactors%id(:)), dim=1))
Expand Down
13 changes: 0 additions & 13 deletions src/collision/collision_regime.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -32,10 +32,6 @@ module subroutine collision_regime_impactors(self, nbody_system, param)
select type(param)
class is (swiftest_parameters)

mtot = sum(impactors%mass(:))
impactors%rbcom(:) = (impactors%mass(1) * impactors%rb(:,1) + impactors%mass(2) * impactors%rb(:,2)) / mtot
impactors%vbcom(:) = (impactors%mass(1) * impactors%vb(:,1) + impactors%mass(2) * impactors%vb(:,2)) / mtot

select case(param%collision_model)
case("MERGE")
impactors%regime = COLLRESOLVE_REGIME_MERGE
Expand Down Expand Up @@ -105,15 +101,6 @@ subroutine collision_regime_LS12(impactors, nbody_system, param)
impactors%mass_dist(2) = min(max(mslr, 0.0_DP), mtot)
impactors%mass_dist(3) = min(max(mtot - mlr - mslr, 0.0_DP), mtot)

! Find the center of mass of the collisional system
mtot = sum(impactors%mass(:))
impactors%rbcom(:) = (impactors%mass(1) * impactors%rb(:,1) + impactors%mass(2) * impactors%rb(:,2)) / mtot
impactors%vbcom(:) = (impactors%mass(1) * impactors%vb(:,1) + impactors%mass(2) * impactors%vb(:,2)) / mtot

! Find the point of impact between the two bodies
runit(:) = impactors%rb(:,2) - impactors%rb(:,1)
runit(:) = runit(:) / (.mag. runit(:))
impactors%rbimp(:) = impactors%rb(:,1) + impactors%radius(1) * runit(:)

! Convert quantities back to the system units and save them into the fragment system
impactors%mass_dist(:) = (impactors%mass_dist(:) / param%MU2KG)
Expand Down
10 changes: 9 additions & 1 deletion src/collision/collision_util.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -387,7 +387,7 @@ module subroutine collision_util_set_coordinate_system(self)
! Internals
integer(I4B) :: i
real(DP), dimension(NDIM) :: delta_r, delta_v, Ltot
real(DP) :: L_mag
real(DP) :: L_mag, mtot
real(DP), dimension(NDIM, self%fragments%nbody) :: L_sigma

associate(fragments => self%fragments, impactors => self%impactors, nfrag => self%fragments%nbody)
Expand All @@ -413,6 +413,14 @@ module subroutine collision_util_set_coordinate_system(self)
impactors%x_unit(:) = impactors%y_unit(:) .cross. impactors%z_unit(:)
impactors%v_unit(:) = .unit.delta_v(:)

! Find the center of mass of the collisional system
mtot = sum(impactors%mass(:))
impactors%rbcom(:) = (impactors%mass(1) * impactors%rb(:,1) + impactors%mass(2) * impactors%rb(:,2)) / mtot
impactors%vbcom(:) = (impactors%mass(1) * impactors%vb(:,1) + impactors%mass(2) * impactors%vb(:,2)) / mtot

! Find the point of impact between the two bodies
impactors%rbimp(:) = impactors%rb(:,1) + impactors%radius(1) * impactors%y_unit(:)

! The "bounce" unit vector is the projection of
impactors%vbimp(:) = dot_product(delta_v(:),impactors%x_unit(:)) * impactors%x_unit(:)

Expand Down
124 changes: 81 additions & 43 deletions src/fraggle/fraggle_generate.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -118,7 +118,7 @@ module subroutine fraggle_generate_hitandrun(collider, nbody_system, param, t)
class is (swiftest_nbody_system)
select type(after => collider%after)
class is (swiftest_nbody_system)
associate(impactors => collider%impactors, fragments => collider%fragments, pl => nbody_system%pl)
associate(impactors => collider%impactors,pl => nbody_system%pl)
message = "Hit and run between"
call collision_io_collider_message(nbody_system%pl, impactors%id, message)
call swiftest_io_log_one_message(COLLISION_LOG_OUT, trim(adjustl(message)))
Expand Down Expand Up @@ -150,7 +150,7 @@ module subroutine fraggle_generate_hitandrun(collider, nbody_system, param, t)
call swiftest_io_log_one_message(COLLISION_LOG_OUT, "Should have been a pure hit and run instead")
nfrag = 0
else
nfrag = fragments%nbody
nfrag = collider%fragments%nbody
write(message, *) nfrag
call swiftest_io_log_one_message(COLLISION_LOG_OUT, "Generating " // trim(adjustl(message)) // " fragments")
end if
Expand All @@ -163,9 +163,9 @@ module subroutine fraggle_generate_hitandrun(collider, nbody_system, param, t)
allocate(after%pl, source=before%pl) ! Be sure to save the pl so that snapshots still work
else
ibiggest = impactors%id(maxloc(pl%Gmass(impactors%id(:)), dim=1))
fragments%id(1) = pl%id(ibiggest)
fragments%id(2:nfrag) = [(i, i = param%maxid + 1, param%maxid + nfrag - 1)]
param%maxid = fragments%id(nfrag)
collider%fragments%id(1) = pl%id(ibiggest)
collider%fragments%id(2:nfrag) = [(i, i = param%maxid + 1, param%maxid + nfrag - 1)]
param%maxid = collider%fragments%id(nfrag)
status = HIT_AND_RUN_DISRUPT
call collision_resolve_mergeaddsub(nbody_system, param, t, status)
end if
Expand Down Expand Up @@ -293,44 +293,46 @@ module subroutine fraggle_generate_fragments(collider, nbody_system, param, lfai
call collider%get_energy_and_momentum(nbody_system, param, lbefore=.false.)
call collider%set_budgets()

call fraggle_generate_tan_vel(collider, lfailure)
if (lfailure) then
call swiftest_io_log_one_message(COLLISION_LOG_OUT, "Fraggle failed to find tangential velocities")
cycle
end if

call fraggle_generate_rad_vel(collider, lfailure)
if (lfailure) then
call swiftest_io_log_one_message(COLLISION_LOG_OUT, "Fraggle failed to find radial velocities")
cycle
end if

call fraggle_generate_spins(collider, f_spin, lfailure)
if (lfailure) then
call swiftest_io_log_one_message(COLLISION_LOG_OUT, "Fraggle failed to find spins")
cycle
end if

call collider%get_energy_and_momentum(nbody_system, param, lbefore=.false.)
dEtot = collider%Etot(2) - collider%Etot(1)
dLmag = .mag. (collider%Ltot(:,2) - collider%Ltot(:,1))
exit

lfailure = ((abs(dEtot + impactors%Qloss) > FRAGGLE_ETOL) .or. (dEtot > 0.0_DP))
if (lfailure) then
write(message, *) dEtot, abs(dEtot + impactors%Qloss) / FRAGGLE_ETOL
call swiftest_io_log_one_message(COLLISION_LOG_OUT, "Fraggle failed due to high energy error: " // &
trim(adjustl(message)))
cycle
end if

lfailure = ((abs(dLmag) / (.mag.collider%Ltot(:,1))) > FRAGGLE_LTOL)
if (lfailure) then
write(message,*) dLmag / (.mag.collider%Ltot(:,1))
call swiftest_io_log_one_message(COLLISION_LOG_OUT, "Fraggle failed due to high angular momentum error: " // &
trim(adjustl(message)))
cycle
end if
call fraggle_generate_vel_vec_guess(collider)

! call fraggle_generate_tan_vel(collider, lfailure)
! if (lfailure) then
! call swiftest_io_log_one_message(COLLISION_LOG_OUT, "Fraggle failed to find tangential velocities")
! cycle
! end if

! call fraggle_generate_rad_vel(collider, lfailure)
! if (lfailure) then
! call swiftest_io_log_one_message(COLLISION_LOG_OUT, "Fraggle failed to find radial velocities")
! cycle
! end if

! call fraggle_generate_spins(collider, f_spin, lfailure)
! if (lfailure) then
! call swiftest_io_log_one_message(COLLISION_LOG_OUT, "Fraggle failed to find spins")
! cycle
! end if

! call collider%get_energy_and_momentum(nbody_system, param, lbefore=.false.)
! dEtot = collider%Etot(2) - collider%Etot(1)
! dLmag = .mag. (collider%Ltot(:,2) - collider%Ltot(:,1))
! exit

! lfailure = ((abs(dEtot + impactors%Qloss) > FRAGGLE_ETOL) .or. (dEtot > 0.0_DP))
! if (lfailure) then
! write(message, *) dEtot, abs(dEtot + impactors%Qloss) / FRAGGLE_ETOL
! call swiftest_io_log_one_message(COLLISION_LOG_OUT, "Fraggle failed due to high energy error: " // &
! trim(adjustl(message)))
! cycle
! end if

! lfailure = ((abs(dLmag) / (.mag.collider%Ltot(:,1))) > FRAGGLE_LTOL)
! if (lfailure) then
! write(message,*) dLmag / (.mag.collider%Ltot(:,1))
! call swiftest_io_log_one_message(COLLISION_LOG_OUT, "Fraggle failed due to high angular momentum error: " // &
! trim(adjustl(message)))
! cycle
! end if

! Check if any of the usual floating point exceptions happened, and fail the try if so
call ieee_get_flag(ieee_usual, fpe_flag)
Expand Down Expand Up @@ -442,6 +444,42 @@ subroutine fraggle_generate_pos_vec(collider, r_max_start)
end subroutine fraggle_generate_pos_vec


subroutine fraggle_generate_vel_vec_guess(collider)
!! Author: David A. Minton
!!
!! Generates an initial "guess" for the velocitity vectors
implicit none
! Arguments
class(collision_fraggle), intent(inout) :: collider !! Fraggle collision system object
! Internals
integer(I4B) :: i, j
logical :: lcat
real(DP), dimension(NDIM) :: vimp_unit
real(DP), dimension(NDIM,collider%fragments%nbody) :: vnoise
real(DP), parameter :: VNOISE_MAG = 0.10_DP
real(DP) :: vmag

associate(fragments => collider%fragments, impactors => collider%impactors, nfrag => collider%fragments%nbody)
lcat = (impactors%regime == COLLRESOLVE_REGIME_SUPERCATASTROPHIC)

! "Bounce" the first two bodies
do i = 1,2
fragments%vc(:,i) = impactors%vb(:,i) - impactors%vbcom(:) - 2 * impactors%vbimp(:)
end do

vmag = .mag.impactors%vbcom(:)
vimp_unit(:) = .unit. impactors%vbimp(:)
call random_number(vnoise)
vnoise = (2 * vnoise - 1.0_DP) * vmag
do i = 3, nfrag
vimp_unit(:) = .unit. (fragments%rc(:,i) + impactors%rbcom(:) - impactors%rbimp(:))
fragments%vc(:,i) = vmag * vimp_unit(:) + vnoise(:,i)
end do
end associate
return
end subroutine fraggle_generate_vel_vec_guess


subroutine fraggle_generate_spins(collider, f_spin, lfailure)
!! Author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton
!!
Expand Down
3 changes: 2 additions & 1 deletion src/misc/minimizer_module.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -140,9 +140,10 @@ module subroutine minimize_bfgs(f, N, x0, eps, maxloop, lerr, x1)

call ieee_get_status(original_fpe_status) ! Save the original floating point exception status
call ieee_set_flag(ieee_all, .false.) ! Set all flags to quiet
allocate(fpe_flag(size(ieee_usual)))
if (.not.allocated(fpe_flag)) allocate(fpe_flag(size(ieee_usual)))

lerr = .false.
if (allocated(x1)) deallocate(x1)
allocate(x1, source=x0)
! Initialize approximate Hessian with the identity matrix (i.e. begin with method of steepest descent)
! Get initial gradient and initialize arrays for updated values of gradient and x
Expand Down
15 changes: 12 additions & 3 deletions src/swiftest/swiftest_util.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -235,11 +235,13 @@ module subroutine swiftest_util_append_body(self, source, lsource_mask)
nsrc = source%nbody
nnew = count(lsource_mask(1:nsrc))

call swiftest_util_append(self%info, source%info, nold, nsrc, lsource_mask)
call swiftest_util_append(self%id, source%id, nold, nsrc, lsource_mask)
call swiftest_util_append(self%info, source%info, nold, nsrc, lsource_mask)
call swiftest_util_append(self%lmask, source%lmask, nold, nsrc, lsource_mask)
call swiftest_util_append(self%status, source%status, nold, nsrc, lsource_mask)
call swiftest_util_append(self%ldiscard, source%ldiscard, nold, nsrc, lsource_mask)
call swiftest_util_append(self%lmask, source%lmask, nold, nsrc, lsource_mask)
call swiftest_util_append(self%lencounter, source%lencounter, nold, nsrc, lsource_mask)
call swiftest_util_append(self%lcollision, source%lcollision, nold, nsrc, lsource_mask)
call swiftest_util_append(self%mu, source%mu, nold, nsrc, lsource_mask)
call swiftest_util_append(self%rh, source%rh, nold, nsrc, lsource_mask)
call swiftest_util_append(self%vh, source%vh, nold, nsrc, lsource_mask)
Expand All @@ -250,6 +252,9 @@ module subroutine swiftest_util_append_body(self, source, lsource_mask)
call swiftest_util_append(self%atide, source%atide, nold, nsrc, lsource_mask)
call swiftest_util_append(self%agr, source%agr, nold, nsrc, lsource_mask)
call swiftest_util_append(self%ir3h, source%ir3h, nold, nsrc, lsource_mask)
call swiftest_util_append(self%isperi, source%isperi, nold, nsrc, lsource_mask)
call swiftest_util_append(self%peri, source%peri, nold, nsrc, lsource_mask)
call swiftest_util_append(self%atp, source%atp, nold, nsrc, lsource_mask)
call swiftest_util_append(self%a, source%a, nold, nsrc, lsource_mask)
call swiftest_util_append(self%e, source%e, nold, nsrc, lsource_mask)
call swiftest_util_append(self%inc, source%inc, nold, nsrc, lsource_mask)
Expand Down Expand Up @@ -282,15 +287,19 @@ module subroutine swiftest_util_append_pl(self, source, lsource_mask)
call swiftest_util_append(self%rhill, source%rhill, nold, nsrc, lsource_mask)
call swiftest_util_append(self%renc, source%renc, nold, nsrc, lsource_mask)
call swiftest_util_append(self%radius, source%radius, nold, nsrc, lsource_mask)
call swiftest_util_append(self%density, source%density, nold, nsrc, lsource_mask)
call swiftest_util_append(self%rbeg, source%rbeg, nold, nsrc, lsource_mask)
call swiftest_util_append(self%rend, source%rend, nold, nsrc, lsource_mask)
call swiftest_util_append(self%vbeg, source%vbeg, nold, nsrc, lsource_mask)
call swiftest_util_append(self%density, source%density, nold, nsrc, lsource_mask)
call swiftest_util_append(self%Ip, source%Ip, nold, nsrc, lsource_mask)
call swiftest_util_append(self%rot, source%rot, nold, nsrc, lsource_mask)
call swiftest_util_append(self%k2, source%k2, nold, nsrc, lsource_mask)
call swiftest_util_append(self%Q, source%Q, nold, nsrc, lsource_mask)
call swiftest_util_append(self%tlag, source%tlag, nold, nsrc, lsource_mask)
call swiftest_util_append(self%kin, source%kin, nold, nsrc, lsource_mask)
call swiftest_util_append(self%lmtiny, source%lmtiny, nold, nsrc, lsource_mask)
call swiftest_util_append(self%nplenc, source%nplenc, nold, nsrc, lsource_mask)
call swiftest_util_append(self%ntpenc, source%ntpenc, nold, nsrc, lsource_mask)

if (allocated(self%k_plpl)) deallocate(self%k_plpl)

Expand Down

0 comments on commit 9a3723f

Please sign in to comment.