From 785d8d8477bb490b81aca4ce6b96465fad049557 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Thu, 5 Jan 2023 09:31:36 -0500 Subject: [PATCH] Fixed a bunch of stuff related to taking collision snapshots, including consolidation of collision NetCDF into a single file. --- src/base/base_module.f90 | 14 +- src/collision/collision_generate.f90 | 36 +++-- src/collision/collision_io.f90 | 166 ++++++++++++++++------ src/collision/collision_module.f90 | 60 ++++---- src/collision/collision_resolve.f90 | 61 ++++---- src/collision/collision_util.f90 | 161 +++++++-------------- src/fraggle/fraggle_generate.f90 | 20 +-- src/fraggle/fraggle_util.f90 | 20 +-- src/netcdf_io/netcdf_io_module.f90 | 18 +-- src/swiftest/swiftest_io.f90 | 202 +++++++++++++-------------- src/swiftest/swiftest_module.f90 | 34 ++--- src/swiftest/swiftest_util.f90 | 26 ++-- src/symba/symba_discard.f90 | 28 ++-- 13 files changed, 442 insertions(+), 404 deletions(-) diff --git a/src/base/base_module.f90 b/src/base/base_module.f90 index fe2f0710b..cc8d06034 100644 --- a/src/base/base_module.f90 +++ b/src/base/base_module.f90 @@ -85,15 +85,15 @@ module base logical :: ltides = .false. !! Include tidal dissipation ! 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) :: E_orbit_orig = 0.0_DP !! Initial orbital energy 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) :: Lescape = 0.0_DP !! Angular momentum of bodies that escaped the system (used for bookeeping) + real(DP), dimension(NDIM) :: L_total_orig = 0.0_DP !! Initial total angular momentum vector + real(DP), dimension(NDIM) :: L_orbit_orig = 0.0_DP !! Initial orbital angular momentum + real(DP), dimension(NDIM) :: L_spin_orig = 0.0_DP !! Initial spin angular momentum vector + real(DP), dimension(NDIM) :: L_escape = 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) :: E_collisions = 0.0_DP !! Energy lost from system due to collisions + real(DP) :: E_untracked = 0.0_DP !! Energy gained from system due to escaped bodies logical :: lfirstenergy = .true. !! This is the first time computing energe logical :: lfirstkick = .true. !! Initiate the first kick in a symplectic step logical :: lrestart = .false. !! Indicates whether or not this is a restarted run diff --git a/src/collision/collision_generate.f90 b/src/collision/collision_generate.f90 index b946a387e..5cb2e68aa 100644 --- a/src/collision/collision_generate.f90 +++ b/src/collision/collision_generate.f90 @@ -55,6 +55,13 @@ module subroutine collision_generate_bounce(self, nbody_system, param, t) associate(impactors => nbody_system%collider%impactors, fragments => nbody_system%collider%fragments) select case (impactors%regime) case (COLLRESOLVE_REGIME_DISRUPTION, COLLRESOLVE_REGIME_SUPERCATASTROPHIC) + + ! Manually save the before/after snapshots because this case doesn't use the mergeaddsub procedure + select type(before => self%before) + class is (swiftest_nbody_system) + allocate(before%pl, source=pl) + end select + nfrag = size(impactors%id(:)) do i = 1, nfrag j = impactors%id(i) @@ -68,13 +75,12 @@ module subroutine collision_generate_bounce(self, nbody_system, param, t) pl%ldiscard(j) = .false. pl%lcollision(j) = .false. end do - select type(before => self%before) - class is (swiftest_nbody_system) + select type(after => self%after) class is (swiftest_nbody_system) - allocate(after%pl, source=before%pl) ! Be sure to save the pl so that snapshots still work - end select + allocate(after%pl, source=pl) end select + case (COLLRESOLVE_REGIME_HIT_AND_RUN) call self%hitandrun(nbody_system, param, t) case (COLLRESOLVE_REGIME_MERGE, COLLRESOLVE_REGIME_GRAZE_AND_MERGE) @@ -107,7 +113,6 @@ module subroutine collision_generate_hitandrun(self, nbody_system, param, t) ! Internals character(len=STRMAX) :: message - select type(nbody_system) class is (swiftest_nbody_system) select type(pl => nbody_system%pl) @@ -117,18 +122,21 @@ module subroutine collision_generate_hitandrun(self, nbody_system, param, t) call collision_io_collider_message(nbody_system%pl, impactors%id, message) call swiftest_io_log_one_message(COLLISION_LOG_OUT, trim(adjustl(message))) + ! Manually save the before/after snapshots because this case doesn't use the mergeaddsub procedure + select type(before => self%before) + class is (swiftest_nbody_system) + allocate(before%pl, source=pl) + end select + status = HIT_AND_RUN_PURE pl%status(impactors%id(:)) = ACTIVE pl%ldiscard(impactors%id(:)) = .false. pl%lcollision(impactors%id(:)) = .false. - ! Be sure to save the pl so that snapshots still work - select type(before => self%before) - class is (swiftest_nbody_system) + select type(after => self%after) class is (swiftest_nbody_system) - allocate(after%pl, source=before%pl) + allocate(after%pl, source=pl) end select - end select end associate end select @@ -154,7 +162,7 @@ module subroutine collision_generate_merge(self, nbody_system, param, t) real(DP), intent(in) :: t !! The time of the collision ! Internals integer(I4B) :: i, j, k, ibiggest - real(DP), dimension(NDIM) :: Lspin_new + real(DP), dimension(NDIM) :: L_spin_new real(DP) :: volume, G character(len=STRMAX) :: message @@ -193,12 +201,12 @@ module subroutine collision_generate_merge(self, nbody_system, param, t) if (param%lrotation) then do concurrent(i = 1:NDIM) fragments%Ip(i,1) = sum(impactors%mass(:) * impactors%Ip(i,:)) - Lspin_new(i) = sum(impactors%Lorbit(i,:) + impactors%Lorbit(i,:)) + L_spin_new(i) = sum(impactors%L_orbit(i,:) + impactors%L_orbit(i,:)) end do fragments%Ip(:,1) = fragments%Ip(:,1) / fragments%mass(1) - fragments%rot(:,1) = Lspin_new(:) / (fragments%Ip(3,1) * fragments%mass(1) * fragments%radius(1)**2) + fragments%rot(:,1) = L_spin_new(:) / (fragments%Ip(3,1) * fragments%mass(1) * fragments%radius(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 - nbody_system%Lescape(:) = nbody_system%Lescape(:) + impactors%Lorbit(:,1) + impactors%Lorbit(:,2) + nbody_system%L_escape(:) = nbody_system%L_escape(:) + impactors%L_orbit(:,1) + impactors%L_orbit(:,2) end if ! The fragment trajectory will be the barycentric trajectory diff --git a/src/collision/collision_io.f90 b/src/collision/collision_io.f90 index 820f69af1..583b19bb9 100644 --- a/src/collision/collision_io.f90 +++ b/src/collision/collision_io.f90 @@ -98,15 +98,12 @@ module subroutine collision_io_netcdf_dump(self, param) select type(nc => self%nc) class is (collision_netcdf_parameters) select type(param) - class is (swiftest_parameters) + class is (base_parameters) if (self%iframe > 0) then nc%file_number = nc%file_number + 1 call self%make_index_map() - nc%event_dimsize = self%nt - nc%name_dimsize = self%nid - write(nc%file_name, '("collision_",I0.6,".nc")') nc%file_number - call nc%initialize(param) + call nc%open(param) do i = 1, self%nframes if (allocated(self%frame(i)%item)) then @@ -175,70 +172,74 @@ module subroutine collision_io_netcdf_initialize_output(self, param) nc%lfile_is_open = .true. ! Dimensions - call netcdf_io_check( nf90_def_dim(nc%id, nc%event_dimname, nc%event_dimsize, nc%event_dimid), "collision_io_netcdf_initialize_output nf90_def_dim event_dimid" ) ! Dimension to store individual collision events - call netcdf_io_check( nf90_def_dim(nc%id, nc%space_dimname, NDIM, nc%space_dimid), "collision_io_netcdf_initialize_output nf90_def_dim space_dimid" ) ! 3D space dimension - call netcdf_io_check( nf90_def_dim(nc%id, nc%name_dimname, nc%name_dimsize, nc%name_dimid), "collision_io_netcdf_initialize_output nf90_def_dim name_dimid" ) ! Dimension to store particle id numbers - call netcdf_io_check( nf90_def_dim(nc%id, nc%str_dimname, NAMELEN, nc%str_dimid), "collision_io_netcdf_initialize_output nf90_def_dim str_dimid" ) ! Dimension for string variables (aka character arrays) - call netcdf_io_check( nf90_def_dim(nc%id, nc%stage_dimname, 2, nc%stage_dimid), "collision_io_netcdf_initialize_output nf90_def_dim stage_dimid" ) ! Dimension for stage variables (aka "before" vs. "after" + call netcdf_io_check( nf90_def_dim(nc%id, nc%collision_id_varname, NF90_UNLIMITED, nc%collision_id_dimid), "collision_io_netcdf_initialize_output nf90_def_dim collision_id_dimid" ) ! Dimension to store individual collision events + call netcdf_io_check( nf90_def_dim(nc%id, nc%space_dimname, NDIM, nc%space_dimid), "collision_io_netcdf_initialize_output nf90_def_dim space_dimid" ) ! 3D space dimension + call netcdf_io_check( nf90_def_dim(nc%id, nc%name_dimname, NF90_UNLIMITED, nc%name_dimid), "collision_io_netcdf_initialize_output nf90_def_dim name_dimid" ) ! Dimension to store particle id numbers + call netcdf_io_check( nf90_def_dim(nc%id, nc%str_dimname, NAMELEN, nc%str_dimid), "collision_io_netcdf_initialize_output nf90_def_dim str_dimid" ) ! Dimension for string variables (aka character arrays) + call netcdf_io_check( nf90_def_dim(nc%id, nc%stage_dimname, 2, nc%stage_dimid), "collision_io_netcdf_initialize_output nf90_def_dim stage_dimid" ) ! Dimension for stage variables (aka "before" vs. "after" ! Dimension coordinates - call netcdf_io_check( nf90_def_var(nc%id, nc%space_dimname, NF90_CHAR, nc%space_dimid, nc%space_varid), "collision_io_netcdf_initialize_output nf90_def_var space_varid" ) - call netcdf_io_check( nf90_def_var(nc%id, nc%name_dimname, NF90_CHAR, [nc%str_dimid, nc%name_dimid], nc%name_varid), "collision_io_netcdf_initialize_output nf90_def_var name_varid") - call netcdf_io_check( nf90_def_var(nc%id, nc%stage_dimname, NF90_CHAR, [nc%str_dimid, nc%stage_dimid], nc%stage_varid), "collision_io_netcdf_initialize_output nf90_def_var stage_varid" ) + call netcdf_io_check( nf90_def_var(nc%id, nc%collision_id_varname, NF90_INT, [nc%collision_id_dimid], nc%collision_id_varid), "collision_io_netcdf_initialize_output nf90_def_var collision_id_varid") + call netcdf_io_check( nf90_def_var(nc%id, nc%space_dimname, NF90_CHAR, nc%space_dimid, nc%space_varid), "collision_io_netcdf_initialize_output nf90_def_var space_varid" ) + call netcdf_io_check( nf90_def_var(nc%id, nc%name_dimname, NF90_CHAR, [nc%str_dimid, nc%name_dimid], nc%name_varid), "collision_io_netcdf_initialize_output nf90_def_var name_varid") + call netcdf_io_check( nf90_def_var(nc%id, nc%stage_dimname, NF90_CHAR, [nc%str_dimid, nc%stage_dimid], nc%stage_varid), "collision_io_netcdf_initialize_output nf90_def_var stage_varid" ) ! Variables call netcdf_io_check( nf90_def_var(nc%id, nc%id_varname, NF90_INT, nc%name_dimid, nc%id_varid), "collision_io_netcdf_initialize_output nf90_def_var id_varid" ) call netcdf_io_check( nf90_def_var(nc%id, nc%time_dimname, nc%out_type, & - nc%event_dimid, nc%time_varid), "collision_io_netcdf_initialize_output nf90_def_var time_varid" ) + nc%collision_id_dimid, nc%time_varid), "collision_io_netcdf_initialize_output nf90_def_var time_varid" ) call netcdf_io_check( nf90_def_var(nc%id, nc%regime_varname, NF90_CHAR, & - [nc%str_dimid, nc%event_dimid], nc%regime_varid), "collision_io_netcdf_initialize_output nf90_def_var regime_varid") + [nc%str_dimid, nc%collision_id_dimid], nc%regime_varid), "collision_io_netcdf_initialize_output nf90_def_var regime_varid") call netcdf_io_check( nf90_def_var(nc%id, nc%Qloss_varname, nc%out_type, & - [ nc%event_dimid], nc%Qloss_varid), "collision_io_netcdf_initialize_output nf90_def_var Qloss_varid") + [ nc%collision_id_dimid], nc%Qloss_varid), "collision_io_netcdf_initialize_output nf90_def_var Qloss_varid") call netcdf_io_check( nf90_def_var(nc%id, nc%ptype_varname, NF90_CHAR, & - [nc%str_dimid, nc%name_dimid, nc%stage_dimid, nc%event_dimid], nc%ptype_varid), "collision_io_netcdf_initialize_output nf90_def_var ptype_varid") + [nc%str_dimid, nc%name_dimid, nc%stage_dimid, nc%collision_id_dimid], nc%ptype_varid), "collision_io_netcdf_initialize_output nf90_def_var ptype_varid") call netcdf_io_check( nf90_def_var(nc%id, nc%loop_varname, NF90_INT, & - [ nc%event_dimid], nc%loop_varid), "collision_io_netcdf_initialize_output nf90_def_var loop_varid") + [ nc%collision_id_dimid], nc%loop_varid), "collision_io_netcdf_initialize_output nf90_def_var loop_varid") call netcdf_io_check( nf90_def_var(nc%id, nc%rh_varname, nc%out_type,& - [ nc%space_dimid, nc%name_dimid, nc%stage_dimid, nc%event_dimid], nc%rh_varid), "collision_io_netcdf_initialize_output nf90_def_var rh_varid") + [ nc%space_dimid, nc%name_dimid, nc%stage_dimid, nc%collision_id_dimid], nc%rh_varid), "collision_io_netcdf_initialize_output nf90_def_var rh_varid") call netcdf_io_check( nf90_def_var(nc%id, nc%vh_varname, nc%out_type,& - [ nc%space_dimid, nc%name_dimid, nc%stage_dimid, nc%event_dimid], nc%vh_varid), "collision_io_netcdf_initialize_output nf90_def_var vh_varid") + [ nc%space_dimid, nc%name_dimid, nc%stage_dimid, nc%collision_id_dimid], nc%vh_varid), "collision_io_netcdf_initialize_output nf90_def_var vh_varid") call netcdf_io_check( nf90_def_var(nc%id, nc%Gmass_varname, nc%out_type,& - [ nc%name_dimid, nc%stage_dimid, nc%event_dimid], nc%Gmass_varid), "collision_io_netcdf_initialize_output nf90_def_var Gmass_varid") + [ nc%name_dimid, nc%stage_dimid, nc%collision_id_dimid], nc%Gmass_varid), "collision_io_netcdf_initialize_output nf90_def_var Gmass_varid") call netcdf_io_check( nf90_def_var(nc%id, nc%radius_varname, nc%out_type,& - [ nc%name_dimid, nc%stage_dimid, nc%event_dimid], nc%radius_varid), "collision_io_netcdf_initialize_output nf90_def_var radius_varid") + [ nc%name_dimid, nc%stage_dimid, nc%collision_id_dimid], nc%radius_varid), "collision_io_netcdf_initialize_output nf90_def_var radius_varid") - call netcdf_io_check( nf90_def_var(nc%id, nc%Ip_varname, nc%out_type,& - [ nc%space_dimid, nc%name_dimid, nc%stage_dimid, nc%event_dimid], nc%Ip_varid), "collision_io_netcdf_initialize_output nf90_def_var Ip_varid") + if (param%lrotation) then + call netcdf_io_check( nf90_def_var(nc%id, nc%Ip_varname, nc%out_type,& + [ nc%space_dimid, nc%name_dimid, nc%stage_dimid, nc%collision_id_dimid], nc%Ip_varid), "collision_io_netcdf_initialize_output nf90_def_var Ip_varid") - call netcdf_io_check( nf90_def_var(nc%id, nc%rot_varname, nc%out_type,& - [ nc%space_dimid, nc%name_dimid, nc%stage_dimid, nc%event_dimid], nc%rot_varid), "collision_io_netcdf_initialize_output nf90_def_var rot_varid") - - if (param%lenergy) then + call netcdf_io_check( nf90_def_var(nc%id, nc%rot_varname, nc%out_type,& + [ nc%space_dimid, nc%name_dimid, nc%stage_dimid, nc%collision_id_dimid], nc%rot_varid), "collision_io_netcdf_initialize_output nf90_def_var rot_varid") + end if + if (param%lenergy) then call netcdf_io_check( nf90_def_var(nc%id, nc%ke_orb_varname, nc%out_type,& - [ nc%stage_dimid, nc%event_dimid], nc%KE_orb_varid), "collision_io_netcdf_initialize_output nf90_def_var KE_orb_varid") + [ nc%stage_dimid, nc%collision_id_dimid], nc%KE_orb_varid), "collision_io_netcdf_initialize_output nf90_def_var KE_orb_varid") call netcdf_io_check( nf90_def_var(nc%id, nc%ke_spin_varname, nc%out_type,& - [ nc%stage_dimid, nc%event_dimid], nc%KE_spin_varid), "collision_io_netcdf_initialize_output nf90_def_var KE_spin_varid" ) + [ nc%stage_dimid, nc%collision_id_dimid], nc%KE_spin_varid), "collision_io_netcdf_initialize_output nf90_def_var KE_spin_varid" ) call netcdf_io_check( nf90_def_var(nc%id, nc%pe_varname, nc%out_type,& - [ nc%stage_dimid, nc%event_dimid], nc%PE_varid), "collision_io_netcdf_initialize_output nf90_def_var PE_varid" ) + [ nc%stage_dimid, nc%collision_id_dimid], nc%PE_varid), "collision_io_netcdf_initialize_output nf90_def_var PE_varid" ) call netcdf_io_check( nf90_def_var(nc%id, nc%be_varname, nc%out_type,& - [ nc%stage_dimid, nc%event_dimid], nc%BE_varid), "collision_io_netcdf_initialize_output nf90_def_var BE_varid" ) + [ nc%stage_dimid, nc%collision_id_dimid], nc%BE_varid), "collision_io_netcdf_initialize_output nf90_def_var BE_varid" ) - call netcdf_io_check( nf90_def_var(nc%id, nc%L_orb_varname, nc%out_type, & - [ nc%space_dimid, nc%stage_dimid, nc%event_dimid], nc%L_orb_varid), "collision_io_netcdf_initialize_output nf90_def_var L_orb_varid" ) + if (param%lrotation) then + call netcdf_io_check( nf90_def_var(nc%id, nc%L_orbit_varname, nc%out_type, & + [ nc%space_dimid, nc%stage_dimid, nc%collision_id_dimid], nc%L_orbit_varid), "collision_io_netcdf_initialize_output nf90_def_var L_orbit_varid" ) - call netcdf_io_check( nf90_def_var(nc%id, nc%Lspin_varname, nc%out_type,& - [ nc%space_dimid, nc%stage_dimid, nc%event_dimid], nc%Lspin_varid), "collision_io_netcdf_initialize_output nf90_def_var Lspin_varid" ) + call netcdf_io_check( nf90_def_var(nc%id, nc%L_spin_varname, nc%out_type,& + [ nc%space_dimid, nc%stage_dimid, nc%collision_id_dimid], nc%L_spin_varid), "collision_io_netcdf_initialize_output nf90_def_var L_spin_varid" ) + end if end if call netcdf_io_check( nf90_inquire(nc%id, nVariables=nvar), "collision_io_netcdf_initialize_output nf90_inquire nVariables" ) @@ -274,6 +275,88 @@ module subroutine collision_io_netcdf_initialize_output(self, param) end subroutine collision_io_netcdf_initialize_output + module subroutine collision_io_netcdf_open(self, param, readonly) + !! author: Carlisle A. Wishard, Dana Singh, and David A. Minton + !! + !! Opens a NetCDF file and does the variable inquiries to activate variable ids + use netcdf + implicit none + ! Arguments + class(collision_netcdf_parameters), intent(inout) :: self !! Parameters used to identify a particular NetCDF dataset + class(base_parameters), intent(in) :: param !! Current run configuration parameters + logical, optional, intent(in) :: readonly !! Logical flag indicating that this should be open read only + ! Internals + integer(I4B) :: mode + character(len=STRMAX) :: errmsg + logical fileExists + + mode = NF90_WRITE + if (present(readonly)) then + if (readonly) mode = NF90_NOWRITE + end if + + select type(param) + class is (base_parameters) + associate(nc => self) + + inquire(file=nc%file_name, exist=fileExists) + if (.not.fileExists) then + call nc%initialize(param) + return + end if + + write(errmsg,*) "collision_io_netcdf_opennf90_open ",trim(adjustl(nc%file_name)) + call netcdf_io_check( nf90_open(nc%file_name, mode, nc%id), errmsg) + self%lfile_is_open = .true. + + ! Dimensions + call netcdf_io_check( nf90_inq_dimid(nc%id, nc%collision_id_varname, nc%time_dimid), "collision_io_netcdf_opennf90_inq_dimid time_dimid" ) + call netcdf_io_check( nf90_inq_dimid(nc%id, nc%space_dimname, nc%space_dimid), "collision_io_netcdf_opennf90_inq_dimid space_dimid" ) + call netcdf_io_check( nf90_inq_dimid(nc%id, nc%name_dimname, nc%name_dimid), "collision_io_netcdf_opennf90_inq_dimid name_dimid" ) + call netcdf_io_check( nf90_inq_dimid(nc%id, nc%str_dimname, nc%str_dimid), "collision_io_netcdf_opennf90_inq_dimid str_dimid" ) + call netcdf_io_check( nf90_inq_dimid(nc%id, nc%stage_dimname, nc%stage_dimid), "collision_io_netcdf_opennf90_inq_dimid stage_dimid" ) + + ! Dimension coordinates + call netcdf_io_check( nf90_inq_varid(nc%id, nc%collision_id_varname, nc%time_varid), "collision_io_netcdf_opennf90_inq_varid time_varid" ) + call netcdf_io_check( nf90_inq_varid(nc%id, nc%space_dimname, nc%space_varid), "collision_io_netcdf_opennf90_inq_varid space_varid" ) + call netcdf_io_check( nf90_inq_varid(nc%id, nc%name_dimname, nc%name_varid), "collision_io_netcdf_opennf90_inq_varid name_varid" ) + call netcdf_io_check( nf90_inq_varid(nc%id, nc%stage_dimname, nc%stage_varid), "collision_io_netcdf_opennf90_inq_varid stage_varid" ) + + ! Required Variables + call netcdf_io_check( nf90_inq_varid(nc%id, nc%id_varname, nc%id_varid), "collision_io_netcdf_opennf90_inq_varid name_varid" ) + call netcdf_io_check( nf90_inq_varid(nc%id, nc%time_dimname, nc%time_varid), "collision_io_netcdf_opennf90_inq_varid time_varid" ) + call netcdf_io_check( nf90_inq_varid(nc%id, nc%regime_varname, nc%regime_varid), "collision_io_netcdf_opennf90_inq_varid regime_varid" ) + call netcdf_io_check( nf90_inq_varid(nc%id, nc%Qloss_varname, nc%Qloss_varid), "collision_io_netcdf_opennf90_inq_varid Qloss_varid" ) + call netcdf_io_check( nf90_inq_varid(nc%id, nc%ptype_varname, nc%ptype_varid), "collision_io_netcdf_opennf90_inq_varid ptype_varid" ) + call netcdf_io_check( nf90_inq_varid(nc%id, nc%loop_varname, nc%loop_varid), "collision_io_netcdf_opennf90_inq_varid loop_varid" ) + call netcdf_io_check( nf90_inq_varid(nc%id, nc%rh_varname, nc%rh_varid), "collision_io_netcdf_opennf90_inq_varid rh_varid" ) + call netcdf_io_check( nf90_inq_varid(nc%id, nc%vh_varname, nc%vh_varid), "collision_io_netcdf_opennf90_inq_varid vh_varid" ) + call netcdf_io_check( nf90_inq_varid(nc%id, nc%Gmass_varname, nc%Gmass_varid), "collision_io_netcdf_opennf90_inq_varid Gmass_varid" ) + call netcdf_io_check( nf90_inq_varid(nc%id, nc%radius_varname, nc%radius_varid), "collision_io_netcdf_opennf90_inq_varid radius_varid" ) + + if (param%lrotation) then + call netcdf_io_check( nf90_inq_varid(nc%id, nc%Ip_varname, nc%Ip_varid), "collision_io_netcdf_opennf90_inq_varid Ip_varid" ) + call netcdf_io_check( nf90_inq_varid(nc%id, nc%rot_varname, nc%rot_varid), "collision_io_netcdf_opennf90_inq_varid rot_varid" ) + end if + + if (param%lenergy) then + call netcdf_io_check( nf90_inq_varid(nc%id, nc%ke_orb_varname, nc%ke_orb_varid), "collision_io_netcdf_opennf90_inq_varid ke_orb_varid" ) + call netcdf_io_check( nf90_inq_varid(nc%id, nc%pe_varname, nc%pe_varid), "collision_io_netcdf_opennf90_inq_varid pe_varid" ) + call netcdf_io_check( nf90_inq_varid(nc%id, nc%be_varname, nc%be_varid), "collision_io_netcdf_opennf90_inq_varid be_varid" ) + call netcdf_io_check( nf90_inq_varid(nc%id, nc%L_orbit_varname, nc%L_orbit_varid), "collision_io_netcdf_opennf90_inq_varid L_orbit_varid" ) + if (param%lrotation) then + call netcdf_io_check( nf90_inq_varid(nc%id, nc%ke_spin_varname, nc%ke_spin_varid), "collision_io_netcdf_opennf90_inq_varid ke_spin_varid" ) + call netcdf_io_check( nf90_inq_varid(nc%id, nc%L_spin_varname, nc%L_spin_varid), "collision_io_netcdf_opennf90_inq_varid L_spin_varid" ) + end if + end if + + end associate + end select + + return + end subroutine collision_io_netcdf_open + + module subroutine collision_io_netcdf_write_frame_snapshot(self, history, param) !! author: David A. Minton !! @@ -325,19 +408,22 @@ module subroutine collision_io_netcdf_write_frame_snapshot(self, history, param) call netcdf_io_check( nf90_put_var(nc%id, nc%vh_varid, pl%vh(:,i), start=[1, idslot, stage, eslot], count=[NDIM,1,1,1]), "collision_io_netcdf_write_frame_snapshot nf90_put_var vh_varid" ) call netcdf_io_check( nf90_put_var(nc%id, nc%Gmass_varid, pl%Gmass(i), start=[ idslot, stage, eslot]), "collision_io_netcdf_write_frame_snapshot nf90_put_var Gmass_varid" ) call netcdf_io_check( nf90_put_var(nc%id, nc%radius_varid, pl%radius(i), start=[ idslot, stage, eslot]), "collision_io_netcdf_write_frame_snapshot nf90_put_var radius_varid" ) - call netcdf_io_check( nf90_put_var(nc%id, nc%Ip_varid, pl%Ip(:,i), start=[1, idslot, stage, eslot], count=[NDIM,1,1,1]), "collision_io_netcdf_write_frame_snapshot nf90_put_var Ip_varid" ) - call netcdf_io_check( nf90_put_var(nc%id, nc%rot_varid, pl%rot(:,i), start=[1, idslot, stage, eslot], count=[NDIM,1,1,1]), "collision_io_netcdf_write_frame_snapshot nf90_put_var rotx_varid" ) + if (param%lrotation) then + call netcdf_io_check( nf90_put_var(nc%id, nc%Ip_varid, pl%Ip(:,i), start=[1, idslot, stage, eslot], count=[NDIM,1,1,1]), "collision_io_netcdf_write_frame_snapshot nf90_put_var Ip_varid" ) + call netcdf_io_check( nf90_put_var(nc%id, nc%rot_varid, pl%rot(:,i), start=[1, idslot, stage, eslot], count=[NDIM,1,1,1]), "collision_io_netcdf_write_frame_snapshot nf90_put_var rotx_varid" ) + end if end do end do end select end select + if (param%lenergy) then call netcdf_io_check( nf90_put_var(nc%id, nc%ke_orb_varid, collider%ke_orbit(:), start=[ 1, eslot], count=[ 2, 1]), "collision_io_netcdf_write_frame_snapshot nf90_put_var ke_orb_varid before" ) call netcdf_io_check( nf90_put_var(nc%id, nc%ke_spin_varid, collider%ke_spin(:), start=[ 1, eslot], count=[ 2, 1]), "collision_io_netcdf_write_frame_snapshot nf90_put_var ke_spin_varid before" ) call netcdf_io_check( nf90_put_var(nc%id, nc%pe_varid, collider%pe(:), start=[ 1, eslot], count=[ 2, 1]), "collision_io_netcdf_write_frame_snapshot nf90_put_var pe_varid before" ) call netcdf_io_check( nf90_put_var(nc%id, nc%be_varid, collider%be(:), start=[ 1, eslot], count=[ 2, 1]), "collision_io_netcdf_write_frame_snapshot nf90_put_var pe_varid before" ) - call netcdf_io_check( nf90_put_var(nc%id, nc%L_orb_varid, collider%Lorbit(:,:), start=[1, 1, eslot], count=[NDIM, 2, 1]), "collision_io_netcdf_write_frame_snapshot nf90_put_var L_orb_varid before" ) - call netcdf_io_check( nf90_put_var(nc%id, nc%Lspin_varid, collider%Lspin(:,:), start=[1, 1, eslot], count=[NDIM, 2, 1]), "collision_io_netcdf_write_frame_snapshot nf90_put_var Lspin_varid before" ) + call netcdf_io_check( nf90_put_var(nc%id, nc%L_orbit_varid, collider%L_orbit(:,:), start=[1, 1, eslot], count=[NDIM, 2, 1]), "collision_io_netcdf_write_frame_snapshot nf90_put_var L_orbit_varid before" ) + call netcdf_io_check( nf90_put_var(nc%id, nc%L_spin_varid, collider%L_spin(:,:), start=[1, 1, eslot], count=[NDIM, 2, 1]), "collision_io_netcdf_write_frame_snapshot nf90_put_var L_spin_varid before" ) end if call netcdf_io_check( nf90_set_fill(nc%id, old_mode, old_mode) ) diff --git a/src/collision/collision_module.f90 b/src/collision/collision_module.f90 index 3fa91a220..73e7b4243 100644 --- a/src/collision/collision_module.f90 +++ b/src/collision/collision_module.f90 @@ -18,7 +18,8 @@ module collision implicit none public - character(len=*), parameter :: COLLISION_LOG_OUT = "collision.log" !! Name of log file for collision diagnostic information + character(len=*), parameter :: COLLISION_OUTFILE = 'collision.nc' !! Name of NetCDF output file for collision information + character(len=*), parameter :: COLLISION_LOG_OUT = "collision.log" !! Name of log file for collision diagnostic information !>Symbolic names for collisional outcomes from collresolve_resolve: integer(I4B), parameter :: COLLRESOLVE_REGIME_MERGE = 1 @@ -57,8 +58,8 @@ module collision real(DP), dimension(NDIM,2) :: rc !! Two-body equivalent position vectors of the collider bodies prior to collision in collision center of mass coordinates real(DP), dimension(NDIM,2) :: vc !! Two-body equivalent velocity vectors of the collider bodies prior to collision in collision center of mass coordinates real(DP), dimension(NDIM,2) :: rot !! Two-body equivalent principal axes moments of inertia the collider bodies prior to collision - real(DP), dimension(NDIM,2) :: Lspin !! Two-body equivalent spin angular momentum vectors of the collider bodies prior to collision - real(DP), dimension(NDIM,2) :: Lorbit !! Two-body equivalent orbital angular momentum vectors of the collider bodies prior to collision + real(DP), dimension(NDIM,2) :: L_spin !! Two-body equivalent spin angular momentum vectors of the collider bodies prior to collision + real(DP), dimension(NDIM,2) :: L_orbit !! Two-body equivalent orbital angular momentum vectors of the collider bodies prior to collision real(DP), dimension(NDIM,2) :: Ip !! Two-body equivalent principal axes moments of inertia the collider bodies prior to collision real(DP), dimension(2) :: Gmass !! Two-body equivalent G*mass of the collider bodies prior to the collision real(DP), dimension(2) :: mass !! Two-body equivalent mass of the collider bodies prior to the collision @@ -112,10 +113,10 @@ module collision real(DP), dimension(NDIM,nbody) :: t_unit !! Array of tangential direction unit vectors of individual fragments in the collisional coordinate frame real(DP), dimension(NDIM,nbody) :: n_unit !! Array of normal direction unit vectors of individual fragments in the collisional coordinate frame integer(I1B), dimension(nbody) :: origin_body !! Array of indices indicating which impactor body (1 or 2) the fragment originates from - real(DP), dimension(NDIM) :: Lorbit_tot !! Orbital angular momentum vector of all fragments - real(DP), dimension(NDIM) :: Lspin_tot !! Spin angular momentum vector of all fragments - real(DP), dimension(NDIM,nbody) :: Lorbit !! Orbital angular momentum vector of each individual fragment - real(DP), dimension(NDIM,nbody) :: Lspin !! Spin angular momentum vector of each individual fragment + real(DP), dimension(NDIM) :: L_orbit_tot !! Orbital angular momentum vector of all fragments + real(DP), dimension(NDIM) :: L_spin_tot !! Spin angular momentum vector of all fragments + real(DP), dimension(NDIM,nbody) :: L_orbit !! Orbital angular momentum vector of each individual fragment + real(DP), dimension(NDIM,nbody) :: L_spin !! Spin angular momentum vector of each individual fragment real(DP) :: ke_orbit_tot !! Orbital kinetic energy of all fragments real(DP) :: ke_spin_tot !! Spin kinetic energy of all fragments real(DP) :: pe !! Potential energy of all fragments @@ -138,27 +139,27 @@ module collision !! to resolve collision by defining extended types of encounters_impactors and/or encounetr_fragments !! !! The generate method for this class is the merge model. This allows any extended type to have access to the merge procedure by selecting the collision_basic parent class - class(collision_fragments(:)), allocatable :: fragments !! Object containing information on the pre-collision system - class(collision_impactors), allocatable :: impactors !! Object containing information on the post-collision system - class(base_nbody_system), allocatable :: before !! A snapshot of the subset of the nbody_system involved in the collision - class(base_nbody_system), allocatable :: after !! A snapshot of the subset of the nbody_system containing products of the collision - integer(I4B) :: status !! Status flag to pass to the collision list once the collision has been resolved + class(collision_fragments(:)), allocatable :: fragments !! Object containing information on the pre-collision system + class(collision_impactors), allocatable :: impactors !! Object containing information on the post-collision system + class(base_nbody_system), allocatable :: before !! A snapshot of the subset of the nbody_system involved in the collision + class(base_nbody_system), allocatable :: after !! A snapshot of the subset of the nbody_system containing products of the collision + integer(I4B) :: status !! Status flag to pass to the collision list once the collision has been resolved + integer(I4B) :: collision_id !! ID number of this collision event ! For the following variables, index 1 refers to the *entire* n-body nbody_system in its pre-collisional state and index 2 refers to the nbody_system in its post-collisional state - real(DP), dimension(NDIM,2) :: Lorbit !! Before/after orbital angular momentum - real(DP), dimension(NDIM,2) :: Lspin !! Before/after spin angular momentum - real(DP), dimension(NDIM,2) :: Ltot !! Before/after total nbody_system angular momentum + real(DP), dimension(NDIM,2) :: L_orbit !! Before/after orbital angular momentum + real(DP), dimension(NDIM,2) :: L_spin !! Before/after spin angular momentum + real(DP), dimension(NDIM,2) :: L_total !! Before/after total nbody_system angular momentum real(DP), dimension(2) :: ke_orbit !! Before/after orbital kinetic energy real(DP), dimension(2) :: ke_spin !! Before/after spin kinetic energy real(DP), dimension(2) :: pe !! Before/after potential energy real(DP), dimension(2) :: be !! Before/after binding energy - real(DP), dimension(2) :: Etot !! Before/after total nbody_system energy + real(DP), dimension(2) :: te !! Before/after total system energy contains procedure :: setup => collision_util_setup_collider !! Initializer for the encounter collision system and the before/after snapshots procedure :: setup_impactors => collision_util_setup_impactors_collider !! Initializer for the impactors for the encounter collision system. Deallocates old impactors before creating new ones procedure :: setup_fragments => collision_util_setup_fragments_collider !! Initializer for the fragments of the collision system. procedure :: add_fragments => collision_util_add_fragments_to_collider !! Add fragments to nbody_system - procedure :: construct_temporary_system => collision_util_construct_temporary_system !! Constructs temporary n-body nbody_system in order to compute pre- and post-impact energy and momentum procedure :: get_energy_and_momentum => collision_util_get_energy_momentum !! Calculates total nbody_system energy in either the pre-collision outcome state (lbefore = .true.) or the post-collision outcome state (lbefore = .false.) procedure :: reset => collision_util_reset_system !! Deallocates all allocatables procedure :: set_budgets => collision_util_set_budgets !! Sets the energy and momentum budgets of the fragments based on the collider value @@ -168,6 +169,7 @@ module collision procedure :: merge => collision_generate_merge !! Merges the impactors to make a single final body end type collision_basic + type, extends(collision_basic) :: collision_bounce contains procedure :: generate => collision_generate_bounce !! If a collision would result in a disruption, "bounce" the bodies instead. @@ -175,8 +177,6 @@ module collision end type collision_bounce - - !! NetCDF dimension and variable names for the enounter save object type, extends(encounter_netcdf_parameters) :: collision_netcdf_parameters integer(I4B) :: stage_dimid !! ID for the stage dimension @@ -184,10 +184,7 @@ module collision character(NAMELEN) :: stage_dimname = "stage" !! name of the stage dimension (before/after) character(len=6), dimension(2) :: stage_coords = ["before", "after"] !! The stage coordinate labels - character(NAMELEN) :: event_dimname = "collision" !! Name of collision event dimension - integer(I4B) :: event_dimid !! ID for the collision event dimension - integer(I4B) :: event_varid !! ID for the collision event variable - integer(I4B) :: event_dimsize = 0 !! Number of events + integer(I4B) :: collision_id_dimid !! ID for the collision event dimension character(NAMELEN) :: Qloss_varname = "Qloss" !! name of the energy loss variable integer(I4B) :: Qloss_varid !! ID for the energy loss variable @@ -195,6 +192,7 @@ module collision integer(I4B) :: regime_varid !! ID for the collision regime variable contains procedure :: initialize => collision_io_netcdf_initialize_output !! Initialize a set of parameters used to identify a NetCDF output object + procedure :: open => collision_io_netcdf_open !! Opens an old file final :: collision_final_netcdf_parameters !! Finalizer closes the NetCDF file end type collision_netcdf_parameters @@ -276,6 +274,13 @@ module subroutine collision_io_netcdf_initialize_output(self, param) class(base_parameters), intent(in) :: param !! Current run configuration parameters end subroutine collision_io_netcdf_initialize_output + module subroutine collision_io_netcdf_open(self, param, readonly) + implicit none + class(collision_netcdf_parameters), intent(inout) :: self !! Parameters used to identify a particular NetCDF dataset + class(base_parameters), intent(in) :: param !! Current run configuration parameters + logical, optional, intent(in) :: readonly !! Logical flag indicating that this should be open read only + end subroutine collision_io_netcdf_open + module subroutine collision_io_netcdf_write_frame_snapshot(self, history, param) implicit none class(collision_snapshot), intent(in) :: self !! Swiftest encounter structure @@ -375,15 +380,6 @@ module subroutine collision_util_add_fragments_to_collider(self, nbody_system, p class(base_parameters), intent(in) :: param !! Current swiftest run configuration parameters end subroutine collision_util_add_fragments_to_collider - module subroutine collision_util_construct_temporary_system(self, nbody_system, param, tmpsys, tmpparam) - implicit none - class(collision_basic), intent(inout) :: self !! Collision system object - class(base_nbody_system), intent(in) :: nbody_system !! Original swiftest nbody system object - class(base_parameters), intent(in) :: param !! Current swiftest run configuration parameters - class(base_nbody_system), allocatable, intent(out) :: tmpsys !! Output temporary swiftest nbody system object - class(base_parameters), allocatable, intent(out) :: tmpparam !! Output temporary configuration run parameters - end subroutine collision_util_construct_temporary_system - module subroutine collision_util_get_angular_momentum(self) implicit none class(collision_fragments(*)), intent(inout) :: self !! Fragment system object diff --git a/src/collision/collision_resolve.f90 b/src/collision/collision_resolve.f90 index a745f7b3e..e33f9b3b3 100644 --- a/src/collision/collision_resolve.f90 +++ b/src/collision/collision_resolve.f90 @@ -16,7 +16,7 @@ module subroutine collision_resolve_consolidate_impactors(self, nbody_system, pa !! author: David A. Minton !! !! Loops through the pl-pl collision list and groups families together by index. Outputs the indices of all impactors%id members, - !! and pairs of quantities (x and v vectors, mass, radius, Lspin, and Ip) that can be used to resolve the collisional outcome. + !! and pairs of quantities (x and v vectors, mass, radius, L_spin, and Ip) that can be used to resolve the collisional outcome. implicit none ! Arguments class(collision_impactors), intent(out) :: self !! Collision impactors object @@ -83,7 +83,7 @@ module subroutine collision_resolve_consolidate_impactors(self, nbody_system, pa impactors%ncoll = count(pl%lcollision(impactors%id(:))) impactors%id = pack(impactors%id(:), pl%lcollision(impactors%id(:))) - impactors%Lspin(:,:) = 0.0_DP + impactors%L_spin(:,:) = 0.0_DP impactors%Ip(:,:) = 0.0_DP ! Find the barycenter of each body along with its children, if it has any @@ -93,7 +93,7 @@ module subroutine collision_resolve_consolidate_impactors(self, nbody_system, pa ! Assume principal axis rotation about axis corresponding to highest moment of inertia (3rd Ip) if (param%lrotation) then impactors%Ip(:, j) = impactors%mass(j) * pl%Ip(:, idx_parent(j)) - impactors%Lspin(:, j) = impactors%Ip(3, j) * impactors%radius(j)**2 * pl%rot(:, idx_parent(j)) + impactors%L_spin(:, j) = impactors%Ip(3, j) * impactors%radius(j)**2 * pl%rot(:, idx_parent(j)) end if if (nchild(j) > 0) then @@ -113,14 +113,14 @@ module subroutine collision_resolve_consolidate_impactors(self, nbody_system, pa xc(:) = impactors%rb(:, j) - xcom(:) vc(:) = impactors%vb(:, j) - vcom(:) xcrossv(:) = xc(:) .cross. vc(:) - impactors%Lspin(:, j) = impactors%Lspin(:, j) + impactors%mass(j) * xcrossv(:) + impactors%L_spin(:, j) = impactors%L_spin(:, j) + impactors%mass(j) * xcrossv(:) xc(:) = xchild(:) - xcom(:) vc(:) = vchild(:) - vcom(:) xcrossv(:) = xc(:) .cross. vc(:) - impactors%Lspin(:, j) = impactors%Lspin(:, j) + mchild * xcrossv(:) + impactors%L_spin(:, j) = impactors%L_spin(:, j) + mchild * xcrossv(:) - impactors%Lspin(:, j) = impactors%Lspin(:, j) + mchild * pl%Ip(3, idx_child) & + impactors%L_spin(:, j) = impactors%L_spin(:, j) + mchild * pl%Ip(3, idx_child) & * pl%radius(idx_child)**2 & * pl%rot(:, idx_child) impactors%Ip(:, j) = impactors%Ip(:, j) + mchild * pl%Ip(:, idx_child) @@ -144,7 +144,7 @@ module subroutine collision_resolve_consolidate_impactors(self, nbody_system, pa mxc(:, 2) = impactors%mass(2) * (impactors%rb(:, 2) - xcom(:)) vcc(:, 1) = impactors%vb(:, 1) - vcom(:) vcc(:, 2) = impactors%vb(:, 2) - vcom(:) - impactors%Lorbit(:,:) = mxc(:,:) .cross. vcc(:,:) + impactors%L_orbit(:,:) = mxc(:,:) .cross. vcc(:,:) ! Destroy the kinship relationships for all members of this impactors%id call pl%reset_kinship(impactors%id(:)) @@ -334,9 +334,6 @@ module subroutine collision_resolve_mergeaddsub(nbody_system, param, t, status) nimpactors = impactors%ncoll nfrag = fragments%nbody - param%maxid_collision = max(param%maxid_collision, maxval(nbody_system%pl%info(:)%collision_id)) - param%maxid_collision = param%maxid_collision + 1 - ! Setup new bodies allocate(plnew, mold=pl) call plnew%setup(nfrag, param) @@ -440,12 +437,6 @@ module subroutine collision_resolve_mergeaddsub(nbody_system, param, t, status) plnew%info(1:nfrag)%particle_type = PL_TYPE_NAME end where - ! Log the properties of the new bodies - select type(after => collider%after) - class is (swiftest_nbody_system) - allocate(after%pl, source=plnew) - end select - ! Append the new merged body to the list call pl_adds%append(plnew, lsource_mask=[(.true., i=1, nfrag)]) @@ -457,16 +448,22 @@ module subroutine collision_resolve_mergeaddsub(nbody_system, param, t, status) lmask(:) = .false. lmask(impactors%id(:)) = .true. - call plnew%setup(0, param) - deallocate(plnew) - allocate(plsub, mold=pl) call pl%spill(plsub, lmask, ldestructive=.false.) call pl_discards%append(plsub, lsource_mask=[(.true., i = 1, nimpactors)]) - call plsub%setup(0, param) - deallocate(plsub) + ! Log the properties of the old and new bodies + select type(before => collider%before) + class is (swiftest_nbody_system) + call move_alloc(plsub, before%pl) + end select + + select type(after => collider%after) + class is (swiftest_nbody_system) + call move_alloc(plnew, after%pl) + end select + end associate end select @@ -490,14 +487,13 @@ module subroutine collision_resolve_plpl(self, nbody_system, param, t, dt, irec) real(DP), intent(in) :: dt !! Current simulation step size integer(I4B), intent(in) :: irec !! Current recursion level ! Internals - real(DP) :: Eorbit_before, Eorbit_after + real(DP) :: E_orbit_before, E_orbit_after logical :: lplpl_collision character(len=STRMAX) :: timestr integer(I4B), dimension(2) :: idx_parent !! Index of the two bodies considered the "parents" of the collision logical :: lgoodcollision integer(I4B) :: i, loop, ncollisions integer(I4B), parameter :: MAXCASCADE = 1000 - real(DP) :: dpe select type (nbody_system) class is (swiftest_nbody_system) @@ -518,7 +514,7 @@ module subroutine collision_resolve_plpl(self, nbody_system, param, t, dt, irec) ! Get the energy before the collision is resolved if (param%lenergy) then call nbody_system%get_energy_and_momentum(param) - Eorbit_before = nbody_system%te + E_orbit_before = nbody_system%te end if do loop = 1, MAXCASCADE @@ -538,17 +534,20 @@ module subroutine collision_resolve_plpl(self, nbody_system, param, t, dt, irec) idx_parent(2) = pl%kin(idx2(i))%parent call impactors%consolidate(nbody_system, param, idx_parent, lgoodcollision) if ((.not. lgoodcollision) .or. any(pl%status(idx_parent(:)) /= COLLIDED)) cycle + + ! Advance the collision id number and save it + param%maxid_collision = max(param%maxid_collision, maxval(nbody_system%pl%info(:)%collision_id)) + param%maxid_collision = param%maxid_collision + 1 + collider%collision_id = param%maxid_collision + + ! Get the collision regime call impactors%get_regime(nbody_system, param) call collision_history%take_snapshot(param,nbody_system, t, "before") - call nbody_system%get_energy_and_momentum(param) - collider%pe(1) = nbody_system%pe - + ! Generate the new bodies resulting from the collision call collider%generate(nbody_system, param, t) - call nbody_system%get_energy_and_momentum(param) - call collision_history%take_snapshot(param,nbody_system, t, "after") plpl_collision%status(i) = collider%status @@ -585,8 +584,8 @@ module subroutine collision_resolve_plpl(self, nbody_system, param, t, dt, irec) if (param%lenergy) then call nbody_system%get_energy_and_momentum(param) - Eorbit_after = nbody_system%te - nbody_system%Ecollisions = nbody_system%Ecollisions + (Eorbit_after - Eorbit_before) + E_orbit_after = nbody_system%te + nbody_system%E_collisions = nbody_system%E_collisions + (E_orbit_after - E_orbit_before) end if end associate diff --git a/src/collision/collision_util.f90 b/src/collision/collision_util.f90 index 7f073eecf..882bb8f29 100644 --- a/src/collision/collision_util.f90 +++ b/src/collision/collision_util.f90 @@ -61,67 +61,6 @@ module subroutine collision_util_add_fragments_to_collider(self, nbody_system, p end subroutine collision_util_add_fragments_to_collider - module subroutine collision_util_construct_temporary_system(self, nbody_system, param, tmpsys, tmpparam) - !! Author: David A. Minton - !! - !! Constructs a temporary internal system consisting of active bodies and additional fragments. This internal temporary system is used to calculate system energy with and without fragments - implicit none - ! Arguments - class(collision_basic), intent(inout) :: self !! Fraggle collision system object - class(base_nbody_system), intent(in) :: nbody_system !! Original swiftest nbody system object - class(base_parameters), intent(in) :: param !! Current swiftest run configuration parameters - class(base_nbody_system), allocatable, intent(out) :: tmpsys !! Output temporary swiftest nbody system object - class(base_parameters), allocatable, intent(out) :: tmpparam !! Output temporary configuration run parameters - ! Internals - logical, dimension(:), allocatable :: linclude - integer(I4B) :: npl_tot - ! The following are needed in order to deal with typing requirements - class(swiftest_nbody_system), allocatable :: tmpsys_local - class(swiftest_parameters), allocatable :: tmpparam_local - - select type(nbody_system) - class is (swiftest_nbody_system) - select type(param) - class is (swiftest_parameters) - associate(fragments => self%fragments, nfrag => self%fragments%nbody, pl => nbody_system%pl, npl => nbody_system%pl%nbody, cb => nbody_system%cb) - ! Set up a new system based on the original - if (allocated(tmpparam)) deallocate(tmpparam) - if (allocated(tmpsys)) deallocate(tmpsys) - allocate(tmpparam_local, source=param) - select type(tmpparam_local) - class is (swiftest_parameters) - tmpparam_local%system_history%nc%lfile_is_open = .false. - end select - call swiftest_util_setup_construct_system(tmpsys_local, tmpparam_local) - - ! No test particles necessary for energy/momentum calcs - call tmpsys_local%tp%setup(0, tmpparam_local) - - ! Replace the empty central body object with a copy of the original - deallocate(tmpsys_local%cb) - allocate(tmpsys_local%cb, source=cb) - - ! Make space for the fragments - npl_tot = npl + nfrag - call tmpsys_local%pl%setup(npl_tot, tmpparam_local) - allocate(linclude(npl_tot)) - - ! Fill up the temporary system with all of the original bodies, leaving the spaces for fragments empty until we add them in later - linclude(1:npl) = .true. - linclude(npl+1:npl_tot) = .false. - call tmpsys_local%pl%fill(pl, linclude) - - call move_alloc(tmpsys_local, tmpsys) - call move_alloc(tmpparam_local, tmpparam) - - end associate - end select - end select - - return - end subroutine collision_util_construct_temporary_system - - module subroutine collision_util_get_idvalues_snapshot(self, idvals) !! author: David A. Minton !! @@ -188,12 +127,12 @@ module subroutine collision_util_get_angular_momentum(self) associate(fragments => self, nfrag => self%nbody) do i = 1, nfrag - fragments%Lorbit(:,i) = fragments%mass(i) * (fragments%rc(:,i) .cross. fragments%vc(:, i)) - fragments%Lspin(:,i) = fragments%mass(i) * fragments%radius(i)**2 * fragments%Ip(:,i) * fragments%rot(:,i) + fragments%L_orbit(:,i) = fragments%mass(i) * (fragments%rc(:,i) .cross. fragments%vc(:, i)) + fragments%L_spin(:,i) = fragments%mass(i) * fragments%radius(i)**2 * fragments%Ip(:,i) * fragments%rot(:,i) end do - fragments%Lorbit_tot(:) = sum(fragments%Lorbit, dim=2) - fragments%Lspin_tot(:) = sum(fragments%Lspin, dim=2) + fragments%L_orbit_tot(:) = sum(fragments%L_orbit, dim=2) + fragments%L_spin_tot(:) = sum(fragments%L_spin, dim=2) end associate return @@ -253,9 +192,8 @@ module subroutine collision_util_get_energy_momentum(self, nbody_system, param, logical, intent(in) :: lbefore !! Flag indicating that this the "before" state of the nbody_system, with impactors included and fragments excluded or vice versa ! Internals integer(I4B) :: stage,i,j - real(DP), dimension(NDIM) :: Lorbit, Lspin + real(DP), dimension(NDIM) :: L_orbit, L_spin real(DP) :: ke_orbit, ke_spin, pe, be - real(DP), dimension(self%fragments%nbody) :: pepl select type(nbody_system) class is (swiftest_nbody_system) @@ -265,8 +203,8 @@ module subroutine collision_util_get_energy_momentum(self, nbody_system, param, if (lbefore) then - Lorbit(:) = sum(impactors%Lorbit(:,:), dim=2) - Lspin(:) = sum(impactors%Lspin(:,:), dim=2) + L_orbit(:) = sum(impactors%L_orbit(:,:), dim=2) + L_spin(:) = sum(impactors%L_spin(:,:), dim=2) ke_orbit = 0.0_DP ke_spin = 0.0_DP do concurrent(i = 1:2) @@ -281,8 +219,8 @@ module subroutine collision_util_get_energy_momentum(self, nbody_system, param, else call fragments%get_angular_momentum() - Lorbit(:) = fragments%Lorbit_tot(:) - Lspin(:) = fragments%Lspin_tot(:) + L_orbit(:) = fragments%L_orbit_tot(:) + L_spin(:) = fragments%L_spin_tot(:) call fragments%get_energy() ke_orbit = fragments%ke_orbit_tot @@ -297,14 +235,14 @@ module subroutine collision_util_get_energy_momentum(self, nbody_system, param, else stage = 2 end if - self%Lorbit(:,stage) = Lorbit(:) - self%Lspin(:,stage) = Lspin(:) - self%Ltot(:,stage) = Lorbit(:) + Lspin(:) + self%L_orbit(:,stage) = L_orbit(:) + self%L_spin(:,stage) = L_spin(:) + self%L_total(:,stage) = L_orbit(:) + L_spin(:) self%ke_orbit(stage) = ke_orbit self%ke_spin(stage) = ke_spin self%pe(stage) = pe self%be(stage) = be - self%Etot(stage) = ke_orbit + ke_spin + pe + be + self%te(stage) = ke_orbit + ke_spin + pe + be end associate end select end select @@ -351,8 +289,8 @@ module subroutine collision_util_reset_impactors(self) self%rb(:,:) = 0.0_DP self%vb(:,:) = 0.0_DP self%rot(:,:) = 0.0_DP - self%Lspin(:,:) = 0.0_DP - self%Lorbit(:,:) = 0.0_DP + self%L_spin(:,:) = 0.0_DP + self%L_orbit(:,:) = 0.0_DP self%Ip(:,:) = 0.0_DP self%mass(:) = 0.0_DP self%radius(:) = 0.0_DP @@ -420,13 +358,13 @@ module subroutine collision_util_reset_system(self) if (allocated(after%tp)) deallocate(after%tp) end select - self%Lorbit(:,:) = 0.0_DP - self%Lspin(:,:) = 0.0_DP - self%Ltot(:,:) = 0.0_DP + self%L_orbit(:,:) = 0.0_DP + self%L_spin(:,:) = 0.0_DP + self%L_total(:,:) = 0.0_DP self%ke_orbit(:) = 0.0_DP self%ke_spin(:) = 0.0_DP self%pe(:) = 0.0_DP - self%Etot(:) = 0.0_DP + self%te(:) = 0.0_DP if (allocated(self%impactors)) call self%impactors%reset() if (allocated(self%fragments)) deallocate(self%fragments) @@ -445,8 +383,8 @@ module subroutine collision_util_set_budgets(self) associate(impactors => self%impactors, fragments => self%fragments) - fragments%L_budget(:) = self%Ltot(:,1) - fragments%E_budget = self%Etot(1) - impactors%Qloss + fragments%L_budget(:) = self%L_total(:,1) + fragments%E_budget = self%te(1) - impactors%Qloss end associate @@ -509,7 +447,7 @@ module subroutine collision_util_set_coordinate_impactors(self) ! Arguments class(collision_impactors), intent(inout) :: self !! Collisional nbody_system ! Internals - real(DP), dimension(NDIM) :: delta_r, delta_v, Ltot + real(DP), dimension(NDIM) :: delta_r, delta_v, L_total real(DP) :: L_mag, mtot associate(impactors => self) @@ -521,11 +459,11 @@ module subroutine collision_util_set_coordinate_impactors(self) ! y-axis is the separation distance impactors%y_unit(:) = .unit.delta_r(:) - Ltot = impactors%Lorbit(:,1) + impactors%Lorbit(:,2) + impactors%Lspin(:,1) + impactors%Lspin(:,2) + L_total = impactors%L_orbit(:,1) + impactors%L_orbit(:,2) + impactors%L_spin(:,1) + impactors%L_spin(:,2) - L_mag = .mag.Ltot(:) + L_mag = .mag.L_total(:) if (L_mag > sqrt(tiny(L_mag))) then - impactors%z_unit(:) = .unit.Ltot(:) + impactors%z_unit(:) = .unit.L_total(:) else ! Not enough angular momentum to determine a z-axis direction. We'll just pick a random direction call random_number(impactors%z_unit(:)) impactors%z_unit(:) = .unit.impactors%z_unit(:) @@ -628,8 +566,8 @@ module subroutine collision_util_setup_fragments_collider(self, nfrag) self%fragments%density(:) = 0.0_DP self%fragments%rmag(:) = 0.0_DP self%fragments%vmag(:) = 0.0_DP - self%fragments%Lorbit_tot(:) = 0.0_DP - self%fragments%Lspin_tot(:) = 0.0_DP + self%fragments%L_orbit_tot(:) = 0.0_DP + self%fragments%L_spin_tot(:) = 0.0_DP self%fragments%L_budget(:) = 0.0_DP self%fragments%ke_orbit_tot = 0.0_DP self%fragments%ke_spin_tot = 0.0_DP @@ -727,8 +665,7 @@ module subroutine collision_util_snapshot(self, param, nbody_system, t, arg) real(DP), intent(in), optional :: t !! Time of snapshot if different from nbody_system time character(*), intent(in), optional :: arg !! "before": takes a snapshot just before the collision. "after" takes the snapshot just after the collision. ! Arguments - class(collision_snapshot), allocatable :: snapshot - class(swiftest_pl), allocatable :: pl + class(collision_snapshot), allocatable, save :: snapshot character(len=:), allocatable :: stage if (present(arg)) then @@ -744,27 +681,35 @@ module subroutine collision_util_snapshot(self, param, nbody_system, t, arg) select case(stage) case("before") ! Saves the states of the bodies involved in the collision before the collision is resolved - associate (idx => nbody_system%collider%impactors%id, ncoll => nbody_system%collider%impactors%ncoll) - allocate(pl, mold=nbody_system%pl) - call pl%setup(ncoll, param) - pl%id(:) = nbody_system%pl%id(idx(:)) - pl%Gmass(:) = nbody_system%pl%Gmass(idx(:)) - pl%radius(:) = nbody_system%pl%radius(idx(:)) - pl%rot(:,:) = nbody_system%pl%rot(:,idx(:)) - pl%Ip(:,:) = nbody_system%pl%Ip(:,idx(:)) - pl%rh(:,:) = nbody_system%pl%rh(:,idx(:)) - pl%vh(:,:) = nbody_system%pl%vh(:,idx(:)) - pl%info(:) = nbody_system%pl%info(idx(:)) - select type (before => nbody_system%collider%before) - class is (swiftest_nbody_system) - call move_alloc(pl, before%pl) - end select - end associate - case("after") allocate(collision_snapshot :: snapshot) allocate(snapshot%collider, source=nbody_system%collider) snapshot%t = t + + ! Get and record the energy of the system before the collision + call nbody_system%get_energy_and_momentum(param) + snapshot%collider%L_orbit(:,1) = nbody_system%L_orbit(:) + snapshot%collider%L_spin(:,1) = nbody_system%L_spin(:) + snapshot%collider%L_total(:,1) = nbody_system%L_total(:) + snapshot%collider%ke_orbit(1) = nbody_system%ke_orbit + snapshot%collider%ke_spin(1) = nbody_system%ke_spin + snapshot%collider%pe(1) = nbody_system%pe + snapshot%collider%be(1) = nbody_system%be + snapshot%collider%te(1) = nbody_system%te + case("after") + ! Get record the energy of the sytem after the collision + call nbody_system%get_energy_and_momentum(param) + snapshot%collider%L_orbit(:,2) = nbody_system%L_orbit(:) + snapshot%collider%L_spin(:,2) = nbody_system%L_spin(:) + snapshot%collider%L_total(:,2) = nbody_system%L_total(:) + snapshot%collider%ke_orbit(2) = nbody_system%ke_orbit + snapshot%collider%ke_spin(2) = nbody_system%ke_spin + snapshot%collider%pe(2) = nbody_system%pe + snapshot%collider%be(2) = nbody_system%be + snapshot%collider%te(2) = nbody_system%te + + ! Save the snapshot for posterity call collision_util_save_snapshot(nbody_system%collision_history,snapshot) + deallocate(snapshot) case default write(*,*) "collision_util_snapshot requies either 'before' or 'after' passed to 'arg'" end select diff --git a/src/fraggle/fraggle_generate.f90 b/src/fraggle/fraggle_generate.f90 index d8dbeb0de..c72254ca9 100644 --- a/src/fraggle/fraggle_generate.f90 +++ b/src/fraggle/fraggle_generate.f90 @@ -136,8 +136,8 @@ module subroutine fraggle_generate_disrupt(self, nbody_system, param, t, lfailur call self%get_energy_and_momentum(nbody_system, param, lbefore=.false.) call self%set_original_scale() - dE = self%Etot(2) - self%Etot(1) - dL = .mag.(self%Ltot(:,2) - self%Ltot(:,1)) + dE = self%te(2) - self%te(1) + dL = .mag.(self%L_total(:,2) - self%L_total(:,1)) write(message,*) dE call swiftest_io_log_one_message(COLLISION_LOG_OUT, "Estimated energy change: " // trim(adjustl(message))) @@ -326,7 +326,7 @@ module subroutine fraggle_generate_rot_vec(collider) ! Arguments class(collision_fraggle), intent(inout) :: collider !! Fraggle collision system object ! Internals - real(DP), dimension(NDIM) :: Lbefore, Lafter, Lspin, rotdir + real(DP), dimension(NDIM) :: Lbefore, Lafter, L_spin, rotdir real(DP) :: v_init, v_final, mass_init, mass_final, rotmag integer(I4B) :: i @@ -341,8 +341,8 @@ module subroutine fraggle_generate_rot_vec(collider) Lbefore(:) = mass_init * (impactors%rb(:,2) - impactors%rb(:,1)) .cross. (impactors%vb(:,2) - impactors%vb(:,1)) Lafter(:) = mass_final * (impactors%rb(:,2) - impactors%rb(:,1)) .cross. (v_final * impactors%bounce_unit(:)) - Lspin(:) = impactors%Lspin(:,1) + (Lbefore(:) - Lafter(:)) - fragments%rot(:,1) = Lspin(:) / (fragments%mass(1) * fragments%radius(1)**2 * fragments%Ip(3,1)) + L_spin(:) = impactors%L_spin(:,1) + (Lbefore(:) - Lafter(:)) + fragments%rot(:,1) = L_spin(:) / (fragments%mass(1) * fragments%radius(1)**2 * fragments%Ip(3,1)) ! Add in some random spin noise. The magnitude will be scaled by the before-after amount and the direction will be random do concurrent(i = 2:nfrag) @@ -488,7 +488,7 @@ module subroutine fraggle_generate_vel_vec(collider, lfailure) if (lsupercat) then ! Put some of the residual angular momentum into velocity shear. Not too much, or we get some weird trajectories call fragments%get_angular_momentum() - Lresidual(:) = fragments%L_budget(:) - (fragments%Lorbit_tot(:) + fragments%Lspin_tot(:)) + Lresidual(:) = fragments%L_budget(:) - (fragments%L_orbit_tot(:) + fragments%L_spin_tot(:)) do concurrent(i = istart:nfrag) vunit(:) = .unit. (Lresidual(:) .cross. fragments%r_unit(:,i)) vshear(:) = vunit(:) * (.mag.Lresidual(:) / ((nfrag-istart+1)*fragments%mass(i) * fragments%rmag(i))) @@ -497,14 +497,14 @@ module subroutine fraggle_generate_vel_vec(collider, lfailure) end if ! Check for any residual angular momentum, and if there is any, put it into spin call fragments%get_angular_momentum() - Lresidual(:) = fragments%L_budget(:) - (fragments%Lorbit_tot(:) + fragments%Lspin_tot(:)) + Lresidual(:) = fragments%L_budget(:) - (fragments%L_orbit_tot(:) + fragments%L_spin_tot(:)) do concurrent(i = 1:nfrag) - fragments%Lspin(:,i) = fragments%Lspin(:,i) + Lresidual(:) / nfrag - fragments%rot(:,i) = fragments%Lspin(:,i) / (fragments%mass(i) * fragments%radius(i)**2 * fragments%Ip(:,i)) + fragments%L_spin(:,i) = fragments%L_spin(:,i) + Lresidual(:) / nfrag + fragments%rot(:,i) = fragments%L_spin(:,i) / (fragments%mass(i) * fragments%radius(i)**2 * fragments%Ip(:,i)) end do call fragments%get_angular_momentum() - Lresidual(:) = fragments%L_budget(:) - (fragments%Lorbit_tot(:) + fragments%Lspin_tot(:)) + Lresidual(:) = fragments%L_budget(:) - (fragments%L_orbit_tot(:) + fragments%L_spin_tot(:)) end do ! We didn't converge. Try another configuration and see if we get a better result diff --git a/src/fraggle/fraggle_util.f90 b/src/fraggle/fraggle_util.f90 index b69a7aea3..caf72fbd3 100644 --- a/src/fraggle/fraggle_util.f90 +++ b/src/fraggle/fraggle_util.f90 @@ -260,11 +260,11 @@ module subroutine fraggle_util_set_natural_scale_factors(self) impactors%Gmass(:) = impactors%Gmass(:) / (collider%dscale**3/collider%tscale**2) impactors%Mcb = impactors%Mcb / collider%mscale impactors%radius(:) = impactors%radius(:) / collider%dscale - impactors%Lspin(:,:) = impactors%Lspin(:,:) / collider%Lscale - impactors%Lorbit(:,:) = impactors%Lorbit(:,:) / collider%Lscale + impactors%L_spin(:,:) = impactors%L_spin(:,:) / collider%Lscale + impactors%L_orbit(:,:) = impactors%L_orbit(:,:) / collider%Lscale do concurrent(i = 1:2) - impactors%rot(:,i) = impactors%Lspin(:,i) / (impactors%mass(i) * impactors%radius(i)**2 * impactors%Ip(3,i)) + impactors%rot(:,i) = impactors%L_spin(:,i) / (impactors%mass(i) * impactors%radius(i)**2 * impactors%Ip(3,i)) end do fragments%mtot = fragments%mtot / collider%mscale @@ -309,10 +309,10 @@ module subroutine fraggle_util_set_original_scale_factors(self) impactors%vb = impactors%vb * collider%vscale impactors%rc = impactors%rc * collider%dscale impactors%vc = impactors%vc * collider%vscale - impactors%Lspin = impactors%Lspin * collider%Lscale - impactors%Lorbit = impactors%Lorbit * collider%Lscale + impactors%L_spin = impactors%L_spin * collider%Lscale + impactors%L_orbit = impactors%L_orbit * collider%Lscale do concurrent(i = 1:2) - impactors%rot(:,i) = impactors%Lspin(:,i) * (impactors%mass(i) * impactors%radius(i)**2 * impactors%Ip(3,i)) + impactors%rot(:,i) = impactors%L_spin(:,i) * (impactors%mass(i) * impactors%radius(i)**2 * impactors%Ip(3,i)) end do fragments%mtot = fragments%mtot * collider%mscale @@ -327,14 +327,14 @@ module subroutine fraggle_util_set_original_scale_factors(self) impactors%Qloss = impactors%Qloss * collider%Escale - collider%Lorbit(:,:) = collider%Lorbit(:,:) * collider%Lscale - collider%Lspin(:,:) = collider%Lspin(:,:) * collider%Lscale - collider%Ltot(:,:) = collider%Ltot(:,:) * collider%Lscale + collider%L_orbit(:,:) = collider%L_orbit(:,:) * collider%Lscale + collider%L_spin(:,:) = collider%L_spin(:,:) * collider%Lscale + collider%L_total(:,:) = collider%L_total(:,:) * collider%Lscale collider%ke_orbit(:) = collider%ke_orbit(:) * collider%Escale collider%ke_spin(:) = collider%ke_spin(:) * collider%Escale collider%pe(:) = collider%pe(:) * collider%Escale collider%be(:) = collider%be(:) * collider%Escale - collider%Etot(:) = collider%Etot(:) * collider%Escale + collider%te(:) = collider%te(:) * collider%Escale collider%mscale = 1.0_DP collider%dscale = 1.0_DP diff --git a/src/netcdf_io/netcdf_io_module.f90 b/src/netcdf_io/netcdf_io_module.f90 index 3926d2eea..f49ad81cb 100644 --- a/src/netcdf_io/netcdf_io_module.f90 +++ b/src/netcdf_io/netcdf_io_module.f90 @@ -78,7 +78,7 @@ module netcdf_io integer(I4B) :: vh_varid !! ID for the heliocentric velocity vector variable character(NAMELEN) :: gr_pseudo_vh_varname = "gr_pseudo_vh" !! name of the heliocentric pseudovelocity vector variable (used in GR only) integer(I4B) :: gr_pseudo_vh_varid !! ID for the heliocentric pseudovelocity vector variable (used in GR) - character(NAMELEN) :: gmass_varname = "Gmass" !! name of the mass variable + character(NAMELEN) :: Gmass_varname = "Gmass" !! name of the mass variable integer(I4B) :: Gmass_varid !! ID for the mass variable character(NAMELEN) :: rhill_varname = "rhill" !! name of the hill radius variable integer(I4B) :: rhill_varid !! ID for the hill radius variable @@ -104,16 +104,16 @@ module netcdf_io integer(I4B) :: PE_varid !! ID for the system potential energy variable character(NAMELEN) :: be_varname = "BE" !! name of the system binding energy variable integer(I4B) :: BE_varid !! ID for the system binding energy variable - character(NAMELEN) :: L_orb_varname = "L_orb" !! name of the orbital angular momentum vector variable - integer(I4B) :: L_orb_varid !! ID for the system orbital angular momentum vector variable - character(NAMELEN) :: Lspin_varname = "Lspin" !! name of the spin angular momentum vector variable - integer(I4B) :: Lspin_varid !! ID for the system spin angular momentum vector variable + character(NAMELEN) :: L_orbit_varname = "L_orbit" !! name of the orbital angular momentum vector variable + integer(I4B) :: L_orbit_varid !! ID for the system orbital angular momentum vector variable + character(NAMELEN) :: L_spin_varname = "L_spin" !! name of the spin angular momentum vector variable + integer(I4B) :: L_spin_varid !! ID for the system spin angular momentum vector variable character(NAMELEN) :: L_escape_varname = "L_escape" !! name of the escaped angular momentum vector variable integer(I4B) :: L_escape_varid !! ID for the escaped angular momentum vector variable - character(NAMELEN) :: Ecollisions_varname = "Ecollisions" !! name of the escaped angular momentum y variable - integer(I4B) :: Ecollisions_varid !! ID for the energy lost in collisions variable - character(NAMELEN) :: Euntracked_varname = "Euntracked" !! name of the energy that is untracked due to loss (untracked potential energy due to mergers and body energy for escaped bodies) - integer(I4B) :: Euntracked_varid !! ID for the energy that is untracked due to loss (untracked potential energy due to mergers and body energy for escaped bodies) + character(NAMELEN) :: E_collisions_varname = "E_collisions" !! name of the escaped angular momentum y variable + integer(I4B) :: E_collisions_varid !! ID for the energy lost in collisions variable + character(NAMELEN) :: E_untracked_varname = "E_untracked" !! name of the energy that is untracked due to loss (untracked potential energy due to mergers and body energy for escaped bodies) + integer(I4B) :: E_untracked_varid !! ID for the energy that is untracked due to loss (untracked potential energy due to mergers and body energy for escaped bodies) character(NAMELEN) :: GMescape_varname = "GMescape" !! name of the G*Mass of bodies that escape the system integer(I4B) :: GMescape_varid !! ID for the G*Mass of bodies that escape the system character(NAMELEN) :: origin_type_varname = "origin_type" !! name of the origin type variable (Initial Conditions, Disruption, etc.) diff --git a/src/swiftest/swiftest_io.f90 b/src/swiftest/swiftest_io.f90 index ba53fde8f..063a637ad 100644 --- a/src/swiftest/swiftest_io.f90 +++ b/src/swiftest/swiftest_io.f90 @@ -41,11 +41,11 @@ module subroutine swiftest_io_compact_output(self, param, timer) formatted_output = formatted_output // fmt("NPLM",pl%nplm) end select if (param%lenergy) then - formatted_output = formatted_output // fmt("LTOTERR",self%Ltot_error) // fmt("ETOTERR",self%Etot_error) // fmt("MTOTERR",self%Mtot_error) & - // fmt("KEOERR",self%ke_orbit_error) // fmt("PEERR",self%pe_error) // fmt("EORBERR",self%Eorbit_error) & - // fmt("EUNTRERR",self%Euntracked_error) // fmt("LESCERR",self%Lescape_error) // fmt("MESCERR",self%Mescape_error) + formatted_output = formatted_output // fmt("LTOTERR",self%L_total_error) // fmt("ETOTERR",self%te_error) // fmt("MTOTERR",self%Mtot_error) & + // fmt("KEOERR",self%ke_orbit_error) // fmt("PEERR",self%pe_error) // fmt("EORBERR",self%E_orbit_error) & + // fmt("EUNTRERR",self%E_untracked_error) // fmt("LESCERR",self%L_escape_error) // fmt("MESCERR",self%Mescape_error) if (param%lclose) formatted_output = formatted_output // fmt("ECOLLERR",self%Ecoll_error) - if (param%lrotation) formatted_output = formatted_output // fmt("KESPINERR",self%ke_spin_error) // fmt("LSPINERR",self%Lspin_error) + if (param%lrotation) formatted_output = formatted_output // fmt("KESPINERR",self%ke_spin_error) // fmt("LSPINERR",self%L_spin_error) end if if (.not. timer%main_is_started) then ! This is the start of a new run @@ -120,14 +120,14 @@ module subroutine swiftest_io_conservation_report(self, param, lterminal) class(swiftest_parameters), intent(inout) :: param !! Input colleciton of user-defined parameters logical, intent(in) :: lterminal !! Indicates whether to output information to the terminal screen ! Internals - real(DP), dimension(NDIM) :: Ltot_now, Lorbit_now, Lspin_now - real(DP) :: ke_orbit_now, ke_spin_now, pe_now, Eorbit_now, be_now + real(DP), dimension(NDIM) :: L_total_now, L_orbit_now, L_spin_now + real(DP) :: ke_orbit_now, ke_spin_now, pe_now, E_orbit_now, be_now real(DP) :: GMtot_now character(len=STRMAX) :: errmsg integer(I4B), parameter :: EGYIU = 72 character(len=*), parameter :: EGYTERMFMT = '(" DL/L0 = ", ES12.5 & - "; DEcollisions/|E0| = ", ES12.5, & - "; D(Eorbit+Ecollisions)/|E0| = ", ES12.5, & + "; DE_collisions/|E0| = ", ES12.5, & + "; D(E_orbit+E_collisions)/|E0| = ", ES12.5, & "; DM/M0 = ", ES12.5)' associate(nbody_system => self, pl => self%pl, cb => self%cb, npl => self%pl%nbody, display_unit => param%display_unit, nc => param%system_history%nc) @@ -140,10 +140,10 @@ module subroutine swiftest_io_conservation_report(self, param, lterminal) ke_spin_now = nbody_system%ke_spin pe_now = nbody_system%pe be_now = nbody_system%be - Lorbit_now(:) = nbody_system%Lorbit(:) - Lspin_now(:) = nbody_system%Lspin(:) - Eorbit_now = ke_orbit_now + ke_spin_now + pe_now + be_now - Ltot_now(:) = nbody_system%Ltot(:) + nbody_system%Lescape(:) + L_orbit_now(:) = nbody_system%L_orbit(:) + L_spin_now(:) = nbody_system%L_spin(:) + E_orbit_now = ke_orbit_now + ke_spin_now + pe_now + be_now + L_total_now(:) = nbody_system%L_total(:) + nbody_system%L_escape(:) GMtot_now = nbody_system%GMtot + nbody_system%GMescape if (param%lfirstenergy) then @@ -151,31 +151,31 @@ module subroutine swiftest_io_conservation_report(self, param, lterminal) nbody_system%ke_spin_orig = ke_spin_now nbody_system%pe_orig = pe_now nbody_system%be_orig = be_now - nbody_system%Eorbit_orig = Eorbit_now + nbody_system%E_orbit_orig = E_orbit_now nbody_system%GMtot_orig = GMtot_now - nbody_system%Lorbit_orig(:) = Lorbit_now(:) - nbody_system%Lspin_orig(:) = Lspin_now(:) - nbody_system%Ltot_orig(:) = Ltot_now(:) + nbody_system%L_orbit_orig(:) = L_orbit_now(:) + nbody_system%L_spin_orig(:) = L_spin_now(:) + nbody_system%L_total_orig(:) = L_total_now(:) param%lfirstenergy = .false. end if if (.not.param%lfirstenergy) then - nbody_system%ke_orbit_error = (ke_orbit_now - nbody_system%ke_orbit_orig) / abs(nbody_system%Eorbit_orig) - nbody_system%ke_spin_error = (ke_spin_now - nbody_system%ke_spin_orig) / abs(nbody_system%Eorbit_orig) - nbody_system%pe_error = (pe_now - nbody_system%pe_orig) / abs(nbody_system%Eorbit_orig) - nbody_system%be_error = (be_now - nbody_system%be_orig) / abs(nbody_system%Eorbit_orig) - nbody_system%Eorbit_error = (Eorbit_now - nbody_system%Eorbit_orig) / abs(nbody_system%Eorbit_orig) - nbody_system%Ecoll_error = nbody_system%Ecollisions / abs(nbody_system%Eorbit_orig) - nbody_system%Euntracked_error = nbody_system%Euntracked / abs(nbody_system%Eorbit_orig) - nbody_system%Etot_error = (Eorbit_now - nbody_system%Ecollisions - nbody_system%Eorbit_orig - nbody_system%Euntracked) / abs(nbody_system%Eorbit_orig) - - nbody_system%Lorbit_error = norm2(Lorbit_now(:) - nbody_system%Lorbit_orig(:)) / norm2(nbody_system%Ltot_orig(:)) - nbody_system%Lspin_error = norm2(Lspin_now(:) - nbody_system%Lspin_orig(:)) / norm2(nbody_system%Ltot_orig(:)) - nbody_system%Lescape_error = norm2(nbody_system%Lescape(:)) / norm2(nbody_system%Ltot_orig(:)) - nbody_system%Ltot_error = norm2(Ltot_now(:) - nbody_system%Ltot_orig(:)) / norm2(nbody_system%Ltot_orig(:)) + nbody_system%ke_orbit_error = (ke_orbit_now - nbody_system%ke_orbit_orig) / abs(nbody_system%E_orbit_orig) + nbody_system%ke_spin_error = (ke_spin_now - nbody_system%ke_spin_orig) / abs(nbody_system%E_orbit_orig) + nbody_system%pe_error = (pe_now - nbody_system%pe_orig) / abs(nbody_system%E_orbit_orig) + nbody_system%be_error = (be_now - nbody_system%be_orig) / abs(nbody_system%E_orbit_orig) + nbody_system%E_orbit_error = (E_orbit_now - nbody_system%E_orbit_orig) / abs(nbody_system%E_orbit_orig) + nbody_system%Ecoll_error = nbody_system%E_collisions / abs(nbody_system%E_orbit_orig) + nbody_system%E_untracked_error = nbody_system%E_untracked / abs(nbody_system%E_orbit_orig) + nbody_system%te_error = (E_orbit_now - nbody_system%E_collisions - nbody_system%E_orbit_orig - nbody_system%E_untracked) / abs(nbody_system%E_orbit_orig) + + nbody_system%L_orbit_error = norm2(L_orbit_now(:) - nbody_system%L_orbit_orig(:)) / norm2(nbody_system%L_total_orig(:)) + nbody_system%L_spin_error = norm2(L_spin_now(:) - nbody_system%L_spin_orig(:)) / norm2(nbody_system%L_total_orig(:)) + nbody_system%L_escape_error = norm2(nbody_system%L_escape(:)) / norm2(nbody_system%L_total_orig(:)) + nbody_system%L_total_error = norm2(L_total_now(:) - nbody_system%L_total_orig(:)) / norm2(nbody_system%L_total_orig(:)) nbody_system%Mescape_error = nbody_system%GMescape / nbody_system%GMtot_orig nbody_system%Mtot_error = (GMtot_now - nbody_system%GMtot_orig) / nbody_system%GMtot_orig - if (lterminal) write(display_unit, EGYTERMFMT) nbody_system%Ltot_error, nbody_system%Ecoll_error, nbody_system%Etot_error,nbody_system%Mtot_error + if (lterminal) write(display_unit, EGYTERMFMT) nbody_system%L_total_error, nbody_system%Ecoll_error, nbody_system%te_error,nbody_system%Mtot_error if (abs(nbody_system%Mtot_error) > 100 * epsilon(nbody_system%Mtot_error)) then write(*,*) "Severe error! Mass not conserved! Halting!" ! Save the frame of data to the bin file in the slot just after the present one for diagnostics @@ -523,16 +523,16 @@ module subroutine swiftest_io_netcdf_get_t0_values_system(self, param) call netcdf_io_check( nf90_get_var(nc%id, nc%BE_varid, rtemp, start=[1], count=[1]), "netcdf_io_get_t0_values_system BE_varid" ) BE_orig = rtemp(1) - call netcdf_io_check( nf90_get_var(nc%id, nc%Ecollisions_varid, self%Ecollisions, start=[1]), "netcdf_io_get_t0_values_system Ecollisions_varid" ) - call netcdf_io_check( nf90_get_var(nc%id, nc%Euntracked_varid, self%Euntracked, start=[1]), "netcdf_io_get_t0_values_system Euntracked_varid" ) + call netcdf_io_check( nf90_get_var(nc%id, nc%E_collisions_varid, self%E_collisions, start=[1]), "netcdf_io_get_t0_values_system E_collisions_varid" ) + call netcdf_io_check( nf90_get_var(nc%id, nc%E_untracked_varid, self%E_untracked, start=[1]), "netcdf_io_get_t0_values_system E_untracked_varid" ) - self%Eorbit_orig = KE_orb_orig + KE_spin_orig + PE_orig + BE_orig + self%Ecollisions + self%Euntracked + self%E_orbit_orig = KE_orb_orig + KE_spin_orig + PE_orig + BE_orig + self%E_collisions + self%E_untracked - call netcdf_io_check( nf90_get_var(nc%id, nc%L_orb_varid, self%Lorbit_orig(:), start=[1,1], count=[NDIM,1]), "netcdf_io_get_t0_values_system L_orb_varid" ) - call netcdf_io_check( nf90_get_var(nc%id, nc%Lspin_varid, self%Lspin_orig(:), start=[1,1], count=[NDIM,1]), "netcdf_io_get_t0_values_system Lspin_varid" ) - call netcdf_io_check( nf90_get_var(nc%id, nc%L_escape_varid, self%Lescape(:), start=[1,1], count=[NDIM,1]), "netcdf_io_get_t0_values_system L_escape_varid" ) + call netcdf_io_check( nf90_get_var(nc%id, nc%L_orbit_varid, self%L_orbit_orig(:), start=[1,1], count=[NDIM,1]), "netcdf_io_get_t0_values_system L_orbit_varid" ) + call netcdf_io_check( nf90_get_var(nc%id, nc%L_spin_varid, self%L_spin_orig(:), start=[1,1], count=[NDIM,1]), "netcdf_io_get_t0_values_system L_spin_varid" ) + call netcdf_io_check( nf90_get_var(nc%id, nc%L_escape_varid, self%L_escape(:), start=[1,1], count=[NDIM,1]), "netcdf_io_get_t0_values_system L_escape_varid" ) - self%Ltot_orig(:) = self%Lorbit_orig(:) + self%Lspin_orig(:) + self%Lescape(:) + self%L_total_orig(:) = self%L_orbit_orig(:) + self%L_spin_orig(:) + self%L_escape(:) call netcdf_io_check( nf90_get_var(nc%id, nc%Gmass_varid, vals, start=[1,1], count=[idmax,1]), "netcdf_io_get_t0_values_system Gmass_varid" ) call netcdf_io_check( nf90_get_var(nc%id, nc%GMescape_varid, self%GMescape, start=[1]), "netcdf_io_get_t0_values_system GMescape_varid" ) @@ -652,7 +652,7 @@ module subroutine swiftest_io_netcdf_initialize_output(self, param) call netcdf_io_check( nf90_def_var(nc%id, nc%cape_varname, nc%out_type, [nc%name_dimid, nc%time_dimid], nc%cape_varid), "netcdf_io_initialize_output nf90_def_var cape_varid" ) end if - call netcdf_io_check( nf90_def_var(nc%id, nc%gmass_varname, nc%out_type, [nc%name_dimid, nc%time_dimid], nc%Gmass_varid), "netcdf_io_initialize_output nf90_def_var Gmass_varid" ) + call netcdf_io_check( nf90_def_var(nc%id, nc%Gmass_varname, nc%out_type, [nc%name_dimid, nc%time_dimid], nc%Gmass_varid), "netcdf_io_initialize_output nf90_def_var Gmass_varid" ) if (param%lrhill_present) then call netcdf_io_check( nf90_def_var(nc%id, nc%rhill_varname, nc%out_type, [nc%name_dimid, nc%time_dimid], nc%rhill_varid), "netcdf_io_initialize_output nf90_def_var rhill_varid" ) @@ -689,11 +689,11 @@ module subroutine swiftest_io_netcdf_initialize_output(self, param) call netcdf_io_check( nf90_def_var(nc%id, nc%ke_spin_varname, nc%out_type, nc%time_dimid, nc%KE_spin_varid), "netcdf_io_initialize_output nf90_def_var KE_spin_varid" ) call netcdf_io_check( nf90_def_var(nc%id, nc%pe_varname, nc%out_type, nc%time_dimid, nc%PE_varid), "netcdf_io_initialize_output nf90_def_var PE_varid" ) call netcdf_io_check( nf90_def_var(nc%id, nc%be_varname, nc%out_type, nc%time_dimid, nc%BE_varid), "netcdf_io_initialize_output nf90_def_var PE_varid" ) - call netcdf_io_check( nf90_def_var(nc%id, nc%L_orb_varname, nc%out_type, [nc%space_dimid, nc%time_dimid], nc%L_orb_varid), "netcdf_io_initialize_output nf90_def_var L_orb_varid" ) - call netcdf_io_check( nf90_def_var(nc%id, nc%Lspin_varname, nc%out_type, [nc%space_dimid, nc%time_dimid], nc%Lspin_varid), "netcdf_io_initialize_output nf90_def_var Lspin_varid" ) + call netcdf_io_check( nf90_def_var(nc%id, nc%L_orbit_varname, nc%out_type, [nc%space_dimid, nc%time_dimid], nc%L_orbit_varid), "netcdf_io_initialize_output nf90_def_var L_orbit_varid" ) + call netcdf_io_check( nf90_def_var(nc%id, nc%L_spin_varname, nc%out_type, [nc%space_dimid, nc%time_dimid], nc%L_spin_varid), "netcdf_io_initialize_output nf90_def_var L_spin_varid" ) call netcdf_io_check( nf90_def_var(nc%id, nc%L_escape_varname, nc%out_type, [nc%space_dimid, nc%time_dimid], nc%L_escape_varid), "netcdf_io_initialize_output nf90_def_var L_escape_varid" ) - call netcdf_io_check( nf90_def_var(nc%id, nc%Ecollisions_varname, nc%out_type, nc%time_dimid, nc%Ecollisions_varid), "netcdf_io_initialize_output nf90_def_var Ecollisions_varid" ) - call netcdf_io_check( nf90_def_var(nc%id, nc%Euntracked_varname, nc%out_type, nc%time_dimid, nc%Euntracked_varid), "netcdf_io_initialize_output nf90_def_var Euntracked_varid" ) + call netcdf_io_check( nf90_def_var(nc%id, nc%E_collisions_varname, nc%out_type, nc%time_dimid, nc%E_collisions_varid), "netcdf_io_initialize_output nf90_def_var E_collisions_varid" ) + call netcdf_io_check( nf90_def_var(nc%id, nc%E_untracked_varname, nc%out_type, nc%time_dimid, nc%E_untracked_varid), "netcdf_io_initialize_output nf90_def_var E_untracked_varid" ) call netcdf_io_check( nf90_def_var(nc%id, nc%GMescape_varname, nc%out_type, nc%time_dimid, nc%GMescape_varid), "netcdf_io_initialize_output nf90_def_var GMescape_varid" ) end if @@ -760,28 +760,28 @@ module subroutine swiftest_io_netcdf_open(self, param, readonly) associate(nc => self) - write(errmsg,*) "netcdf_io_open nf90_open ",trim(adjustl(nc%file_name)) + write(errmsg,*) "swiftest_io_netcdf_opennf90_open ",trim(adjustl(nc%file_name)) call netcdf_io_check( nf90_open(nc%file_name, mode, nc%id), errmsg) self%lfile_is_open = .true. ! Dimensions - call netcdf_io_check( nf90_inq_dimid(nc%id, nc%time_dimname, nc%time_dimid), "netcdf_io_open nf90_inq_dimid time_dimid" ) - call netcdf_io_check( nf90_inq_dimid(nc%id, nc%space_dimname, nc%space_dimid), "netcdf_io_open nf90_inq_dimid space_dimid" ) - call netcdf_io_check( nf90_inq_dimid(nc%id, nc%name_dimname, nc%name_dimid), "netcdf_io_open nf90_inq_dimid name_dimid" ) - call netcdf_io_check( nf90_inq_dimid(nc%id, nc%str_dimname, nc%str_dimid), "netcdf_io_open nf90_inq_dimid str_dimid" ) + call netcdf_io_check( nf90_inq_dimid(nc%id, nc%time_dimname, nc%time_dimid), "swiftest_io_netcdf_opennf90_inq_dimid time_dimid" ) + call netcdf_io_check( nf90_inq_dimid(nc%id, nc%space_dimname, nc%space_dimid), "swiftest_io_netcdf_opennf90_inq_dimid space_dimid" ) + call netcdf_io_check( nf90_inq_dimid(nc%id, nc%name_dimname, nc%name_dimid), "swiftest_io_netcdf_opennf90_inq_dimid name_dimid" ) + call netcdf_io_check( nf90_inq_dimid(nc%id, nc%str_dimname, nc%str_dimid), "swiftest_io_netcdf_opennf90_inq_dimid str_dimid" ) ! Dimension coordinates - call netcdf_io_check( nf90_inq_varid(nc%id, nc%time_dimname, nc%time_varid), "netcdf_io_open nf90_inq_varid time_varid" ) - call netcdf_io_check( nf90_inq_varid(nc%id, nc%space_dimname, nc%space_varid), "netcdf_io_open nf90_inq_varid space_varid" ) - call netcdf_io_check( nf90_inq_varid(nc%id, nc%name_dimname, nc%name_varid), "netcdf_io_open nf90_inq_varid name_varid" ) + call netcdf_io_check( nf90_inq_varid(nc%id, nc%time_dimname, nc%time_varid), "swiftest_io_netcdf_opennf90_inq_varid time_varid" ) + call netcdf_io_check( nf90_inq_varid(nc%id, nc%space_dimname, nc%space_varid), "swiftest_io_netcdf_opennf90_inq_varid space_varid" ) + call netcdf_io_check( nf90_inq_varid(nc%id, nc%name_dimname, nc%name_varid), "swiftest_io_netcdf_opennf90_inq_varid name_varid" ) ! Required Variables - call netcdf_io_check( nf90_inq_varid(nc%id, nc%id_varname, nc%id_varid), "netcdf_io_open nf90_inq_varid name_varid" ) - call netcdf_io_check( nf90_inq_varid(nc%id, nc%gmass_varname, nc%Gmass_varid), "netcdf_io_open nf90_inq_varid Gmass_varid" ) + call netcdf_io_check( nf90_inq_varid(nc%id, nc%id_varname, nc%id_varid), "swiftest_io_netcdf_opennf90_inq_varid name_varid" ) + call netcdf_io_check( nf90_inq_varid(nc%id, nc%Gmass_varname, nc%Gmass_varid), "swiftest_io_netcdf_opennf90_inq_varid Gmass_varid" ) if ((param%out_form == "XV") .or. (param%out_form == "XVEL")) then - call netcdf_io_check( nf90_inq_varid(nc%id, nc%rh_varname, nc%rh_varid), "netcdf_io_open nf90_inq_varid rh_varid" ) - call netcdf_io_check( nf90_inq_varid(nc%id, nc%vh_varname, nc%vh_varid), "netcdf_io_open nf90_inq_varid vh_varid" ) + call netcdf_io_check( nf90_inq_varid(nc%id, nc%rh_varname, nc%rh_varid), "swiftest_io_netcdf_opennf90_inq_varid rh_varid" ) + call netcdf_io_check( nf90_inq_varid(nc%id, nc%vh_varname, nc%vh_varid), "swiftest_io_netcdf_opennf90_inq_varid vh_varid" ) if (param%lgr) then !! check if pseudovelocity vectors exist in this file. If they are, set the correct flag so we know whe should not do the conversion. @@ -795,26 +795,26 @@ module subroutine swiftest_io_netcdf_open(self, param, readonly) end if if ((param%out_form == "EL") .or. (param%out_form == "XVEL")) then - call netcdf_io_check( nf90_inq_varid(nc%id, nc%a_varname, nc%a_varid), "netcdf_io_open nf90_inq_varid a_varid" ) - call netcdf_io_check( nf90_inq_varid(nc%id, nc%e_varname, nc%e_varid), "netcdf_io_open nf90_inq_varid e_varid" ) - call netcdf_io_check( nf90_inq_varid(nc%id, nc%inc_varname, nc%inc_varid), "netcdf_io_open nf90_inq_varid inc_varid" ) - call netcdf_io_check( nf90_inq_varid(nc%id, nc%capom_varname, nc%capom_varid), "netcdf_io_open nf90_inq_varid capom_varid" ) - call netcdf_io_check( nf90_inq_varid(nc%id, nc%omega_varname, nc%omega_varid), "netcdf_io_open nf90_inq_varid omega_varid" ) - call netcdf_io_check( nf90_inq_varid(nc%id, nc%capm_varname, nc%capm_varid), "netcdf_io_open nf90_inq_varid capm_varid" ) + call netcdf_io_check( nf90_inq_varid(nc%id, nc%a_varname, nc%a_varid), "swiftest_io_netcdf_opennf90_inq_varid a_varid" ) + call netcdf_io_check( nf90_inq_varid(nc%id, nc%e_varname, nc%e_varid), "swiftest_io_netcdf_opennf90_inq_varid e_varid" ) + call netcdf_io_check( nf90_inq_varid(nc%id, nc%inc_varname, nc%inc_varid), "swiftest_io_netcdf_opennf90_inq_varid inc_varid" ) + call netcdf_io_check( nf90_inq_varid(nc%id, nc%capom_varname, nc%capom_varid), "swiftest_io_netcdf_opennf90_inq_varid capom_varid" ) + call netcdf_io_check( nf90_inq_varid(nc%id, nc%omega_varname, nc%omega_varid), "swiftest_io_netcdf_opennf90_inq_varid omega_varid" ) + call netcdf_io_check( nf90_inq_varid(nc%id, nc%capm_varname, nc%capm_varid), "swiftest_io_netcdf_opennf90_inq_varid capm_varid" ) end if if (param%lclose) then - call netcdf_io_check( nf90_inq_varid(nc%id, nc%radius_varname, nc%radius_varid), "netcdf_io_open nf90_inq_varid radius_varid" ) + call netcdf_io_check( nf90_inq_varid(nc%id, nc%radius_varname, nc%radius_varid), "swiftest_io_netcdf_opennf90_inq_varid radius_varid" ) end if if (param%lrotation) then - call netcdf_io_check( nf90_inq_varid(nc%id, nc%Ip_varname, nc%Ip_varid), "netcdf_io_open nf90_inq_varid Ip_varid" ) - call netcdf_io_check( nf90_inq_varid(nc%id, nc%rot_varname, nc%rot_varid), "netcdf_io_open nf90_inq_varid rot_varid" ) + call netcdf_io_check( nf90_inq_varid(nc%id, nc%Ip_varname, nc%Ip_varid), "swiftest_io_netcdf_opennf90_inq_varid Ip_varid" ) + call netcdf_io_check( nf90_inq_varid(nc%id, nc%rot_varname, nc%rot_varid), "swiftest_io_netcdf_opennf90_inq_varid rot_varid" ) end if ! if (param%ltides) then - ! call netcdf_io_check( nf90_inq_varid(nc%id, nc%k2_varname, nc%k2_varid), "netcdf_io_open nf90_inq_varid k2_varid" ) - ! call netcdf_io_check( nf90_inq_varid(nc%id, nc%q_varname, nc%Q_varid), "netcdf_io_open nf90_inq_varid Q_varid" ) + ! call netcdf_io_check( nf90_inq_varid(nc%id, nc%k2_varname, nc%k2_varid), "swiftest_io_netcdf_opennf90_inq_varid k2_varid" ) + ! call netcdf_io_check( nf90_inq_varid(nc%id, nc%q_varname, nc%Q_varid), "swiftest_io_netcdf_opennf90_inq_varid Q_varid" ) ! end if ! Optional Variables @@ -853,11 +853,11 @@ module subroutine swiftest_io_netcdf_open(self, param, readonly) status = nf90_inq_varid(nc%id, nc%ke_spin_varname, nc%KE_spin_varid) status = nf90_inq_varid(nc%id, nc%pe_varname, nc%PE_varid) status = nf90_inq_varid(nc%id, nc%be_varname, nc%BE_varid) - status = nf90_inq_varid(nc%id, nc%L_orb_varname, nc%L_orb_varid) - status = nf90_inq_varid(nc%id, nc%Lspin_varname, nc%Lspin_varid) + status = nf90_inq_varid(nc%id, nc%L_orbit_varname, nc%L_orbit_varid) + status = nf90_inq_varid(nc%id, nc%L_spin_varname, nc%L_spin_varid) status = nf90_inq_varid(nc%id, nc%L_escape_varname, nc%L_escape_varid) - status = nf90_inq_varid(nc%id, nc%Ecollisions_varname, nc%Ecollisions_varid) - status = nf90_inq_varid(nc%id, nc%Euntracked_varname, nc%Euntracked_varid) + status = nf90_inq_varid(nc%id, nc%E_collisions_varname, nc%E_collisions_varid) + status = nf90_inq_varid(nc%id, nc%E_untracked_varname, nc%E_untracked_varid) status = nf90_inq_varid(nc%id, nc%GMescape_varname, nc%GMescape_varid) end if @@ -1176,16 +1176,16 @@ module subroutine swiftest_io_netcdf_read_hdr_system(self, nc, param) if (status == nf90_noerr) call netcdf_io_check( nf90_get_var(nc%id, nc%PE_varid, self%pe, start=[tslot]), "netcdf_io_read_hdr_system nf90_getvar PE_varid" ) status = nf90_inq_varid(nc%id, nc%be_varname, nc%BE_varid) if (status == nf90_noerr) call netcdf_io_check( nf90_get_var(nc%id, nc%BE_varid, self%be, start=[tslot]), "netcdf_io_read_hdr_system nf90_getvar BE_varid" ) - status = nf90_inq_varid(nc%id, nc%L_orb_varname, nc%L_orb_varid) - if (status == nf90_noerr) call netcdf_io_check( nf90_get_var(nc%id, nc%L_orb_varid, self%Lorbit(:), start=[1,tslot], count=[NDIM,1]), "netcdf_io_read_hdr_system nf90_getvar L_orb_varid" ) - status = nf90_inq_varid(nc%id, nc%Lspin_varname, nc%Lspin_varid) - if (status == nf90_noerr) call netcdf_io_check( nf90_get_var(nc%id, nc%Lspin_varid, self%Lspin(:), start=[1,tslot], count=[NDIM,1]), "netcdf_io_read_hdr_system nf90_getvar Lspin_varid" ) + status = nf90_inq_varid(nc%id, nc%L_orbit_varname, nc%L_orbit_varid) + if (status == nf90_noerr) call netcdf_io_check( nf90_get_var(nc%id, nc%L_orbit_varid, self%L_orbit(:), start=[1,tslot], count=[NDIM,1]), "netcdf_io_read_hdr_system nf90_getvar L_orbit_varid" ) + status = nf90_inq_varid(nc%id, nc%L_spin_varname, nc%L_spin_varid) + if (status == nf90_noerr) call netcdf_io_check( nf90_get_var(nc%id, nc%L_spin_varid, self%L_spin(:), start=[1,tslot], count=[NDIM,1]), "netcdf_io_read_hdr_system nf90_getvar L_spin_varid" ) status = nf90_inq_varid(nc%id, nc%L_escape_varname, nc%L_escape_varid) - if (status == nf90_noerr) call netcdf_io_check( nf90_get_var(nc%id, nc%L_escape_varid, self%Lescape(:), start=[1, tslot], count=[NDIM,1]), "netcdf_io_read_hdr_system nf90_getvar L_escape_varid" ) - status = nf90_inq_varid(nc%id, nc%Ecollisions_varname, nc%Ecollisions_varid) - if (status == nf90_noerr) call netcdf_io_check( nf90_get_var(nc%id, nc%Ecollisions_varid, self%Ecollisions, start=[tslot]), "netcdf_io_read_hdr_system nf90_getvar Ecollisions_varid" ) - status = nf90_inq_varid(nc%id, nc%Euntracked_varname, nc%Euntracked_varid) - if (status == nf90_noerr) call netcdf_io_check( nf90_get_var(nc%id, nc%Euntracked_varid, self%Euntracked, start=[tslot]), "netcdf_io_read_hdr_system nf90_getvar Euntracked_varid" ) + if (status == nf90_noerr) call netcdf_io_check( nf90_get_var(nc%id, nc%L_escape_varid, self%L_escape(:), start=[1, tslot], count=[NDIM,1]), "netcdf_io_read_hdr_system nf90_getvar L_escape_varid" ) + status = nf90_inq_varid(nc%id, nc%E_collisions_varname, nc%E_collisions_varid) + if (status == nf90_noerr) call netcdf_io_check( nf90_get_var(nc%id, nc%E_collisions_varid, self%E_collisions, start=[tslot]), "netcdf_io_read_hdr_system nf90_getvar E_collisions_varid" ) + status = nf90_inq_varid(nc%id, nc%E_untracked_varname, nc%E_untracked_varid) + if (status == nf90_noerr) call netcdf_io_check( nf90_get_var(nc%id, nc%E_untracked_varid, self%E_untracked, start=[tslot]), "netcdf_io_read_hdr_system nf90_getvar E_untracked_varid" ) status = nf90_inq_varid(nc%id, nc%GMescape_varname, nc%GMescape_varid) if (status == nf90_noerr) call netcdf_io_check( nf90_get_var(nc%id, nc%GMescape_varid, self%GMescape, start=[tslot]), "netcdf_io_read_hdr_system nf90_getvar GMescape_varid" ) end if @@ -1584,11 +1584,11 @@ module subroutine swiftest_io_netcdf_write_hdr_system(self, nc, param) call netcdf_io_check( nf90_put_var(nc%id, nc%KE_spin_varid, self%ke_spin, start=[tslot]), "netcdf_io_write_hdr_system nf90_put_var KE_spin_varid" ) call netcdf_io_check( nf90_put_var(nc%id, nc%PE_varid, self%pe, start=[tslot]), "netcdf_io_write_hdr_system nf90_put_var PE_varid" ) call netcdf_io_check( nf90_put_var(nc%id, nc%BE_varid, self%be, start=[tslot]), "netcdf_io_write_hdr_system nf90_put_var BE_varid" ) - call netcdf_io_check( nf90_put_var(nc%id, nc%L_orb_varid, self%Lorbit(:), start=[1,tslot], count=[NDIM,1]), "netcdf_io_write_hdr_system nf90_put_var L_orb_varid" ) - call netcdf_io_check( nf90_put_var(nc%id, nc%Lspin_varid, self%Lspin(:), start=[1,tslot], count=[NDIM,1]), "netcdf_io_write_hdr_system nf90_put_var Lspin_varid" ) - call netcdf_io_check( nf90_put_var(nc%id, nc%L_escape_varid, self%Lescape(:), start=[1,tslot], count=[NDIM,1]), "netcdf_io_write_hdr_system nf90_put_var L_escape_varid" ) - call netcdf_io_check( nf90_put_var(nc%id, nc%Ecollisions_varid, self%Ecollisions, start=[tslot]), "netcdf_io_write_hdr_system nf90_put_var Ecollisions_varid" ) - call netcdf_io_check( nf90_put_var(nc%id, nc%Euntracked_varid, self%Euntracked, start=[tslot]), "netcdf_io_write_hdr_system nf90_put_var Euntracked_varid" ) + call netcdf_io_check( nf90_put_var(nc%id, nc%L_orbit_varid, self%L_orbit(:), start=[1,tslot], count=[NDIM,1]), "netcdf_io_write_hdr_system nf90_put_var L_orbit_varid" ) + call netcdf_io_check( nf90_put_var(nc%id, nc%L_spin_varid, self%L_spin(:), start=[1,tslot], count=[NDIM,1]), "netcdf_io_write_hdr_system nf90_put_var L_spin_varid" ) + call netcdf_io_check( nf90_put_var(nc%id, nc%L_escape_varid, self%L_escape(:), start=[1,tslot], count=[NDIM,1]), "netcdf_io_write_hdr_system nf90_put_var L_escape_varid" ) + call netcdf_io_check( nf90_put_var(nc%id, nc%E_collisions_varid, self%E_collisions, start=[tslot]), "netcdf_io_write_hdr_system nf90_put_var E_collisions_varid" ) + call netcdf_io_check( nf90_put_var(nc%id, nc%E_untracked_varid, self%E_untracked, start=[tslot]), "netcdf_io_write_hdr_system nf90_put_var E_untracked_varid" ) call netcdf_io_check( nf90_put_var(nc%id, nc%GMescape_varid, self%GMescape, start=[tslot]), "netcdf_io_write_hdr_system nf90_put_var GMescape_varid" ) end if @@ -1852,43 +1852,43 @@ module subroutine swiftest_io_param_reader(self, unit, iotype, v_list, iostat, i call swiftest_io_toupper(param_value) if (param_value == "NO" .or. param_value == 'F') param%lfirstenergy = .false. case("EORBIT_ORIG") - read(param_value, *, err = 667, iomsg = iomsg) param%Eorbit_orig + read(param_value, *, err = 667, iomsg = iomsg) param%E_orbit_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) + read(param_value, *, err = 667, iomsg = iomsg) param%L_total_orig(1) do i = 2, NDIM ifirst = ilast + 2 param_value = swiftest_io_get_token(line, ifirst, ilast, iostat) - read(param_value, *, err = 667, iomsg = iomsg) param%Ltot_orig(i) + read(param_value, *, err = 667, iomsg = iomsg) param%L_total_orig(i) end do case("LORBIT_ORIG") - read(param_value, *, err = 667, iomsg = iomsg) param%Lorbit_orig(1) + read(param_value, *, err = 667, iomsg = iomsg) param%L_orbit_orig(1) do i = 2, NDIM ifirst = ilast + 2 param_value = swiftest_io_get_token(line, ifirst, ilast, iostat) - read(param_value, *, err = 667, iomsg = iomsg) param%Lorbit_orig(i) + read(param_value, *, err = 667, iomsg = iomsg) param%L_orbit_orig(i) end do case("LSPIN_ORIG") - read(param_value, *, err = 667, iomsg = iomsg) param%Lspin_orig(1) + read(param_value, *, err = 667, iomsg = iomsg) param%L_spin_orig(1) do i = 2, NDIM ifirst = ilast + 2 param_value = swiftest_io_get_token(line, ifirst, ilast, iostat) - read(param_value, *, err = 667, iomsg = iomsg) param%Lspin_orig(i) + read(param_value, *, err = 667, iomsg = iomsg) param%L_spin_orig(i) end do case("LESCAPE") - read(param_value, *, err = 667, iomsg = iomsg) param%Lescape(1) + read(param_value, *, err = 667, iomsg = iomsg) param%L_escape(1) do i = 2, NDIM ifirst = ilast + 2 param_value = swiftest_io_get_token(line, ifirst, ilast, iostat) - read(param_value, *, err = 667, iomsg = iomsg) param%Lescape(i) + read(param_value, *, err = 667, iomsg = iomsg) param%L_escape(i) end do case("GMESCAPE") read(param_value, *, err = 667, iomsg = iomsg) param%GMescape case("ECOLLISIONS") - read(param_value, *, err = 667, iomsg = iomsg) param%Ecollisions + read(param_value, *, err = 667, iomsg = iomsg) param%E_collisions case("EUNTRACKED") - read(param_value, *, err = 667, iomsg = iomsg) param%Euntracked + read(param_value, *, err = 667, iomsg = iomsg) param%E_untracked case ("MAXID") read(param_value, *, err = 667, iomsg = iomsg) param%maxid case ("MAXID_COLLISION") @@ -2586,14 +2586,14 @@ module subroutine swiftest_io_read_in_system(self, param) call self%pl%read_in(param) call self%tp%read_in(param) ! Copy over param file variable inputs - self%Eorbit_orig = param%Eorbit_orig + self%E_orbit_orig = param%E_orbit_orig self%GMtot_orig = param%GMtot_orig - self%Ltot_orig(:) = param%Ltot_orig(:) - self%Lorbit_orig(:) = param%Lorbit_orig(:) - self%Lspin_orig(:) = param%Lspin_orig(:) - self%Lescape(:) = param%Lescape(:) - self%Ecollisions = param%Ecollisions - self%Euntracked = param%Euntracked + self%L_total_orig(:) = param%L_total_orig(:) + self%L_orbit_orig(:) = param%L_orbit_orig(:) + self%L_spin_orig(:) = param%L_spin_orig(:) + self%L_escape(:) = param%L_escape(:) + self%E_collisions = param%E_collisions + self%E_untracked = param%E_untracked else allocate(tmp_param, source=param) tmp_param%system_history%nc%file_name = param%nc_in diff --git a/src/swiftest/swiftest_module.f90 b/src/swiftest/swiftest_module.f90 index beb883abe..c38470b69 100644 --- a/src/swiftest/swiftest_module.f90 +++ b/src/swiftest/swiftest_module.f90 @@ -333,36 +333,36 @@ module swiftest real(DP) :: be = 0.0_DP !! nbody_system binding energy of all bodies real(DP) :: te = 0.0_DP !! nbody_system total energy real(DP) :: oblpot = 0.0_DP !! nbody_system potential energy due to oblateness of the central body - real(DP), dimension(NDIM) :: Lorbit = 0.0_DP !! nbody_system orbital angular momentum vector - real(DP), dimension(NDIM) :: Lspin = 0.0_DP !! nbody_system spin angular momentum vector - real(DP), dimension(NDIM) :: Ltot = 0.0_DP !! nbody_system angular momentum vector + real(DP), dimension(NDIM) :: L_orbit = 0.0_DP !! nbody_system orbital angular momentum vector + real(DP), dimension(NDIM) :: L_spin = 0.0_DP !! nbody_system spin angular momentum vector + real(DP), dimension(NDIM) :: L_total = 0.0_DP !! nbody_system angular momentum vector real(DP) :: ke_orbit_orig = 0.0_DP !! Initial orbital kinetic energy real(DP) :: ke_spin_orig = 0.0_DP !! Initial spin kinetic energy real(DP) :: pe_orig = 0.0_DP !! Initial potential energy real(DP) :: be_orig = 0.0_DP !! Initial binding energy - real(DP) :: Eorbit_orig = 0.0_DP !! Initial orbital energy + real(DP) :: E_orbit_orig = 0.0_DP !! Initial orbital energy real(DP) :: GMtot_orig = 0.0_DP !! Initial nbody_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) :: Lescape = 0.0_DP !! Angular momentum of bodies that escaped the nbody_system (used for bookeeping) + real(DP), dimension(NDIM) :: L_total_orig = 0.0_DP !! Initial total angular momentum vector + real(DP), dimension(NDIM) :: L_orbit_orig = 0.0_DP !! Initial orbital angular momentum + real(DP), dimension(NDIM) :: L_spin_orig = 0.0_DP !! Initial spin angular momentum vector + real(DP), dimension(NDIM) :: L_escape = 0.0_DP !! Angular momentum of bodies that escaped the nbody_system (used for bookeeping) real(DP) :: GMescape = 0.0_DP !! Mass of bodies that escaped the nbody_system (used for bookeeping) - real(DP) :: Ecollisions = 0.0_DP !! Energy lost from nbody_system due to collisions - real(DP) :: Euntracked = 0.0_DP !! Energy gained from nbody_system due to escaped bodies + real(DP) :: E_collisions = 0.0_DP !! Energy lost from nbody_system due to collisions + real(DP) :: E_untracked = 0.0_DP !! Energy gained from nbody_system due to escaped bodies ! Energy, momentum, and mass errors (used in error reporting) real(DP) :: ke_orbit_error = 0.0_DP real(DP) :: ke_spin_error = 0.0_DP real(DP) :: pe_error = 0.0_DP real(DP) :: be_error = 0.0_DP - real(DP) :: Eorbit_error = 0.0_DP + real(DP) :: E_orbit_error = 0.0_DP real(DP) :: Ecoll_error = 0.0_DP - real(DP) :: Euntracked_error = 0.0_DP - real(DP) :: Etot_error = 0.0_DP - real(DP) :: Lorbit_error = 0.0_DP - real(DP) :: Lspin_error = 0.0_DP - real(DP) :: Lescape_error = 0.0_DP - real(DP) :: Ltot_error = 0.0_DP + real(DP) :: E_untracked_error = 0.0_DP + real(DP) :: te_error = 0.0_DP + real(DP) :: L_orbit_error = 0.0_DP + real(DP) :: L_spin_error = 0.0_DP + real(DP) :: L_escape_error = 0.0_DP + real(DP) :: L_total_error = 0.0_DP real(DP) :: Mtot_error = 0.0_DP real(DP) :: Mescape_error = 0.0_DP diff --git a/src/swiftest/swiftest_util.f90 b/src/swiftest/swiftest_util.f90 index 9d34b1af5..667ad1627 100644 --- a/src/swiftest/swiftest_util.f90 +++ b/src/swiftest/swiftest_util.f90 @@ -1242,9 +1242,9 @@ module subroutine swiftest_util_get_energy_momentum_system(self, param) real(DP) :: hx, hy, hz associate(nbody_system => self, pl => self%pl, npl => self%pl%nbody, cb => self%cb) - nbody_system%Lorbit(:) = 0.0_DP - nbody_system%Lspin(:) = 0.0_DP - nbody_system%Ltot(:) = 0.0_DP + nbody_system%L_orbit(:) = 0.0_DP + nbody_system%L_spin(:) = 0.0_DP + nbody_system%L_total(:) = 0.0_DP nbody_system%ke_orbit = 0.0_DP nbody_system%ke_spin = 0.0_DP @@ -1312,20 +1312,20 @@ module subroutine swiftest_util_get_energy_momentum_system(self, param) nbody_system%ke_orbit = 0.5_DP * (kecb + sum(kepl(1:npl), pl%lmask(1:npl))) if (param%lrotation) nbody_system%ke_spin = 0.5_DP * (kespincb + sum(kespinpl(1:npl), pl%lmask(1:npl))) - nbody_system%Lorbit(1) = Lcborbit(1) + sum(Lplorbitx(1:npl), pl%lmask(1:npl)) - nbody_system%Lorbit(2) = Lcborbit(2) + sum(Lplorbity(1:npl), pl%lmask(1:npl)) - nbody_system%Lorbit(3) = Lcborbit(3) + sum(Lplorbitz(1:npl), pl%lmask(1:npl)) + nbody_system%L_orbit(1) = Lcborbit(1) + sum(Lplorbitx(1:npl), pl%lmask(1:npl)) + nbody_system%L_orbit(2) = Lcborbit(2) + sum(Lplorbity(1:npl), pl%lmask(1:npl)) + nbody_system%L_orbit(3) = Lcborbit(3) + sum(Lplorbitz(1:npl), pl%lmask(1:npl)) if (param%lrotation) then - nbody_system%Lspin(1) = Lcbspin(1) + sum(Lplspinx(1:npl), pl%lmask(1:npl)) - nbody_system%Lspin(2) = Lcbspin(2) + sum(Lplspiny(1:npl), pl%lmask(1:npl)) - nbody_system%Lspin(3) = Lcbspin(3) + sum(Lplspinz(1:npl), pl%lmask(1:npl)) + nbody_system%L_spin(1) = Lcbspin(1) + sum(Lplspinx(1:npl), pl%lmask(1:npl)) + nbody_system%L_spin(2) = Lcbspin(2) + sum(Lplspiny(1:npl), pl%lmask(1:npl)) + nbody_system%L_spin(3) = Lcbspin(3) + sum(Lplspinz(1:npl), pl%lmask(1:npl)) end if nbody_system%be = sum(-3*pl%Gmass(1:npl)*pl%mass(1:npl)/(5*pl%radius(1:npl)), pl%lmask(1:npl)) nbody_system%te = nbody_system%ke_orbit + nbody_system%ke_spin + nbody_system%pe + nbody_system%be - nbody_system%Ltot(:) = nbody_system%Lorbit(:) + nbody_system%Lspin(:) + nbody_system%L_total(:) = nbody_system%L_orbit(:) + nbody_system%L_spin(:) end associate return @@ -2617,7 +2617,11 @@ module subroutine swiftest_util_setup_construct_system(nbody_system, param) call collision_history%reset() select type(nc => collision_history%nc) class is (collision_netcdf_parameters) - nc%file_number = param%iloop / param%dump_cadence + nc%file_name = COLLISION_OUTFILE + if (.not.param%lrestart) then + call nc%initialize(param) + call nc%close() + end if end select allocate(nbody_system%collision_history, source=collision_history) diff --git a/src/symba/symba_discard.f90 b/src/symba/symba_discard.f90 index 216ab5f28..75abb741f 100644 --- a/src/symba/symba_discard.f90 +++ b/src/symba/symba_discard.f90 @@ -116,7 +116,7 @@ subroutine symba_discard_conserve_mtm(pl, nbody_system, param, ipl, lescape_body integer(I4B), intent(in) :: ipl logical, intent(in) :: lescape_body ! Internals - real(DP), dimension(NDIM) :: Lpl, Ltot, Lcb, xcom, vcom + real(DP), dimension(NDIM) :: Lpl, L_total, Lcb, xcom, vcom real(DP) :: pe, be, ke_orbit, ke_spin integer(I4B) :: i, oldstat @@ -142,12 +142,12 @@ subroutine symba_discard_conserve_mtm(pl, nbody_system, param, ipl, lescape_body pe = pe - pl%Gmass(i) * pl%mass(ipl) / norm2(pl%rb(:, ipl) - pl%rb(:, i)) end do - Ltot(:) = 0.0_DP + L_total(:) = 0.0_DP do i = 1, pl%nbody Lpl(:) = pL%mass(i) * (pl%rb(:,i) .cross. pl%vb(:, i)) - Ltot(:) = Ltot(:) + Lpl(:) + L_total(:) = L_total(:) + Lpl(:) end do - Ltot(:) = Ltot(:) + cb%mass * (cb%rb(:) .cross. cb%vb(:)) + L_total(:) = L_total(:) + cb%mass * (cb%rb(:) .cross. cb%vb(:)) call pl%b2h(cb) oldstat = pl%status(ipl) pl%status(ipl) = INACTIVE @@ -156,11 +156,11 @@ subroutine symba_discard_conserve_mtm(pl, nbody_system, param, ipl, lescape_body do i = 1, pl%nbody if (i == ipl) cycle Lpl(:) = pl%mass(i) * (pl%rb(:,i) .cross. pl%vb(:, i)) - Ltot(:) = Ltot(:) - Lpl(:) + L_total(:) = L_total(:) - Lpl(:) end do - Ltot(:) = Ltot(:) - cb%mass * (cb%rb(:) .cross. cb%vb(:)) - nbody_system%Lescape(:) = nbody_system%Lescape(:) + Ltot(:) - if (param%lrotation) nbody_system%Lescape(:) = nbody_system%Lescape + pl%mass(ipl) * pl%radius(ipl)**2 & + L_total(:) = L_total(:) - cb%mass * (cb%rb(:) .cross. cb%vb(:)) + nbody_system%L_escape(:) = nbody_system%L_escape(:) + L_total(:) + if (param%lrotation) nbody_system%L_escape(:) = nbody_system%L_escape + pl%mass(ipl) * pl%radius(ipl)**2 & * pl%Ip(3, ipl) * pl%rot(:, ipl) else @@ -195,8 +195,8 @@ subroutine symba_discard_conserve_mtm(pl, nbody_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 - nbody_system%Ecollisions = nbody_system%Ecollisions + ke_orbit + ke_spin + pe + be - nbody_system%Euntracked = nbody_system%Euntracked - (ke_orbit + ke_spin + pe + be) + nbody_system%E_collisions = nbody_system%E_collisions + ke_orbit + ke_spin + pe + be + nbody_system%E_untracked = nbody_system%E_untracked - (ke_orbit + ke_spin + pe + be) end if end select @@ -345,7 +345,7 @@ module subroutine symba_discard_pl(self, nbody_system, param) class(swiftest_nbody_system), intent(inout) :: nbody_system !! Swiftest nbody system object class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters ! Internals - real(DP) :: Eorbit_before, Eorbit_after + real(DP) :: E_orbit_before, E_orbit_after select type(nbody_system) class is (symba_nbody_system) @@ -362,7 +362,7 @@ module subroutine symba_discard_pl(self, nbody_system, param) if (param%lenergy) then call nbody_system%get_energy_and_momentum(param) - Eorbit_before = nbody_system%te + E_orbit_before = nbody_system%te end if call symba_discard_nonplpl_conservation(self, nbody_system, param) @@ -374,8 +374,8 @@ module subroutine symba_discard_pl(self, nbody_system, param) if (param%lenergy) then call nbody_system%get_energy_and_momentum(param) - Eorbit_after = nbody_system%te - nbody_system%Ecollisions = nbody_system%Ecollisions + (Eorbit_after - Eorbit_before) + E_orbit_after = nbody_system%te + nbody_system%E_collisions = nbody_system%E_collisions + (E_orbit_after - E_orbit_before) end if end associate