diff --git a/examples/symba_chambers_2013/param.in b/examples/symba_chambers_2013/param.in index 28500d32b..2e263e1b5 100644 --- a/examples/symba_chambers_2013/param.in +++ b/examples/symba_chambers_2013/param.in @@ -8,8 +8,8 @@ CB_IN sun_MsunAUYR.in PL_IN pl_chambers_2013.in TP_IN tp.in IN_TYPE ASCII -ISTEP_OUT 6250 ! output cadence -ISTEP_DUMP 6250 ! system dump cadence +ISTEP_OUT 625 ! output cadence +ISTEP_DUMP 625 ! system dump cadence BIN_OUT bin.dat PARTICLE_OUT particle.dat OUT_TYPE REAL8 ! double precision real output @@ -18,7 +18,6 @@ OUT_STAT REPLACE CHK_CLOSE yes ! check for planetary close encounters CHK_RMAX 100000.0 ! discard outside of EXTRA_FORCE no ! no extra user-defined forces -DISCARD_OUT discard.out BIG_DISCARD no ! output all planets if anything discarded RHILL_PRESENT yes ! Hill's sphere radii in input file MU2KG 1.98847e30 ! (M_sun-> kg) @@ -29,4 +28,5 @@ ENERGY yes ENERGY_OUT energy.dat ROTATION yes FRAGMENTATION yes +DISCARD_OUT discard.out SEED 8 12261555 871132 92734722 21132443 36344777 4334443 219291656 3848566 diff --git a/examples/symba_clement_2018/param.in b/examples/symba_clement_2018/param.in index 79890595e..f4c163875 100644 --- a/examples/symba_clement_2018/param.in +++ b/examples/symba_clement_2018/param.in @@ -27,6 +27,7 @@ TU2S 3.15569259747e7 ! time unit to seconds (years --> seconds) GMTINY 1e-10 ENERGY yes ENERGY_OUT energy.dat +ENERGY_OUT energy.dat ROTATION yes FRAGMENTATION yes SEED 8 12261555 871132 92734722 21132443 36344777 4334443 219291656 3848566 diff --git a/examples/symba_energy_momentum/param.escape.in b/examples/symba_energy_momentum/param.escape.in index 90d118017..5f8ef2d0d 100644 --- a/examples/symba_energy_momentum/param.escape.in +++ b/examples/symba_energy_momentum/param.escape.in @@ -9,20 +9,16 @@ ISTEP_OUT 1 ! output cadence every year BIN_OUT bin.escape.dat PARTICLE_OUT particle.escape.dat OUT_TYPE REAL8 ! double precision real output -OUT_FORM XV ! osculating element output +OUT_FORM EL ! osculating element output OUT_STAT REPLACE ISTEP_DUMP 1 ! system dump cadence -J2 0.0 ! no J2 term -J4 0.0 ! no J4 term CHK_CLOSE yes ! check for planetary close encounters -CHK_RMIN 0.00465047 ! check for close solar encounters in AU CHK_RMAX 10000.0 ! discard outside of -CHK_EJECT -1.0 ! ignore this check -CHK_QMIN -1.0 ! ignore this check ENC_OUT enc.escape.dat +DISCARD_OUT discard.escape.dat EXTRA_FORCE no ! no extra user-defined forces BIG_DISCARD no ! output all planets if anything discarded -RHILL_PRESENT no ! Hill's sphere radii in input file +RHILL_PRESENT yes ! Hill's sphere radii in input file GMTINY 1.0e-16 FRAGMENTATION yes MU2KG 1.98908e30 diff --git a/examples/symba_energy_momentum/param.sun.in b/examples/symba_energy_momentum/param.sun.in index a7748b19c..351d8adb2 100644 --- a/examples/symba_energy_momentum/param.sun.in +++ b/examples/symba_energy_momentum/param.sun.in @@ -19,15 +19,17 @@ CHK_RMIN 0.005 CHK_RMAX 1e2 CHK_EJECT -1.0 ! ignore this check CHK_QMIN -1.0 ! ignore this check -ENC_OUT enc.escape.dat +ENC_OUT enc.sun.dat +DISCARD_OUT discard.sun.dat EXTRA_FORCE no ! no extra user-defined forces BIG_DISCARD no ! output all planets if anything discarded RHILL_PRESENT no ! Hill's sphere radii in input file GMTINY 1.0e-16 -FRAGMENTATION yes +FRAGMENTATION no MU2KG 1.98908e30 TU2S 3.1556925e7 DU2M 1.49598e11 ENERGY yes +ENERGY_OUT energy.sun.out ROTATION yes SEED 8 1230834 2346113 123409874 -123121105 -767545 -534058022 343309814 -12535638 diff --git a/src/discard/discard.f90 b/src/discard/discard.f90 index 2998e8804..5c6f704e1 100644 --- a/src/discard/discard.f90 +++ b/src/discard/discard.f90 @@ -112,8 +112,9 @@ subroutine discard_cb_tp(tp, system, param) ! Internals integer(I4B) :: i real(DP) :: energy, vb2, rb2, rh2, rmin2, rmax2, rmaxu2 + character(len=STRMAX) :: idstr, timestr - associate(ntp => tp%nbody, cb => system%cb, t => param%t, Gmtot => system%Gmtot) + associate(ntp => tp%nbody, cb => system%cb, Gmtot => system%Gmtot) rmin2 = max(param%rmin * param%rmin, cb%radius * cb%radius) rmax2 = param%rmax**2 rmaxu2 = param%rmaxu**2 @@ -122,12 +123,16 @@ subroutine discard_cb_tp(tp, system, param) rh2 = dot_product(tp%xh(:, i), tp%xh(:, i)) if ((param%rmax >= 0.0_DP) .and. (rh2 > rmax2)) then tp%status(i) = DISCARDED_RMAX - write(*, *) "Particle ", tp%id(i), " too far from sun at t = ", t + write(idstr, *) tp%id(i) + write(timestr, *) param%t + write(*, *) "Particle " // trim(adjustl(idstr)) // " too far from the central body at t = " // trim(adjustl(timestr)) tp%ldiscard(i) = .true. tp%lmask(i) = .false. else if ((param%rmin >= 0.0_DP) .and. (rh2 < rmin2)) then tp%status(i) = DISCARDED_RMIN - write(*, *) "Particle ", tp%id(i), " too close to sun at t = ", t + write(idstr, *) tp%id(i) + write(timestr, *) param%t + write(*, *) "Particle " // trim(adjustl(idstr)) // " too close to the central body at t = " // trim(adjustl(timestr)) tp%ldiscard(i) = .true. tp%lmask(i) = .false. else if (param%rmaxu >= 0.0_DP) then @@ -136,7 +141,9 @@ subroutine discard_cb_tp(tp, system, param) energy = 0.5_DP * vb2 - Gmtot / sqrt(rb2) if ((energy > 0.0_DP) .and. (rb2 > rmaxu2)) then tp%status(i) = DISCARDED_RMAXU - write(*, *) "Particle ", tp%id(i), " is unbound and too far from barycenter at t = ", t + write(idstr, *) tp%id(i) + write(timestr, *) param%t + write(*, *) "Particle " // trim(adjustl(idstr)) // " is unbound and too far from barycenter at t = " // trim(adjustl(timestr)) tp%ldiscard(i) = .true. tp%lmask(i) = .false. end if diff --git a/src/io/io.f90 b/src/io/io.f90 index f7532e7b9..670d98db8 100644 --- a/src/io/io.f90 +++ b/src/io/io.f90 @@ -14,11 +14,10 @@ module subroutine io_conservation_report(self, param, lterminal) ! Internals real(DP), dimension(NDIM) :: Ltot_now, Lorbit_now, Lspin_now real(DP), dimension(NDIM), save :: Ltot_last, Lorbit_last, Lspin_last - 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 + real(DP) :: GMtot_now + real(DP) :: Lerror, Merror character(len=STRMAX) :: errmsg character(len=*), parameter :: EGYFMT = '(ES23.16,10(",",ES23.16,:))' ! Format code for all simulation output character(len=*), parameter :: EGYHEADER = '("t,Eorbit,Ecollisions,Lx,Ly,Lz,Mtot")' @@ -28,12 +27,9 @@ 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, & - Ltot_orig => param%Ltot_orig(:), Lmag_orig => param%Lmag_orig, Lorbit_orig => param%Lorbit_orig(:), Lspin_orig => param%Lspin_orig(:), & - lfirst => param%lfirstenergy) + associate(system => self, pl => self%pl, cb => self%cb, npl => self%pl%nbody) if (param%energy_out /= "") then - if (lfirst .and. (param%out_stat /= "OLD")) then + if (param%lfirstenergy .and. (param%out_stat /= "OLD")) then open(unit = EGYIU, file = param%energy_out, form = "formatted", status = "replace", action = "write", err = 667, iomsg = errmsg) write(EGYIU,EGYHEADER, err = 667, iomsg = errmsg) else @@ -46,48 +42,34 @@ module subroutine io_conservation_report(self, param, lterminal) ke_orbit_now = system%ke_orbit ke_spin_now = system%ke_spin pe_now = system%pe - Lorbit_now = system%Lorbit - Lspin_now = system%Lspin + Lorbit_now(:) = system%Lorbit(:) + 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 - if (lfirst) then - 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. + Ltot_now(:) = system%Ltot(:) + param%Lescape(:) + GMtot_now = system%GMtot + param%GMescape + + if (param%lfirstenergy) then + param%Eorbit_orig = Eorbit_now + param%GMtot_orig = GMtot_now + param%Lorbit_orig(:) = Lorbit_now(:) + param%Lspin_orig(:) = Lspin_now(:) + param%Ltot_orig(:) = Ltot_now(:) + param%lfirstenergy = .false. 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, param%Ecollisions, Ltot_now, GMtot_now close(EGYIU, err = 667, iomsg = errmsg) end if - if (.not.lfirst .and. lterminal) then - Lmag_now = norm2(Ltot_now) - Lerror = norm2(Ltot_now - Ltot_orig) / Lmag_orig - 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 + + if (.not.param%lfirstenergy .and. lterminal) then + Lerror = norm2(Ltot_now(:) - param%Ltot_orig(:)) / norm2(param%Ltot_orig(:)) + Eorbit_error = (Eorbit_now - param%Eorbit_orig) / abs(param%Eorbit_orig) + Ecoll_error = param%Ecollisions / abs(param%Eorbit_orig) + Etotal_error = (Eorbit_now - param%Ecollisions - param%Eorbit_orig - param%Euntracked) / abs(param%Eorbit_orig) + Merror = (GMtot_now - param%GMtot_orig) / param%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!' - ! 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_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 @@ -212,6 +194,7 @@ module subroutine io_dump_system(self, param) dump_param%out_stat = 'APPEND' dump_param%in_type = REAL8_TYPE dump_param%T0 = param%t + call dump_param%dump(param_file_name) call self%cb%dump(dump_param) @@ -527,8 +510,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 @@ -536,7 +519,6 @@ 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%Ltot_orig(i) end do - param%Lmag_orig = norm2(param%Ltot_orig(:)) case("LORBIT_ORIG") read(param_value, *, err = 667, iomsg = iomsg) param%Lorbit_orig(1) do i = 2, NDIM @@ -558,8 +540,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 +794,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 +919,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..f08e8eb05 100644 --- a/src/modules/swiftest_classes.f90 +++ b/src/modules/swiftest_classes.f90 @@ -60,14 +60,12 @@ 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) :: Lmag_orig = 0.0_DP !! Initial total angular momentum magnitude + real(DP) :: GMtot_orig = 0.0_DP !! Initial system mass 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 @@ -283,17 +281,14 @@ module swiftest_classes class(swiftest_tp), allocatable :: tp !! Test particle data structure class(swiftest_tp), allocatable :: tp_discards !! Discarded test particle data structure class(swiftest_pl), allocatable :: pl_discards !! Discarded massive body particle data structure - real(DP) :: Gmtot = 0.0_DP !! Total system mass - used for barycentric coordinate conversion + real(DP) :: GMtot = 0.0_DP !! Total system mass - used for barycentric coordinate conversion real(DP) :: ke_orbit = 0.0_DP !! System orbital kinetic energy real(DP) :: ke_spin = 0.0_DP !! System spin kinetic energy real(DP) :: pe = 0.0_DP !! System potential energy real(DP) :: te = 0.0_DP !! System total energy 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) :: 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 + real(DP), dimension(NDIM) :: Ltot = 0.0_DP !! System angular momentum vector logical :: lbeg !! True if this is the beginning of a step. This is used so that test particle steps can be calculated !! separately from massive bodies. Massive body variables are saved at half steps, and passed to !! the test particles diff --git a/src/modules/symba_classes.f90 b/src/modules/symba_classes.f90 index 056e6f390..f2b9134d8 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 @@ -215,14 +215,14 @@ module subroutine symba_collision_resolve_fragmentations(self, system, param) implicit none class(symba_plplenc), intent(inout) :: self !! SyMBA pl-pl encounter list class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object - class(symba_parameters), intent(in) :: param !! Current run configuration parameters with SyMBA additions + class(symba_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA additions end subroutine symba_collision_resolve_fragmentations module subroutine symba_collision_resolve_mergers(self, system, param) implicit none class(symba_plplenc), intent(inout) :: self !! SyMBA pl-pl encounter list class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object - class(symba_parameters), intent(in) :: param !! Current run configuration parameters with SyMBA additions + class(symba_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA additions end subroutine symba_collision_resolve_mergers module subroutine symba_collision_resolve_plplenc(self, system, param, t, dt, irec) @@ -309,7 +309,7 @@ end function symba_encounter_check_tp module function symba_collision_casedisruption(system, param, family, x, v, mass, radius, L_spin, Ip, mass_res, Qloss) result(status) implicit none class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object - class(symba_parameters), intent(in) :: param !! Current run configuration parameters with SyMBA additions + class(symba_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA additions integer(I4B), dimension(:), intent(in) :: family !! List of indices of all bodies inovlved in the collision real(DP), dimension(:,:), intent(inout) :: x, v, L_spin, Ip !! Input values that represent a 2-body equivalent of a possibly 2+ body collision real(DP), dimension(:), intent(inout) :: mass, radius !! Input values that represent a 2-body equivalent of a possibly 2+ body collision @@ -321,7 +321,7 @@ end function symba_collision_casedisruption module function symba_collision_casehitandrun(system, param, family, x, v, mass, radius, L_spin, Ip, mass_res, Qloss) result(status) implicit none class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object - class(symba_parameters), intent(in) :: param !! Current run configuration parameters with SyMBA additions + class(symba_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA additions integer(I4B), dimension(:), intent(in) :: family !! List of indices of all bodies inovlved in the collision real(DP), dimension(:,:), intent(inout) :: x, v, L_spin, Ip !! Input values that represent a 2-body equivalent of a possibly 2+ body collision real(DP), dimension(:), intent(inout) :: mass, radius !! Input values that represent a 2-body equivalent of a possibly 2+ body collision @@ -333,7 +333,7 @@ end function symba_collision_casehitandrun module function symba_collision_casemerge(system, param, family, x, v, mass, radius, L_spin, Ip) result(status) implicit none class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object - class(symba_parameters), intent(in) :: param !! Current run configuration parameters with SyMBA additions + class(symba_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA additions integer(I4B), dimension(:), intent(in) :: family !! List of indices of all bodies inovlved in the collision real(DP), dimension(:,:), intent(in) :: x, v, L_spin, Ip !! Input values that represent a 2-body equivalent of a possibly 2+ body collision real(DP), dimension(:), intent(in) :: mass, radius !! Input values that represent a 2-body equivalent of a possibly 2+ body collisio @@ -343,7 +343,7 @@ end function symba_collision_casemerge module function symba_collision_casesupercatastrophic(system, param, family, x, v, mass, radius, L_spin, Ip, mass_res, Qloss) result(status) implicit none class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object - class(symba_parameters), intent(in) :: param !! Current run configuration parameters with SyMBA additions + class(symba_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA additions integer(I4B), dimension(:), intent(in) :: family !! List of indices of all bodies inovlved in the collision real(DP), dimension(:,:), intent(inout) :: x, v, L_spin, Ip !! Input values that represent a 2-body equivalent of a possibly 2+ body collision real(DP), dimension(:), intent(inout) :: mass, radius !! Input values that represent a 2-body equivalent of a possibly 2+ body collision diff --git a/src/symba/symba_collision.f90 b/src/symba/symba_collision.f90 index b00a4a8e6..5ba08558c 100644 --- a/src/symba/symba_collision.f90 +++ b/src/symba/symba_collision.f90 @@ -13,7 +13,7 @@ module function symba_collision_casedisruption(system, param, family, x, v, mass implicit none ! Arguments class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object - class(symba_parameters), intent(in) :: param !! Current run configuration parameters with SyMBA additions + class(symba_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA additions integer(I4B), dimension(:), intent(in) :: family !! List of indices of all bodies inovlved in the collision real(DP), dimension(:,:), intent(inout) :: x, v, L_spin, Ip !! Input values that represent a 2-body equivalent of a possibly 2+ body collision real(DP), dimension(:), intent(inout) :: mass, radius !! Input values that represent a 2-body equivalent of a possibly 2+ body collision @@ -106,7 +106,7 @@ module function symba_collision_casehitandrun(system, param, family, x, v, mass, implicit none ! Arguments class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object - class(symba_parameters), intent(in) :: param !! Current run configuration parameters with SyMBA additions + class(symba_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA additions integer(I4B), dimension(:), intent(in) :: family !! List of indices of all bodies inovlved in the collision real(DP), dimension(:,:), intent(inout) :: x, v, L_spin, Ip !! Input values that represent a 2-body equivalent of a possibly 2+ body collision real(DP), dimension(:), intent(inout) :: mass, radius !! Input values that represent a 2-body equivalent of a possibly 2+ body collision @@ -213,7 +213,7 @@ module function symba_collision_casemerge(system, param, family, x, v, mass, rad implicit none ! Arguments class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object - class(symba_parameters), intent(in) :: param !! Current run configuration parameters with SyMBA additions + class(symba_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA additions integer(I4B), dimension(:), intent(in) :: family !! List of indices of all bodies inovlved in the collision real(DP), dimension(:,:), intent(in) :: x, v, L_spin, Ip !! Input values that represent a 2-body equivalent of a possibly 2+ body collision real(DP), dimension(:), intent(in) :: mass, radius !! Input values that represent a 2-body equivalent of a possibly 2+ body collision @@ -269,7 +269,7 @@ module function symba_collision_casemerge(system, param, family, x, v, mass, rad ! Assume prinicpal axis rotation on 3rd Ip axis rot_frag(:,1) = L_spin_new(:) / (Ip_frag(3,1) * m_frag(1) * rad_frag(1)**2) else ! If spin is not enabled, we will consider the lost pre-collision angular momentum as "escaped" and add it to our bookkeeping variable - system%Lescape(:) = system%Lescape(:) + L_orb_old(:) + param%Lescape(:) = param%Lescape(:) + L_orb_old(:) end if ! Keep track of the component of potential energy due to the pre-impact family for book-keeping @@ -280,8 +280,8 @@ module function symba_collision_casemerge(system, param, family, x, v, mass, rad pe = pe - pl%Gmass(i) * pl%mass(j) / norm2(pl%xb(:, i) - pl%xb(:, j)) end do end do - system%Ecollisions = system%Ecollisions + pe - system%Euntracked = system%Euntracked - pe + param%Ecollisions = param%Ecollisions + pe + param%Euntracked = param%Euntracked - pe status = MERGED call symba_collision_mergeaddsub(system, param, family, id_frag, Ip_frag, m_frag, rad_frag, xb_frag, vb_frag, rot_frag, status) @@ -300,7 +300,7 @@ module function symba_collision_casesupercatastrophic(system, param, family, x, implicit none ! Arguments class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object - class(symba_parameters), intent(in) :: param !! Current run configuration parameters with SyMBA additions + class(symba_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA additions integer(I4B), dimension(:), intent(in) :: family !! List of indices of all bodies inovlved in the collision real(DP), dimension(:,:), intent(inout) :: x, v, L_spin, Ip !! Input values that represent a 2-body equivalent of a possibly 2+ body collision real(DP), dimension(:), intent(inout) :: mass, radius !! Input values that represent a 2-body equivalent of a possibly 2+ body collision @@ -930,7 +930,7 @@ module subroutine symba_collision_resolve_fragmentations(self, system, param) ! Arguments class(symba_plplenc), intent(inout) :: self !! SyMBA pl-pl encounter list class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object - class(symba_parameters), intent(in) :: param !! Current run configuration parameters with SyMBA additions + class(symba_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA additions ! Internals ! Internals integer(I4B), dimension(:), allocatable :: family !! List of indices of all bodies inovlved in the collision @@ -1023,7 +1023,7 @@ module subroutine symba_collision_resolve_mergers(self, system, param) ! Arguments class(symba_plplenc), intent(inout) :: self !! SyMBA pl-pl encounter list class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object - class(symba_parameters), intent(in) :: param !! Current run configuration parameters with SyMBA additions + class(symba_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA additions ! Internals integer(I4B), dimension(:), allocatable :: family !! List of indices of all bodies inovlved in the collision integer(I4B), dimension(2) :: idx_parent !! Index of the two bodies considered the "parents" of the collision @@ -1119,7 +1119,7 @@ module subroutine symba_collision_resolve_plplenc(self, system, param, t, dt, ir if (param%lenergy) then call system%get_energy_and_momentum(param) Eorbit_after = system%te - system%Ecollisions = system%Ecollisions + (Eorbit_after - Eorbit_before) + param%Ecollisions = param%Ecollisions + (Eorbit_after - Eorbit_before) end if end select diff --git a/src/symba/symba_discard.f90 b/src/symba/symba_discard.f90 index 2bf7b2a90..193461bc7 100644 --- a/src/symba/symba_discard.f90 +++ b/src/symba/symba_discard.f90 @@ -20,6 +20,7 @@ subroutine symba_discard_cb_pl(pl, system, param) ! Internals integer(I4B) :: i, j real(DP) :: energy, vb2, rb2, rh2, rmin2, rmax2, rmaxu2 + character(len=STRMAX) :: idstr, timestr associate(npl => pl%nbody, cb => system%cb) call system%set_msys() @@ -33,12 +34,16 @@ subroutine symba_discard_cb_pl(pl, system, param) pl%ldiscard(i) = .true. pl%lcollision(i) = .false. pl%status(i) = DISCARDED_RMAX - write(*, *) "Massive body ", pl%id(i), " too far from the central body at t = ", param%t + write(idstr, *) pl%id(i) + write(timestr, *) param%t + write(*, *) "Massive body " // trim(adjustl(idstr)) // " too far from the central body at t = " // trim(adjustl(timestr)) else if ((param%rmin >= 0.0_DP) .and. (rh2 < rmin2)) then pl%ldiscard(i) = .true. pl%lcollision(i) = .false. pl%status(i) = DISCARDED_RMIN - write(*, *) "Massive body ", pl%id(i), " too close to the central body at t = ", param%t + write(idstr, *) pl%id(i) + write(timestr, *) param%t + write(*, *) "Massive body " // trim(adjustl(idstr)) // " too close to the central body at t = " // trim(adjustl(timestr)) else if (param%rmaxu >= 0.0_DP) then rb2 = dot_product(pl%xb(:,i), pl%xb(:,i)) vb2 = dot_product(pl%vb(:,i), pl%vb(:,i)) @@ -47,7 +52,9 @@ subroutine symba_discard_cb_pl(pl, system, param) pl%ldiscard(i) = .true. pl%lcollision(i) = .false. pl%status(i) = DISCARDED_RMAXU - write(*, *) "Massive body ", pl%id(i), " is unbound and too far from barycenter at t = ", param%t + write(idstr, *) pl%id(i) + write(timestr, *) param%t + write(*, *) "Massive body " // trim(adjustl(idstr)) // " is unbound and too far from barycenter at t = " // trim(adjustl(timestr)) end if end if end if @@ -78,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)) @@ -89,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) + param%GMescape = param%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)) @@ -112,8 +119,8 @@ subroutine symba_discard_conserve_mtm(pl, system, param, ipl, lescape_body) Ltot(:) = Ltot(:) - Lpl(:) end do Ltot(:) = Ltot(:) - cb%mass * (cb%xb(:) .cross. cb%vb(:)) - system%Lescape(:) = system%Lescape(:) + Ltot(:) - if (param%lrotation) system%Lescape(:) = system%Lescape + pl%mass(ipl) * pl%radius(ipl)**2 * pl%Ip(3, ipl) * pl%rot(:, ipl) + param%Lescape(:) = param%Lescape(:) + Ltot(:) + if (param%lrotation) param%Lescape(:) = param%Lescape + pl%mass(ipl) * pl%radius(ipl)**2 * pl%Ip(3, ipl) * pl%rot(:, ipl) else xcom(:) = (pl%mass(ipl) * pl%xb(:, ipl) + cb%mass * cb%xb(:)) / (cb%mass + pl%mass(ipl)) @@ -126,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 @@ -147,11 +154,11 @@ subroutine symba_discard_conserve_mtm(pl, system, param, ipl, lescape_body) ! We must do this for proper book-keeping, since we can no longer track this body's contribution to energy directly if (lescape_body) then - system%Ecollisions = system%Ecollisions + ke_orbit + ke_spin + pe - system%Euntracked = system%Euntracked - (ke_orbit + ke_spin + pe) + param%Ecollisions = param%Ecollisions + ke_orbit + ke_spin + pe + param%Euntracked = param%Euntracked - (ke_orbit + ke_spin + pe) else - system%Ecollisions = system%Ecollisions + pe - system%Euntracked = system%Euntracked - pe + param%Ecollisions = param%Ecollisions + pe + param%Euntracked = param%Euntracked - pe end if end select @@ -306,7 +313,7 @@ module subroutine symba_discard_pl(self, system, param) if (param%lenergy) then call system%get_energy_and_momentum(param) Eorbit_after = system%te - system%Ecollisions = system%Ecollisions + (Eorbit_after - Eorbit_before) + param%Ecollisions = param%Ecollisions + (Eorbit_after - Eorbit_before) end if end associate diff --git a/src/symba/symba_io.f90 b/src/symba/symba_io.f90 index 06472aae9..1c0981232 100644 --- a/src/symba/symba_io.f90 +++ b/src/symba/symba_io.f90 @@ -278,7 +278,7 @@ module subroutine symba_io_read_particle(system, param) lmatch = .false. read(LUN, err = 667, iomsg = errmsg, end = 333) id - if (idx == cb%id) then + if (id == cb%id) then read(LUN, err = 667, iomsg = errmsg) cb%info lmatch = .true. else diff --git a/src/util/util_get_energy_momentum.f90 b/src/util/util_get_energy_momentum.f90 index dffdc29ae..4243070ac 100644 --- a/src/util/util_get_energy_momentum.f90 +++ b/src/util/util_get_energy_momentum.f90 @@ -28,6 +28,7 @@ module subroutine util_get_energy_momentum_system(self, param) associate(system => self, pl => self%pl, npl => self%pl%nbody, cb => self%cb) system%Lorbit(:) = 0.0_DP system%Lspin(:) = 0.0_DP + system%Ltot(:) = 0.0_DP system%ke_orbit = 0.0_DP system%ke_spin = 0.0_DP @@ -38,8 +39,10 @@ module subroutine util_get_energy_momentum_system(self, param) Lplspinx(:) = 0.0_DP Lplspiny(:) = 0.0_DP Lplspinz(:) = 0.0_DP + lstatus(1:npl) = pl%status(1:npl) /= INACTIVE + system%GMtot = cb%Gmass + sum(pl%Gmass(1:npl), lstatus(1:npl)) kecb = cb%mass * dot_product(cb%vb(:), cb%vb(:)) Lcborbit(:) = cb%mass * (cb%xb(:) .cross. cb%vb(:)) @@ -121,6 +124,7 @@ module subroutine util_get_energy_momentum_system(self, param) end if system%te = system%ke_orbit + system%ke_spin + system%pe + system%Ltot(:) = system%Lorbit(:) + system%Lspin(:) end associate return