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