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

Commit

Permalink
Added more careful bookkeeping of energy by including the solar gravi…
Browse files Browse the repository at this point in the history 
…tational binding energy
  • Loading branch information
@daminton
daminton committed Jan 10, 2023
1 parent 10f6f9e commit 79c2542
Show file tree
Hide file tree
Showing 7 changed files with 96 additions and 76 deletions.
4 changes: 4 additions & 0 deletions src/collision/collision_io.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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, &
Expand Down Expand Up @@ -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" )
Expand Down Expand Up @@ -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
Expand Down
54 changes: 24 additions & 30 deletions src/collision/collision_util.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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))
Expand Down Expand Up @@ -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
Expand All @@ -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)
Expand All @@ -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

Expand Down
8 changes: 5 additions & 3 deletions src/netcdf_io/netcdf_io_module.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
Loading

0 comments on commit 79c2542

Please sign in to comment.