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

Commit

Permalink
Merge branch '12-fix-problem-causing-the-grav-harmonics-dimensions-to…
Browse files Browse the repository at this point in the history 
…-be-created-in-the-output-file-even-when-not-being-used' into oblate_coord_rot
  • Loading branch information
@daminton
daminton committed Mar 1, 2024
2 parents 787017c + 849e2a7 commit deac0c8
Show file tree
Hide file tree
Showing 2 changed files with 59 additions and 38 deletions.
78 changes: 52 additions & 26 deletions src/swiftest/swiftest_gr.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -23,9 +23,12 @@ pure module subroutine swiftest_gr_kick_getaccb_ns_body(self, nbody_system, para
!! Adapted from David A. Minton's Swifter routine routine gr_kick_getaccb_ns.f90
implicit none
! Arguments
class(swiftest_body), intent(inout) :: self !! Swiftest generic body object
class(swiftest_nbody_system), intent(inout) :: nbody_system !! Swiftest nbody system object
class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters
class(swiftest_body), intent(inout) :: self
!! Swiftest generic body object
class(swiftest_nbody_system), intent(inout) :: nbody_system
!! Swiftest nbody system object
class(swiftest_parameters), intent(in) :: param
!! Current run configuration parameters
! Internals
real(DP) :: rmag, rdotv, vmag2
integer(I4B) :: i
Expand Down Expand Up @@ -62,12 +65,18 @@ pure module subroutine swiftest_gr_kick_getacch(mu, x, lmask, n, inv_c2, agr)
!! Adapted from David A. Minton's Swifter routine routine gr_whm_kick_getacch.f90
implicit none
! Arguments
real(DP), dimension(:), intent(in) :: mu !! Gravitational constant
real(DP), dimension(:,:), intent(in) :: x !! Position vectors
logical, dimension(:), intent(in) :: lmask !! Logical mask indicating which bodies to compute
integer(I4B), intent(in) :: n !! Total number of bodies
real(DP), intent(in) :: inv_c2 !! Inverse speed of light squared: 1 / c**2
real(DP), dimension(:,:), intent(out) :: agr !! Accelerations
real(DP), dimension(:), intent(in) :: mu
!! Gravitational constant
real(DP), dimension(:,:), intent(in) :: x
!! Position vectors
logical, dimension(:), intent(in) :: lmask
!! Logical mask indicating which bodies to compute
integer(I4B), intent(in) :: n
!! Total number of bodies
real(DP), intent(in) :: inv_c2
!! Inverse speed of light squared: 1 / c**2
real(DP), dimension(:,:), intent(out) :: agr
!! Accelerations
! Internals
integer(I4B) :: i
real(DP) :: beta, rjmag4
Expand Down Expand Up @@ -100,10 +109,14 @@ pure elemental module subroutine swiftest_gr_p4_pos_kick(inv_c2, rx, ry, rz, vx,
!! Adapted from David A. Minton's Swifter routine gr_whm_p4.f90
implicit none
! Arguments
real(DP), intent(in) :: inv_c2 !! One over speed of light squared (1/c**2)
real(DP), intent(inout) :: rx, ry, rz !! Position vector
real(DP), intent(in) :: vx, vy, vz !! Velocity vector
real(DP), intent(in) :: dt !! Step size
real(DP), intent(in) :: inv_c2
!! One over speed of light squared (1/c**2)
real(DP), intent(inout) :: rx, ry, rz
!! Position vector
real(DP), intent(in) :: vx, vy, vz
!! Velocity vector
real(DP), intent(in) :: dt
!! Step size
! Internals
real(DP) :: drx, dry, drz
real(DP) :: vmag2
Expand Down Expand Up @@ -133,11 +146,16 @@ pure module subroutine swiftest_gr_pseudovel2vel(param, mu, rh, pv, vh)
!! Adapted from David A. Minton's Swifter routine gr_pseudovel2vel.f90
implicit none
! Arguments
class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters
real(DP), intent(in) :: mu !! G * (Mcb + m), G = gravitational constant, Mcb = mass of central body, m = mass of body
real(DP), dimension(:), intent(in) :: rh !! Heliocentric position vector
real(DP), dimension(:), intent(in) :: pv !! Pseudovelocity velocity vector - see Saha & Tremain (1994), eq. (32)
real(DP), dimension(:), intent(out) :: vh !! Heliocentric velocity vector
class(swiftest_parameters), intent(in) :: param
!! Current run configuration parameters
real(DP), intent(in) :: mu
!! G * (Mcb + m), G = gravitational constant, Mcb = mass of central body, m = mass of body
real(DP), dimension(:), intent(in) :: rh
!! Heliocentric position vector
real(DP), dimension(:), intent(in) :: pv
!! Pseudovelocity velocity vector - see Saha & Tremain (1994), eq. (32)
real(DP), dimension(:), intent(out) :: vh
!! Heliocentric velocity vector
! Internals
real(DP) :: vmag2, rmag, grterm

Expand Down Expand Up @@ -191,11 +209,16 @@ pure module subroutine swiftest_gr_vel2pseudovel(param, mu, rh, vh, pv)
!! Adapted from David A. Minton's Swifter routine gr_vel2pseudovel.f90
implicit none
! Arguments
class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters
real(DP), intent(in) :: mu !! G * (Mcb + m), G = gravitational constant, Mcb = mass of central body, m = mass of body
real(DP), dimension(:), intent(in) :: rh !! Heliocentric position vector
real(DP), dimension(:), intent(in) :: vh !! Heliocentric velocity vector
real(DP), dimension(:), intent(out) :: pv !! Pseudovelocity vector - see Saha & Tremain (1994), eq. (32)
class(swiftest_parameters), intent(in) :: param
!! Current run configuration parameters
real(DP), intent(in) :: mu
!! G * (Mcb + m), G = gravitational constant, Mcb = mass of central body, m = mass of body
real(DP), dimension(:), intent(in) :: rh
!! Heliocentric position vector
real(DP), dimension(:), intent(in) :: vh
!! Heliocentric velocity vector
real(DP), dimension(:), intent(out) :: pv
!! Pseudovelocity vector - see Saha & Tremain (1994), eq. (32)
! Internals
real(DP) :: v2, G, pv2, rterm, det
real(DP), dimension(NDIM,NDIM) :: J,Jinv
Expand Down Expand Up @@ -264,11 +287,14 @@ pure module subroutine swiftest_gr_vh2pv_body(self, param)
!! Wrapper function that converts from heliocentric velocity to pseudovelocity for Swiftest bodies
implicit none
! Arguments
class(swiftest_body), intent(inout) :: self !! Swiftest particle object
class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters
class(swiftest_body), intent(inout) :: self
!! Swiftest particle object
class(swiftest_parameters), intent(in) :: param
!! Current run configuration parameters
! Internals
integer(I4B) :: i
real(DP), dimension(:,:), allocatable :: pv !! Temporary holder of pseudovelocity for in-place conversion
real(DP), dimension(:,:), allocatable :: pv
!! Temporary holder of pseudovelocity for in-place conversion

associate(n => self%nbody)
if (n == 0) return
Expand Down
19 changes: 7 additions & 12 deletions src/swiftest/swiftest_io.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1072,14 +1072,13 @@ module subroutine swiftest_io_netcdf_initialize_output(self, param)
"netcdf_io_initialize_output nf90_put_var space" )

if (param%lrotation .and. nc%lc_lm_exists) then

! Populate coordinate values for l and m and export to hdf file
call netcdf_io_check( nf90_put_var(nc%id, nc%l_varid, [(i, i=0, nc%l_dim_max-1)]), &
"netcdf_io_write_frame_cb nf90_put_var l_varid")
call netcdf_io_check( nf90_put_var(nc%id, nc%m_varid, [(i, i=0, nc%m_dim_max-1)]), &
"netcdf_io_write_frame_cb nf90_put_var m_varid")
call netcdf_io_check( nf90_put_var(nc%id, nc%sign_varid, [1,-1]), &
"netcdf_io_write_frame_cb nf90_put_var sign_varid")
! Populate coordinate values for l, m, and sign
call netcdf_io_check( nf90_put_var(nc%id, nc%l_varid, [(i, i=0, nc%l_dim_max-1)], start=[1], count=[nc%l_dim_max]), &
"netcdf_io_netcdf_initialize_output nf90_put_var l_varid")
call netcdf_io_check( nf90_put_var(nc%id, nc%m_varid, [(i, i=0, nc%m_dim_max-1)], start=[1], count=[nc%m_dim_max]), &
"netcdf_io_netcdf_initialize_output nf90_put_var m_varid")
call netcdf_io_check( nf90_put_var(nc%id, nc%sign_varid, nc%sign_coords, start=[1], count=[2] ), &
"netcdf_io_netcdf_initialize_output nf90_put_var sign_varid")
end if

end associate
Expand Down Expand Up @@ -1617,10 +1616,6 @@ module function swiftest_io_netcdf_read_frame_system(self, nc, param) result(ier
! if (npl > 0) pl%Q(:) = pack(rtemp, plmask)
! end if





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

if (param%in_form == "EL") then
Expand Down

0 comments on commit deac0c8

Please sign in to comment.