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 Feb 14, 2023
2 parents 4021d3f + e220a6f commit 91ffbe2
Show file tree
Hide file tree
Showing 3 changed files with 33 additions and 23 deletions.
12 changes: 11 additions & 1 deletion examples/Fragmentation/Fragmentation_Movie.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -124,6 +124,16 @@
"merge_spinner" : 5.0e-3,
}

nfrag_reduction = {"disruption_headon" : 10.0,
"disruption_off_axis" : 10.0,
"supercatastrophic_headon" : 10.0,
"supercatastrophic_off_axis" : 10.0,
"hitandrun_disrupt" : 1.0,
"hitandrun_pure" : 1.0,
"merge" : 1.0,
"merge_spinner" : 1.0,
}

density = 3000 * swiftest.AU2M**3 / swiftest.MSun
GU = swiftest.GMSun * swiftest.YR2S**2 / swiftest.AU2M**3
body_radius = body_Gmass.copy()
Expand Down Expand Up @@ -348,7 +358,7 @@ def vec_props(self, c):
# Set fragmentation parameters
minimum_fragment_gmass = 0.01 * body_Gmass[style][1]
gmtiny = 0.10 * body_Gmass[style][1]
sim.set_parameter(collision_model="fraggle", encounter_save="both", gmtiny=gmtiny, minimum_fragment_gmass=minimum_fragment_gmass, nfrag_reduction=10.0)
sim.set_parameter(collision_model="fraggle", encounter_save="both", gmtiny=gmtiny, minimum_fragment_gmass=minimum_fragment_gmass, nfrag_reduction=nfrag_reduction[style])
sim.run(dt=5e-4, tstop=tstop[style], istep_out=1, dump_cadence=0)

