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

Commit

Permalink
Fixed a whole mess of bugs and improved Fraggle's ability to converge…
Browse files Browse the repository at this point in the history 
… on the energy constraint
  • Loading branch information
@daminton
daminton committed Dec 30, 2022
1 parent 7b1b9cd commit 83d5490
Show file tree
Hide file tree
Showing 3 changed files with 72 additions and 58 deletions.
3 changes: 1 addition & 2 deletions examples/Fragmentation/Fragmentation_Movie.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -152,8 +152,7 @@ def __init__(self, sim, animfile, title, style, nskip=1):
self.fig, self.ax = self.setup_plot()

# Then setup FuncAnimation.
self.ani = animation.FuncAnimation(self.fig, self.update_plot, interval=1, frames=range(0,nframes,nskip),
blit=True)
self.ani = animation.FuncAnimation(self.fig, self.update_plot, interval=1, frames=range(0,nframes,nskip), blit=True)
self.ani.save(animfile, fps=60, dpi=300, extra_args=['-vcodec', 'libx264'])
print(f"Finished writing {animfile}")

Expand Down
64 changes: 39 additions & 25 deletions src/fraggle/fraggle_generate.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -133,7 +133,7 @@ module subroutine fraggle_generate_disrupt(self, nbody_system, param, t, lfailur
end if
call ieee_set_flag(ieee_all, .false.) ! Set all fpe flags to quiet

!call self%set_natural_scale()
call self%set_natural_scale()

call fragments%reset()

Expand All @@ -157,11 +157,10 @@ module subroutine fraggle_generate_disrupt(self, nbody_system, param, t, lfailur
call fraggle_generate_rot_vec(self)
call fraggle_generate_vel_vec(self)
call self%get_energy_and_momentum(nbody_system, param, lbefore=.false.)
exit
dEtot = self%Etot(2) - self%Etot(1)
dLmag = .mag. (self%Ltot(:,2) - self%Ltot(:,1))

lfailure_local = (dEtot > FRAGGLE_ETOL)
lfailure_local = (dEtot > -impactors%Qloss)
if (lfailure_local) then
write(message, *) "dEtot: ",dEtot, "dEtot/Qloss", dEtot / impactors%Qloss
call swiftest_io_log_one_message(COLLISION_LOG_OUT, "Fraggle failed due to energy gain: " // &
Expand Down Expand Up @@ -195,7 +194,7 @@ module subroutine fraggle_generate_disrupt(self, nbody_system, param, t, lfailur
trim(adjustl(message)) // " tries")
end if

!call self%set_original_scale()
call self%set_original_scale()

! Restore the big array
if (lk_plpl) call pl%flatten(param)
Expand Down Expand Up @@ -234,7 +233,7 @@ module subroutine fraggle_generate_hitandrun(self, nbody_system, param, t)
class is (swiftest_nbody_system)
select type(pl => nbody_system%pl)
class is (swiftest_pl)
associate(impactors => self%impactors)
associate(impactors => self%impactors, nfrag => self%fragments%nbody)
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 @@ -310,7 +309,7 @@ module subroutine fraggle_generate_pos_vec(collider)
logical :: lcat, lhitandrun
integer(I4B), parameter :: MAXLOOP = 10000
real(DP) :: rdistance
real(DP), parameter :: fail_scale = 1.1_DP ! Scale factor to apply to cloud radius and distance if cloud generation fails
real(DP), parameter :: fail_scale = 1.01_DP ! Scale factor to apply to cloud radius and distance if cloud generation fails


associate(fragments => collider%fragments, impactors => collider%impactors, nfrag => collider%fragments%nbody)
Expand Down Expand Up @@ -431,13 +430,14 @@ module subroutine fraggle_generate_vel_vec(collider)
! Arguments
class(collision_fraggle), intent(inout) :: collider !! Fraggle collision system object
! Internals
integer(I4B) :: i,j, loop, istart
integer(I4B) :: i, j, loop, istart, n, ndof
logical :: lhitandrun, lnoncat
real(DP), dimension(NDIM) :: vimp_unit, rimp, vrot, Lresidual
real(DP) :: vmag, vesc, rotmag, ke_residual, ke_per_dof
real(DP) :: vmag, vesc, rotmag, ke_residual, ke_per_dof, ke_tot
integer(I4B), dimension(collider%fragments%nbody) :: vsign
real(DP), dimension(collider%fragments%nbody) :: vscale, mass_vscale, ke_avail
integer(I4B), parameter :: MAXLOOP = 100
real(DP), parameter :: TOL = 1e-4

associate(fragments => collider%fragments, impactors => collider%impactors, nfrag => collider%fragments%nbody)
lhitandrun = (impactors%regime == COLLRESOLVE_REGIME_HIT_AND_RUN)
Expand Down Expand Up @@ -501,30 +501,44 @@ module subroutine fraggle_generate_vel_vec(collider)
call collider%set_coordinate_system()
call fragments%get_kinetic_energy()
ke_residual = fragments%ke_budget - (fragments%ke_orbit + fragments%ke_spin)
if (abs(ke_residual) <= fragments%ke_budget * TOL) exit
! Make sure we don't take away too much orbital kinetic energy, otherwise the fragment can't escape
ke_avail(:) = fragments%ke_orbit_frag(:) - impactors%Gmass(1)*impactors%mass(2)/fragments%vmag(:)
ke_per_dof = -ke_residual/((nfrag - istart + 1))
ke_avail(:) = fragments%ke_orbit_frag(:) - impactors%Gmass(1)*impactors%mass(2)/fragments%rmag(:)
ke_tot = 0.0_DP
ke_per_dof = -ke_residual
do i = 1, 2*(nfrag - istart + 1)
n = count(ke_avail(istart:nfrag) > -ke_residual/i) + count(fragments%ke_spin_frag(istart:nfrag) > -ke_residual/i)
if (abs(n * ke_per_dof) > ke_tot) then
ke_per_dof = -ke_residual/i
ke_tot = n * ke_per_dof
ndof = i
if (abs(ke_tot) > abs(ke_residual)) then
ke_tot = -ke_residual
ke_per_dof = ke_tot/n
exit
end if
end if
end do
do concurrent(i = istart:nfrag, ke_avail(i) > ke_per_dof)
fragments%vmag(i) = sqrt(2 * (fragments%ke_orbit_frag(i) - ke_per_dof)/fragments%mass(i))
fragments%vc(:,i) = fragments%vmag(i) * fragments%v_unit(:,i)
vmag = max(fragments%vmag(i)**2 - 2*ke_per_dof/fragments%mass(i),vesc**2)
fragments%vmag(i) = sqrt(vmag)
fragments%vc(:,i) = fragments%vmag(i) * .unit.fragments%vc(:,i)
end do
do concurrent(i = istart:nfrag, fragments%ke_spin_frag(i) > ke_per_dof)
rotmag = fragments%rotmag(i)**2 - 2*ke_per_dof/(fragments%mass(i) * fragments%radius(i)**2 * fragments%Ip(3,i))
rotmag = max(rotmag, 0.0_DP)
fragments%rotmag(i) = sqrt(rotmag)
fragments%rot(:,i) = rotmag * .unit.fragments%rot(:,i)
end do
! do concurrent(i = istart:nfrag, fragments%ke_spin_frag(i) > ke_per_dof)
! fragments%rotmag(i) = sqrt(2 * (fragments%ke_spin_frag(i) - ke_per_dof)/(fragments%mass(i) * fragments%radius(i)**2 * fragments%Ip(3,i)))
! fragments%rot(:,i) = fragments%rotmag(i) * .unit.fragments%rot(:,i)
! end do
call fragments%get_kinetic_energy()
ke_residual = fragments%ke_budget - (fragments%ke_orbit + fragments%ke_spin)

! Check for any residual angular momentum, and if there is any, put it into spin of the biggest body
! Check for any residual angular momentum, and if there is any, put it into spin
call collider%set_coordinate_system()
call fragments%get_angular_momentum()
Lresidual(:) = fragments%L_budget(:) - (fragments%Lorbit(:) + fragments%Lspin(:))
rotmag = .mag. fragments%rot(:,1)
fragments%rot(:,1) = fragments%rot(:,1) + Lresidual(:) / (fragments%mass(1) * fragments%radius(1)**2 * fragments%Ip(:,1))
rotmag = .mag. fragments%rot(:,1)

if (ke_residual >= 0.0_DP) exit

do concurrent(i = 1:nfrag)
rotmag = .mag. fragments%rot(:,i)
fragments%rot(:,i) = fragments%rot(:,i) + Lresidual(:) / (nfrag*fragments%mass(i) * fragments%radius(i)**2 * fragments%Ip(:,i))
end do
end do

do concurrent(i = 1:nfrag)
Expand Down
63 changes: 32 additions & 31 deletions src/fraggle/fraggle_util.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -231,39 +231,40 @@ module subroutine fraggle_util_set_natural_scale_factors(self)
integer(I4B) :: i

associate(collider => self, fragments => self%fragments, impactors => self%impactors)
! Set scale factors
collider%Escale = 0.5_DP * ( impactors%mass(1) * dot_product(impactors%vb(:,1), impactors%vb(:,1)) &
+ impactors%mass(2) * dot_product(impactors%vb(:,2), impactors%vb(:,2)))
collider%dscale = sum(impactors%radius(:))
! Set primary scale factors (mass, length, and time) based on the impactor properties at the time of collision
collider%mscale = fragments%mtot
collider%vscale = sqrt(collider%Escale / collider%mscale)
collider%tscale = collider%dscale / collider%vscale
collider%dscale = sum(impactors%radius(:))
collider%tscale = collider%dscale / (.mag.(impactors%vc(:,2) - impactors%vc(:,1)))

! Set secondary scale factors for convenience
collider%vscale = collider%dscale / collider%tscale
collider%Escale = collider%mscale * collider%vscale**2
collider%Lscale = collider%mscale * collider%dscale * collider%vscale

! Scale all dimensioned quantities of impactors and fragments
impactors%rbcom(:) = impactors%rbcom(:) / collider%dscale
impactors%vbcom(:) = impactors%vbcom(:) / collider%vscale
impactors%rbimp(:) = impactors%rbimp(:) / collider%dscale
impactors%bounce_unit(:) = impactors%bounce_unit(:) / collider%vscale
impactors%rb(:,:) = impactors%rb(:,:) / collider%dscale
impactors%vb(:,:) = impactors%vb(:,:) / collider%vscale
impactors%rc(:,:) = impactors%rc(:,:) / collider%dscale
impactors%vc(:,:) = impactors%vc(:,:) / collider%vscale
impactors%mass(:) = impactors%mass(:) / collider%mscale
impactors%Gmass(:) = impactors%Gmass(:) / 0.5_DP * collider%vscale**2 * collider%dscale
impactors%Mcb = impactors%Mcb / collider%mscale
impactors%radius(:) = impactors%radius(:) / collider%dscale
impactors%Lspin(:,:) = impactors%Lspin(:,:) / collider%Lscale
impactors%Lorbit(:,:) = impactors%Lorbit(:,:) / collider%Lscale
impactors%rbcom(:) = impactors%rbcom(:) / collider%dscale
impactors%vbcom(:) = impactors%vbcom(:) / collider%vscale
impactors%rbimp(:) = impactors%rbimp(:) / collider%dscale
impactors%rb(:,:) = impactors%rb(:,:) / collider%dscale
impactors%vb(:,:) = impactors%vb(:,:) / collider%vscale
impactors%rc(:,:) = impactors%rc(:,:) / collider%dscale
impactors%vc(:,:) = impactors%vc(:,:) / collider%vscale
impactors%mass(:) = impactors%mass(:) / collider%mscale
impactors%Gmass(:) = impactors%Gmass(:) / (collider%dscale**3/collider%tscale**2)
impactors%Mcb = impactors%Mcb / collider%mscale
impactors%radius(:) = impactors%radius(:) / collider%dscale
impactors%Lspin(:,:) = impactors%Lspin(:,:) / collider%Lscale
impactors%Lorbit(:,:) = impactors%Lorbit(:,:) / collider%Lscale
impactors%bounce_unit(:) = impactors%bounce_unit(:) / collider%vscale

do i = 1, 2
impactors%rot(:,i) = impactors%Lspin(:,i) / (impactors%mass(i) * impactors%radius(i)**2 * impactors%Ip(3, i))
impactors%rot(:,i) = impactors%Lspin(:,i) / (impactors%mass(i) * impactors%radius(i)**2 * impactors%Ip(3,i))
end do

fragments%mtot = fragments%mtot / collider%mscale
fragments%mass = fragments%mass / collider%mscale
fragments%radius = fragments%radius / collider%dscale
impactors%Qloss = impactors%Qloss / collider%Escale
fragments%mtot = fragments%mtot / collider%mscale
fragments%mass = fragments%mass / collider%mscale
fragments%radius = fragments%radius / collider%dscale
impactors%Qloss = impactors%Qloss / collider%Escale
end associate

return
Expand All @@ -288,13 +289,13 @@ module subroutine fraggle_util_set_original_scale_factors(self)
associate(collider => self, fragments => self%fragments, impactors => self%impactors)

! Restore scale factors
impactors%rbcom(:) = impactors%rbcom(:) * collider%dscale
impactors%vbcom(:) = impactors%vbcom(:) * collider%vscale
impactors%rbimp(:) = impactors%rbimp(:) * collider%dscale
impactors%rbcom(:) = impactors%rbcom(:) * collider%dscale
impactors%vbcom(:) = impactors%vbcom(:) * collider%vscale
impactors%rbimp(:) = impactors%rbimp(:) * collider%dscale
impactors%bounce_unit(:) = impactors%bounce_unit(:) * collider%vscale

impactors%mass = impactors%mass * collider%mscale
impactors%Gmass(:) = impactors%Gmass(:) * 0.5_DP * collider%vscale**2 * collider%dscale
impactors%Gmass(:) = impactors%Gmass(:) * (collider%dscale**3/collider%tscale**2)
impactors%Mcb = impactors%Mcb * collider%mscale
impactors%mass_dist = impactors%mass_dist * collider%mscale
impactors%radius = impactors%radius * collider%dscale
Expand All @@ -305,7 +306,7 @@ module subroutine fraggle_util_set_original_scale_factors(self)
impactors%Lspin = impactors%Lspin * collider%Lscale
impactors%Lorbit = impactors%Lorbit * collider%Lscale
do i = 1, 2
impactors%rot(:,i) = impactors%Lspin(:,i) * (impactors%mass(i) * impactors%radius(i)**2 * impactors%Ip(3, i))
impactors%rot(:,i) = impactors%Lspin(:,i) * (impactors%mass(i) * impactors%radius(i)**2 * impactors%Ip(3,i))
end do

fragments%mtot = fragments%mtot * collider%mscale
Expand Down

0 comments on commit 83d5490

Please sign in to comment.