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 
Revamped energy diagnostic reporting
  • Loading branch information
@daminton
daminton committed May 20, 2021
1 parent 4cfeaba commit ded6b0d
Showing 1 changed file with 43 additions and 24 deletions.
67 changes: 43 additions & 24 deletions src/symba/symba_frag_pos.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -113,7 +113,7 @@ subroutine symba_frag_pos(param, symba_plA, family, x, v, L_spin, Ip, mass, radi
logical, intent(out) :: lmerge ! Answers the question: Should this have been a merger instead?
real(DP), intent(inout) :: Qloss
! Internals
real(DP), parameter :: f_spin = 0.00_DP !! Fraction of pre-impact orbital angular momentum that is converted to fragment spin
real(DP), parameter :: f_spin = 0.20_DP !! Fraction of pre-impact orbital angular momentum that is converted to fragment spin
real(DP) :: mscale = 1.0_DP, rscale = 1.0_DP, vscale = 1.0_DP, tscale = 1.0_DP, Lscale = 1.0_DP, Escale = 1.0_DP! Scale factors that reduce quantities to O(~1) in the collisional system
real(DP) :: mtot
real(DP), dimension(NDIM) :: xcom, vcom
Expand All @@ -129,7 +129,7 @@ subroutine symba_frag_pos(param, symba_plA, family, x, v, L_spin, Ip, mass, radi
real(DP), dimension(NDIM) :: L_frag_spin, L_frag_tot, L_frag_orb, L_offset
real(DP) :: ke_family, ke_target, ke_frag, ke_offset
real(DP), dimension(NDIM) :: x_col_unit, y_col_unit, z_col_unit
character(len=*), parameter :: fmtlabel = "(A14,10(ES9.2,1X,:))"
character(len=*), parameter :: fmtlabel = "(A14,10(ES11.4,1X,:))"
integer(I4B), dimension(:), allocatable :: seed
integer(I4B) :: nseed

Expand All @@ -156,15 +156,13 @@ subroutine symba_frag_pos(param, symba_plA, family, x, v, L_spin, Ip, mass, radi
Lmag_before = norm2(Ltot_before(:))
Etot_before = ke_orb_before + ke_spin_before + pe_before

write(*, "(' ---------------------------------------------------------------------------')")
write(*, "(' -------------------------------------------------------------------------------------')")
write(*, "(' Energy normalized by |Etot_before|')")
write(*, "(' ---------------------------------------------------------------------------')")
write(*,fmtlabel) ' Qloss |',-Qloss / abs(Etot_before)

call define_coordinate_system()
call define_pre_collisional_family()

write(*, "(' ---------------------------------------------------------------------------')")
write(*, "(' -------------------------------------------------------------------------------------')")
write(*, "(' First pass to get angular momentum ')")

call set_fragment_position_vectors()
Expand All @@ -173,39 +171,60 @@ subroutine symba_frag_pos(param, symba_plA, family, x, v, L_spin, Ip, mass, radi
Etot_after = ke_orb_after + ke_spin_after + pe_after
Lmag_after = norm2(Ltot_after(:))

write(*, "(' ---------------------------------------------------------------------------')")
write(*, "(' | T_orb T_spin T pe Etot Ltot')")
write(*, "(' ---------------------------------------------------------------------------')")
write(*, "(' -------------------------------------------------------------------------------------')")
write(*, "(' | T_orb T_spin T pe Etot Ltot')")
write(*, "(' -------------------------------------------------------------------------------------')")
write(*,fmtlabel) ' change |',(ke_orb_after - ke_orb_before) / abs(Etot_before), &
(ke_spin_after - ke_spin_before)/ abs(Etot_before), &
(ke_orb_after + ke_spin_after - ke_orb_before - ke_spin_before)/ abs(Etot_before), &
(pe_after - pe_before) / abs(Etot_before), &
(Etot_after - Etot_before) / abs(Etot_before), &
norm2(Ltot_after - Ltot_before) / Lmag_before

write(*, "(' ---------------------------------------------------------------------------')")
write(*, "(' Second pass to get energy ')")
write(*, "(' ---------------------------------------------------------------------------')")
write(*, "(' -------------------------------------------------------------------------------------')")
write(*, "(' If T_offset > 0, this will be a merge')")
write(*,fmtlabel) ' T_offset |',ke_offset / abs(Etot_before)


call set_fragment_radial_velocities(lmerge)

call calculate_system_energy(ke_orb_after, ke_spin_after, pe_after, Ltot_after, linclude_fragments=.true.)
Etot_after = ke_orb_after + ke_spin_after + pe_after
Lmag_after = norm2(Ltot_after(:))

write(*,fmtlabel) ' T_frag calc |',ke_frag / abs(Etot_before)
write(*,fmtlabel) ' residual |',(ke_frag - ke_target) / abs(Etot_before)
write(*, "(' ---------------------------------------------------------------------------')")
write(*, "(' | T_orb T_spin T pe Etot Ltot')")
write(*, "(' ---------------------------------------------------------------------------')")
write(*,fmtlabel) ' change |',(ke_orb_after - ke_orb_before) / abs(Etot_before), &
(ke_spin_after - ke_spin_before)/ abs(Etot_before), &
(ke_orb_after + ke_spin_after - ke_orb_before - ke_spin_before)/ abs(Etot_before), &
(pe_after - pe_before) / abs(Etot_before), &
(Etot_after - Etot_before) / abs(Etot_before), &
norm2(Ltot_after - Ltot_before) / Lmag_before
if (.not.lmerge) then
write(*, "(' -------------------------------------------------------------------------------------')")
write(*, "(' Second pass to get energy ')")
write(*, "(' -------------------------------------------------------------------------------------')")
write(*,fmtlabel) ' T_frag targ |',ke_target / abs(Etot_before)
write(*,fmtlabel) ' T_frag calc |',ke_frag / abs(Etot_before)
write(*, "(' T_offset should now be small')")
write(*,fmtlabel) ' T_offset |',ke_offset / abs(Etot_before)
write(*, "(' -------------------------------------------------------------------------------------')")
write(*, "(' | T_orb T_spin T pe Etot Ltot')")
write(*, "(' -------------------------------------------------------------------------------------')")
write(*,fmtlabel) ' change |',(ke_orb_after - ke_orb_before) / abs(Etot_before), &
(ke_spin_after - ke_spin_before)/ abs(Etot_before), &
(ke_orb_after + ke_spin_after - ke_orb_before - ke_spin_before)/ abs(Etot_before), &
(pe_after - pe_before) / abs(Etot_before), &
(Etot_after - Etot_before) / abs(Etot_before), &
norm2(Ltot_after - Ltot_before) / Lmag_before
end if

lmerge = lmerge .or. ((Etot_after - Etot_before) / abs(Etot_before) > 0._DP)
write(*, "(' -------------------------------------------------------------------------------------')")
write(*, "(' Final diagnostic')")
write(*, "(' -------------------------------------------------------------------------------------')")
write(*, "(' dL_tot: Should be very small' )")
write(*,fmtlabel) ' dL_tot |', norm2(Ltot_after - Ltot_before) / Lmag_before
write(*, "(' dE_tot: If not a merge, should be negative and equal to Qloss' )")
write(*,fmtlabel) ' dE_tot |', (Etot_after - Etot_before) / abs(Etot_before)
write(*,fmtlabel) ' Qloss |', -Qloss / abs(Etot_before)
write(*,fmtlabel) ' dE - Qloss |', (Etot_after - Etot_before + Qloss) / abs(Etot_before)
write(*, "(' -------------------------------------------------------------------------------------')")
write(*,*) "Flag as merge?", lmerge
write(*, "(' -------------------------------------------------------------------------------------')")

call restore_scale_factors()

return
Expand Down Expand Up @@ -518,7 +537,7 @@ subroutine set_fragment_position_vectors()
end do
loverlap(:) = .false.
do j = 1, nfrag
do i = j+1, nfrag
do i = j + 1, nfrag
dis = norm2(x_frag(:,j) - x_frag(:,i))
loverlap(i) = loverlap(i) .or. (dis <= (rad_frag(i) + rad_frag(j)))
end do
Expand Down

0 comments on commit ded6b0d

Please sign in to comment.