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

Commit

Permalink
When GR is enabled and NetCDF is the output or input file type, the c…
Browse files Browse the repository at this point in the history 
…onversion between velocity and pseudovelocity is now handled with greater care. Pseudovelocity values are now stored in the NetCDF file so that bit-identical runs can be started from previous output files and dump files.
  • Loading branch information
@daminton
daminton committed Sep 15, 2022
1 parent 1dfcb8e commit b80bb42
Show file tree
Hide file tree
Showing 6 changed files with 131 additions and 35 deletions.
22 changes: 12 additions & 10 deletions src/io/io.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2184,24 +2184,26 @@ module subroutine io_write_frame_system(self, param)
end if
end if

if (param%lgr) then
call pl%pv2v(param)
call tp%pv2v(param)
end if

if ((param%out_form == EL) .or. (param%out_form == XVEL)) then ! Do an orbital element conversion prior to writing out the frame, as we have access to the central body here
call pl%xv2el(cb)
call tp%xv2el(cb)
end if

! Write out each data type frame
if ((param%out_type == REAL4_TYPE) .or. (param%out_type == REAL8_TYPE)) then
! For these data types, do these conversion here before writing the output.
if (param%lgr) then
call pl%pv2v(param)
call tp%pv2v(param)
end if

if ((param%out_form == EL) .or. (param%out_form == XVEL)) then ! Do an orbital element conversion prior to writing out the frame, as we have access to the central body here
call pl%xv2el(cb)
call tp%xv2el(cb)
end if

