From 1b61f96c23e8a557cb51ff82ec9647efdca0f254 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Wed, 21 Dec 2022 09:49:00 -0500 Subject: [PATCH] More cleanup --- src/CMakeLists.txt | 4 +- src/base/base_module.f90 | 3 + src/collision/collision_util.f90 | 40 ++++---- src/swiftest/swiftest_driver.f90 | 4 +- src/swiftest/swiftest_io.f90 | 4 +- src/swiftest/swiftest_module.f90 | 23 +---- src/swiftest/swiftest_setup.f90 | 4 +- .../swiftest_user.f90} | 16 ++-- src/swiftest/swiftest_util.f90 | 22 +++++ src/tides/tides_getacch_pl.f90 | 64 +++++++------ src/tides/tides_module.f90 | 92 +++++++++++++++++++ src/tides/tides_spin_step.f90 | 58 ++++-------- 12 files changed, 208 insertions(+), 126 deletions(-) rename src/{user/user_getacch.f90 => swiftest/swiftest_user.f90} (59%) create mode 100644 src/tides/tides_module.f90 diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index aaec591d1..33a8f5536 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -18,6 +18,7 @@ SET(STRICT_MATH_FILES ${SRC}/rmvs/rmvs_kick.f90 ${SRC}/symba/symba_kick.f90 ${SRC}/whm/whm_kick.f90 + ${SRC}/swiftest/swiftest_user.f90 ) SET(FAST_MATH_FILES @@ -82,9 +83,6 @@ SET(FAST_MATH_FILES ${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}/walltime/walltime_implementations.f90 ${SRC}/whm/whm_coord.f90 ${SRC}/whm/whm_drift.f90 diff --git a/src/base/base_module.f90 b/src/base/base_module.f90 index ea28a5e18..96e529c9d 100644 --- a/src/base/base_module.f90 +++ b/src/base/base_module.f90 @@ -292,6 +292,7 @@ end subroutine abstract_io_netcdf_initialize_output final :: final_storage_frame end type + type, abstract :: base_storage(nframes) !! An class that establishes the pattern for various storage objects integer(I4B), len :: nframes = 4096 !! Total number of frames that can be stored @@ -321,11 +322,13 @@ end subroutine abstract_io_netcdf_initialize_output type, abstract :: base_object end type base_object + type, abstract :: base_multibody(nbody) integer(I4B), len :: nbody integer(I4B), dimension(nbody) :: id end type base_multibody + !> Class definition for the kinship relationships used in bookkeeping multiple collisions bodies in a single time step. type, abstract :: base_kinship end type base_kinship diff --git a/src/collision/collision_util.f90 b/src/collision/collision_util.f90 index 02eaafdd3..a7506c5db 100644 --- a/src/collision/collision_util.f90 +++ b/src/collision/collision_util.f90 @@ -75,6 +75,9 @@ module subroutine collision_util_construct_temporary_system(self, nbody_system, ! Internals logical, dimension(:), allocatable :: linclude integer(I4B) :: npl_tot + ! The following are needed in order to deal with typing requirements + class(swiftest_nbody_system), allocatable :: tmpsys_local + class(swiftest_parameters), allocatable :: tmpparam_local select type(nbody_system) class is (swiftest_nbody_system) @@ -85,30 +88,27 @@ module subroutine collision_util_construct_temporary_system(self, nbody_system, if (allocated(tmpparam)) deallocate(tmpparam) if (allocated(tmpsys)) deallocate(tmpsys) allocate(tmpparam, source=param) - call swiftest_setup_construct_system(tmpsys, tmpparam) - select type(tmpsys) - class is (swiftest_nbody_system) - select type(tmpparam) - class is (swiftest_parameters) + call swiftest_setup_construct_system(tmpsys_local, tmpparam_local) - ! No test particles necessary for energy/momentum calcs - call tmpsys%tp%setup(0, param) + ! No test particles necessary for energy/momentum calcs + call tmpsys_local%tp%setup(0, tmpparam_local) - ! Replace the empty central body object with a copy of the original - deallocate(tmpsys%cb) - allocate(tmpsys%cb, source=cb) + ! Replace the empty central body object with a copy of the original + deallocate(tmpsys_local%cb) + allocate(tmpsys_local%cb, source=cb) - ! Make space for the fragments - npl_tot = npl + nfrag - call tmpsys%pl%setup(npl_tot, tmpparam) - allocate(linclude(npl_tot)) + ! Make space for the fragments + npl_tot = npl + nfrag + call tmpsys_local%pl%setup(npl_tot, tmpparam_local) + allocate(linclude(npl_tot)) - ! Fill up the temporary system with all of the original bodies, leaving the spaces for fragments empty until we add them in later - linclude(1:npl) = .true. - linclude(npl+1:npl_tot) = .false. - call tmpsys%pl%fill(pl, linclude) - end select - end select + ! Fill up the temporary system with all of the original bodies, leaving the spaces for fragments empty until we add them in later + linclude(1:npl) = .true. + linclude(npl+1:npl_tot) = .false. + call tmpsys_local%pl%fill(pl, linclude) + + call move_alloc(tmpsys_local, tmpsys) + call move_alloc(tmpparam_local, tmpparam) end associate end select diff --git a/src/swiftest/swiftest_driver.f90 b/src/swiftest/swiftest_driver.f90 index 6b64fcdb2..b581cd666 100644 --- a/src/swiftest/swiftest_driver.f90 +++ b/src/swiftest/swiftest_driver.f90 @@ -41,7 +41,7 @@ program swiftest_driver character(*), parameter :: symbacompactfmt = '(";NPLM",ES22.15,$)' !type(base_storage(nframes=:)), allocatable :: system_history - call io_get_args(integrator, param_file_name, display_style) + call swiftest_io_get_args(integrator, param_file_name, display_style) !> Read in the user-defined parameters file and the initial conditions of the system allocate(swiftest_parameters :: param) @@ -71,7 +71,7 @@ program swiftest_driver if (dump_cadence == 0) dump_cadence = ceiling(nloops / (1.0_DP * istep_out), kind=I8B) ! Construct the main n-body system using the user-input integrator to choose the type of system - call setup_construct_system(system, param) + call swiftest_setup_construct_system(system, param) !> Define the maximum number of threads nthreads = 1 ! In the *serial* case diff --git a/src/swiftest/swiftest_io.f90 b/src/swiftest/swiftest_io.f90 index 581e7e001..cf2fe02d6 100644 --- a/src/swiftest/swiftest_io.f90 +++ b/src/swiftest/swiftest_io.f90 @@ -310,7 +310,7 @@ module subroutine swiftest_io_dump_storage(self, param) end subroutine swiftest_io_dump_storage - module subroutine swiftest_io_get_args(integrator, param_file_name, display_style) + module subroutine swiftest_swiftest_io_get_args(integrator, param_file_name, display_style) !! author: David A. Minton !! !! Reads in the name of the parameter file from command line arguments. @@ -381,7 +381,7 @@ module subroutine swiftest_io_get_args(integrator, param_file_name, display_styl end if return - end subroutine swiftest_io_get_args + end subroutine swiftest_swiftest_io_get_args module function swiftest_io_get_token(buffer, ifirst, ilast, ierr) result(token) diff --git a/src/swiftest/swiftest_module.f90 b/src/swiftest/swiftest_module.f90 index e7925944a..a7bae5336 100644 --- a/src/swiftest/swiftest_module.f90 +++ b/src/swiftest/swiftest_module.f90 @@ -585,12 +585,12 @@ module subroutine swiftest_io_dump_storage(self, param) class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters end subroutine swiftest_io_dump_storage - module subroutine swiftest_io_get_args(integrator, param_file_name, display_style) + module subroutine swiftest_swiftest_io_get_args(integrator, param_file_name, display_style) implicit none character(len=:), allocatable, intent(inout) :: integrator !! Symbolic code of the requested integrator character(len=:), allocatable, intent(inout) :: param_file_name !! Name of the input parameters file character(len=:), allocatable, intent(inout) :: display_style !! Style of the output display {"STANDARD", "COMPACT"}). Default is "STANDARD" - end subroutine swiftest_io_get_args + end subroutine swiftest_swiftest_io_get_args module function swiftest_io_get_token(buffer, ifirst, ilast, ierr) result(token) implicit none @@ -1005,8 +1005,8 @@ end subroutine swiftest_setup_body module subroutine swiftest_setup_construct_system(system, param) implicit none - class(base_nbody_system), allocatable, intent(inout) :: system !! Swiftest system object - class(base_parameters), intent(inout) :: param !! Current run configuration parameters + class(swiftest_nbody_system), allocatable, intent(inout) :: system !! Swiftest system object + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters end subroutine swiftest_setup_construct_system module subroutine swiftest_setup_initialize_particle_info_system(self, param) @@ -1035,21 +1035,6 @@ module subroutine swiftest_setup_tp(self, n, param) class(swiftest_parameters), intent(in) :: param !! Current run configuration parametersr end subroutine swiftest_setup_tp - ! TODO: Implement the tides model - module subroutine swiftest_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 swiftest_tides_kick_getacch_pl - - ! module subroutine swiftest_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 swiftest_tides_step_spin_system - module subroutine swiftest_user_kick_getacch_body(self, system, param, t, lbeg) implicit none class(swiftest_body), intent(inout) :: self !! Swiftest massive body particle data structure diff --git a/src/swiftest/swiftest_setup.f90 b/src/swiftest/swiftest_setup.f90 index a69d62299..5429681c6 100644 --- a/src/swiftest/swiftest_setup.f90 +++ b/src/swiftest/swiftest_setup.f90 @@ -21,8 +21,8 @@ module subroutine swiftest_setup_construct_system(system, param) !! implicit none ! Arguments - class(base_nbody_system), allocatable, intent(inout) :: system !! Swiftest system object - class(base_parameters), intent(inout) :: param !! Swiftest parameters + class(swiftest_nbody_system), allocatable, intent(inout) :: system !! Swiftest system object + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters ! Internals type(encounter_storage) :: encounter_history type(collision_storage) :: collision_history diff --git a/src/user/user_getacch.f90 b/src/swiftest/swiftest_user.f90 similarity index 59% rename from src/user/user_getacch.f90 rename to src/swiftest/swiftest_user.f90 index deb54e93f..9d9c30783 100644 --- a/src/user/user_getacch.f90 +++ b/src/swiftest/swiftest_user.f90 @@ -7,10 +7,10 @@ !! You should have received a copy of the GNU General Public License along with Swiftest. !! If not, see: https://www.gnu.org/licenses. -submodule(base) s_user_kick_getacch +submodule(swiftest) s_user_kick_getacch use swiftest contains - module subroutine user_kick_getacch_body(self, system, param, t, lbeg) + module subroutine swiftest_user_kick_getacch_body(self, system, param, t, lbeg) !! author: David A. Minton !! !! Add user-supplied heliocentric accelerations to planets. @@ -18,13 +18,13 @@ module subroutine user_kick_getacch_body(self, system, param, t, lbeg) !! Adapted from David E. Kaufmann's Swifter routine whm_user_kick_getacch.f90 implicit none ! Arguments - class(swiftest_body), intent(inout) :: self !! Swiftest massive body particle data structure - class(swiftest_system), intent(inout) :: system !! Swiftest nbody_system_object - class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters user parameters - real(DP), intent(in) :: t !! Current time - logical, intent(in) :: lbeg !! Logical flag that determines whether or not this is the beginning or end of the ste + class(swiftest_body), intent(inout) :: self !! Swiftest massive body particle data structure + class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody_system_object + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters user parameters + real(DP), intent(in) :: t !! Current time + logical, intent(in) :: lbeg !! Logical flag that determines whether or not this is the beginning or end of the ste return - end subroutine user_kick_getacch_body + end subroutine swiftest_user_kick_getacch_body end submodule s_user_kick_getacch diff --git a/src/swiftest/swiftest_util.f90 b/src/swiftest/swiftest_util.f90 index 64868db02..7e303b214 100644 --- a/src/swiftest/swiftest_util.f90 +++ b/src/swiftest/swiftest_util.f90 @@ -2480,6 +2480,28 @@ module subroutine swiftest_util_rescale_system(self, param, mscale, dscale, tsca end subroutine swiftest_util_rescale_system + module subroutine swiftest_util_reset_kinship_pl(self, idx) + !! author: David A. Minton + !! + !! Resets the kinship status of bodies. + !! + implicit none + class(swiftest_pl), intent(inout) :: self !! SyMBA massive body object + integer(I4B), dimension(:), intent(in) :: idx !! Index array of bodies to reset + ! Internals + integer(I4B) :: i, j + + self%kin(idx(:))%parent = idx(:) + self%kin(idx(:))%nchild = 0 + do j = 1, size(idx(:)) + i = idx(j) + if (allocated(self%kin(i)%child)) deallocate(self%kin(i)%child) + end do + + return + end subroutine swiftest_util_reset_kinship_pl + + module subroutine swiftest_util_resize_arr_char_string(arr, nnew) !! author: David A. Minton !! diff --git a/src/tides/tides_getacch_pl.f90 b/src/tides/tides_getacch_pl.f90 index aed941e8f..680c22704 100644 --- a/src/tides/tides_getacch_pl.f90 +++ b/src/tides/tides_getacch_pl.f90 @@ -1,4 +1,4 @@ -submodule(base) s_tides_kick_getacch +submodule(tides) s_tides_kick_getacch use swiftest contains @@ -25,39 +25,45 @@ module subroutine tides_kick_getacch_pl(self, system) 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%rh(:,i)) - vmag = norm2(pl%vh(:,i)) - r_unit(:) = pl%rh(:,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(:)) + select type(pl => self) + class is (swiftest_pl) + select type(system) + class is (swiftest_nbody_system) + associate(npl => pl%nbody, cb => system%cb) + pl%atide(:,:) = 0.0_DP + cb%atide(:) = 0.0_DP + do i = 1, npl + rmag = norm2(pl%rh(:,i)) + vmag = norm2(pl%vh(:,i)) + r_unit(:) = pl%rh(:,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 + ! 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 + 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) + 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 + 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 + do i = 1, npl + pl%ah(:,i) = pl%ah(:,i) + pl%atide(:,i) + cb%atide(:) + end do + end associate + end select + end select return end subroutine tides_kick_getacch_pl diff --git a/src/tides/tides_module.f90 b/src/tides/tides_module.f90 new file mode 100644 index 000000000..6d9c5f465 --- /dev/null +++ b/src/tides/tides_module.f90 @@ -0,0 +1,92 @@ +!! Copyright 2022 - David Minton, Carlisle Wishard, Jennifer Pouplin, Jake Elliott, & Dana Singh +!! This file is part of Swiftest. +!! Swiftest is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License +!! as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version. +!! Swiftest is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty +!! of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. +!! You should have received a copy of the GNU General Public License along with Swiftest. +!! If not, see: https://www.gnu.org/licenses. + +module tides + !! author: The Purdue Swiftest Team - David A. Minton, Carlisle A. Wishard, Jennifer L.L. Pouplin, and Jacob R. Elliott + !! + !! Definition of classes and methods used to determine close encounters + use base + use lambda_function + implicit none + public + + + 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 :: rbeg + real(DP), dimension(:,:), allocatable :: rend + 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, rbeg, rend) result(y) + ! Template for a 0 argument function + import DP, base_nbody_system + real(DP), dimension(:), intent(in) :: x + real(DP), intent(in) :: t + real(DP), intent(in) :: dt + real(DP), dimension(:,:), intent(in) :: rbeg + real(DP), dimension(:,:), intent(in) :: rend + real(DP), dimension(:), allocatable :: y + end function + end interface + + + interface + module subroutine tides_kick_getacch_pl(self, system) + implicit none + class(base_object), intent(inout) :: self !! Swiftest massive body object + class(base_nbody_system), intent(inout) :: system !! Swiftest nbody system object + end subroutine tides_kick_getacch_pl + + module function tides_derivs_init(lambda, dt, rbeg, rend) result(f) + implicit none + procedure(tidederiv) :: lambda + real(DP), intent(in) :: dt + real(DP), dimension(:,:), intent(in) :: rbeg + real(DP), dimension(:,:), intent(in) :: rend + type(tides_derivs_func) :: f + end function tides_derivs_init + + module function tides_derivs_eval(self, x, t) result(y) + class(tides_derivs_func), intent(inout) :: self + real(DP), dimension(:), intent(in) :: x + real(DP), intent(in) :: t + real(DP), dimension(:), allocatable :: y + end function tides_derivs_eval + + module function tides_spin_derivs(rot_pl_cb, t, dt, rbeg, rend) result(drot) !! Need to add more arguments so we can pull in mass, radius, Ip, J2, etc... + 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) :: rbeg + real(DP), dimension(:,:), intent(in) :: rend + real(DP), dimension(:,:), allocatable :: drot + end function tides_spin_derivs + + module subroutine tides_step_spin_system(self, param, t, dt) + implicit none + class(base_nbody_system), intent(inout) :: self !! Swiftest nbody system object + class(base_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 + + end interface + + +end module \ No newline at end of file diff --git a/src/tides/tides_spin_step.f90 b/src/tides/tides_spin_step.f90 index fd74e4a21..adf602b77 100644 --- a/src/tides/tides_spin_step.f90 +++ b/src/tides/tides_spin_step.f90 @@ -1,33 +1,6 @@ -submodule(base) s_tides_step_spin +submodule(tides) 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 :: rbeg - real(DP), dimension(:,:), allocatable :: rend - 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, rbeg, rend) result(y) - ! Template for a 0 argument function - import DP, base_nbody_system - real(DP), dimension(:), intent(in) :: x - real(DP), intent(in) :: t - real(DP), intent(in) :: dt - real(DP), dimension(:,:), intent(in) :: rbeg - real(DP), dimension(:,:), intent(in) :: rend - real(DP), dimension(:), allocatable :: y - end function - end interface - contains module subroutine tides_step_spin_system(self, param, t, dt) @@ -46,22 +19,25 @@ module subroutine tides_step_spin_system(self, param, t, dt) 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(:) = swiftest_util_solve_rkf45(lambda_obj(tides_spin_derivs, subdt, pl%rbeg, pl%rend), rot0, dt, subdt tol) - ! Recover with unpack - !pl%rot(:,1:npl) = unpack(rot1... - !cb%rot(:) = unpack(rot1... - end associate + select type(self) + class is (swiftest_nbody_system) + 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(:) = swiftest_util_solve_rkf45(lambda_obj(tides_spin_derivs, subdt, pl%rbeg, pl%rend), rot0, dt, subdt,tol) + ! ! Recover with unpack + ! pl%rot(:,1:npl) = unpack(rot1... + ! cb%rot(:) = unpack(rot1... + end associate + end select return end subroutine tides_step_spin_system - function tides_spin_derivs(rot_pl_cb, t, dt, rbeg, rend) result(drot) !! Need to add more arguments so we can pull in mass, radius, Ip, J2, etc... + module function tides_spin_derivs(rot_pl_cb, t, dt, rbeg, rend) 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 @@ -95,7 +71,7 @@ function tides_spin_derivs(rot_pl_cb, t, dt, rbeg, rend) result(drot) !! Need to return end function tides_spin_derivs - function tides_derivs_eval(self, x, t) result(y) + module function tides_derivs_eval(self, x, t) result(y) implicit none ! Arguments class(tides_derivs_func), intent(inout) :: self @@ -112,7 +88,7 @@ function tides_derivs_eval(self, x, t) result(y) return end function tides_derivs_eval - function tides_derivs_init(lambda, dt, rbeg, rend) result(f) + module function tides_derivs_init(lambda, dt, rbeg, rend) result(f) implicit none ! Arguments procedure(tidederiv) :: lambda