More cleanup

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

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

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

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

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)

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

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

src/user/user_getacch.f90 → src/swiftest/swiftest_user.f90

@@ -7,24 +7,24 @@
 !! 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.
       !!
       !! 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

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
       !!

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