diff --git a/src/collision/collision_io.f90 b/src/collision/collision_io.f90 index 476700a52..28870ba85 100644 --- a/src/collision/collision_io.f90 +++ b/src/collision/collision_io.f90 @@ -226,6 +226,8 @@ module subroutine collision_io_netcdf_initialize_output(self, param) call netcdf_io_check( nf90_def_var(nc%id, nc%be_varname, nc%out_type,& [ 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%te_varname, nc%out_type,& + [ nc%stage_dimid, nc%collision_id_dimid], nc%TE_varid), "collision_io_netcdf_initialize_output nf90_def_var TE_varid" ) if (param%lrotation) then call netcdf_io_check( nf90_def_var(nc%id, nc%L_orbit_varname, nc%out_type, & @@ -337,6 +339,7 @@ module subroutine collision_io_netcdf_open(self, param, readonly) call netcdf_io_check( nf90_inq_varid(nc%id, nc%ke_orb_varname, nc%ke_orb_varid), "collision_io_netcdf_open nf90_inq_varid ke_orb_varid" ) call netcdf_io_check( nf90_inq_varid(nc%id, nc%pe_varname, nc%pe_varid), "collision_io_netcdf_open nf90_inq_varid pe_varid" ) call netcdf_io_check( nf90_inq_varid(nc%id, nc%be_varname, nc%be_varid), "collision_io_netcdf_open nf90_inq_varid be_varid" ) + call netcdf_io_check( nf90_inq_varid(nc%id, nc%te_varname, nc%te_varid), "collision_io_netcdf_open nf90_inq_varid te_varid" ) call netcdf_io_check( nf90_inq_varid(nc%id, nc%L_orbit_varname, nc%L_orbit_varid), "collision_io_netcdf_open nf90_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_open nf90_inq_varid ke_spin_varid" ) @@ -416,6 +419,7 @@ module subroutine collision_io_netcdf_write_frame_snapshot(self, history, param) 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%te_varid, collider%te(:), start=[ 1, eslot], count=[ 2, 1]), "collision_io_netcdf_write_frame_snapshot nf90_put_var pe_varid tefore" ) 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 diff --git a/src/collision/collision_util.f90 b/src/collision/collision_util.f90 index 6e2627126..acf85f893 100644 --- a/src/collision/collision_util.f90 +++ b/src/collision/collision_util.f90 @@ -231,7 +231,6 @@ module subroutine collision_util_get_energy_and_momentum(self, nbody_system, par self%pe(phase_val) = constraint_system%pe / self%Escale self%be(phase_val) = constraint_system%be / self%Escale self%te(phase_val) = constraint_system%te / self%Escale - if (phase_val == 2) then do concurrent(i = 1:nfrag) fragments%ke_orbit(i) = 0.5_DP * fragments%mass(i) * dot_product(fragments%vc(:,i), fragments%vc(:,i)) @@ -683,6 +682,7 @@ module subroutine collision_util_snapshot(self, param, nbody_system, t, arg) ! Arguments class(collision_snapshot), allocatable, save :: snapshot character(len=:), allocatable :: stage + integer(I4B) :: phase_val if (present(arg)) then stage = arg @@ -694,36 +694,32 @@ module subroutine collision_util_snapshot(self, param, nbody_system, t, arg) class is (swiftest_nbody_system) select type(param) class is (swiftest_parameters) - select case(stage) - case("before") - ! Saves the states of the bodies involved in the collision before the collision is resolved + + select case (stage) + case ("before") + phase_val = 1 allocate(collision_snapshot :: snapshot) allocate(snapshot%collider, source=nbody_system%collider) snapshot%t = t + case ("after") + phase_val = 2 + case default + write(*,*) "collision_util_snapshot requies either 'before' or 'after' passed to 'arg'" + return + end select - ! 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 - + ! Get and record the energy of the system before the collision + call nbody_system%get_energy_and_momentum(param) + snapshot%collider%L_orbit(:,phase_val) = nbody_system%L_orbit(:) + snapshot%collider%L_spin(:,phase_val) = nbody_system%L_spin(:) + snapshot%collider%L_total(:,phase_val) = nbody_system%L_total(:) + snapshot%collider%ke_orbit(phase_val) = nbody_system%ke_orbit + snapshot%collider%ke_spin(phase_val) = nbody_system%ke_spin + snapshot%collider%pe(phase_val) = nbody_system%pe + snapshot%collider%be(phase_val) = nbody_system%be + snapshot%collider%te(phase_val) = nbody_system%te + + if (stage == "after") then select type(before_snap => snapshot%collider%before ) class is (swiftest_nbody_system) select type(before_orig => nbody_system%collider%before) @@ -743,9 +739,7 @@ module subroutine collision_util_snapshot(self, param, nbody_system, t, arg) ! Save the snapshot for posterity call self%save(snapshot) deallocate(snapshot) - case default - write(*,*) "collision_util_snapshot requies either 'before' or 'after' passed to 'arg'" - end select + end if end select end select diff --git a/src/netcdf_io/netcdf_io_module.f90 b/src/netcdf_io/netcdf_io_module.f90 index 74325f809..b28b2080e 100644 --- a/src/netcdf_io/netcdf_io_module.f90 +++ b/src/netcdf_io/netcdf_io_module.f90 @@ -106,15 +106,17 @@ 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_orbit_varname = "L_orbit" !! name of the orbital angular momentum vector variable + character(NAMELEN) :: te_varname = "TE" !! name of the system binding energy variable + integer(I4B) :: TE_varid !! ID for the system binding energy 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) :: E_collisions_varname = "E_collisions" !! name of the escaped angular momentum y variable + 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) + 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 diff --git a/src/swiftest/swiftest_io.f90 b/src/swiftest/swiftest_io.f90 index 7bf7c1117..b9294867b 100644 --- a/src/swiftest/swiftest_io.f90 +++ b/src/swiftest/swiftest_io.f90 @@ -121,11 +121,11 @@ module subroutine swiftest_io_conservation_report(self, param, lterminal) logical, intent(in) :: lterminal !! Indicates whether to output information to the terminal screen ! Internals 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) :: ke_orbit_now, ke_spin_now, pe_now, E_orbit_now, be_now, be_cb_now, be_cb_orig, te_now real(DP) :: GMtot_now character(len=STRMAX) :: errmsg integer(I4B), parameter :: EGYIU = 72 - character(len=*), parameter :: EGYTERMFMT = '(" DL/L0 = ", ES12.5, "; DE/|E0| = ", ES12.5, "; DM/M0 = ", ES12.5)' + character(len=*), parameter :: EGYTERMFMT = '(" DL/L0 = ", ES12.5, "; DE_orbit/|E0| = ", ES12.5, "; DE_total/|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 => self%system_history%nc) @@ -137,9 +137,11 @@ 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 + be_cb_now = nbody_system%be_cb + te_now = nbody_system%te 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 + E_orbit_now = ke_orbit_now + pe_now L_total_now(:) = nbody_system%L_total(:) + nbody_system%L_escape(:) GMtot_now = nbody_system%GMtot + nbody_system%GMescape @@ -148,6 +150,7 @@ 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%te_orig = te_now nbody_system%E_orbit_orig = E_orbit_now nbody_system%GMtot_orig = GMtot_now nbody_system%L_orbit_orig(:) = L_orbit_now(:) @@ -157,22 +160,28 @@ module subroutine swiftest_io_conservation_report(self, param, lterminal) end if if (.not.param%lfirstenergy) then + 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) + + be_cb_orig = -(3 * cb%GM0**2 / param%GU) / (5 * cb%R0) + nbody_system%be_error = (be_now - nbody_system%be_orig) / abs(nbody_system%te_orig) + (be_cb_now - be_cb_orig) / abs(nbody_system%te_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%Ecoll_error = nbody_system%E_collisions / abs(nbody_system%te_orig) + nbody_system%E_untracked_error = nbody_system%E_untracked / abs(nbody_system%te_orig) + nbody_system%te_error = (nbody_system%te - nbody_system%te_orig - nbody_system%E_collisions - nbody_system%E_untracked) / abs(nbody_system%te_orig) + (be_cb_now - be_cb_orig) / abs(nbody_system%te_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%L_total_error, nbody_system%te_error,nbody_system%Mtot_error + + if (lterminal) write(display_unit, EGYTERMFMT) nbody_system%L_total_error, nbody_system%E_orbit_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 @@ -567,7 +576,7 @@ module subroutine swiftest_io_netcdf_get_t0_values_system(self, param) real(DP), dimension(:), allocatable :: vals real(DP), dimension(1) :: rtemp real(DP), dimension(NDIM) :: rot0, Ip0 - real(DP) :: KE_orb_orig, KE_spin_orig, PE_orig, BE_orig, mass0 + real(DP) :: mass0 associate (nc => self%system_history%nc, cb => self%cb) call nc%open(param, readonly=.true.) @@ -580,18 +589,21 @@ module subroutine swiftest_io_netcdf_get_t0_values_system(self, param) if (param%lenergy) then call netcdf_io_check( nf90_get_var(nc%id, nc%KE_orb_varid, rtemp, start=[tslot], count=[1]), "netcdf_io_get_t0_values_system KE_orb_varid" ) - KE_orb_orig = rtemp(1) + self%ke_orbit_orig = rtemp(1) call netcdf_io_check( nf90_get_var(nc%id, nc%KE_spin_varid, rtemp, start=[tslot], count=[1]), "netcdf_io_get_t0_values_system KE_spin_varid" ) - KE_spin_orig = rtemp(1) + self%ke_spin_orig = rtemp(1) call netcdf_io_check( nf90_get_var(nc%id, nc%PE_varid, rtemp, start=[tslot], count=[1]), "netcdf_io_get_t0_values_system PE_varid" ) - PE_orig = rtemp(1) + self%pe_orig = rtemp(1) call netcdf_io_check( nf90_get_var(nc%id, nc%BE_varid, rtemp, start=[tslot], count=[1]), "netcdf_io_get_t0_values_system BE_varid" ) - BE_orig = rtemp(1) + self%be_orig = rtemp(1) + + call netcdf_io_check( nf90_get_var(nc%id, nc%TE_varid, rtemp, start=[tslot], count=[1]), "netcdf_io_get_t0_values_system TE_varid" ) + self%te_orig = rtemp(1) - self%E_orbit_orig = KE_orb_orig + KE_spin_orig + PE_orig + BE_orig + self%E_orbit_orig = self%ke_orbit_orig + self%pe_orig call netcdf_io_check( nf90_get_var(nc%id, nc%L_orbit_varid, self%L_orbit_orig(:), start=[1,tslot], 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,tslot], count=[NDIM,1]), "netcdf_io_get_t0_values_system L_spin_varid" ) @@ -603,6 +615,7 @@ module subroutine swiftest_io_netcdf_get_t0_values_system(self, param) cb%GM0 = vals(1) cb%dGM = cb%Gmass - cb%GM0 + mass0 = cb%GM0 / param%GU call netcdf_io_check( nf90_get_var(nc%id, nc%radius_varid, rtemp, start=[1,tslot], count=[1,1]), "netcdf_io_get_t0_values_system radius_varid" ) cb%R0 = rtemp(1) @@ -610,7 +623,6 @@ module subroutine swiftest_io_netcdf_get_t0_values_system(self, param) if (param%lrotation) then call netcdf_io_check( nf90_get_var(nc%id, nc%rot_varid, rot0, start=[1,1,tslot], count=[NDIM,1,1]), "netcdf_io_get_t0_values_system rot_varid" ) call netcdf_io_check( nf90_get_var(nc%id, nc%Ip_varid, Ip0, start=[1,1,tslot], count=[NDIM,1,1]), "netcdf_io_get_t0_values_system Ip_varid" ) - mass0 = cb%GM0 / param%GU cb%L0(:) = Ip0(3) * mass0 * cb%R0**2 * rot0(:) end if @@ -753,7 +765,8 @@ module subroutine swiftest_io_netcdf_initialize_output(self, param) call netcdf_io_check( nf90_def_var(nc%id, nc%ke_orb_varname, nc%out_type, nc%time_dimid, nc%KE_orb_varid), "netcdf_io_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%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%be_varname, nc%out_type, nc%time_dimid, nc%BE_varid), "netcdf_io_initialize_output nf90_def_var BE_varid" ) + call netcdf_io_check( nf90_def_var(nc%id, nc%te_varname, nc%out_type, nc%time_dimid, nc%TE_varid), "netcdf_io_initialize_output nf90_def_var TE_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" ) @@ -920,6 +933,7 @@ 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%te_varname, nc%TE_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) @@ -1257,6 +1271,8 @@ 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%te_varname, nc%TE_varid) + if (status == NF90_NOERR) call netcdf_io_check( nf90_get_var(nc%id, nc%TE_varid, self%te, start=[tslot]), "netcdf_io_read_hdr_system nf90_getvar TE_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) @@ -1667,6 +1683,7 @@ 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%TE_varid, self%te, start=[tslot]), "netcdf_io_write_hdr_system nf90_put_var TE_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" ) diff --git a/src/swiftest/swiftest_module.f90 b/src/swiftest/swiftest_module.f90 index 940948278..211e048f3 100644 --- a/src/swiftest/swiftest_module.f90 +++ b/src/swiftest/swiftest_module.f90 @@ -334,22 +334,24 @@ 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) :: 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), 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) :: E_orbit_orig = 0.0_DP !! Initial orbital energy + real(DP) :: be_orig = 0.0_DP !! Initial gravitational binding energy + real(DP) :: te_orig = 0.0_DP !! Initial total energy (sum of all sources of energy tracked) + real(DP) :: be_cb = 0.0_DP !! Binding energy of central body (usually orders of magnitude larger than the rest of the system, and therefore tracked seperately) + 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) :: 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), 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) :: 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 + 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 diff --git a/src/swiftest/swiftest_util.f90 b/src/swiftest/swiftest_util.f90 index b83ffc8e1..891038a6e 100644 --- a/src/swiftest/swiftest_util.f90 +++ b/src/swiftest/swiftest_util.f90 @@ -1390,6 +1390,7 @@ module subroutine swiftest_util_get_energy_and_momentum_system(self, param) nbody_system%GMtot = cb%Gmass + sum(pl%Gmass(1:npl), pl%lmask(1:npl)) kecb = cb%mass * dot_product(cb%vb(:), cb%vb(:)) + nbody_system%be_cb = -3*cb%Gmass * cb%mass / (5 * cb%radius) Lcborbit(:) = cb%mass * (cb%rb(:) .cross. cb%vb(:)) do concurrent (i = 1:npl, pl%lmask(i)) @@ -1445,8 +1446,8 @@ module subroutine swiftest_util_get_energy_and_momentum_system(self, param) nbody_system%L_orbit(j) = Lcborbit(j) + sum(Lplorbit(j,1:npl), pl%lmask(1:npl)) end do - 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%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%L_total(:) = nbody_system%L_orbit(:) + nbody_system%L_spin(:) end associate diff --git a/src/symba/symba_discard.f90 b/src/symba/symba_discard.f90 index 674772768..471fe230c 100644 --- a/src/symba/symba_discard.f90 +++ b/src/symba/symba_discard.f90 @@ -104,7 +104,7 @@ subroutine symba_discard_cb_pl(pl, nbody_system, param) end subroutine symba_discard_cb_pl - subroutine symba_discard_conserve_mtm(pl, nbody_system, param, ipl, lescape_body) + subroutine symba_discard_conserve_energy_and_momentum(pl, nbody_system, param, ipl, lescape_body) !! author: David A. Minton !! !! Conserves nbody_system momentum when a body is lost from the nbody_system or collides with central body @@ -117,7 +117,7 @@ subroutine symba_discard_conserve_mtm(pl, nbody_system, param, ipl, lescape_body logical, intent(in) :: lescape_body ! Internals real(DP), dimension(NDIM) :: Lpl, L_total, Lcb, xcom, vcom, drot0, drot1 - real(DP) :: pe, be, ke_orbit, ke_spin + real(DP) :: pe, be, ke_orbit, ke_spin, becb0, becb1 integer(I4B) :: i, oldstat select type(cb => nbody_system%cb) @@ -133,14 +133,15 @@ subroutine symba_discard_conserve_mtm(pl, nbody_system, param, ipl, lescape_body ke_spin = 0.0_DP end if - ! Add the pre-collision ke of the central body to the records - ! Add planet mass to central body accumulator if (lescape_body) then nbody_system%GMescape = nbody_system%GMescape + pl%Gmass(ipl) do i = 1, pl%nbody if (i == ipl) cycle pe = pe - pl%Gmass(i) * pl%mass(ipl) / norm2(pl%rb(:, ipl) - pl%rb(:, i)) end do + + 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) L_total(:) = 0.0_DP do i = 1, pl%nbody @@ -175,36 +176,35 @@ subroutine symba_discard_conserve_mtm(pl, nbody_system, param, ipl, lescape_body ke_orbit = ke_orbit + 0.5_DP * cb%mass * dot_product(cb%vb(:), cb%vb(:)) if (param%lrotation) ke_spin = ke_spin + 0.5_DP * cb%mass * cb%radius**2 * cb%Ip(3) * dot_product(cb%rot(:), cb%rot(:)) ! Update mass of central body to be consistent with its total mass + becb0 = -(3 * cb%Gmass * cb%mass) / (5 * cb%radius) cb%dGM = cb%dGM + pl%Gmass(ipl) cb%dR = cb%dR + 1.0_DP / 3.0_DP * (pl%radius(ipl) / cb%radius)**3 - 2.0_DP / 9.0_DP * (pl%radius(ipl) / cb%radius)**6 cb%Gmass = cb%GM0 + cb%dGM cb%mass = cb%Gmass / param%GU cb%radius = cb%R0 + cb%dR param%rmin = cb%radius + becb1 = -(3 * cb%Gmass * cb%mass) / (5 * cb%radius) + ! Add planet angular momentum to central body accumulator cb%dL(:) = Lpl(:) + Lcb(:) + cb%dL(:) ! Update rotation of central body to by consistent with its angular momentum if (param%lrotation) then drot0(:) = cb%L0(:)/ (cb%Ip(3) * cb%mass * cb%radius**2) drot1(:) = cb%dL(:) / (cb%Ip(3) * cb%mass * cb%radius**2) - cb%rot(:) = drot0(:) + drot1(:) ke_spin = ke_spin - 0.5_DP * cb%mass * cb%radius**2 * cb%Ip(3) * dot_product(cb%rot(:), cb%rot(:)) end if cb%rb(:) = xcom(:) cb%vb(:) = vcom(:) ke_orbit = ke_orbit - 0.5_DP * cb%mass * dot_product(cb%vb(:), cb%vb(:)) + + ! Add the change in central body binding energy to the collision energy tracker + nbody_system%E_collisions = nbody_system%E_collisions + (becb1 - becb0) end if call pl%b2h(cb) - ! 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%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 return - end subroutine symba_discard_conserve_mtm + end subroutine symba_discard_conserve_energy_and_momentum subroutine symba_discard_nonplpl(pl, nbody_system, param) @@ -280,7 +280,7 @@ subroutine symba_discard_nonplpl_conservation(pl, nbody_system, param) cycle end if ! Conserve all the quantities - call symba_discard_conserve_mtm(pl, nbody_system, param, discard_index_list(i), lescape) + call symba_discard_conserve_energy_and_momentum(pl, nbody_system, param, discard_index_list(i), lescape) end do end associate