Fixed bugs related to restarts. Now save the maxid value to the param…

…eter dump so that duplicate ids don't get created. Also rearranged the dump output format so that runs that output to EL will still output to EL on restart instead of reverting to XV and screwing up the bin.dat.
MintonGroup · Aug 23, 2021 · ce023d3 · ce023d3
1 parent c3022d3
commit ce023d3
Show file tree

Hide file tree

Showing 4 changed files with 15 additions and 16 deletions.
diff --git a/src/io/io.f90 b/src/io/io.f90
@@ -189,14 +189,14 @@ module subroutine io_dump_system(self, param)
          dump_param%incbfile = trim(adjustl(DUMP_CB_FILE(idx))) 
          dump_param%inplfile = trim(adjustl(DUMP_PL_FILE(idx))) 
          dump_param%intpfile = trim(adjustl(DUMP_TP_FILE(idx)))
-         dump_param%out_form = XV
          dump_param%in_form  = XV
          dump_param%out_stat = 'APPEND'
          dump_param%in_type = REAL8_TYPE
          dump_param%T0 = param%t
 
          call dump_param%dump(param_file_name)
 
+         dump_param%out_form = XV
          call self%cb%dump(dump_param)
          call self%pl%dump(dump_param)
          call self%tp%dump(dump_param)
@@ -454,7 +454,7 @@ module subroutine io_param_reader(self, unit, iotype, v_list, iostat, iomsg)
                   call io_toupper(param_value)
                   param%out_stat = param_value
                case ("ISTEP_DUMP")
-                  read(param_value, *) param%istep_dump
+                  read(param_value, *, err = 667, iomsg = iomsg) param%istep_dump
                case ("CHK_CLOSE")
                   call io_toupper(param_value)
                   if (param_value == "YES" .or. param_value == 'T') param%lclose = .true.
@@ -551,6 +551,8 @@ module subroutine io_param_reader(self, unit, iotype, v_list, iostat, iomsg)
                   read(param_value, *, err = 667, iomsg = iomsg) param%Ecollisions
                case("EUNTRACKED")
                   read(param_value, *, err = 667, iomsg = iomsg) param%Euntracked
+               case ("MAXID")
+                  read(param_value, *, err = 667, iomsg = iomsg) param%maxid 
                case ("NPLMAX", "NTPMAX", "GMTINY", "MIN_GMFRAG", "PARTICLE_OUT", "FRAGMENTATION", "SEED", "YARKOVSKY", "YORP") ! Ignore SyMBA-specific, not-yet-implemented, or obsolete input parameters
                case default
                   write(iomsg,*) "Unknown parameter -> ",param_name
@@ -760,6 +762,7 @@ module subroutine io_param_writer(self, unit, iotype, v_list, iostat, iomsg)
          write(param_name, Afmt) "TP_IN"; write(param_value, Afmt) trim(adjustl(param%intpfile)); write(unit, Afmt, err = 667, iomsg = iomsg) adjustl(param_name), adjustl(param_value)
          write(param_name, Afmt) "IN_TYPE"; write(param_value, Afmt) trim(adjustl(param%in_type)); write(unit, Afmt, err = 667, iomsg = iomsg) adjustl(param_name), adjustl(param_value)
          write(param_name, Afmt) "IN_FORM"; write(param_value, Afmt) trim(adjustl(param%in_form)); write(unit, Afmt, err = 667, iomsg = iomsg) adjustl(param_name), adjustl(param_value)
+         if (param%istep_dump > 0) write(param_name, Afmt) "ISTEP_DUMP"; write(param_value, Ifmt) param%istep_dump; write(unit, Afmt, err = 667, iomsg = iomsg) adjustl(param_name), adjustl(param_value)
          if (param%istep_out > 0) then
             write(param_name, Afmt) "ISTEP_OUT"; write(param_value, Ifmt) param%istep_out; write(unit, Afmt, err = 667, iomsg = iomsg) adjustl(param_name), adjustl(param_value)
             write(param_name, Afmt) "BIN_OUT"; write(param_value, Afmt) trim(adjustl(param%outfile)); write(unit, Afmt, err = 667, iomsg = iomsg) adjustl(param_name), adjustl(param_value)
