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

Commit

Permalink
Browse files Browse the repository at this point in the history 
Added back in output info to help troubleshoot problems keeping angular momentum and energy under control in fragmentations
  • Loading branch information
@daminton
daminton committed May 10, 2021
1 parent 8c09120 commit b88dc07
Show file tree
Hide file tree
Showing 5 changed files with 57 additions and 43 deletions.
2 changes: 1 addition & 1 deletion examples/symba_mars_disk/param.in
View file
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
!Parameter file for the SyMBA-RINGMOONS test
T0 0.0
TSTOP 31536000000.0
TSTOP 3e3
DT 500.0
PL_IN mars.in
TP_IN tp.in
Expand Down
1 change: 1 addition & 0 deletions src/main/swiftest_symba.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -70,6 +70,7 @@ program swiftest_symba
ntp => symba_tpA%helio%swiftest%nbody)

call util_version
call random_seed()

call get_command_argument(1, inparfile, status = ierr)
if (ierr /= 0) then
Expand Down
2 changes: 1 addition & 1 deletion src/symba/symba_casedisruption.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -32,7 +32,7 @@ function symba_casedisruption (symba_plA, family, nmergeadd, mergeadd_list, x, v
logical :: lmerge

! Collisional fragments will be uniformly distributed around the pre-impact barycenter
nfrag = 12 ! This value is set for testing. This needs to be updated such that it is calculated or set by the user
nfrag = 22 ! This value is set for testing. This needs to be updated such that it is calculated or set by the user
allocate(m_frag(nfrag))
allocate(rad_frag(nfrag))
allocate(xb_frag(NDIM, nfrag))
Expand Down
2 changes: 1 addition & 1 deletion src/symba/symba_casesupercatastrophic.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -32,7 +32,7 @@ function symba_casesupercatastrophic (symba_plA, family, nmergeadd, mergeadd_lis
logical :: lmerge

! Collisional fragments will be uniformly distributed around the pre-impact barycenter
nfrag = 12 ! This value is set for testing. This needs to be updated such that it is calculated or set by the user
nfrag = 42 ! This value is set for testing. This needs to be updated such that it is calculated or set by the user
allocate(m_frag(nfrag))
allocate(rad_frag(nfrag))
allocate(xb_frag(NDIM, nfrag))
Expand Down
93 changes: 53 additions & 40 deletions src/symba/symba_frag_pos.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -92,32 +92,32 @@ subroutine symba_frag_pos (param, symba_plA, family, x, v, L_spin, Ip, mass, rad
Etot_after = ke_after + ke_spin_after + pe_after
Lmag_after = norm2(Ltot_after(:))

!write(*, "(' ---------------------------------------------------------------------------')")
!write(*, "(' Energy normalized by |Etot_before|')")
!write(*, "(' | T_orb T_spin T pe Etot Ltot')")
!write(*, "(' ---------------------------------------------------------------------------')")
!write(*, "(' ---------------------------------------------------------------------------')")
!write(*, "(' First pass to get angular momentum ')")
!write(*, "(' ---------------------------------------------------------------------------')")
!write(*,fmtlabel) ' change |',(ke_after - ke_before) / abs(Etot_before), &
! (ke_spin_after - ke_spin_before)/ abs(Etot_before), &
! (ke_after + ke_spin_after - ke_before - ke_spin_before)/ abs(Etot_before), &
! (pe_after - pe_before) / abs(Etot_before), &
! (Etot_after - Etot_before) / abs(Etot_before), &
! (Lmag_after - Lmag_before) / Lmag_before
!write(*, "(' ---------------------------------------------------------------------------')")
!write(*, "(' Second pass to get energy ')")
!write(*, "(' ---------------------------------------------------------------------------')")
!write(*,fmtlabel) ' Q_loss |',-Qloss / abs(Etot_before)
!write(*, "(' ---------------------------------------------------------------------------')")
write(*, "(' ---------------------------------------------------------------------------')")
write(*, "(' Energy normalized by |Etot_before|')")
write(*, "(' | T_orb T_spin T pe Etot Ltot')")
write(*, "(' ---------------------------------------------------------------------------')")
write(*, "(' ---------------------------------------------------------------------------')")
write(*, "(' First pass to get angular momentum ')")
write(*, "(' ---------------------------------------------------------------------------')")
write(*,fmtlabel) ' change |',(ke_after - ke_before) / abs(Etot_before), &
(ke_spin_after - ke_spin_before)/ abs(Etot_before), &
(ke_after + ke_spin_after - ke_before - ke_spin_before)/ abs(Etot_before), &
(pe_after - pe_before) / abs(Etot_before), &
(Etot_after - Etot_before) / abs(Etot_before), &
(Lmag_after - Lmag_before) / Lmag_before
write(*, "(' ---------------------------------------------------------------------------')")
write(*, "(' Second pass to get energy ')")
write(*, "(' ---------------------------------------------------------------------------')")
write(*,fmtlabel) ' Q_loss |',-Qloss / abs(Etot_before)
write(*, "(' ---------------------------------------------------------------------------')")


! Set the "target" ke_after (the value of the orbital kinetic energy that the fragments ought to have)
ke_target = ke_family + (ke_spin_before - ke_spin_after) + (pe_before - pe_after) - Qloss
call symba_frag_pos_kinetic_energy(xcom, vcom, L_orb, m_frag, x_frag, v_frag, ke_target, lmerge)

!write(*, "(' ---------------------------------------------------------------------------')")
!write(*,fmtlabel) ' T_frag targ |',ke_target / abs(Etot_before)
write(*, "(' ---------------------------------------------------------------------------')")
write(*,fmtlabel) ' T_frag targ |',ke_target / abs(Etot_before)

! Shift the fragments into the system barycenter frame
do i = 1, nfrag
Expand All @@ -130,6 +130,15 @@ subroutine symba_frag_pos (param, symba_plA, family, x, v, L_spin, Ip, mass, rad
nfrag=nfrag, Ip_frag=Ip_frag, m_frag=m_frag, rad_frag=rad_frag, xb_frag=xb_frag, vb_frag=vb_frag, rot_frag=rot_frag)
Etot_after = ke_after + ke_spin_after + pe_after
Lmag_after = norm2(Ltot_after(:))


write(*, "(' ---------------------------------------------------------------------------')")
write(*,fmtlabel) ' change |',(ke_after - ke_before) / abs(Etot_before), &
(ke_spin_after - ke_spin_before)/ abs(Etot_before), &
(ke_after + ke_spin_after - ke_before - ke_spin_before)/ abs(Etot_before), &
(pe_after - pe_before) / abs(Etot_before), &
(Etot_after - Etot_before) / abs(Etot_before), &
(Lmag_after - Lmag_before) / Lmag_before

lmerge = lmerge .or. ((Etot_after - Etot_before) / abs(Etot_before) > 0._DP)
if (.not.lmerge) then
Expand All @@ -139,17 +148,23 @@ subroutine symba_frag_pos (param, symba_plA, family, x, v, L_spin, Ip, mass, rad
rot_frag(:,i) = rot_frag(:,i) + L_spin_frag(:) / (Ip_frag(3, i) * m_frag(i) * rad_frag(i)**2)
end do
end if
call symba_frag_pos_energy_calc(npl, symba_plA, lexclude, ke_after, ke_spin_after, pe_after, Ltot_after,&
nfrag=nfrag, Ip_frag=Ip_frag, m_frag=m_frag, rad_frag=rad_frag, xb_frag=xb_frag, vb_frag=vb_frag, rot_frag=rot_frag)
Etot_after = ke_after + ke_spin_after + pe_after
Lmag_after = norm2(Ltot_after(:))

! write(*, "(' ---------------------------------------------------------------------------')")
! write(*, "(' | T_orb T_spin T pe Etot Ltot')")
! write(*, "(' ---------------------------------------------------------------------------')")
! write(*,fmtlabel) ' change |',(ke_after - ke_before) / abs(Etot_before), &
! (ke_spin_after - ke_spin_before)/ abs(Etot_before), &
! (ke_after + ke_spin_after - ke_before - ke_spin_before)/ abs(Etot_before), &
! (pe_after - pe_before) / abs(Etot_before), &
! (Etot_after - Etot_before) / abs(Etot_before), &
! (Lmag_after - Lma_before) / Lmag_before
! write(*, "(' ---------------------------------------------------------------------------')")
write(*, "(' ---------------------------------------------------------------------------')")
write(*, "(' Third pass for correcting any residual angular momentum ')")
write(*, "(' ---------------------------------------------------------------------------')")
!write(*, "(' | T_orb T_spin T pe Etot Ltot')")
!write(*, "(' ---------------------------------------------------------------------------')")
write(*,fmtlabel) ' change |',(ke_after - ke_before) / abs(Etot_before), &
(ke_spin_after - ke_spin_before)/ abs(Etot_before), &
(ke_after + ke_spin_after - ke_before - ke_spin_before)/ abs(Etot_before), &
(pe_after - pe_before) / abs(Etot_before), &
(Etot_after - Etot_before) / abs(Etot_before), &
(Lmag_after - Lmag_before) / Lmag_before
write(*, "(' ---------------------------------------------------------------------------')")
!****************************************************************************************************************

end associate
Expand Down Expand Up @@ -213,15 +228,13 @@ subroutine symba_frag_pos_initialize_fragments(nfrag, xcom, vcom, x, v, L_orb, L

! Re-normalize position and velocity vectors by the fragment number so that for our initial guess we weight each
! fragment position by the mass and assume equipartition of energy for the velocity
r_col_norm = 1.5_DP * 4 * sum(rad_frag(:)) / theta / (nfrag - 2)
r_col_norm = 1.5_DP * 4 * sum(rad_frag(:)) / theta / (nfrag - 1)

! We will treat the first two fragments of the list as special cases. They get placed on the "poles" of the
! fragment distribution (aligned with the collisional system z-axis)
x_frag(:, 1) = - z_col_unit(:)
x_frag(:, 2) = z_col_unit(:)
call random_number(x_frag(:,3:nfrag))

x_frag(:,:) = x_frag(:,:) * r_col_norm
! We will treat the first fragment of the list as a special case.
x_frag(:, 1) = -z_col_unit(:)
call random_number(x_frag(:,2:nfrag))

x_frag(:, :) = x_frag(:, :) * r_col_norm
! The rest of the fragments will go in a ring on the x-y plane
!do i = 3, nfrag
! frag_vec(:) = [b2a * cos(theta * (i - 1)), sin(theta * (i - 1))]
Expand Down Expand Up @@ -324,7 +337,7 @@ subroutine symba_frag_pos_kinetic_energy(xcom, vcom, L_orb, m_frag, x_frag, v_fr
C = 2 * ke_target
A = 0._DP
do i = 1, nfrag
A = A + m_frag(i) * (mtot / m_frag(i))**2
A = A + m_frag(i)
C = C - m_frag(i) * (dot_product(vcom(:), vcom(:)) + dot_product(v_t(:,i), v_t(:,i)) + dot_product(vcom(:), v_t(:, i)))
end do

Expand All @@ -339,7 +352,7 @@ subroutine symba_frag_pos_kinetic_energy(xcom, vcom, L_orb, m_frag, x_frag, v_fr

! Shift the fragments into the system barycenter frame
do i = 1, nfrag
v_frag(:,i) = v_corrected * mtot / m_frag(i) * v_r_unit(:, i) + v_t(:, i)
v_frag(:,i) = v_corrected * v_r_unit(:, i) + v_t(:, i)
end do

call symba_frag_pos_com_adjust(vcom, m_frag, v_frag)
Expand Down

0 comments on commit b88dc07

Please sign in to comment.