From a18b6aba19bffce5cc6d7833d0b6d570b364666e Mon Sep 17 00:00:00 2001 From: Carlisle Wishard Date: Fri, 28 Oct 2022 17:28:25 -0400 Subject: [PATCH 1/3] Revert "deleted tides directory" This reverts commit f9db199bd3efa5bd3e538aaad255f497d7597811. --- src/tides/tides_getacch_pl.f90 | 65 ++++++++++++++++ src/tides/tides_spin_step.f90 | 132 +++++++++++++++++++++++++++++++++ 2 files changed, 197 insertions(+) create mode 100644 src/tides/tides_getacch_pl.f90 create mode 100644 src/tides/tides_spin_step.f90 diff --git a/src/tides/tides_getacch_pl.f90 b/src/tides/tides_getacch_pl.f90 new file mode 100644 index 000000000..4feb76221 --- /dev/null +++ b/src/tides/tides_getacch_pl.f90 @@ -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 \ No newline at end of file diff --git a/src/tides/tides_spin_step.f90 b/src/tides/tides_spin_step.f90 new file mode 100644 index 000000000..576aff8d7 --- /dev/null +++ b/src/tides/tides_spin_step.f90 @@ -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 \ No newline at end of file From 9e667febd42c467740ec8291b2c207a18a4cef27 Mon Sep 17 00:00:00 2001 From: Carlisle Wishard Date: Mon, 31 Oct 2022 15:06:02 -0400 Subject: [PATCH 2/3] Revert "deleted tides files from CMakeLists.txt" This reverts commit 87dd80c7e822bfbef2195ea9945421891abd0ed7. --- src/CMakeLists.txt | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 858018a52..5ea24cc0d 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -61,6 +61,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 @@ -126,4 +128,4 @@ IF(WIN32) ELSE() SET(CMAKE_INSTALL_PREFIX /usr/local) ENDIF(WIN32) -INSTALL(TARGETS ${FOOEXE} RUNTIME DESTINATION bin) +INSTALL(TARGETS ${FOOEXE} RUNTIME DESTINATION bin) \ No newline at end of file From a0eb0cc088c247f7eb58f474e869c83adc874e86 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Fri, 4 Nov 2022 16:53:54 -0400 Subject: [PATCH 3/3] Added back the interfaces for the tide subroutines --- src/modules/swiftest_classes.f90 | 26 +++++++++++++------------- 1 file changed, 13 insertions(+), 13 deletions(-) diff --git a/src/modules/swiftest_classes.f90 b/src/modules/swiftest_classes.f90 index aace128bf..6c0ac2be3 100644 --- a/src/modules/swiftest_classes.f90 +++ b/src/modules/swiftest_classes.f90 @@ -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