diff --git a/src/io/io.f90 b/src/io/io.f90 index f7532e7b9..37a69a6ec 100644 --- a/src/io/io.f90 +++ b/src/io/io.f90 @@ -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 @@ -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 @@ -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(:) @@ -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 @@ -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!' @@ -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 @@ -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") @@ -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 @@ -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 diff --git a/src/modules/swiftest_classes.f90 b/src/modules/swiftest_classes.f90 index 9a674c830..465170151 100644 --- a/src/modules/swiftest_classes.f90 +++ b/src/modules/swiftest_classes.f90 @@ -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 @@ -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 diff --git a/src/modules/symba_classes.f90 b/src/modules/symba_classes.f90 index 056e6f390..9384a11a6 100644 --- a/src/modules/symba_classes.f90 +++ b/src/modules/symba_classes.f90 @@ -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 diff --git a/src/symba/symba_discard.f90 b/src/symba/symba_discard.f90 index 7f3ed3b21..32879a9a5 100644 --- a/src/symba/symba_discard.f90 +++ b/src/symba/symba_discard.f90 @@ -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)) @@ -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)) @@ -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