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

Commit

Permalink
Restructured the code so that the spherical harmonics code is its own…
Browse files Browse the repository at this point in the history 
… small module. Also refactored loblatecb to lnon_spherical_cb and now select which accelration model to use based on whether cb%c_lm is allocated or not.
  • Loading branch information
@daminton
daminton committed Feb 27, 2024
1 parent 46433e7 commit dfa48e1
Show file tree
Hide file tree
Showing 12 changed files with 123 additions and 127 deletions.
2 changes: 1 addition & 1 deletion src/CMakeLists.txt
View file
Original file line number Diff line number Diff line change
Expand Up @@ -33,14 +33,14 @@ SET(STRICT_MATH_FILES
${SRC}/operator/operator_unit.f90
${SRC}/rmvs/rmvs_kick.f90
${SRC}/rmvs/rmvs_step.f90
${SRC}/shgrav/shgrav_accel.f90
${SRC}/swiftest/swiftest_drift.f90
${SRC}/swiftest/swiftest_gr.f90
${SRC}/swiftest/swiftest_io.f90
${SRC}/swiftest/swiftest_kick.f90
${SRC}/swiftest/swiftest_user.f90
${SRC}/swiftest/swiftest_obl.f90
${SRC}/swiftest/swiftest_orbel.f90
${SRC}/swiftest/swiftest_sph.f90
${SRC}/symba/symba_drift.f90
${SRC}/symba/symba_gr.f90
${SRC}/symba/symba_kick.f90
Expand Down
4 changes: 2 additions & 2 deletions src/base/base_module.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -152,7 +152,7 @@ module base
!! Turn on close encounters
logical :: lenergy = .false.
!! Track the total energy of the system
logical :: loblatecb = .false.
logical :: lnon_spherical_cb = .false.
!! Calculate acceleration from oblate central body (automatically turns true if nonzero J2, J4, or c_lm is input)
logical :: lrotation = .false.
!! Include rotation states of big bodies
Expand Down Expand Up @@ -2489,7 +2489,7 @@ subroutine base_coclone_param(self)
call coclone(self%lbig_discard)
call coclone(self%lclose)
call coclone(self%lenergy)
call coclone(self%loblatecb)
call coclone(self%lnon_spherical_cb)
call coclone(self%lrotation)
call coclone(self%ltides)
call coclone(self%E_orbit_orig)
Expand Down
6 changes: 3 additions & 3 deletions src/helio/helio_kick.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -30,8 +30,8 @@ module subroutine helio_kick_getacch_pl(self, nbody_system, param, t, lbeg)

associate(cb => nbody_system%cb, pl => self, npl => self%nbody)
call pl%accel_int(param)
if (param%loblatecb) then
call pl%accel_obl(nbody_system)
if (param%lnon_spherical_cb) then
call pl%accel_non_spherical_cb(nbody_system)
if (lbeg) then
cb%aoblbeg = cb%aobl
else
Expand Down Expand Up @@ -79,7 +79,7 @@ module subroutine helio_kick_getacch_tp(self, nbody_system, param, t, lbeg)
else
call tp%accel_int(param, pl%Gmass(1:npl), pl%rend(:,1:npl), npl)
end if
if (param%loblatecb) call tp%accel_obl(nbody_system)
if (param%lnon_spherical_cb) call tp%accel_non_spherical_cb(nbody_system)
if (param%lextra_force) call tp%accel_user(nbody_system, param, t, lbeg)
if (param%lgr) call tp%accel_gr(param)
end associate
Expand Down
5 changes: 2 additions & 3 deletions src/rmvs/rmvs_kick.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -52,8 +52,7 @@ module subroutine rmvs_kick_getacch_tp(self, nbody_system, param, t, lbeg)

! Temporarily turn off the heliocentric-dependent acceleration terms during an inner encounter using a copy of the parameter list with all of the heliocentric-specific acceleration terms turned off
allocate(param_planetocen, source=param)
param_planetocen%loblatecb = .false.
param_planetocen%lshgrav = .false.
param_planetocen%lnon_spherical_cb = .false.
param_planetocen%lextra_force = .false.
param_planetocen%lgr = .false.

Expand Down Expand Up @@ -91,7 +90,7 @@ module subroutine rmvs_kick_getacch_tp(self, nbody_system, param, t, lbeg)
cb%Gmass = tp%cb_heliocentric%Gmass

! If the heliocentric-specifc acceleration terms are requested, compute those now
if (param%loblatecb) call tp%accel_obl(system_planetocen)
if (param%lnon_spherical_cb) call tp%accel_non_spherical_cb(system_planetocen)
if (param%lextra_force) call tp%accel_user(system_planetocen, param, t, lbeg)
if (param%lgr) call tp%accel_gr(param)

Expand Down
20 changes: 10 additions & 10 deletions src/rmvs/rmvs_step.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -202,7 +202,7 @@ subroutine rmvs_step_out(cb, pl, tp, nbody_system, param, t, dt)
call tp%step(nbody_system, param, outer_time, dto)
tp%lfirst = lfirsttp
else
if (param%loblatecb) then
if (param%lnon_spherical_cb) then
call swiftest_obl_acc(npl, cb%Gmass, cb%j2rp2, cb%j4rp4, pl%rbeg, pl%lmask, pl%outer(outer_index-1)%aobl, cb%rot,&
pl%Gmass, cb%aoblbeg)
call swiftest_obl_acc(npl, cb%Gmass, cb%j2rp2, cb%j4rp4, pl%rend, pl%lmask, pl%outer(outer_index)%aobl, cb%rot, &
Expand Down Expand Up @@ -267,14 +267,14 @@ subroutine rmvs_interp_in(cb, pl, nbody_system, param, dt, outer_index)
xtmp(:, 1:npl) = pl%inner(0)%x(:, 1:npl)
vtmp(:, 1:npl) = pl%inner(0)%v(:, 1:npl)

if ((param%loblatecb) .or. (param%ltides)) then
if ((param%lnon_spherical_cb) .or. (param%ltides)) then
allocate(rh_original, source=pl%rh)
allocate(ah_original, source=pl%ah)
pl%rh(:, 1:npl) = xtmp(:, 1:npl) ! Temporarily replace heliocentric position with inner substep values to calculate the
! oblateness terms
end if
if (param%loblatecb) then
call pl%accel_obl(nbody_system)
if (param%lnon_spherical_cb) then
call pl%accel_non_spherical_cb(nbody_system)
pl%inner(0)%aobl(:, 1:npl) = pl%aobl(:, 1:npl) ! Save the oblateness acceleration on the planet for this substep
end if
! TODO: Implement tides
Expand Down Expand Up @@ -328,9 +328,9 @@ subroutine rmvs_interp_in(cb, pl, nbody_system, param, dt, outer_index)
pl%inner(inner_index)%x(:, 1:npl) = pl%inner(inner_index)%x(:, 1:npl) + frac * xtmp(:, 1:npl)
pl%inner(inner_index)%v(:, 1:npl) = pl%inner(inner_index)%v(:, 1:npl) + frac * vtmp(:, 1:npl)

if (param%loblatecb) then
if (param%lnon_spherical_cb) then
pl%rh(:,1:npl) = pl%inner(inner_index)%x(:, 1:npl)
call pl%accel_obl(nbody_system)
call pl%accel_non_spherical_cb(nbody_system)
pl%inner(inner_index)%aobl(:, 1:npl) = pl%aobl(:, 1:npl)
end if
! TODO: Implement tides
Expand All @@ -339,10 +339,10 @@ subroutine rmvs_interp_in(cb, pl, nbody_system, param, dt, outer_index)
! pl%inner(inner_index)%atide(:, 1:npl) = pl%atide(:, 1:npl)
! end if
end do
if (param%loblatecb) then
if (param%lnon_spherical_cb) then
! Calculate the final value of oblateness accelerations at the final inner substep
pl%rh(:, 1:npl) = pl%inner(NTPHENC)%x(:, 1:npl)
call pl%accel_obl(nbody_system)
call pl%accel_non_spherical_cb(nbody_system)
pl%inner(NTPHENC)%aobl(:, 1:npl) = pl%aobl(:, 1:npl)
end if
! TODO: Implement tides
Expand Down Expand Up @@ -405,7 +405,7 @@ subroutine rmvs_step_in(cb, pl, tp, param, outer_time, dto)
call plenci%set_beg_end(rbeg = plenci%inner(inner_index - 1)%x, &
rend = plenci%inner(inner_index)%x)

if (param%loblatecb) then
if (param%lnon_spherical_cb) then
cbenci%aoblbeg = cbenci%inner(inner_index - 1)%aobl(:, 1)
cbenci%aoblend = cbenci%inner(inner_index )%aobl(:, 1)
end if
Expand Down Expand Up @@ -495,7 +495,7 @@ subroutine rmvs_make_planetocentric(param, cb, pl, tp)
plenci%inner(inner_index)%x(:,1) = -cbenci%inner(inner_index)%x(:,1)
plenci%inner(inner_index)%v(:,1) = -cbenci%inner(inner_index)%v(:,1)

if (param%loblatecb) then
if (param%lnon_spherical_cb) then
allocate(plenci%inner(inner_index)%aobl, mold=pl%inner(inner_index)%aobl)
allocate(cbenci%inner(inner_index)%aobl(NDIM,1))
cbenci%inner(inner_index)%aobl(:,1) = pl%inner(inner_index)%aobl(:, i)
Expand Down
148 changes: 69 additions & 79 deletions src/swiftest/swiftest_sph.f90 → src/shgrav/shgrav_accel.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -10,41 +10,58 @@

!! Swiftest submodule to calculate higher order terms for gravitational acceleration given spherical harmonic coefficients (c_lm)

submodule (swiftest) s_swiftest_sph
use operators
submodule (shgrav) s_shgrav_accel
use swiftest
use SHTOOLS

contains

module subroutine swiftest_sph_g_acc_one(GMcb, r_0, phi_cb, rh, c_lm, g_sph, GMpl, aoblcb)
subroutine shgrav_g_acc_one(GMcb, r_0, phi_cb, rh, c_lm, g_sph, GMpl, aoblcb)
!! author: Kaustub P. Anand
!!
!! Calculate the acceleration terms for one pair of bodies given c_lm, theta, phi, r
!!

implicit none
! Arguments
real(DP), intent(in) :: GMcb !! GMass of the central body
real(DP), intent(in) :: r_0 !! radius of the central body
real(DP), intent(in) :: phi_cb !! rotation phase angle of the central body
real(DP), intent(in), dimension(:) :: rh !! distance vector of body
real(DP), intent(in), dimension(:, :, :) :: c_lm !! Spherical Harmonic coefficients
real(DP), intent(out), dimension(NDIM) :: g_sph !! acceleration vector
real(DP), intent(in), optional :: GMpl !! Mass of input body if it is not a test particle
real(DP), dimension(:), intent(inout), optional :: aoblcb!! Barycentric acceleration of central body (only for massive input bodies)
real(DP), intent(in) :: GMcb
!! GMass of the central body
real(DP), intent(in) :: r_0
!! radius of the central body
real(DP), intent(in) :: phi_cb
!! rotation phase angle of the central body
real(DP), intent(in), dimension(:) :: rh
!! distance vector of body
real(DP), intent(in), dimension(:, :, :) :: c_lm
!! Spherical Harmonic coefficients
real(DP), intent(out), dimension(NDIM) :: g_sph
!! acceleration vector
real(DP), intent(in), optional :: GMpl
!! Mass of input body if it is not a test particle
real(DP), dimension(:), intent(inout), optional :: aoblcb
!! Barycentric acceleration of central body (only for massive input bodies)

! Internals
integer :: l, m !! SPH coefficients
integer :: l_max !! max Spherical Harmonic l order value
integer(I4B) :: N, lmindex !! Length of Legendre polynomials and index at a given l, m
real(DP) :: r_mag !! magnitude of rh
real(DP) :: phi, phi_bar !! Azimuthal/Phase angle (radians) wrt coordinate axes, and central body rotation phase
real(DP) :: theta !! Inclination/Zenith angle (radians)
real(DP) :: plm, plm1 !! Associated Legendre polynomials at a given l, m
real(DP) :: ccss, cssc !! See definition in source code
real(DP) :: cos_theta, sin_theta !! cos(theta) and sin(theta)
real(DP), dimension(:), allocatable :: p !! Associated Lengendre Polynomials at a given cos(theta)
real(DP) :: fac1, fac2, r_fac !! calculation factors
integer :: l, m
!! SPH coefficients
integer :: l_max
!! max Spherical Harmonic l order value
integer(I4B) :: N, lmindex
!! Length of Legendre polynomials and index at a given l, m
real(DP) :: r_mag
!! magnitude of rh
real(DP) :: phi, phi_bar
!! Azimuthal/Phase angle (radians) wrt coordinate axes, and central body rotation phase
real(DP) :: theta
!! Inclination/Zenith angle (radians)
real(DP) :: plm, plm1
!! Associated Legendre polynomials at a given l, m
real(DP) :: ccss, cssc
!! See definition in source code
real(DP) :: cos_theta, sin_theta
!! cos(theta) and sin(theta)
real(DP), dimension(:), allocatable :: p
!! Associated Lengendre Polynomials at a given cos(theta)
real(DP) :: fac1, fac2, r_fac
!! calculation factors

g_sph(:) = 0.0_DP
theta = atan2(sqrt(rh(1)**2 + rh(2)**2), rh(3))
Expand Down Expand Up @@ -117,70 +134,43 @@ module subroutine swiftest_sph_g_acc_one(GMcb, r_0, phi_cb, rh, c_lm, g_sph, GMp
end if

return
end subroutine swiftest_sph_g_acc_one
end subroutine shgrav_g_acc_one

module subroutine swiftest_sph_g_acc_pl_all(self, nbody_system)
module subroutine shgrav_acc(body, nbody_system)
!! author: Kaustub P. Anand
!!
!! Calculate the acceleration terms for all massive bodies given c_lm
!! Calculate the acceleration terms for bodies given c_lm values for the central body
!!
implicit none
! Arguments
class(swiftest_pl), intent(inout) :: self !! Swiftest massive body object
class(swiftest_nbody_system), intent(inout) :: nbody_system !! Swiftest nbody system object
class(swiftest_body), intent(inout) :: body
!! Swiftest body object
class(swiftest_nbody_system), intent(inout) :: nbody_system
!! Swiftest nbody system object
! Internals
integer(I4B) :: i = 1
real(DP), dimension(NDIM) :: g_sph !! Gravitational terms from Spherical Harmonics
integer(I4B) :: i
real(DP), dimension(NDIM) :: g_sph
!! Gravitational terms from Spherical Harmonics

associate(pl => self, npl => self%nbody, cb => nbody_system%cb, rh => self%rh)
associate(cb => nbody_system%cb)
cb%aobl(:) = 0.0_DP

do i = 1, npl
if (pl%lmask(i)) then
! theta = atan2(sqrt(rh(1,i)**2 + rh(2,i)**2), rh(3,i))
! phi = atan2(rh(2,i), rh(1,i)) ! - cb%rotphase

call swiftest_sph_g_acc_one(cb%Gmass, cb%radius, cb%rotphase, rh(:,i), cb%c_lm, g_sph, pl%Gmass(i), cb%aobl)
pl%ah(:, i) = pl%ah(:, i) + g_sph(:) - cb%aobl(:)
pl%aobl(:, i) = g_sph(:)
end if
end do
select type(body)
class is (swiftest_pl)
do i = 1, body%nbody
if (body%lmask(i)) then
call shgrav_g_acc_one(cb%Gmass, cb%radius, cb%rotphase, body%rh(:,i), cb%c_lm, body%aobl, &
GMpl=body%Gmass(i), aoblcb=cb%aobl)
end if
end do
class is (swiftest_tp)
do i = 1, body%nbody
if (body%lmask(i)) then
call shgrav_g_acc_one(cb%Gmass, cb%radius, cb%rotphase, body%rh(:,i), cb%c_lm, body%aobl)
end if
end do
end select
end associate
return
end subroutine swiftest_sph_g_acc_pl_all

module subroutine swiftest_sph_g_acc_tp_all(self, nbody_system)
!! author: Kaustub P. Anand
!!
!! Calculate the acceleration terms for all test particles given c_lm
!!
implicit none
! Arguments
class(swiftest_tp), intent(inout) :: self !! Swiftest test particle object
class(swiftest_nbody_system), intent(inout) :: nbody_system !! Swiftest nbody system object
! Internals
integer(I4B) :: i = 1
real(DP), dimension(NDIM) :: g_sph !! Gravitational terms from Spherical Harmonics
real(DP), dimension(NDIM) :: aoblcb !! Temporary variable for central body oblateness acceleration

associate(tp => self, ntp => self%nbody, cb => nbody_system%cb, rh => self%rh)

if (nbody_system%lbeg) then
aoblcb = cb%aoblbeg
else
aoblcb = cb%aoblend
end if

do i = 1, ntp
if (tp%lmask(i)) then

call swiftest_sph_g_acc_one(cb%Gmass, cb%radius, cb%rotphase, rh(:,i), cb%c_lm, g_sph)
tp%ah(:, i) = tp%ah(:, i) + g_sph(:) - aoblcb(:)
tp%aobl(:, i) = g_sph(:)
end if
end do
end associate
return
end subroutine swiftest_sph_g_acc_tp_all
end subroutine shgrav_acc

end submodule s_swiftest_sph
end submodule s_shgrav_accel
11 changes: 11 additions & 0 deletions src/shgrav/shgrav_module.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -16,5 +16,16 @@ module shgrav
implicit none
public

interface
module subroutine shgrav_acc(body, nbody_system)
implicit none
class(swiftest_body), intent(inout) :: body
!! Swiftest body object
class(swiftest_nbody_system), intent(inout) :: nbody_system
!! Swiftest nbody system object
end subroutine shgrav_acc

end interface

end module shgrav

6 changes: 2 additions & 4 deletions src/swiftest/swiftest_io.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3370,10 +3370,8 @@ module subroutine swiftest_io_read_in_system(self, nc, param)
if (ierr /=0) call base_util_exit(FAILURE,param%display_unit)
end if

param%lshgrav = allocated(self%cb%c_lm)

param%loblatecb = ((self%cb%j2rp2 /= 0.0_DP) .or. (self%cb%j4rp4 /= 0.0_DP)) .and. (.not. param%lshgrav)
if (.not.param%loblatecb .and. .not.param%lshgrav) then
param%lnon_spherical_cb = (self%cb%j2rp2 /= 0.0_DP) .or. (self%cb%j4rp4 /= 0.0_DP) .or. allocated(self%cb%c_lm)
if (.not.param%lnon_spherical_cb) then
if (allocated(self%pl%aobl)) deallocate(self%pl%aobl)
if (allocated(self%tp%aobl)) deallocate(self%tp%aobl)
else
Expand Down
Loading

0 comments on commit dfa48e1

Please sign in to comment.