call self%write_hdr(iu, param)
call cb%write_frame(iu, param)
call pl%write_frame(iu, param)
call tp%write_frame(iu, param)
close(iu, err = 667, iomsg = errmsg)
else if ((param%out_type == NETCDF_FLOAT_TYPE) .or. (param%out_type == NETCDF_DOUBLE_TYPE)) then
! For NetCDF output, because we want to store the pseudovelocity separately from the true velocity, we need to do the orbital element conversion internally
call self%write_hdr(param%nciu, param)
call cb%write_frame(param%nciu, param)
call pl%write_frame(param%nciu, param)
Expand Down
17 changes: 17 additions & 0 deletions src/modules/swiftest_classes.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -32,6 +32,9 @@ module swiftest_classes
integer(I4B) :: vhx_varid !! NetCDF ID for the heliocentric velocity x variable
integer(I4B) :: vhy_varid !! NetCDF ID for the heliocentric velocity y variable
integer(I4B) :: vhz_varid !! NetCDF ID for the heliocentric velocity z variable
integer(I4B) :: gr_pseudo_vhx_varid !! NetCDF ID for the heliocentric pseudovelocity x variable (used in GR)
integer(I4B) :: gr_pseudo_vhy_varid !! NetCDF ID for the heliocentric pseudovelocity y variable (used in GR)
integer(I4B) :: gr_pseudo_vhz_varid !! NetCDF ID for the heliocentric psuedovelocity z variable (used in GR)
integer(I4B) :: Gmass_varid !! NetCDF ID for the mass variable
integer(I4B) :: rhill_varid !! NetCDF ID for the hill radius variable
integer(I4B) :: radius_varid !! NetCDF ID for the radius variable
Expand Down Expand Up @@ -80,6 +83,7 @@ module swiftest_classes
integer(I4B) :: discard_body_id_varid !! NetCDF ID for the id of the other body involved in the discard
integer(I4B) :: id_chunk !! Chunk size for the id dimension variables
integer(I4B) :: time_chunk !! Chunk size for the time dimension variables
logical :: lpseudo_vel_exists = .false. !! Logical flag to indicate whether or not the pseudovelocity vectors were present in an old file.
contains
procedure :: close => netcdf_close !! Closes an open NetCDF file
procedure :: flush => netcdf_flush !! Flushes the current buffer to disk by closing and re-opening the file.
Expand Down Expand Up @@ -1108,6 +1112,19 @@ module pure subroutine orbel_xv2aqt(mu, px, py, pz, vx, vy, vz, a, q, capm, tper
real(DP), intent(out) :: tperi !! time of pericenter passage
end subroutine orbel_xv2aqt

module pure subroutine orbel_xv2el(mu, px, py, pz, vx, vy, vz, a, e, inc, capom, omega, capm)
implicit none
real(DP), intent(in) :: mu !! Gravitational constant
real(DP), intent(in) :: px,py,pz !! Position vector
real(DP), intent(in) :: vx,vy,vz !! Velocity vector
real(DP), intent(out) :: a !! semimajor axis
real(DP), intent(out) :: e !! eccentricity
real(DP), intent(out) :: inc !! inclination
real(DP), intent(out) :: capom !! longitude of ascending node
real(DP), intent(out) :: omega !! argument of periapsis
real(DP), intent(out) :: capm !! mean anomaly
end subroutine orbel_xv2el

module subroutine orbel_xv2el_vec(self, cb)
implicit none
class(swiftest_body), intent(inout) :: self !! Swiftest body object
Expand Down
3 changes: 3 additions & 0 deletions src/modules/swiftest_globals.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -157,6 +157,9 @@ module swiftest_globals
character(*), parameter :: VHX_VARNAME = "vhx" !! NetCDF name of the heliocentric velocity x variable
character(*), parameter :: VHY_VARNAME = "vhy" !! NetCDF name of the heliocentric velocity y variable
character(*), parameter :: VHZ_VARNAME = "vhz" !! NetCDF name of the heliocentric velocity z variable
character(*), parameter :: GR_PSEUDO_VHX_VARNAME = "gr_pseudo_vhx" !! NetCDF name of the heliocentric pseudovelocity x variable (used in GR only)
character(*), parameter :: GR_PSEUDO_VHY_VARNAME = "gr_pseudo_vhy" !! NetCDF name of the heliocentric pseudovelocity y variable (used in GR only)
character(*), parameter :: GR_PSEUDO_VHZ_VARNAME = "gr_pseudo_vhz" !! NetCDF name of the heliocentric pseudovelocity z variable (used in GR only)
character(*), parameter :: GMASS_VARNAME = "Gmass" !! NetCDF name of the mass variable
character(*), parameter :: RHILL_VARNAME = "rhill" !! NetCDF name of the hill radius variable
character(*), parameter :: RADIUS_VARNAME = "radius" !! NetCDF name of the radius variable
Expand Down
120 changes: 97 additions & 23 deletions src/netcdf/netcdf.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -213,6 +213,16 @@ module subroutine netcdf_initialize_output(self, param)
call check( nf90_def_var(self%ncid, VHX_VARNAME, self%out_type, [self%id_dimid, self%time_dimid], self%vhx_varid) )
call check( nf90_def_var(self%ncid, VHY_VARNAME, self%out_type, [self%id_dimid, self%time_dimid], self%vhy_varid) )
call check( nf90_def_var(self%ncid, VHZ_VARNAME, self%out_type, [self%id_dimid, self%time_dimid], self%vhz_varid) )

!! When GR is enabled, we need to save the pseudovelocity vectors in addition to the true heliocentric velocity vectors, otherwise
!! we cannnot expect bit-identical runs from restarted runs with GR enabled due to floating point errors during the conversion.
if (param%lgr) then
call check( nf90_def_var(self%ncid, GR_PSEUDO_VHX_VARNAME, self%out_type, [self%id_dimid, self%time_dimid], self%gr_pseudo_vhx_varid) )
call check( nf90_def_var(self%ncid, GR_PSEUDO_VHY_VARNAME, self%out_type, [self%id_dimid, self%time_dimid], self%gr_pseudo_vhy_varid) )
call check( nf90_def_var(self%ncid, GR_PSEUDO_VHZ_VARNAME, self%out_type, [self%id_dimid, self%time_dimid], self%gr_pseudo_vhz_varid) )
self%lpseudo_vel_exists = .true.
end if

end if

if ((param%out_form == EL) .or. (param%out_form == XVEL)) then
Expand Down Expand Up @@ -329,7 +339,7 @@ module subroutine netcdf_open(self, param, readonly)
class(swiftest_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
integer(I4B) :: mode, status
character(len=NF90_MAX_NAME) :: str_dim_name

mode = NF90_WRITE
Expand Down Expand Up @@ -359,6 +369,25 @@ module subroutine netcdf_open(self, param, readonly)
call check( nf90_inq_varid(self%ncid, VHX_VARNAME, self%vhx_varid))
call check( nf90_inq_varid(self%ncid, VHY_VARNAME, self%vhy_varid))
call check( nf90_inq_varid(self%ncid, VHZ_VARNAME, self%vhz_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.
status = nf90_inq_varid(self%ncid, GR_PSEUDO_VHX_VARNAME, self%gr_pseudo_vhx_varid)
self%lpseudo_vel_exists = (status == nf90_noerr)
if (self%lpseudo_vel_exists) then
status = nf90_inq_varid(self%ncid, GR_PSEUDO_VHY_VARNAME, self%gr_pseudo_vhy_varid)
self%lpseudo_vel_exists = (status == nf90_noerr)
if (self%lpseudo_vel_exists) then
status = nf90_inq_varid(self%ncid, GR_PSEUDO_VHZ_VARNAME, self%gr_pseudo_vhz_varid)
self%lpseudo_vel_exists = (status == nf90_noerr)
end if
end if
if (.not.self%lpseudo_vel_exists) then
write(*,*) "Warning! Pseudovelocity not found in input file for GR enabled run."
end if

end if

end if

if ((param%out_form == EL) .or. (param%out_form == XVEL)) then
Expand Down Expand Up @@ -430,8 +459,6 @@ module subroutine netcdf_open(self, param, readonly)
call check( nf90_inq_varid(self%ncid, J2RP2_VARNAME, self%j2rp2_varid) )
call check( nf90_inq_varid(self%ncid, J4RP4_VARNAME, self%j4rp4_varid) )

if (param%lgr) then
end if

return
end subroutine netcdf_open
Expand Down Expand Up @@ -516,17 +543,31 @@ module function netcdf_read_frame_system(self, iu, param) result(ierr)
if (npl > 0) pl%xh(3,:) = pack(rtemp, plmask)
if (ntp > 0) tp%xh(3,:) = pack(rtemp, tpmask)

call check( nf90_get_var(iu%ncid, iu%vhx_varid, rtemp, start=[1, tslot]) )
if (npl > 0) pl%vh(1,:) = pack(rtemp, plmask)
if (ntp > 0) tp%vh(1,:) = pack(rtemp, tpmask)

call check( nf90_get_var(iu%ncid, iu%vhy_varid, rtemp, start=[1, tslot]) )
if (npl > 0) pl%vh(2,:) = pack(rtemp, plmask)
if (ntp > 0) tp%vh(2,:) = pack(rtemp, tpmask)

call check( nf90_get_var(iu%ncid, iu%vhz_varid, rtemp, start=[1, tslot]) )
if (npl > 0) pl%vh(3,:) = pack(rtemp, plmask)
if (ntp > 0) tp%vh(3,:) = pack(rtemp, tpmask)
if (param%lgr .and. iu%lpseudo_vel_exists) then
call check( nf90_get_var(iu%ncid, iu%gr_pseudo_vhx_varid, rtemp, start=[1, tslot]) )
if (npl > 0) pl%vh(1,:) = pack(rtemp, plmask)
if (ntp > 0) tp%vh(1,:) = pack(rtemp, tpmask)

call check( nf90_get_var(iu%ncid, iu%gr_pseudo_vhy_varid, rtemp, start=[1, tslot]) )
if (npl > 0) pl%vh(2,:) = pack(rtemp, plmask)
if (ntp > 0) tp%vh(2,:) = pack(rtemp, tpmask)

call check( nf90_get_var(iu%ncid, iu%gr_pseudo_vhz_varid, rtemp, start=[1, tslot]) )
if (npl > 0) pl%vh(3,:) = pack(rtemp, plmask)
if (ntp > 0) tp%vh(3,:) = pack(rtemp, tpmask)
else
call check( nf90_get_var(iu%ncid, iu%vhx_varid, rtemp, start=[1, tslot]) )
if (npl > 0) pl%vh(1,:) = pack(rtemp, plmask)
if (ntp > 0) tp%vh(1,:) = pack(rtemp, tpmask)

call check( nf90_get_var(iu%ncid, iu%vhy_varid, rtemp, start=[1, tslot]) )
if (npl > 0) pl%vh(2,:) = pack(rtemp, plmask)
if (ntp > 0) tp%vh(2,:) = pack(rtemp, tpmask)

call check( nf90_get_var(iu%ncid, iu%vhz_varid, rtemp, start=[1, tslot]) )
if (npl > 0) pl%vh(3,:) = pack(rtemp, plmask)
if (ntp > 0) tp%vh(3,:) = pack(rtemp, tpmask)
end if
end if

if ((param%in_form == EL) .or. (param%in_form == XVEL)) then
Expand Down Expand Up @@ -620,6 +661,15 @@ module function netcdf_read_frame_system(self, iu, param) result(ierr)
call check( nf90_get_var(iu%ncid, iu%j4rp4_varid, cb%j4rp4, start=[tslot]) )

call self%read_particle_info(iu, param, plmask, tpmask)

! if this is a GR-enabled run, check to see if we got the pseudovelocities in. Otherwise, we'll need to generate them.
if (param%lgr .and. .not.(iu%lpseudo_vel_exists)) then
call pl%set_mu(cb)
call tp%set_mu(cb)
call pl%v2pv(param)
call tp%v2pv(param)
end if

end associate

call iu%close()
Expand Down Expand Up @@ -865,6 +915,8 @@ module subroutine netcdf_write_frame_base(self, iu, param)
! Internals
integer(I4B) :: i, j, tslot, idslot, old_mode
integer(I4B), dimension(:), allocatable :: ind
real(DP), dimension(NDIM) :: vh !! Temporary variable to store heliocentric velocity values when converting from pseudovelocity in GR-enabled runs
real(DP) :: a, e, inc, omega, capom, capm

call self%write_particle_info(iu, param)

Expand All @@ -882,22 +934,44 @@ module subroutine netcdf_write_frame_base(self, iu, param)
j = ind(i)
idslot = self%id(j) + 1

!! Convert from pseudovelocity to heliocentric without replacing the current value of pseudovelocity
if (param%lgr) call gr_pseudovel2vel(param, self%mu(j), self%xh(:, j), self%vh(:, j), vh(:))

if ((param%out_form == XV) .or. (param%out_form == XVEL)) then
call check( nf90_put_var(iu%ncid, iu%xhx_varid, self%xh(1, j), start=[idslot, tslot]) )
call check( nf90_put_var(iu%ncid, iu%xhy_varid, self%xh(2, j), start=[idslot, tslot]) )
call check( nf90_put_var(iu%ncid, iu%xhz_varid, self%xh(3, j), start=[idslot, tslot]) )
call check( nf90_put_var(iu%ncid, iu%vhx_varid, self%vh(1, j), start=[idslot, tslot]) )
call check( nf90_put_var(iu%ncid, iu%vhy_varid, self%vh(2, j), start=[idslot, tslot]) )
call check( nf90_put_var(iu%ncid, iu%vhz_varid, self%vh(3, j), start=[idslot, tslot]) )
if (param%lgr) then !! Convert from pseudovelocity to heliocentric without replacing the current value of pseudovelocity
call check( nf90_put_var(iu%ncid, iu%vhx_varid, vh(1), start=[idslot, tslot]) )
call check( nf90_put_var(iu%ncid, iu%vhy_varid, vh(2), start=[idslot, tslot]) )
call check( nf90_put_var(iu%ncid, iu%vhz_varid, vh(3), start=[idslot, tslot]) )
call check( nf90_put_var(iu%ncid, iu%gr_pseudo_vhx_varid, self%vh(1, j), start=[idslot, tslot]) )
call check( nf90_put_var(iu%ncid, iu%gr_pseudo_vhy_varid, self%vh(2, j), start=[idslot, tslot]) )
call check( nf90_put_var(iu%ncid, iu%gr_pseudo_vhz_varid, self%vh(3, j), start=[idslot, tslot]) )

else
call check( nf90_put_var(iu%ncid, iu%vhx_varid, self%vh(1, j), start=[idslot, tslot]) )
call check( nf90_put_var(iu%ncid, iu%vhy_varid, self%vh(2, j), start=[idslot, tslot]) )
call check( nf90_put_var(iu%ncid, iu%vhz_varid, self%vh(3, j), start=[idslot, tslot]) )
end if
end if

if ((param%out_form == EL) .or. (param%out_form == XVEL)) then
call check( nf90_put_var(iu%ncid, iu%a_varid, self%a(j), start=[idslot, tslot]) )
call check( nf90_put_var(iu%ncid, iu%e_varid, self%e(j), start=[idslot, tslot]) )
call check( nf90_put_var(iu%ncid, iu%inc_varid, self%inc(j) * RAD2DEG, start=[idslot, tslot]) )
call check( nf90_put_var(iu%ncid, iu%capom_varid, self%capom(j) * RAD2DEG, start=[idslot, tslot]) )
call check( nf90_put_var(iu%ncid, iu%omega_varid, self%omega(j) * RAD2DEG, start=[idslot, tslot]) )
call check( nf90_put_var(iu%ncid, iu%capm_varid, self%capm(j) * RAD2DEG, start=[idslot, tslot]) )
if (param%lgr) then !! For GR-enabled runs, use the true value of velocity computed above
call orbel_xv2el(self%mu(j), self%xh(1,j), self%xh(2,j), self%xh(3,j), &
vh(1), vh(2), vh(3), &
a, e, inc, capom, omega, capm)
else !! For non-GR runs just convert from the velocity we have
call orbel_xv2el(self%mu(j), self%xh(1,j), self%xh(2,j), self%xh(3,j), &
self%vh(1,j), self%vh(2,j), self%vh(3,j), &
a, e, inc, capom, omega, capm)
end if
call check( nf90_put_var(iu%ncid, iu%a_varid, a, start=[idslot, tslot]) )
call check( nf90_put_var(iu%ncid, iu%e_varid, e, start=[idslot, tslot]) )
call check( nf90_put_var(iu%ncid, iu%inc_varid, inc * RAD2DEG, start=[idslot, tslot]) )
call check( nf90_put_var(iu%ncid, iu%capom_varid, capom * RAD2DEG, start=[idslot, tslot]) )
call check( nf90_put_var(iu%ncid, iu%omega_varid, omega * RAD2DEG, start=[idslot, tslot]) )
call check( nf90_put_var(iu%ncid, iu%capm_varid, capm * RAD2DEG, start=[idslot, tslot]) )
end if

select type(self)
Expand Down
2 changes: 1 addition & 1 deletion src/orbel/orbel.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -888,7 +888,7 @@ module subroutine orbel_xv2el_vec(self, cb)
end subroutine orbel_xv2el_vec


pure subroutine orbel_xv2el(mu, px, py, pz, vx, vy, vz, a, e, inc, capom, omega, capm)
module pure subroutine orbel_xv2el(mu, px, py, pz, vx, vy, vz, a, e, inc, capom, omega, capm)
!! author: David A. Minton
!!
!! Compute osculating orbital elements from relative Cartesian position and velocity
Expand Down
2 changes: 1 addition & 1 deletion src/whm/whm_setup.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -80,7 +80,7 @@ module subroutine whm_setup_initialize_system(self, param)
call self%tp_discards%setup(0, param)
call self%pl%set_mu(self%cb)
call self%tp%set_mu(self%cb)
if (param%lgr) then
if (param%lgr .and. ((param%in_type == REAL8_TYPE) .or. (param%in_type == REAL4_TYPE))) then !! pseudovelocity conversion for NetCDF input files is handled by NetCDF routines
call self%pl%v2pv(param)
call self%tp%v2pv(param)
end if
Expand Down

0 comments on commit b80bb42

Please sign in to comment.