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

Commit

Permalink
Switched to tracking initial G*M instead of M and G*Mescape instead o…
Browse files Browse the repository at this point in the history 
…f Mescape
  • Loading branch information
@daminton
daminton committed Aug 19, 2021
1 parent fa9932f commit d86888f
Show file tree
Hide file tree
Showing 4 changed files with 25 additions and 24 deletions.
29 changes: 15 additions & 14 deletions src/io/io.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -17,7 +17,7 @@ module subroutine io_conservation_report(self, param, lterminal)
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) :: GMtot_now, Merror
real(DP) :: Lmag_now, Lerror
character(len=STRMAX) :: errmsg
character(len=*), parameter :: EGYFMT = '(ES23.16,10(",",ES23.16,:))' ! Format code for all simulation output
Expand All @@ -28,8 +28,8 @@ module subroutine io_conservation_report(self, param, lterminal)
"; D(Eorbit+Ecollisions)/|E0| = ", ES12.5, &
"; DM/M0 = ", ES12.5)'

associate(system => self, pl => self%pl, cb => self%cb, npl => self%pl%nbody, Ecollisions => self%Ecollisions, Lescape => self%Lescape, Mescape => self%Mescape, &
Euntracked => self%Euntracked, Eorbit_orig => param%Eorbit_orig, Mtot_orig => param%Mtot_orig, &
associate(system => self, pl => self%pl, cb => self%cb, npl => self%pl%nbody, Ecollisions => self%Ecollisions, Lescape => self%Lescape, GMescape => self%GMescape, &
Euntracked => self%Euntracked, Eorbit_orig => param%Eorbit_orig, GMtot_orig => param%GMtot_orig, &
Ltot_orig => param%Ltot_orig(:), Lmag_orig => param%Lmag_orig, Lorbit_orig => param%Lorbit_orig(:), Lspin_orig => param%Lspin_orig(:), &
lfirst => param%lfirstenergy)
if (param%energy_out /= "") then
Expand All @@ -50,10 +50,10 @@ module subroutine io_conservation_report(self, param, lterminal)
Lspin_now = system%Lspin
Eorbit_now = ke_orbit_now + ke_spin_now + pe_now
Ltot_now(:) = Lorbit_now(:) + Lspin_now(:) + Lescape(:)
Mtot_now = cb%mass + sum(pl%mass(1:npl)) + system%Mescape
GMtot_now = cb%Gmass + sum(pl%Gmass(1:npl)) + system%GMescape
if (lfirst) then
Eorbit_orig = Eorbit_now
Mtot_orig = Mtot_now
GMtot_orig = GMtot_now
Lorbit_orig(:) = Lorbit_now(:)
Lspin_orig(:) = Lspin_now(:)
Ltot_orig(:) = Ltot_now(:)
Expand All @@ -62,7 +62,7 @@ module subroutine io_conservation_report(self, param, lterminal)
end if

if (param%energy_out /= "") then
write(EGYIU,EGYFMT, err = 667, iomsg = errmsg) param%t, Eorbit_now, Ecollisions, Ltot_now, Mtot_now
write(EGYIU,EGYFMT, err = 667, iomsg = errmsg) param%t, Eorbit_now, Ecollisions, Ltot_now, GMtot_now
close(EGYIU, err = 667, iomsg = errmsg)
end if
if (.not.lfirst .and. lterminal) then
Expand All @@ -71,7 +71,7 @@ module subroutine io_conservation_report(self, param, lterminal)
Eorbit_error = (Eorbit_now - Eorbit_orig) / abs(Eorbit_orig)
Ecoll_error = Ecollisions / abs(Eorbit_orig)
Etotal_error = (Eorbit_now - Ecollisions - Eorbit_orig - Euntracked) / abs(Eorbit_orig)
Merror = (Mtot_now - Mtot_orig) / Mtot_orig
Merror = (GMtot_now - GMtot_orig) / GMtot_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!'
Expand Down Expand Up @@ -527,8 +527,8 @@ module subroutine io_param_reader(self, unit, iotype, v_list, iostat, iomsg)
if (param_value == "NO" .or. param_value == 'F') param%lfirstenergy = .false.
case("EORBIT_ORIG")
read(param_value, *, err = 667, iomsg = iomsg) param%Eorbit_orig
case("MTOT_ORIG")
read(param_value, *, err = 667, iomsg = iomsg) param%Mtot_orig
case("GMTOT_ORIG")
read(param_value, *, err = 667, iomsg = iomsg) param%GMtot_orig
case("LTOT_ORIG")
read(param_value, *, err = 667, iomsg = iomsg) param%Ltot_orig(1)
do i = 2, NDIM
Expand Down Expand Up @@ -558,8 +558,8 @@ module subroutine io_param_reader(self, unit, iotype, v_list, iostat, iomsg)
param_value = io_get_token(line, ifirst, ilast, iostat)
read(param_value, *, err = 667, iomsg = iomsg) param%Lescape(i)
end do
case("MESCAPE")
read(param_value, *, err = 667, iomsg = iomsg) param%Mescape
case("GMESCAPE")
read(param_value, *, err = 667, iomsg = iomsg) param%GMescape
case("ECOLLISIONS")
read(param_value, *, err = 667, iomsg = iomsg) param%Ecollisions
case("EUNTRACKED")
Expand Down Expand Up @@ -812,13 +812,13 @@ module subroutine io_param_writer(self, unit, iotype, v_list, iostat, iomsg)
if (param%lenergy) then
write(param_name, Afmt) "FIRSTENERGY"; write(param_value, Lfmt) param%lfirstenergy; write(unit, Afmt, err = 667, iomsg = iomsg) adjustl(param_name), adjustl(param_value)
write(param_name, Afmt) "EORBIT_ORIG"; write(param_value, Rfmt) param%Eorbit_orig; write(unit, Afmt, err = 667, iomsg = iomsg) adjustl(param_name), adjustl(param_value)
write(param_name, Afmt) "MTOT_ORIG"; write(param_value, Rfmt) param%Mtot_orig; write(unit, Afmt, err = 667, iomsg = iomsg) adjustl(param_name), adjustl(param_value)
write(param_name, Afmt) "GMTOT_ORIG"; write(param_value, Rfmt) param%GMtot_orig; write(unit, Afmt, err = 667, iomsg = iomsg) adjustl(param_name), adjustl(param_value)
write(unit, '("LTOT_ORIG ",3(1X,ES25.17))') param%Ltot_orig(:)
write(unit, '("LORBIT_ORIG",3(1X,ES25.17))') param%Lorbit_orig(:)
write(unit, '("LSPIN_ORIG ",3(1X,ES25.17))') param%Lspin_orig(:)
write(unit, '("LESCAPE ",3(1X,ES25.17))') param%Lescape(:)

write(param_name, Afmt) "MESCAPE"; write(param_value, Rfmt) param%Mescape; write(unit, Afmt, err = 667, iomsg = iomsg) adjustl(param_name), adjustl(param_value)
write(param_name, Afmt) "GMescape"; write(param_value, Rfmt) param%GMescape; write(unit, Afmt, err = 667, iomsg = iomsg) adjustl(param_name), adjustl(param_value)
write(param_name, Afmt) "ECOLLISIONS"; write(param_value, Rfmt) param%Ecollisions; write(unit, Afmt, err = 667, iomsg = iomsg) adjustl(param_name), adjustl(param_value)
write(param_name, Afmt) "EUNTRACKED"; write(param_value, Rfmt) param%Euntracked; write(unit, Afmt, err = 667, iomsg = iomsg) adjustl(param_name), adjustl(param_value)
end if
Expand Down Expand Up @@ -937,7 +937,8 @@ module subroutine io_read_in_cb(self, param)

select type(cb => self)
class is (symba_cb)
cb%M0 = cb%mass
cb%GM0 = cb%Gmass
cb%dGM = 0.0_DP
cb%R0 = cb%radius
cb%L0(:) = cb%Ip(3) * cb%mass * cb%radius**2 * cb%rot(:)
end select
Expand Down
6 changes: 3 additions & 3 deletions src/modules/swiftest_classes.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -60,14 +60,14 @@ module swiftest_classes

! Initial values to pass to the energy report subroutine (usually only used in the case of a restart, otherwise these will be updated with initial conditions values)
real(DP) :: Eorbit_orig = 0.0_DP !! Initial orbital energy
real(DP) :: Mtot_orig = 0.0_DP !! Initial system mass
real(DP) :: GMtot_orig = 0.0_DP !! Initial system mass
real(DP) :: Lmag_orig = 0.0_DP !! Initial total angular momentum magnitude
real(DP), dimension(NDIM) :: Ltot_orig = 0.0_DP !! Initial total angular momentum vector
real(DP), dimension(NDIM) :: Lorbit_orig = 0.0_DP !! Initial orbital angular momentum
real(DP), dimension(NDIM) :: Lspin_orig = 0.0_DP !! Initial spin angular momentum vector
real(DP), dimension(NDIM) :: Ltot = 0.0_DP !! System angular momentum vector
real(DP), dimension(NDIM) :: Lescape = 0.0_DP !! Angular momentum of bodies that escaped the system (used for bookeeping)
real(DP) :: Mescape = 0.0_DP !! Mass of bodies that escaped the system (used for bookeeping)
real(DP) :: GMescape = 0.0_DP !! Mass of bodies that escaped the system (used for bookeeping)
real(DP) :: Ecollisions = 0.0_DP !! Energy lost from system due to collisions
real(DP) :: Euntracked = 0.0_DP !! Energy gained from system due to escaped bodies
logical :: lfirstenergy = .true. !! This is the first time computing energe
Expand Down Expand Up @@ -291,7 +291,7 @@ module swiftest_classes
real(DP), dimension(NDIM) :: Lorbit = 0.0_DP !! System orbital angular momentum vector
real(DP), dimension(NDIM) :: Lspin = 0.0_DP !! System spin angular momentum vector
real(DP), dimension(NDIM) :: Lescape = 0.0_DP !! Angular momentum of bodies that escaped the system (used for bookeeping)
real(DP) :: Mescape = 0.0_DP !! Mass of bodies that escaped the system (used for bookeeping)
real(DP) :: GMescape = 0.0_DP !! Mass of bodies that escaped the system (used for bookeeping)
real(DP) :: Ecollisions = 0.0_DP !! Energy lost from system due to collisions
real(DP) :: Euntracked = 0.0_DP !! Energy gained from system due to escaped bodies
logical :: lbeg !! True if this is the beginning of a step. This is used so that test particle steps can be calculated
Expand Down
4 changes: 2 additions & 2 deletions src/modules/symba_classes.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -55,8 +55,8 @@ module symba_classes
!*******************************************************************************************************************************
!> SyMBA central body particle class
type, extends(helio_cb) :: symba_cb
real(DP) :: M0 = 0.0_DP !! Initial mass of the central body
real(DP) :: dM = 0.0_DP !! Change in mass of the central body
real(DP) :: GM0 = 0.0_DP !! Initial G*mass of the central body
real(DP) :: dGM = 0.0_DP !! Change in G*mass of the central body
real(DP) :: R0 = 0.0_DP !! Initial radius of the central body
real(DP) :: dR = 0.0_DP !! Change in the radius of the central body
type(symba_particle_info) :: info
Expand Down
10 changes: 5 additions & 5 deletions src/symba/symba_discard.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -85,7 +85,7 @@ subroutine symba_discard_conserve_mtm(pl, system, param, ipl, lescape_body)
class is (symba_cb)

! Add the potential and kinetic energy of the lost body to the records
pe = -cb%mass * pl%mass(ipl) / norm2(pl%xb(:, ipl) - cb%xb(:))
pe = -cb%Gmass * pl%mass(ipl) / norm2(pl%xb(:, ipl) - cb%xb(:))
ke_orbit = 0.5_DP * pl%mass(ipl) * dot_product(pl%vb(:, ipl), pl%vb(:, ipl))
if (param%lrotation) then
ke_spin = 0.5_DP * pl%mass(ipl) * pl%radius(ipl)**2 * pl%Ip(3, ipl) * dot_product(pl%rot(:, ipl), pl%rot(:, ipl))
Expand All @@ -96,7 +96,7 @@ subroutine symba_discard_conserve_mtm(pl, system, param, ipl, lescape_body)
! Add the pre-collision ke of the central body to the records
! Add planet mass to central body accumulator
if (lescape_body) then
system%Mescape = system%Mescape + pl%mass(ipl)
system%GMescape = system%GMescape + pl%Gmass(ipl)
do i = 1, pl%nbody
if (i == ipl) cycle
pe = pe - pl%Gmass(i) * pl%mass(ipl) / norm2(pl%xb(:, ipl) - pl%xb(:, i))
Expand Down Expand Up @@ -133,10 +133,10 @@ subroutine symba_discard_conserve_mtm(pl, system, param, ipl, lescape_body)
ke_orbit = ke_orbit + 0.5_DP * cb%mass * dot_product(cb%vb(:), cb%vb(:))
if (param%lrotation) ke_spin = ke_spin + 0.5_DP * cb%mass * cb%radius**2 * cb%Ip(3) * dot_product(cb%rot(:), cb%rot(:))
! Update mass of central body to be consistent with its total mass
cb%dM = cb%dM + pl%mass(ipl)
cb%dGM = cb%dGM + pl%Gmass(ipl)
cb%dR = cb%dR + 1.0_DP / 3.0_DP * (pl%radius(ipl) / cb%radius)**3 - 2.0_DP / 9.0_DP * (pl%radius(ipl) / cb%radius)**6
cb%mass = cb%M0 + cb%dM
cb%Gmass = param%GU * cb%mass
cb%Gmass = cb%GM0 + cb%dGM
cb%mass = cb%Gmass / param%GU
cb%radius = cb%R0 + cb%dR
param%rmin = cb%radius
! Add planet angular momentum to central body accumulator
Expand Down

0 comments on commit d86888f

Please sign in to comment.