Fixed numerous issues that were prevening SyMBA runs from being prope…

…rly restarted such that the produce identical results to non-restarted runs.

examples/symba_test_restart/original/clement_2018.in

@@ -1,4 +1,4 @@
-1106      ! Solar System in unit system AU, M_sun, and years
+1106    ! Solar System in unit system AU, M_sun, and years
 1 3.94784176e+01
 .0 .0 .0        ! x y z
 .0 .0 .0        ! vx vy vz

examples/symba_test_restart/original/param.in

@@ -7,8 +7,8 @@ DT             0.016              ! stepsize in years
 PL_IN          clement_2018.in
 TP_IN          tp.in
 IN_TYPE        ASCII
-ISTEP_OUT      625               ! output cadence 
-ISTEP_DUMP     625               ! system dump cadence
+ISTEP_OUT      625             ! output cadence 
+ISTEP_DUMP     625             ! system dump cadence
 BIN_OUT        bin.dat
 PARTICLE_FILE  particle.dat
 OUT_TYPE       REAL8              ! double precision real output

src/io/io_dump_pl_symba.f90

@@ -0,0 +1,51 @@
+subroutine io_dump_pl_symba(npl, symba_plA, param)
+   !! Author: David A. Minton
+   !!
+   !! Dumps a SyMBA data structure to file.
+   !! 
+! modules
+   use swiftest
+   use module_symba
+   use module_interfaces, except_this_one => io_dump_pl_symba
+   implicit none
+
+! arguments
+   integer(I4B), intent(in)         :: npl
+   type(symba_pl), intent(inout) :: symba_plA
+   type(user_input_parameters),intent(inout) :: param
+
+! internals
+   integer(I4B)                     :: ierr
+   integer(I4B), save               :: idx = 1
+   integer(I4B),parameter           :: lun = 7
+
+   open(unit = lun, file = dump_pl_file(idx), form = "UNFORMATTED", status = 'REPLACE', iostat = ierr)
+   if (ierr /= 0) then
+      write(*, *) "Swiftest error:"
+      write(*, *) "   Unable to open binary dump file ", trim(dump_pl_file(idx))
+      call util_exit(FAILURE)
+   end if
+   write(lun) npl
+   write(lun) symba_plA%helio%swiftest%id(1:npl)
+   write(lun) symba_plA%helio%swiftest%mass(1:npl)
+   if (param%lrhill_present) write(lun) symba_plA%helio%swiftest%rhill(1:npl) 
+   if (param%lclose) write(lun) symba_plA%helio%swiftest%radius(1:npl) 
+   write(lun) symba_plA%helio%swiftest%xh(:,1:npl)
+   write(lun) symba_plA%helio%swiftest%vh(:,1:npl)
+   if (param%lrotation) then
+      write(lun) symba_plA%helio%swiftest%ip(:,1:npl)
+      write(lun) symba_plA%helio%swiftest%rot(:,1:npl)
+   end if
+   write(lun) symba_plA%helio%swiftest%vb(:,1:npl)
+   write(lun) symba_plA%helio%ah(:,1:npl)
+   write(lun) symba_plA%helio%ahi(:,1:npl)
+   close(lun)
+   if (idx == 1) then
+      idx = 2
+   else
+      idx = 1
+   end if
+
+   return
+
+end subroutine io_dump_pl_symba

src/main/swiftest_symba.f90

@@ -93,7 +93,7 @@ program swiftest_symba
       end if
 
       ! reads in initial conditions of all massive bodies from input file
-      call symba_plA%helio%swiftest%read_from_file(param)
+      call symba_read_pl_in(symba_plA, param) 
       call symba_tpA%helio%swiftest%read_from_file(param)
 
       !Temporary until the argument lists get fixed
@@ -106,7 +106,7 @@ program swiftest_symba
          call symba_tp_allocate(symba_tpA,ntpfile)
          call symba_plA%helio%swiftest%dealloc()
          call symba_tpA%helio%swiftest%dealloc()
