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

Commit

Permalink
Removed all initial properties values from the system object and keep…
Browse files Browse the repository at this point in the history 
… it now in the parameter object. I may move this back (and out of param) later, but for now, this is where they live. Fixed up some bookeeping mismatched resulting from having these variables in two places.
  • Loading branch information
@daminton
daminton committed Aug 19, 2021
1 parent d86888f commit b482943
Show file tree
Hide file tree
Showing 6 changed files with 55 additions and 74 deletions.
68 changes: 25 additions & 43 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) :: GMtot_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, 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)
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(:)
GMtot_now = cb%Gmass + sum(pl%Gmass(1:npl)) + system%GMescape
if (lfirst) then
Eorbit_orig = Eorbit_now
GMtot_orig = GMtot_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, GMtot_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 = (GMtot_now - GMtot_orig) / GMtot_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 @@ -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
Expand Down
9 changes: 2 additions & 7 deletions src/modules/swiftest_classes.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -61,11 +61,9 @@ 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) :: 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) :: 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
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) :: 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
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
12 changes: 6 additions & 6 deletions src/modules/symba_classes.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down Expand Up @@ -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
Expand All @@ -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
Expand All @@ -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
Expand All @@ -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
Expand Down
Loading

0 comments on commit b482943

Please sign in to comment.