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

Commit

Permalink
Merge branch 'debug' into OOPTides
Browse files Browse the repository at this point in the history
  • Loading branch information
@daminton
daminton committed Aug 19, 2021
2 parents 8e590cf + 87239f0 commit 7b1f1cb
Show file tree
Hide file tree
Showing 12 changed files with 104 additions and 109 deletions.
6 changes: 3 additions & 3 deletions examples/symba_chambers_2013/param.in
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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)
Expand All @@ -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
1 change: 1 addition & 0 deletions examples/symba_clement_2018/param.in
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
10 changes: 3 additions & 7 deletions examples/symba_energy_momentum/param.escape.in
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
6 changes: 4 additions & 2 deletions examples/symba_energy_momentum/param.sun.in
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
15 changes: 11 additions & 4 deletions src/discard/discard.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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
Expand All @@ -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
Expand Down
83 changes: 33 additions & 50 deletions src/io/io.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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")'
Expand All @@ -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
Expand All @@ -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
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -527,16 +510,15 @@ 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
ifirst = ilast + 1
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
Expand All @@ -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")
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down
13 changes: 4 additions & 9 deletions src/modules/swiftest_classes.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down
Loading

0 comments on commit 7b1f1cb

Please sign in to comment.