From b80bb426297980fa0f001367f52661fa6e3eeceb Mon Sep 17 00:00:00 2001 From: David A Minton Date: Thu, 15 Sep 2022 14:46:06 -0400 Subject: [PATCH] When GR is enabled and NetCDF is the output or input file type, the conversion 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. --- src/io/io.f90 | 22 +++--- src/modules/swiftest_classes.f90 | 17 +++++ src/modules/swiftest_globals.f90 | 3 + src/netcdf/netcdf.f90 | 120 +++++++++++++++++++++++++------ src/orbel/orbel.f90 | 2 +- src/whm/whm_setup.f90 | 2 +- 6 files changed, 131 insertions(+), 35 deletions(-) diff --git a/src/io/io.f90 b/src/io/io.f90 index 8fd3381a8..2ecac13b0 100644 --- a/src/io/io.f90 +++ b/src/io/io.f90 @@ -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) diff --git a/src/modules/swiftest_classes.f90 b/src/modules/swiftest_classes.f90 index 0564a2326..0d8a1fd57 100644 --- a/src/modules/swiftest_classes.f90 +++ b/src/modules/swiftest_classes.f90 @@ -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 @@ -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. @@ -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 diff --git a/src/modules/swiftest_globals.f90 b/src/modules/swiftest_globals.f90 index 9d13600f2..abe128cd8 100644 --- a/src/modules/swiftest_globals.f90 +++ b/src/modules/swiftest_globals.f90 @@ -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 diff --git a/src/netcdf/netcdf.f90 b/src/netcdf/netcdf.f90 index 0adf0f7c6..e012ac389 100644 --- a/src/netcdf/netcdf.f90 +++ b/src/netcdf/netcdf.f90 @@ -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 @@ -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 @@ -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 @@ -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 @@ -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 @@ -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() @@ -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) @@ -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) diff --git a/src/orbel/orbel.f90 b/src/orbel/orbel.f90 index 428320cb4..e35d1e20f 100644 --- a/src/orbel/orbel.f90 +++ b/src/orbel/orbel.f90 @@ -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 diff --git a/src/whm/whm_setup.f90 b/src/whm/whm_setup.f90 index 5c672abd0..fb48dc1cc 100644 --- a/src/whm/whm_setup.f90 +++ b/src/whm/whm_setup.f90 @@ -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