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

Commit

Permalink
Merge branch 'debug'
Browse files Browse the repository at this point in the history
  • Loading branch information
@daminton
daminton committed Nov 4, 2022
2 parents 740ea4f + a0eb0cc commit 531a35d
Show file tree
Hide file tree
Showing 4 changed files with 213 additions and 14 deletions.
4 changes: 3 additions & 1 deletion src/CMakeLists.txt
View file
Original file line number Diff line number Diff line change
Expand Up @@ -70,6 +70,8 @@ SET(FOO_src
${SRC}/symba/symba_setup.f90
${SRC}/symba/symba_step.f90
${SRC}/symba/symba_util.f90
${SRC}/tides/tides_getacch_pl.f90
${SRC}/tides/tides_spin_step.f90
${SRC}/user/user_getacch.f90
${SRC}/util/util_append.f90
${SRC}/util/util_coord.f90
Expand Down Expand Up @@ -135,4 +137,4 @@ IF(WIN32)
ELSE()
SET(CMAKE_INSTALL_PREFIX /usr/local)
ENDIF(WIN32)
INSTALL(TARGETS ${FOOEXE} RUNTIME DESTINATION bin)
INSTALL(TARGETS ${FOOEXE} RUNTIME DESTINATION bin)
26 changes: 13 additions & 13 deletions src/modules/swiftest_classes.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1110,19 +1110,19 @@ module subroutine setup_tp(self, n, param)
end subroutine setup_tp

! TODO: Implement the tides model
! module subroutine tides_kick_getacch_pl(self, system)
! implicit none
! class(swiftest_pl), intent(inout) :: self !! Swiftest massive body object
! class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object
! end subroutine tides_kick_getacch_pl

! module subroutine tides_step_spin_system(self, param, t, dt)
! implicit none
! class(swiftest_nbody_system), intent(inout) :: self !! Swiftest nbody system object
! class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters
! real(DP), intent(in) :: t !! Simulation time
! real(DP), intent(in) :: dt !! Current stepsize
! end subroutine tides_step_spin_system
module subroutine tides_kick_getacch_pl(self, system)
implicit none
class(swiftest_pl), intent(inout) :: self !! Swiftest massive body object
class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object
end subroutine tides_kick_getacch_pl

module subroutine tides_step_spin_system(self, param, t, dt)
implicit none
class(swiftest_nbody_system), intent(inout) :: self !! Swiftest nbody system object
class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters
real(DP), intent(in) :: t !! Simulation time
real(DP), intent(in) :: dt !! Current stepsize
end subroutine tides_step_spin_system

module subroutine user_kick_getacch_body(self, system, param, t, lbeg)
implicit none
Expand Down
65 changes: 65 additions & 0 deletions src/tides/tides_getacch_pl.f90
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,65 @@
submodule(swiftest_classes) s_tides_kick_getacch
use swiftest
contains

module subroutine tides_kick_getacch_pl(self, system)
!! author: Jennifer L.L. Pouplin, Carlisle A. wishard, and David A. Minton
!!
!! Calculated tidal torques from central body to any planet and from any planet to central body
!! planet - planet interactions are considered negligable.
!! This is a constant time lag model.
!!
!! Adapted from Mercury-T code from Bolmont et al. (2015)
!!
!! Reference:
!! Bolmont, E., Raymond, S.N., Leconte, J., Hersant, F., Correia, A.C.M., 2015.
!! Mercury-T : A new code to study tidally evolving multi-planet systems.
!! Applications to Kepler-62. A&A 583, A116. https://doi.org/10.1051/0004-6361/201525909
implicit none
! Arguments
class(swiftest_pl), intent(inout) :: self !! Swiftest massive body object
class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object
! Internals
integer(I4B) :: i
real(DP) :: rmag, vmag
real(DP), dimension(NDIM) :: r_unit, v_unit, h_unit, theta_unit, theta_dot, F_T
real(DP) :: Ftr, Ptopl, Ptocb, r5cbterm, r5plterm

associate(pl => self, npl => self%nbody, cb => system%cb)
pl%atide(:,:) = 0.0_DP
cb%atide(:) = 0.0_DP
do i = 1, npl
rmag = norm2(pl%xh(:,i))
vmag = norm2(pl%vh(:,i))
r_unit(:) = pl%xh(:,i) / rmag
v_unit(:) = pl%vh(:,i) / vmag
h_unit(:) = r_unit(:) .cross. v_unit(:)
theta_unit(:) = h_unit(:) .cross. r_unit(:)
theta_dot = dot_product(pl%vh(:,i), theta_unit(:))

! First calculate the tangential component of the force vector (eq. 5 & 6 of Bolmont et al. 2015)
! The radial component is already computed in the obl_acc methods
r5cbterm = pl%Gmass(i)**2 * cb%k2 * cb%radius**5
r5plterm = cb%Gmass**2 * pl%k2(i) * pl%radius(i)**5

Ptopl = 3 * r5plterm * pl%tlag(i) / rmag**7
Ptocb = 3 * r5cbterm * cb%tlag / rmag**7

Ftr = -3 / rmag**7 * (r5cbterm + r5plterm) - 3 * vmag / rmag * (Ptocb + Ptopl)

F_T(:) = (Ftr + (Ptocb + Ptopl) * dot_product(v_unit, r_unit) / rmag) * r_unit(:) &
+ Ptopl * ((pl%rot(:,i) - theta_dot(:)) .cross. r_unit(:)) &
+ Ptocb * ((cb%rot(:) - theta_dot(:)) .cross. r_unit(:))
cb%atide(:) = cb%atide(:) + F_T(:) / cb%Gmass
pl%atide(:,i) = F_T(:) / pl%Gmass(i)
end do

do i = 1, npl
pl%ah(:,i) = pl%ah(:,i) + pl%atide(:,i) + cb%atide(:)
end do
end associate

return
end subroutine tides_kick_getacch_pl

end submodule s_tides_kick_getacch
132 changes: 132 additions & 0 deletions src/tides/tides_spin_step.f90
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,132 @@
submodule(swiftest_classes) s_tides_step_spin
use swiftest

type, extends(lambda_obj_tvar) :: tides_derivs_func
!! Base class for an lambda function object. This object takes no additional arguments other than the dependent variable x, an array of real numbers
procedure(tidederiv), pointer, nopass :: lambdaptr_tides_deriv
real(DP), dimension(:,:), allocatable :: xbeg
real(DP), dimension(:,:), allocatable :: xend
real(DP) :: dt
contains
generic :: init => tides_derivs_init
procedure :: evalt => tides_derivs_eval
procedure, nopass :: tides_derivs_init
end type
interface lambda_obj
module procedure tides_derivs_init
end interface
abstract interface
function tidederiv(x, t, dt, xbeg, xend) result(y)
! Template for a 0 argument function
import DP, swiftest_nbody_system
real(DP), dimension(:), intent(in) :: x
real(DP), intent(in) :: t
real(DP), intent(in) :: dt
real(DP), dimension(:,:), intent(in) :: xbeg
real(DP), dimension(:,:), intent(in) :: xend
real(DP), dimension(:), allocatable :: y
end function
end interface

contains

module subroutine tides_step_spin_system(self, param, t, dt)
!! author: Jennifer L.L. Pouplin and David A. Minton
!!
!! Integrates the spin equations for central and massive bodies of the system subjected to tides.
implicit none
! Arguments
class(swiftest_nbody_system), intent(inout) :: self !! Swiftest nbody system object
class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters
real(DP), intent(in) :: t !! Simulation time
real(DP), intent(in) :: dt !! Current stepsize
! Internals
real(DP), dimension(:), allocatable :: rot0, rot1
real(DP) :: subt
real(DP), parameter :: tol=1e-6_DP !! Just a guess at the moment
real(DP) :: subdt

associate(pl => self%pl, npl => self%pl%nbody, cb => self%cb)
allocate(rot0(NDIM*(npl+1)))
rot0 = [pack(pl%rot(:,1:npl),.true.), pack(cb%rot(:),.true.)]
! Use this space call the ode_solver, passing tides_spin_derivs as the function:
subdt = dt / 20._DP
!rot1(:) = util_solve_rkf45(lambda_obj(tides_spin_derivs, subdt, pl%xbeg, pl%xend), rot0, dt, subdt tol)
! Recover with unpack
!pl%rot(:,1:npl) = unpack(rot1...
!cb%rot(:) = unpack(rot1...
end associate

return
end subroutine tides_step_spin_system


function tides_spin_derivs(rot_pl_cb, t, dt, xbeg, xend) result(drot) !! Need to add more arguments so we can pull in mass, radius, Ip, J2, etc...
!! author: Jennifer L.L. Pouplin and David A. Minton
!!
!! function used to calculate the derivatives that are fed to the ODE solver
implicit none
! Arguments
real(DP), dimension(:,:), intent(in) :: rot_pl_cb !! Array of rotations. The last element is the central body, and all others are massive bodies
real(DP), intent(in) :: t !! Current time, which is used to interpolate the massive body positions
real(DP), intent(in) :: dt !! Total step size
real(DP), dimension(:,:), intent(in) :: xbeg
real(DP), dimension(:,:), intent(in) :: xend
! Internals
real(DP), dimension(:,:), allocatable :: drot
real(DP), dimension(:), allocatable :: flatrot
real(DP), dimension(NDIM) :: N_Tcb, N_Rcb, N_Tpl, N_Rpl, xinterp
real(DP) :: C_cb, C_pl, r_dot_rot_cb, r_dot_rot_pl, rmag
integer(I4B) :: i, n


n = size(rot_pl_cb,2)
if (allocated(drot)) deallocate(drot)
allocate(drot, mold=rot_pl_cb)
drot(:,:) = 0.0_DP
do i = 1,n-1
xinterp(:) = xbeg(:,i) + t / dt * (xend(:,i) - xbeg(:,i))
! Calculate Ncb and Npl as a function of xinterp
!drot(:,i) = -Mcb / (Mcb + Mpl(i)) * (N_Tpl + N_Rpl)
!drot(:,n) = drot(:,n) - Mcb / (Mcb + Mpl(i) * (N_Tcb + N_Rcb)
!
end do

return
end function tides_spin_derivs

function tides_derivs_eval(self, x, t) result(y)
implicit none
! Arguments
class(tides_derivs_func), intent(inout) :: self
real(DP), dimension(:), intent(in) :: x
real(DP), intent(in) :: t
! Result
real(DP), dimension(:), allocatable :: y
if (associated(self%lambdaptr_tides_deriv)) then
y = self%lambdaptr_tides_deriv(x, t, self%dt, self%xbeg, self%xend)
else
error stop "Lambda function was not initialized"
end if

return
end function tides_derivs_eval

function tides_derivs_init(lambda, dt, xbeg, xend) result(f)
implicit none
! Arguments
procedure(tidederiv) :: lambda
real(DP), intent(in) :: dt
real(DP), dimension(:,:), intent(in) :: xbeg
real(DP), dimension(:,:), intent(in) :: xend
! Result
type(tides_derivs_func) :: f
f%lambdaptr_tides_deriv => lambda
f%dt = dt
allocate(f%xbeg, source = xbeg)
allocate(f%xend, source = xend)

return
end function tides_derivs_init

end submodule s_tides_step_spin

0 comments on commit 531a35d

Please sign in to comment.