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 
Consolidated energy calculations into the nested subroutine
  • Loading branch information
@daminton
daminton committed May 20, 2021
1 parent 63421c2 commit 8394cf5
Showing 1 changed file with 27 additions and 19 deletions.
46 changes: 27 additions & 19 deletions src/symba/symba_frag_pos.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -125,7 +125,7 @@ subroutine symba_frag_pos(param, symba_plA, family, x, v, L_spin, Ip, mass, radi
real(DP), dimension(NDIM) :: Ltot_before
real(DP), dimension(NDIM) :: Ltot_after
real(DP) :: Etot_before, ke_orb_before, ke_spin_before, pe_before, Lmag_before
real(DP) :: Etot_after, ke_orb_after, ke_spin_after, pe_after, Lmag_after
real(DP) :: Etot_after, ke_orb_after, ke_spin_after, pe_after, Lmag_after, dEtot, dLmag
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
Expand Down Expand Up @@ -166,9 +166,8 @@ subroutine symba_frag_pos(param, symba_plA, family, x, v, L_spin, Ip, mass, radi
end associate
call set_scale_factors()

call calculate_system_energy(ke_orb_before, ke_spin_before, pe_before, Ltot_before, linclude_fragments=.false.)
Lmag_before = norm2(Ltot_before(:))
Etot_before = ke_orb_before + ke_spin_before + pe_before
call calculate_system_energy(linclude_fragments=.false.)


write(*, "(' -------------------------------------------------------------------------------------')")
write(*, "(' Energy normalized by |Etot_before|')")
Expand All @@ -184,9 +183,7 @@ subroutine symba_frag_pos(param, symba_plA, family, x, v, L_spin, Ip, mass, radi

call set_fragment_position_vectors()
call set_fragment_tangential_velocities()
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(:))
call calculate_system_energy(linclude_fragments=.true.)

write(*, "(' -------------------------------------------------------------------------------------')")
write(*, "(' | T_orb T_spin T pe Etot Ltot')")
Expand All @@ -195,18 +192,15 @@ subroutine symba_frag_pos(param, symba_plA, family, x, v, L_spin, Ip, mass, radi
(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
dEtot, dLmag

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(:))
call calculate_system_energy(linclude_fragments=.true.)

if (.not.lmerge) then
write(*, "(' -------------------------------------------------------------------------------------')")
Expand All @@ -223,8 +217,7 @@ subroutine symba_frag_pos(param, symba_plA, family, x, v, L_spin, Ip, mass, radi
(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
dEtot, dLmag
end if

lmerge = lmerge .or. ((Etot_after - Etot_before) / abs(Etot_before) > 0._DP)
Expand All @@ -235,9 +228,9 @@ subroutine symba_frag_pos(param, symba_plA, family, x, v, L_spin, Ip, mass, radi
write(*, "(' Final diagnostic')")
write(*, "(' -------------------------------------------------------------------------------------')")
write(*, "(' dL_tot: Should be very small' )")
write(*,fmtlabel) ' dL_tot |', norm2(Ltot_after - Ltot_before) / Lmag_before
write(*,fmtlabel) ' dL_tot |', dLmag
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) ' dE_tot |', dEtot
write(*,fmtlabel) ' Qloss |', -Qloss / abs(Etot_before)
write(*,fmtlabel) ' dE - Qloss |', (Etot_after - Etot_before + Qloss) / abs(Etot_before)
write(*, "(' -------------------------------------------------------------------------------------')")
Expand Down Expand Up @@ -391,18 +384,18 @@ subroutine define_pre_collisional_family()
return
end subroutine define_pre_collisional_family

subroutine calculate_system_energy(ke_orb, ke_spin, pe, Ltot, linclude_fragments)
subroutine calculate_system_energy(linclude_fragments)
!! Author: David A. Minton
!!
!! Calculates total system energy, including all bodies in the symba_plA list that do not have a corresponding value of the lexclude array that is true
!! and optionally including fragments.
use module_swiftestalloc
implicit none
! Arguments
real(DP), intent(out) :: ke_orb, ke_spin, pe
real(DP), dimension(:), intent(out) :: Ltot
logical, intent(in) :: linclude_fragments
! Internals
real(DP) :: ke_orb, ke_spin, pe
real(DP), dimension(NDIM) :: Ltot
integer(I4B) :: i, npl_new, nplm
logical, dimension(:), allocatable :: ltmp
logical :: lk_plpl
Expand Down Expand Up @@ -471,13 +464,28 @@ subroutine calculate_system_energy(ke_orb, ke_spin, pe, Ltot, linclude_fragments

! Calculate the current fragment energy and momentum balances
if (linclude_fragments) then
Ltot_after(:) = Ltot(:)
Lmag_after = norm2(Ltot(:))
ke_orb_after = ke_orb
ke_spin_after = ke_spin
pe_after = pe
Etot_after = ke_orb_after + ke_spin_after + pe_after
dEtot = (Etot_after - Etot_before) / abs(Etot_before)
dLmag = norm2(Ltot_after(:) - Ltot_before(:)) / Lmag_before
ke_frag = 0._DP
do i = 1, nfrag
ke_frag = ke_frag + 0.5_DP * m_frag(i) * dot_product(vb_frag(:, i), vb_frag(:, i))
end do
ke_target = ke_family + (ke_spin_before - ke_spin) + (pe_before - pe) - Qloss
ke_offset = ke_frag - ke_target
L_offset(:) = Ltot_before(:) - Ltot(:)
else
Ltot_before(:) = Ltot(:)
Lmag_before = norm2(Ltot(:))
ke_orb_before = ke_orb
ke_spin_before = ke_spin
pe_before = pe
Etot_before = ke_orb_before + ke_spin_before + pe_before
end if
end associate
return
Expand Down

0 comments on commit 8394cf5

Please sign in to comment.