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 
Refactoring for clarity and adding in additional angular momentum diagnostics
  • Loading branch information
@daminton
daminton committed May 24, 2021
1 parent be56439 commit 4ba2e65
Show file tree
Hide file tree
Showing 5 changed files with 105 additions and 72 deletions.
48 changes: 32 additions & 16 deletions src/io/io_conservation_report.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -15,10 +15,12 @@ module subroutine io_conservation_report(t, symba_plA, npl, j2rp2, j4rp4, param,
type(user_input_parameters), intent(in) :: param !! Input colleciton of user-defined parameters
logical, intent(in) :: lterminal !! Indicates whether to output information to the terminal screen
! Internals
real(DP), dimension(NDIM), save :: Ltot_orig, Ltot_last
real(DP), save :: Eorbit_orig, Mtot_orig, Lmag_orig, ke_orb_last, ke_spin_last, pe_last, Eorbit_last
real(DP) :: ke_orbit, ke_spin, pe, Eorbit
real(DP), dimension(NDIM) :: Ltot_now
real(DP), dimension(NDIM) :: Ltot_now, Lorbit_now, Lspin_now
real(DP), dimension(NDIM), save :: Ltot_orig, Lorbit_orig, Lspin_orig
real(DP), dimension(NDIM), save :: Ltot_last, Lorbit_last, Lspin_last
real(DP), save :: Eorbit_orig, Mtot_orig, Lmag_orig
real(DP), save :: ke_orbit_last, ke_spin_last, pe_last, Eorbit_last
real(DP) :: ke_orbit_now, ke_spin_now, pe_now, Eorbit_now
real(DP) :: Eorbit_error, Etotal_error, Ecoll_error
real(DP) :: Mtot_now, Merror
real(DP) :: Lmag_now, Lerror
Expand All @@ -42,39 +44,53 @@ module subroutine io_conservation_report(t, symba_plA, npl, j2rp2, j4rp4, param,
write(egyiu,egyheader)
end if
end if
call symba_energy_eucl(npl, symba_plA, j2rp2, j4rp4, ke_orbit, ke_spin, pe, Eorbit, Ltot_now)
call symba_energy_eucl(npl, symba_plA, j2rp2, j4rp4, ke_orbit_now, ke_spin_now, pe_now, Lorbit_now, Lspin_now)
Eorbit_now = ke_orbit_now + ke_spin_now + pe_now
Ltot_now(:) = Lorbit_now(:) + Lspin_now(:) + Lescape(:)
Mtot_now = dMcb + sum(mass(2:npl)) + Mcb_initial + Mescape
Ltot_now(:) = Lescape(:) + Ltot_now(:)
if (lfirst) then
Eorbit_orig = Eorbit
Eorbit_orig = Eorbit_now
Mtot_orig = Mtot_now
Lorbit_orig(:) = Lorbit_now(:)
Lspin_orig(:) = Lspin_now(:)
Ltot_orig(:) = Ltot_now(:)
Lmag_orig = norm2(Ltot_orig(:))
lfirst = .false.
end if

write(egyiu,egyfmt) t, Eorbit, Ecollisions, Ltot_now, Mtot_now
write(egyiu,egyfmt) t, Eorbit_now, Ecollisions, Ltot_now, Mtot_now
flush(egyiu)
if (.not.lfirst .and. lterminal) then
Lmag_now = norm2(Ltot_now)
Lerror = norm2(Ltot_now - Ltot_orig) / Lmag_orig
Eorbit_error = (Eorbit - Eorbit_orig) / abs(Eorbit_orig)
Eorbit_error = (Eorbit_now - Eorbit_orig) / abs(Eorbit_orig)
Ecoll_error = Ecollisions / abs(Eorbit_orig)
Etotal_error = (Eorbit - Ecollisions - Eorbit_orig - Euntracked) / abs(Eorbit_orig)
Etotal_error = (Eorbit_now - Ecollisions - Eorbit_orig - Euntracked) / abs(Eorbit_orig)
Merror = (Mtot_now - Mtot_orig) / Mtot_orig
write(*, egytermfmt) Lerror, Ecoll_error, Etotal_error, Merror
if (Ecoll_error > 0.0_DP) then
write(*,*) 'Something has gone wrong! Collisional energy should not be positive!'
write(*,*) 'dke_orbit: ',(ke_orbit - ke_orb_last) / abs(Eorbit_orig)
write(*,*) 'dke_spin : ',(ke_spin - ke_spin_last) / abs(Eorbit_orig)
write(*,*) 'dpe : ',(pe - pe_last) / abs(Eorbit_orig)
write(*,*) 'dke_orbit: ',(ke_orbit_now - ke_orbit_last) / abs(Eorbit_orig)
write(*,*) 'dke_spin : ',(ke_spin_now - ke_spin_last) / abs(Eorbit_orig)
write(*,*) 'dpe : ',(pe_now - pe_last) / abs(Eorbit_orig)
write(*,*)
end if
end if
ke_orb_last = ke_orbit
ke_spin_last = ke_spin
pe_last = pe
Eorbit_last = Eorbit
write(*,*) 'Angular momentum changes since last time'
write(*,*) ' dLorbit : ',(Lorbit_now(:) - Lorbit_last(:)) / Lmag_orig
write(*,*) '|dLorbit|: ',norm2(Lorbit_now(:) - Lorbit_last(:)) / Lmag_orig
write(*,*) ' dLspin : ',(Lspin_now(:) - Lspin_last(:)) / Lmag_orig
write(*,*) '|dLspin| : ',norm2(Lspin_now(:) - Lspin_last(:)) / Lmag_orig
write(*,*) ' dLtot : ',(Ltot_now(:) - Ltot_last(:)) / Lmag_orig
write(*,*) '|dLtot| : ',norm2(Ltot_now(:) - Ltot_last(:)) / Lmag_orig
ke_orbit_last = ke_orbit_now
ke_spin_last = ke_spin_now
pe_last = pe_now
Eorbit_last = Eorbit_now
Lorbit_last(:) = Lorbit_now(:)
Lspin_last(:) = Lspin_now(:)
Ltot_last(:) = Ltot_now(:)
end associate
return

Expand Down
15 changes: 10 additions & 5 deletions src/main/swiftest_symba.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -21,9 +21,9 @@ program swiftest_symba
logical :: lfrag_add, ldiscard_pl, ldiscard_tp
integer(I4B) :: nplm, ntp, ntp0, nsppl, nsptp, iout, idump, iloop, i
integer(I4B) :: nplplenc, npltpenc, nmergeadd, nmergesub
real(DP) :: t, tfrac, tbase, msys, Eorbit_orig
real(DP) :: t, tfrac, tbase, msys
real(DP) :: Ecollision, Eorbit_before, Eorbit_after, ke_orbit_before, ke_spin_before, ke_orbit_after, ke_spin_after, pe_before, pe_after
real(DP), dimension(NDIM) :: Ltot
real(DP), dimension(NDIM) :: Lorbit, Lspin
character(STRMAX) :: inparfile
type(symba_pl) :: symba_plA
type(symba_tp) :: symba_tpA
Expand Down Expand Up @@ -153,7 +153,8 @@ program swiftest_symba
if (param%lenergy) then
call coord_h2b(npl, symba_plA%helio%swiftest, msys)
call io_conservation_report(t, symba_plA, npl, j2rp2, j4rp4, param, lterminal=.false.)
call symba_energy_eucl(npl, symba_plA, j2rp2, j4rp4, ke_orbit_before, ke_spin_before, pe_before, Eorbit_orig, Ltot)
call symba_energy_eucl(npl, symba_plA, j2rp2, j4rp4, ke_orbit_before, ke_spin_before, pe_before, Lorbit, Lspin)

end if
write(*, *) " *************** Main Loop *************** "

Expand All @@ -177,7 +178,10 @@ program swiftest_symba
ldiscard_pl = ldiscard_pl .or. lfrag_add

if (ldiscard_pl .or. ldiscard_tp) then
if (param%lenergy) call symba_energy_eucl(npl, symba_plA, j2rp2, j4rp4, ke_orbit_before, ke_spin_before, pe_before, Eorbit_before, Ltot)
if (param%lenergy) then
call symba_energy_eucl(npl, symba_plA, j2rp2, j4rp4, ke_orbit_before, ke_spin_before, pe_before, Lorbit, Lspin)
Eorbit_before = ke_orbit_before + ke_spin_before + pe_before
end if
call symba_rearray(t, npl, nplm, ntp, nsppl, nsptp, symba_plA, symba_tpA, nmergeadd, mergeadd_list, discard_plA, &
discard_tpA, ldiscard_pl, ldiscard_tp, mtiny, param, discard_l_pl, discard_stat_list)
call io_discard_write_symba(t, mtiny, npl, nsppl, nsptp, nmergesub, symba_plA, &
Expand All @@ -189,7 +193,8 @@ program swiftest_symba
nplm = count(symba_plA%helio%swiftest%mass(1:npl) > mtiny)

if (param%lenergy) then
call symba_energy_eucl(npl, symba_plA, j2rp2, j4rp4, ke_orbit_after, ke_spin_after, pe_after, Eorbit_after, Ltot)
call symba_energy_eucl(npl, symba_plA, j2rp2, j4rp4, ke_orbit_after, ke_spin_after, pe_after, Lorbit, Lspin)
Eorbit_after = ke_orbit_after + ke_spin_after + pe_after
Ecollision = Eorbit_after - Eorbit_before ! Energy change resulting in this collisional event Total running energy offset from collision in this step
symba_plA%helio%swiftest%Ecollisions = symba_plA%helio%swiftest%Ecollisions + Ecollision
end if
Expand Down
6 changes: 3 additions & 3 deletions src/modules/module_interfaces.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1260,16 +1260,16 @@ END SUBROUTINE symba_user_getacch_tp
END INTERFACE

INTERFACE
SUBROUTINE symba_energy_eucl(npl, symba_plA, j2rp2, j4rp4, ke_orbit, ke_spin, pe, te, Ltot)
SUBROUTINE symba_energy_eucl(npl, symba_plA, j2rp2, j4rp4, ke_orbit, ke_spin, pe, Lorbit, Lspin)
USE swiftest_globals
USE swiftest_data_structures
use module_symba
IMPLICIT NONE
integer(I4B), intent(in) :: npl
type(symba_pl), intent(inout) :: symba_plA
real(DP), intent(in) :: j2rp2, j4rp4
real(DP), intent(out) :: ke_orbit, ke_spin, pe, te
real(DP), dimension(:), intent(out) :: Ltot
real(DP), intent(out) :: ke_orbit, ke_spin, pe
real(DP), dimension(:), intent(out) :: Lorbit, Lspin
END SUBROUTINE symba_energy_eucl
END INTERFACE

Expand Down
51 changes: 32 additions & 19 deletions src/symba/symba_energy_eucl.f90
View file
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
subroutine symba_energy_eucl(npl, symba_plA, j2rp2, j4rp4, ke_orbit, ke_spin, pe, te, Ltot)
subroutine symba_energy_eucl(npl, symba_plA, j2rp2, j4rp4, ke_orbit, ke_spin, pe, Lorbit, Lspin)
!! author: David A. Minton
!!
!! Compute total system angular momentum vector and kinetic, potential and total system energy
Expand All @@ -14,31 +14,36 @@ subroutine symba_energy_eucl(npl, symba_plA, j2rp2, j4rp4, ke_orbit, ke_spin, pe
! arguments
integer(I4B), intent(in) :: npl
real(DP), intent(in) :: j2rp2, j4rp4
real(DP), intent(out) :: ke_orbit, ke_spin, pe, te
real(DP), dimension(:), intent(out) :: Ltot
real(DP), intent(out) :: ke_orbit, ke_spin, pe
real(DP), dimension(:), intent(out) :: Lorbit, Lspin
type(symba_pl), intent(inout) :: symba_plA

! internals
integer(I4B) :: i, j
integer(I8B) :: k
real(DP) :: rmag, v2, rot2, oblpot, hx, hy, hz, hsx, hsy, hsz
real(DP), dimension(npl) :: irh, kepl, kespinpl, pecb
real(DP), dimension(npl) :: Lplx, Lply, Lplz
real(DP), dimension(npl) :: Lplorbitx, Lplorbity, Lplorbitz
real(DP), dimension(npl) :: Lplspinx, Lplspiny, Lplspinz
real(DP), dimension(symba_plA%helio%swiftest%num_plpl_comparisons) :: pepl
logical, dimension(symba_plA%helio%swiftest%num_plpl_comparisons) :: lstatpl
logical, dimension(npl) :: lstatus

! executable code

Ltot = 0.0_DP
Lorbit(:) = 0.0_DP
Lspin(:) = 0.0_DP
ke_orbit = 0.0_DP
ke_spin = 0.0_DP
associate(xb => symba_plA%helio%swiftest%xb, vb => symba_plA%helio%swiftest%vb, mass => symba_plA%helio%swiftest%mass, radius => symba_plA%helio%swiftest%radius, &
Ip => symba_plA%helio%swiftest%Ip, rot => symba_plA%helio%swiftest%rot, xh => symba_plA%helio%swiftest%xh, status => symba_plA%helio%swiftest%status)
kepl(:) = 0.0_DP
Lplx(:) = 0.0_DP
Lply(:) = 0.0_DP
Lplz(:) = 0.0_DP
Lplorbitx(:) = 0.0_DP
Lplorbity(:) = 0.0_DP
Lplorbitz(:) = 0.0_DP
Lplspinx(:) = 0.0_DP
Lplspiny(:) = 0.0_DP
Lplspinz(:) = 0.0_DP
lstatus(1:npl) = status(1:npl) /= INACTIVE
!!$omp simd private(v2, rot2, hx, hy, hz)
do i = 1, npl
Expand All @@ -54,26 +59,24 @@ subroutine symba_energy_eucl(npl, symba_plA, j2rp2, j4rp4, ke_orbit, ke_spin, pe
hsz = Ip(3,i) * radius(i)**2 * rot(3,i)

! Angular momentum from orbit and spin
Lplx(i) = mass(i) * (hsx + hx)
Lply(i) = mass(i) * (hsy + hy)
Lplz(i) = mass(i) * (hsz + hz)
Lplorbitx(i) = mass(i) * hx
Lplorbity(i) = mass(i) * hy
Lplorbitz(i) = mass(i) * hz

Lplspinx(i) = mass(i) * hsx
Lplspiny(i) = mass(i) * hsy
Lplspinz(i) = mass(i) * hsz

! Kinetic energy from orbit and spin
kepl(i) = mass(i) * v2
kespinpl(i) = mass(i) * Ip(3, i) * radius(i)**2 * rot2
end do

ke_orbit = 0.5_DP * sum(kepl(1:npl), lstatus(:))
ke_spin = 0.5_DP * sum(kespinpl(1:npl), lstatus(:))
Ltot(1) = sum(Lplx(1:npl), lstatus(1:npl))
Ltot(2) = sum(Lply(1:npl), lstatus(1:npl))
Ltot(3) = sum(Lplz(1:npl), lstatus(1:npl))

! Do the central body potential energy component first
!$omp simd
do i = 2, npl
associate(px => xh(1,i), py => xh(2,i), pz => xh(3,i))
pecb(i) = - mass(1) * mass(i) / sqrt(px**2 + py**2 + pz**2)
pecb(i) = -mass(1) * mass(i) / sqrt(px**2 + py**2 + pz**2)
end associate
end do

Expand All @@ -85,6 +88,9 @@ subroutine symba_energy_eucl(npl, symba_plA, j2rp2, j4rp4, ke_orbit, ke_spin, pe
end associate
end do

ke_orbit = 0.5_DP * sum(kepl(1:npl), lstatus(:))
ke_spin = 0.5_DP * sum(kespinpl(1:npl), lstatus(:))

pe = sum(pepl(:), lstatpl(:)) + sum(pecb(2:npl), lstatus(2:npl))

! Potential energy from the oblateness term
Expand All @@ -97,7 +103,14 @@ subroutine symba_energy_eucl(npl, symba_plA, j2rp2, j4rp4, ke_orbit, ke_spin, pe
pe = pe + oblpot
end if

te = ke_orbit + ke_spin + pe
Lorbit(1) = sum(Lplorbitx(1:npl), lstatus(1:npl))
Lorbit(2) = sum(Lplorbity(1:npl), lstatus(1:npl))
Lorbit(3) = sum(Lplorbitz(1:npl), lstatus(1:npl))

Lspin(1) = sum(Lplspinx(1:npl), lstatus(1:npl))
Lspin(2) = sum(Lplspiny(1:npl), lstatus(1:npl))
Lspin(3) = sum(Lplspinz(1:npl), lstatus(1:npl))

end associate
return

Expand Down
Loading

0 comments on commit 4ba2e65

Please sign in to comment.