@@ -768,9 +771,6 @@ module subroutine io_param_writer(self, unit, iotype, v_list, iostat, iomsg)
             write(param_name, Afmt) "OUT_STAT"; write(param_value, Afmt) "APPEND"; write(unit, Afmt, err = 667, iomsg = iomsg) adjustl(param_name), adjustl(param_value)
          end if
          write(param_name, Afmt) "ENC_OUT"; write(param_value, Afmt) trim(adjustl(param%enc_out)); write(unit, Afmt, err = 667, iomsg = iomsg) adjustl(param_name), adjustl(param_value)
-         if (param%istep_dump > 0) then
-            write(param_name, Afmt) "ISTEP_DUMP"; write(param_value, Ifmt) param%istep_dump; write(unit, Afmt, err = 667, iomsg = iomsg) adjustl(param_name), adjustl(param_value)
-         end if
          write(param_name, Afmt) "CHK_RMIN"; write(param_value, Rfmt) param%rmin; write(unit, Afmt, err = 667, iomsg = iomsg) adjustl(param_name), adjustl(param_value)
          write(param_name, Afmt) "CHK_RMAX"; write(param_value, Rfmt) param%rmax; write(unit, Afmt, err = 667, iomsg = iomsg) adjustl(param_name), adjustl(param_value)
          write(param_name, Afmt) "CHK_EJECT"; write(param_value, Rfmt) param%rmaxu; write(unit, Afmt, err = 667, iomsg = iomsg) adjustl(param_name), adjustl(param_value)
@@ -810,6 +810,7 @@ module subroutine io_param_writer(self, unit, iotype, v_list, iostat, iomsg)
             write(param_name, Afmt) "EUNTRACKED"; write(param_value, Rfmt) param%Euntracked; write(unit, Afmt, err = 667, iomsg = iomsg) adjustl(param_name), adjustl(param_value)
          end if
          write(param_name, Afmt) "FIRSTKICK"; write(param_value, Lfmt) param%lfirstkick; write(unit, Afmt, err = 667, iomsg = iomsg) adjustl(param_name), adjustl(param_value)
+         write(param_name, Afmt) "MAXID"; write(param_value, Ifmt) param%maxid; write(unit, Afmt, err = 667, iomsg = iomsg) adjustl(param_name), adjustl(param_value)
 
          iostat = 0
          iomsg = "UDIO not implemented"

diff --git a/src/modules/swiftest_classes.f90 b/src/modules/swiftest_classes.f90
@@ -15,8 +15,7 @@ module swiftest_classes
    !>    Each paramter is initialized to a default values. 
    type :: swiftest_parameters
       integer(I4B)         :: integrator     = UNKNOWN_INTEGRATOR !! Symbolic name of the nbody integrator  used
-      integer(I4B)         :: nplmax         = -1                 !! Maximum allowed number of massive bodies
-      integer(I4B)         :: ntpmax         = -1                 !! Maximum allowed number of test particles
+      integer(I4B)         :: maxid = -1           !! The current maximum particle id number 
       real(DP)             :: t0             = -1.0_DP            !! Integration start time
       real(DP)             :: t              = -1.0_DP            !! Integration current time
       real(DP)             :: tstop          = -1.0_DP            !! Integration stop time
@@ -295,7 +294,6 @@ module swiftest_classes
       logical                                    :: lbeg                 !! True if this is the beginning of a step. This is used so that test particle steps can be calculated 
                                                                          !!    separately from massive bodies.  Massive body variables are saved at half steps, and passed to 
                                                                          !!    the test particles
-      integer(I4B)                               :: maxid = -1           !! The current maximum particle id number 
    contains
       !> Each integrator will have its own version of the step
       procedure(abstract_step_system), deferred :: step
@@ -1448,7 +1446,7 @@ end subroutine util_spill_tp
       module subroutine util_valid_id_system(self, param)
          implicit none
          class(swiftest_nbody_system), intent(inout) :: self  !! Swiftest nbody system object
-         class(swiftest_parameters),   intent(in)    :: param !! Current run configuration parameters
+         class(swiftest_parameters),   intent(inout) :: param !! Current run configuration parameters
       end subroutine util_valid_id_system
 
       module subroutine util_version()