-         call symba_plA%helio%swiftest%read_from_file(param)
+         call symba_read_pl_in(symba_plA, param)
          call symba_tpA%helio%swiftest%read_from_file(param)
       !**************************************************
       call io_write_particle_pl(symba_plA%helio%swiftest, [(i, i=1,npl)], param)
@@ -116,9 +116,9 @@ program swiftest_symba
       end if
 
       ! reorder by mass 
-      if (out_stat /= "APPEND") call symba_reorder_pl(npl, symba_plA)
+      if (out_stat /= "APPEND") call symba_reorder_pl(npl, symba_plA)  ! This is a new run, so we will sort the massive body list by mass
+
       call util_valid(npl, ntp, symba_plA%helio%swiftest, symba_tpA%helio%swiftest)
-      call util_hills(npl, symba_plA%helio%swiftest)
       ntp0 = ntp
       t = t0
       tbase = t0
@@ -155,6 +155,7 @@ program swiftest_symba
       start = clock_count / (count_rate * 1.0_DP)
       finish = start
       do while ((t < tstop) .and. ((ntp0 == 0) .or. (ntp > 0)))
+         call util_hills(npl, symba_plA%helio%swiftest)
          call symba_step_eucl(t, dt, param,npl,ntp,symba_plA, symba_tpA, nplplenc, npltpenc,&
                plplenc_list, pltpenc_list, nmergeadd, nmergesub, mergeadd_list, mergesub_list)
 
@@ -223,7 +224,7 @@ program swiftest_symba
                write(*,walltimefmt) finish - start, wallperstep
 
                call param%dump_to_file(t)
-               call io_dump_pl(npl, symba_plA%helio%swiftest, param)
+               call io_dump_pl_symba(npl, symba_plA, param)
                call io_dump_tp(ntp, symba_tpA%helio%swiftest)
                idump = istep_dump
             end if

src/modules/module_interfaces.f90

@@ -497,6 +497,16 @@ SUBROUTINE io_dump_pl(npl, swiftest_plA, param)
                type(user_input_parameters),intent(inout) :: param
 
           END SUBROUTINE io_dump_pl
+
+         subroutine io_dump_pl_symba(npl, symba_plA, param)
+            use swiftest_globals
+            use module_symba
+            use swiftest_data_structures
+            implicit none
+            integer(I4B), intent(in)        :: npl
+            type(symba_pl), intent(inout):: symba_plA
+            type(user_input_parameters),intent(inout) :: param
+         end subroutine io_dump_pl_symba
      END INTERFACE
 
      INTERFACE

src/modules/module_symba.f90

@@ -103,4 +103,11 @@ MODULE module_symba
 
    end type symba_merger 
 
+   interface
+      module subroutine symba_read_pl_in(symba_plA, param) 
+         type(symba_pl),              intent(inout) :: symba_plA  !! Swiftest data structure to store massive body initial conditions
+         type(user_input_parameters), intent(inout) :: param    !! Input collection of user-defined parameters
+      end subroutine symba_read_pl_in
+   end interface
+
 END MODULE module_symba

src/symba/symba_helio_getacch.f90