print("Generating animation")
Expand Down
29 changes: 15 additions & 14 deletions src/collision/collision_regime.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -187,7 +187,7 @@ subroutine collision_regime_LS12_SI(Mcb, m1, m2, rad1, rad2, rh1, rh2, vb1, vb2,
real(DP) :: Qr, Qrd_pstar, Qr_erosion, Qr_supercat
real(DP) :: Vhr, Verosion, Vescp, Vhill, Vimp, Vsupercat
real(DP) :: Mint, Mtot, Mtmp, Mbig, Msmall
real(DP) :: Rp, rhill
real(DP) :: Rc1, rhill
real(DP) :: Mresidual
real(DP) :: U_binding

Expand All @@ -209,14 +209,14 @@ subroutine collision_regime_LS12_SI(Mcb, m1, m2, rad1, rad2, rh1, rh2, vb1, vb2,
alpha = 1.0_DP
Mint = m2
end if
Rp = (3 * (m1 / den1 + alpha * m2 / den2) / (4 * PI))**(1.0_DP/3.0_DP) ! (Mustill et al. 2018)
c_star = calc_c_star(Rp)
Rc1 = (3 * Mtot / (4 * PI * DENSITY1))**THIRD ! Stewart and Leinhardt (2009)
c_star = calc_c_star(Rc1)

!calculate Vescp
Vescp = sqrt(2 * GC * Mtot / Rp) !Mustill et al. 2018 eq 6
Vescp = sqrt(2 * GC * Mtot / Rc1) ! Steward and Leinhardt (2012) e.g. Fig. 7 caption.

!calculate rhill
rhill = a1 * (m1 / 3.0_DP / (Mcb + m1))**(1.0_DP/3.0_DP)
rhill = a1 * (m1 / (3 *(Mcb + m1)))**(1.0_DP/3.0_DP)

!calculate Vhill
if ((rad2 + rad1) < rhill) then
Expand Down Expand Up @@ -245,7 +245,7 @@ subroutine collision_regime_LS12_SI(Mcb, m1, m2, rad1, rad2, rh1, rh2, vb1, vb2,
Vhr = Vescp * (C1 * zeta**2 * theta**(2.5_DP) + C2 * zeta**2 + C3 * theta**(2.5_DP) + C4) ! Kokubo & Genda (2010) eq. (3)
bcrit = rad1 / (rad1 + rad2)
! Specific binding energy
U_binding = (3 * GC * Mtot) / (5 * Rp) ! LS12 eq. 27
U_binding = (3 * GC * Mtot) / (5 * Rc1) ! LS12 eq. 27
ke = 0.5_DP * Vimp**2
pe = - GC * m1 * m2 / (Mtot * norm2(rh2 - rh1))

Expand Down Expand Up @@ -325,17 +325,18 @@ function calc_Qrd_pstar(Mtarg, Mp, alpha, c_star) result(Qrd_pstar)
! Result
real(DP) :: Qrd_pstar
! Internals
real(DP) :: Qrd_star1, mu_alpha, mu, Qrd_star
real(DP) :: Qrd_star1, mu_alpha, mu, Qrd_star, gamma

! calc mu, mu_alpha
mu = (Mtarg * Mp) / (Mtarg + Mp) ! [kg]
mu_alpha = (Mtarg * alpha * Mp) / (Mtarg + alpha * Mp) ! [kg]
gamma = Mp / Mtarg
! calc Qrd_star1
Qrd_star1 = (c_star * 4 * PI * DENSITY1 * GC * Rp**2) / 5.0_DP
Qrd_star1 = (c_star * 4 * PI * DENSITY1 * GC * Rc1**2) / 5.0_DP ! LS12 eq. 28
! calc Qrd_star
Qrd_star = Qrd_star1 * (((Mp / Mtarg + 1.0_DP)**2) / (4 * Mp / Mtarg))**(2.0_DP / (3.0_DP * MU_BAR) - 1.0_DP) !(eq 23)
Qrd_star = Qrd_star1 * (((gamma + 1.0_DP)**2) / (4 * gamma))**(2.0_DP / (3 * MU_BAR) - 1.0_DP) !(eq 23)
! calc Qrd_pstar, v_pstar
Qrd_pstar = ((mu / mu_alpha)**(2.0_DP - 3.0_DP * MU_BAR / 2.0_DP)) * Qrd_star ! (eq 15)
Qrd_pstar = ((mu / mu_alpha)**(2.0_DP - 3 * MU_BAR / 2.0_DP)) * Qrd_star ! (eq 15)

return
end function calc_Qrd_pstar
Expand Down Expand Up @@ -365,8 +366,8 @@ function calc_Qrd_rev(Mp, Mtarg, Mint, den1, den2, Vimp, c_star) result(Mslr)
! calc Qrd_star1, v_star1
Qrd_star1 = (c_star * 4 * PI * mtot_rev * GC ) / Rc1 / 5.0_DP
! calc Qrd_pstar_rev
Qrd_star = Qrd_star1 * (((gamma_rev + 1.0_DP)**2) / (4 * gamma_rev)) ** (2.0_DP / (3.0_DP * MU_BAR) - 1.0_DP) !(eq 52)
Qrd_pstar = Qrd_star * ((mu_rev / mu_alpha_rev)**(2.0_DP - 3.0_DP * MU_BAR / 2.0_DP))
Qrd_star = Qrd_star1 * (((gamma_rev + 1.0_DP)**2) / (4 * gamma_rev)) ** (2.0_DP / (3 * MU_BAR) - 1.0_DP) !(eq 52)
Qrd_pstar = Qrd_star * ((mu_rev / mu_alpha_rev)**(2.0_DP - 3 * MU_BAR / 2.0_DP))
Qrd_pstar_rev = Qrd_pstar * (Vhill / Vescp)**CRUFU !Rufu and Aharaonson eq (3)
!calc Qr_supercat_rev
Qr_supercat_rev = 1.8_DP * Qrd_pstar_rev
Expand Down Expand Up @@ -417,8 +418,8 @@ function calc_c_star(Rc1) result(c_star)
!! Result
real(DP) :: c_star
!! Internals
real(DP), parameter :: loR = 1.0e3_DP ! Lower bound of inteRpolation size (m)
real(DP), parameter :: hiR = 1.0e7_DP ! Upper bound of inteRpolation size (m)
real(DP), parameter :: loR = 1.0e3_DP ! Lower bound of interpolation size (m)
real(DP), parameter :: hiR = 1.0e7_DP ! Upper bound of interpolation size (m)
real(DP), parameter :: loval = 5.0_DP ! Value of C* at lower bound
real(DP), parameter :: hival = 1.9_DP ! Value of C* at upper bound

Expand Down
15 changes: 7 additions & 8 deletions src/fraggle/fraggle_generate.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -533,7 +533,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
real(DP), parameter :: ENERGY_SUCCESS_METRIC = 1.0e-2_DP !! Relative energy error to accept as a success (success also must be energy-losing in addition to being within the metric amount)
real(DP), parameter :: ENERGY_SUCCESS_METRIC = 1.0e-3_DP !! Relative energy error to accept as a success (success also must be energy-losing in addition to being within the metric amount)
real(DP) :: MOMENTUM_SUCCESS_METRIC = 10*epsilon(1.0_DP) !! Relative angular momentum error to accept as a success (should be *much* stricter than energy)
integer(I4B) :: i, j, loop, try, istart, nfrag, nsteps, nsteps_best, posloop
logical :: lhitandrun, lsupercat
Expand All @@ -542,13 +542,12 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu
integer(I4B), dimension(:), allocatable :: vsign
real(DP), dimension(:), allocatable :: vscale
real(DP), dimension(:), allocatable :: dLi_mag
real(DP), parameter :: L_ROT_VEL_RATIO = 0.5_DP ! Ratio of angular momentum to put into rotation relative to velocity shear of fragments
! 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), parameter :: hitandrun_vscale = 0.25_DP
real(DP) :: vmin_guess
real(DP) :: vmax_guess
integer(I4B), parameter :: MAXLOOP = 20
integer(I4B), parameter :: MAXTRY = 10
integer(I4B), parameter :: MAXLOOP = 100
integer(I4B), parameter :: MAXTRY = 100
integer(I4B), parameter :: MAXANGMTM = 1000
class(collision_fraggle), allocatable :: collider_local
character(len=STRMAX) :: message
Expand Down Expand Up @@ -653,12 +652,12 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu

do i = istart,fragments%nbody
dL(:) = -dL1_mag * dLi_mag(i) * L_residual_unit(:)
drot(:) = L_ROT_VEL_RATIO * dL(:) / (fragments%mass(i) * fragments%Ip(3,i) * fragments%radius(i)**2)
drot(:) = dL(:) / (fragments%mass(i) * fragments%Ip(3,i) * fragments%radius(i)**2)
rot_new(:) = fragments%rot(:,i) + drot(:)
if (.mag.rot_new(:) < collider_local%max_rot) then
fragments%rot(:,i) = rot_new(:)
fragments%rotmag(i) = .mag.fragments%rot(:,i)
else ! We would break the spin barrier here. Put less into spin and more into velocity shear.
else ! We would break the spin barrier here. Add a random component of rotation that is less than what would break the limit. The rest will go in velocity shear
call random_number(drot)
call random_number(rn)
drot(:) = (rn * collider_local%max_rot - fragments%rotmag(i)) * 2 * (drot(:) - 0.5_DP)
Expand All @@ -669,7 +668,7 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu
fragments%rot(:,i) = fragments%rotmag(i) * .unit. fragments%rot(:,i)
end if
end if
L_residual(:) = L_residual(:) - drot(:) * fragments%Ip(3,i) * fragments%mass(i) * fragments%radius(i)**2
L_residual(:) = L_residual(:) + drot(:) * fragments%Ip(3,i) * fragments%mass(i) * fragments%radius(i)**2
end do

! Put any remaining residual into velocity shear
Expand All @@ -682,7 +681,7 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu
if (all(dL_metric(:) <= 1.0_DP)) exit angmtm

do i = istart, fragments%nbody
dL(:) = -L_residual(:) * fragments%mass(i) / sum(fragments%mass(istart:fragments%nbody))
dL(:) = -L_residual(:) / (fragments%nbody - istart + 1)
call collision_util_velocity_torque(dL, fragments%mass(i), fragments%rc(:,i), fragments%vc(:,i))
call collision_util_shift_vector_to_origin(fragments%mass, fragments%vc)
fragments%vmag(i) = .mag.fragments%vc(:,i)
Expand Down

0 comments on commit 91ffbe2

Please sign in to comment.