diff --git a/src/symba/symba_collision.f90 b/src/symba/symba_collision.f90
@@ -68,7 +68,7 @@ module function symba_collision_casedisruption(system, param, family, x, v, mass
       ! Distribute any residual mass if there is any and set the radius
       m_frag(nfrag) = m_frag(nfrag) + (mtot - sum(m_frag(:)))
       rad_frag(:) = (3 * m_frag(:) / (4 * PI * avg_dens))**(1.0_DP / 3.0_DP)
-      id_frag(:) = [(i, i = system%maxid + 1, system%maxid + nfrag)]
+      id_frag(:) = [(i, i = param%maxid + 1, param%maxid + nfrag)]
 
       do i = 1, nfrag
          Ip_frag(:, i) = Ip_new(:)
@@ -169,7 +169,7 @@ module function symba_collision_casehitandrun(system, param, family, x, v, mass,
          m_frag(2:nfrag) = (mtot - m_frag(1)) / (nfrag - 1) 
          rad_frag(2:nfrag) = (3 * m_frag(2:nfrag) / (4 * PI * avg_dens))**(1.0_DP / 3.0_DP)
          m_frag(nfrag) = m_frag(nfrag) + (mtot - sum(m_frag(:)))
-         id_frag(1:nfrag) = [(i, i = system%maxid + 1, system%maxid + nfrag)]
+         id_frag(1:nfrag) = [(i, i = param%maxid + 1, param%maxid + nfrag)]
 
          do i = 1, nfrag
             Ip_frag(:, i) = Ip(:, jproj)
@@ -372,7 +372,7 @@ module function symba_collision_casesupercatastrophic(system, param, family, x,
       ! Distribute any residual mass if there is any and set the radius
       m_frag(nfrag) = m_frag(nfrag) + (mtot - sum(m_frag(:)))
       rad_frag(:) = (3 * m_frag(:) / (4 * PI * avg_dens))**(1.0_DP / 3.0_DP)
-      id_frag(:) = [(i, i = system%maxid + 1, system%maxid + nfrag)]
+      id_frag(:) = [(i, i = param%maxid + 1, param%maxid + nfrag)]
 
       do i = 1, nfrag
          Ip_frag(:, i) = Ip_new(:)
@@ -821,7 +821,7 @@ subroutine symba_collision_mergeaddsub(system, param, family, id_frag, Ip_frag,
       implicit none
       ! Arguments
       class(symba_nbody_system),       intent(inout) :: system           !! SyMBA nbody system object
-      class(symba_parameters),         intent(in)    :: param            !! Current run configuration parameters with SyMBA additions
+      class(symba_parameters),         intent(inout) :: param            !! Current run configuration parameters with SyMBA additions
       integer(I4B),    dimension(:),   intent(in)    :: family           !! List of indices of all bodies inovlved in the collision
       integer(I4B),    dimension(:),   intent(in)    :: id_frag          !! List of fragment ids
       real(DP),        dimension(:),   intent(in)    :: m_frag, rad_frag !! Distribution of fragment mass and radii
@@ -861,7 +861,7 @@ subroutine symba_collision_mergeaddsub(system, param, family, id_frag, Ip_frag,
 
                ! Copy over identification, information, and physical properties of the new bodies from the fragment list
                plnew%id(1:nfrag) = id_frag(1:nfrag) 
-               system%maxid = system%maxid + nfrag
+               param%maxid = param%maxid + nfrag
                plnew%xb(:, 1:nfrag) = xb_frag(:, 1:nfrag) 
                plnew%vb(:, 1:nfrag) = vb_frag(:, 1:nfrag)
                do i = 1, nfrag

diff --git a/src/util/util_valid.f90 b/src/util/util_valid.f90
@@ -12,7 +12,7 @@ module subroutine util_valid_id_system(self, param)
       implicit none
       ! Arguments
       class(swiftest_nbody_system), intent(inout) :: self   !! Swiftest nbody system object
-      class(swiftest_parameters),   intent(in)    :: param  !! Current run configuration parameters
+      class(swiftest_parameters),   intent(inout) :: param  !! Current run configuration parameters
       ! Internals
       integer(I4B)                  :: i
       integer(I4B), dimension(:), allocatable :: idarr
@@ -34,7 +34,7 @@ module subroutine util_valid_id_system(self, param)
                call util_exit(FAILURE)
             end if
          end do
-         self%maxid = maxval(idarr)
+         param%maxid = max(param%maxid, maxval(idarr))
       end associate
 
       return