@@ -54,10 +54,8 @@ SUBROUTINE symba_helio_getacch(lflag, lextra_force, t, npl, nplm, helio_plA, j2r
 
 ! Executable code
      IF (lflag) THEN
-          DO i = 2, npl
-               helio_plA%ahi(:,i) = (/ 0.0_DP, 0.0_DP, 0.0_DP /)
-          END DO
-          CALL symba_helio_getacch_int(npl, nplm, helio_plA) 
+         helio_plA%ahi(:,2:npl) = 0.0_DP
+         CALL symba_helio_getacch_int(npl, nplm, helio_plA) 
      END IF
      IF (j2rp2 /= 0.0_DP) THEN
           DO i = 2, npl
@@ -70,9 +68,7 @@ SUBROUTINE symba_helio_getacch(lflag, lextra_force, t, npl, nplm, helio_plA, j2r
                helio_plA%ah(:,i) = helio_plA%ahi(:,i) + aobl(:, i) - aobl(:, 1)
           END DO
      ELSE
-          DO i = 2, npl
-               helio_plA%ah(:,i) = helio_plA%ahi(:,i)
-          END DO
+         helio_plA%ah(:,2:npl) = helio_plA%ahi(:,2:npl)
      END IF
      IF (lextra_force) CALL helio_user_getacch(t, npl, helio_plA) 
 

src/symba/symba_read_pl_in.f90

@@ -0,0 +1,136 @@
+submodule (module_symba) s_symba_read_pl_in
+contains
+   module subroutine symba_read_pl_in(symba_plA, param) 
+      !! author: The Purdue Swiftest Team -  David A. Minton, Carlisle A. Wishard, Jennifer L.L. Pouplin, and Jacob R. Elliott
+      !!
+      !! Read in massive body data 
+      !!
+      !! Adapted from David E. Kaufmann's Swifter routine symba_init_pl.f90
+      !! Adapted from Martin Duncan's Swift routine symba_init_pl.f
+      use swiftest
+      use module_interfaces
+      implicit none
+      ! Arguments
+      type(symba_pl),              intent(inout) :: symba_plA  !! Swiftest data structure to store massive body initial conditions
+      type(user_input_parameters), intent(inout) :: param    !! Input collection of user-defined parameters
+      ! Internals
+      integer(I4B), parameter :: LUN = 7              !! Unit number of input file
+      integer(I4B)            :: i, ierr, npl
+      logical                 :: is_ascii 
+
+      ierr = 0
+      is_ascii = (param%in_type == 'ASCII') 
+      if (is_ascii) then
+         open(unit = LUN, file = param%inplfile, status = 'old', form = 'formatted', iostat = ierr)
+      else
+         open(unit = LUN, file = param%inplfile, status = 'old', form = 'unformatted', iostat = ierr)
+      end if
+      if (ierr /=  0) then
+         write(*,*) 'Error opening massive body initial conditions file ',trim(adjustl(param%inplfile))
+         return
+      end if
+
+      if (is_ascii) then
+         read(LUN, *, iostat = ierr) npl
+      else
+         read(LUN, iostat = ierr) npl
+      end if
+      if (npl <= 0) return
+      call symba_plA%helio%swiftest%alloc(npl)
+
+      if (is_ascii) then
+         read(LUN, *, iostat = ierr) symba_plA%helio%swiftest%id(1), symba_plA%helio%swiftest%mass(1)
+         symba_plA%helio%swiftest%rhill(1) = 0.0_DP
+         symba_plA%helio%swiftest%radius(1) = param%rmin
+         read(LUN, *, iostat = ierr) symba_plA%helio%swiftest%xh(:,1)
+         read(LUN, *, iostat = ierr) symba_plA%helio%swiftest%vh(:,1)
+         if (param%lrotation) THEN
+            read(LUN, *, iostat = ierr) symba_plA%helio%swiftest%Ip(:,1)
+            read(LUN, *, iostat = ierr) symba_plA%helio%swiftest%rot(:,1)
+         end if
+         if (ierr /= 0) then
+            write(*,*) 'Error reading central body values in ',trim(adjustl(param%inplfile))
+            return
+         end if
+         do i = 1, NDIM
+            if ((symba_plA%helio%swiftest%xh(i,1) /= 0.0_DP) .or. (symba_plA%helio%swiftest%vh(i,1) /= 0.0_DP)) then
+               write(*, *) "Swiftest error:"
+               write(*, *) " Input must be in heliocentric coordinates."
+               write(*, *) " position/velocity components of body 1 are"
+               write(*, *) symba_plA%helio%swiftest%xh(:,1)
+               write(*, *) symba_plA%helio%swiftest%vh(:,1)
+            end if
+         end do
+         symba_plA%helio%swiftest%status(1) = ACTIVE
+         do i = 2, symba_plA%helio%swiftest%nbody
+            if (param%lrhill_present) then
+               read(LUN, *, iostat = ierr) symba_plA%helio%swiftest%id(i), symba_plA%helio%swiftest%mass(i), symba_plA%helio%swiftest%rhill(i)
+            else
+               read(LUN, *, iostat = ierr) symba_plA%helio%swiftest%id(i), symba_plA%helio%swiftest%mass(i)
+               symba_plA%helio%swiftest%rhill(i) = 0.0_DP
+            end if
+            if (ierr /= 0 ) exit
+            if (param%lclose) then
+               read(LUN, *, iostat = ierr) symba_plA%helio%swiftest%radius(i)
+               if (ierr /= 0 ) exit
+            else
+               symba_plA%helio%swiftest%radius(i) = 0.0_DP
+            end if
+            read(LUN, *, iostat = ierr) symba_plA%helio%swiftest%xh(:,i)
+            read(LUN, *, iostat = ierr) symba_plA%helio%swiftest%vh(:,i)
+            if (param%lrotation) THEN
+               read(LUN, *, iostat = ierr) symba_plA%helio%swiftest%Ip(:,i)
+               read(LUN, *, iostat = ierr) symba_plA%helio%swiftest%rot(:,i)
+            end if
+            if (ierr /= 0 ) exit
+            symba_plA%helio%swiftest%status(i) = ACTIVE
+         end do
+      else
+         read(LUN, iostat = ierr) symba_plA%helio%swiftest%id(:)
+         read(LUN, iostat = ierr) symba_plA%helio%swiftest%mass(:)
+         if (param%lrhill_present) then
+            read(LUN, iostat = ierr) symba_plA%helio%swiftest%rhill(:)
+         else
+            symba_plA%helio%swiftest%rhill(:) = 0.0_DP
+         end if
+         if (param%lclose) then
+            read(LUN, iostat = ierr) symba_plA%helio%swiftest%radius(:)
+         else
+            symba_plA%helio%swiftest%radius(:) = 0.0_DP
+         end if
+         read(LUN, iostat = ierr) symba_plA%helio%swiftest%xh(:,:)
+         read(LUN, iostat = ierr) symba_plA%helio%swiftest%vh(:,:)
+         if (param%lrotation) THEN
+            read(LUN, iostat = ierr) symba_plA%helio%swiftest%Ip(:,:)
+            read(LUN, iostat = ierr) symba_plA%helio%swiftest%rot(:,:)
+         end if
+         if (param%out_stat == 'APPEND') then
+            read(lun, iostat = ierr) symba_plA%helio%swiftest%vb(:,:)
+            read(lun, iostat = ierr) symba_plA%helio%ah(:,:)
+            read(lun, iostat = ierr) symba_plA%helio%ahi(:,:)
+         end if
+         if (ierr /= 0) then
+            write(*,*) 'An error occurred reading in ',trim(adjustl(param%inplfile))
+            call util_exit(FAILURE)
+         end if
+         symba_plA%helio%swiftest%status(:) = ACTIVE
+      end if
+      close(unit = LUN)
+
+      symba_plA%helio%swiftest%info(1)%origin_type = "Central body"
+      do i = 2, npl
+         symba_plA%helio%swiftest%info(i)%origin_xh(:) = symba_plA%helio%swiftest%xh(:,i)
+         symba_plA%helio%swiftest%info(i)%origin_vh(:) = symba_plA%helio%swiftest%vh(:,i)
+         symba_plA%helio%swiftest%info(i)%origin_time = 0.0_DP
+         symba_plA%helio%swiftest%info(i)%origin_type = "Initial conditions"
+      end do
+
+      ! Give massive bodies a positive id
+      symba_plA%helio%swiftest%id(:) = abs(symba_plA%helio%swiftest%id(:))
+      symba_plA%helio%swiftest%maxid = maxval(symba_plA%helio%swiftest%id(:))
+
+      return
+   end subroutine symba_read_pl_in
+
+end submodule s_symba_read_pl_in
+