Merge branch 'debug' into DocumentationTesting

Makefile.Defines

@@ -65,13 +65,13 @@ GPAR = -fopenmp -ftree-parallelize-loops=4
 GMEM = -fsanitize=undefined -fsanitize=address -fsanitize=leak 
 GWARNINGS = -Wall -Warray-bounds  -Wimplicit-interface -Wextra -Warray-temporaries
 
-FFLAGS  = $(IDEBUG) $(HEAPARR)
-#FFLAGS  = -init=snan,arrays -no-wrap-margin -O3 $(STRICTREAL) $(SIMDVEC) $(PAR) 
+#FFLAGS  = $(IDEBUG) $(HEAPARR)
+FFLAGS  = -init=snan,arrays -no-wrap-margin -O3 $(STRICTREAL) $(SIMDVEC) $(PAR) 
 FORTRAN = ifort
 #AR     = xiar
 
 #FORTRAN = gfortran
-#FFLAGS =  -ffree-line-length-none $(GDEBUG) $(GMEM)
+#FFLAGS =  -ffree-line-length-none $(GDEBUG) #$(GMEM)
 AR     = ar
 
 # DO NOT include in CFLAGS the "-c" option to compile object only

docs/src/fragmentation.f90

@@ -369,7 +369,7 @@ subroutine calculate_system_energy(linclude_fragments)
             class is (symba_pl)
                select type(param)
                class is (symba_parameters)
-                  plwksp%nplm = count(plwksp%Gmass > param%mtiny / mscale)
+                  plwksp%nplm = count(plwksp%Gmass > param%Gmtiny / mscale)
                end select
             end select
             call tmpsys%pl%eucl_index()
@@ -381,7 +381,7 @@ subroutine calculate_system_energy(linclude_fragments)
             class is (symba_pl)
                select type(param)
                class is (symba_parameters)
-                  nplm = count(pl%mass > param%mtiny)
+                  nplm = count(pl%mass > param%Gmtiny)
                end select
             end select
             if (lk_plpl) call pl%eucl_index()
@@ -836,7 +836,7 @@ end subroutine fragmentation_initialize
 
 
 
-   module subroutine fragmentation_regime(Mcb, m1, m2, rad1, rad2, xh1, xh2, vb1, vb2, den1, den2, regime, Mlr, Mslr, mtiny, Qloss)
+   module subroutine fragmentation_regime(Mcb, m1, m2, rad1, rad2, xh1, xh2, vb1, vb2, den1, den2, regime, Mlr, Mslr, Gmtiny, Qloss)
       !! Author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton
       !!
       !! Determine the collisional regime of two colliding bodies. 
@@ -857,7 +857,7 @@ module subroutine fragmentation_regime(Mcb, m1, m2, rad1, rad2, xh1, xh2, vb1, v
       ! Arguments
       integer(I4B), intent(out)         :: regime
       real(DP), intent(out)          :: Mlr, Mslr
-      real(DP), intent(in)           :: Mcb, m1, m2, rad1, rad2, den1, den2, mtiny 
+      real(DP), intent(in)           :: Mcb, m1, m2, rad1, rad2, den1, den2, Gmtiny 
       real(DP), dimension(:), intent(in)   :: xh1, xh2, vb1, vb2
       real(DP), intent(out)          :: Qloss !! The residual energy after the collision 
       ! Constants
@@ -931,7 +931,7 @@ module subroutine fragmentation_regime(Mcb, m1, m2, rad1, rad2, xh1, xh2, vb1, v
       Qloss = 0.0_DP
       U_binding = (3.0_DP * Mtot) / (5.0_DP * Rp) ! LS12 eq. 27
 
-      if ((m1 < mtiny).or.(m2 < mtiny)) then 
+      if ((m1 < Gmtiny).or.(m2 < Gmtiny)) then 
          regime = COLLRESOLVE_REGIME_MERGE !perfect merging regime
          Mlr = Mtot
          Mslr = 0.0_DP

docs/src/io.f90

@@ -493,7 +493,7 @@ module subroutine io_param_reader(self, unit, iotype, v_list, iostat, iomsg)
                   read(param_value, *) param%Ecollisions
                case("EUNTRACKED")
                   read(param_value, *) param%Euntracked
-               case ("NPLMAX", "NTPMAX", "MTINY", "PARTICLE_FILE", "FRAGMENTATION", "SEED", "YARKOVSKY", "YORP") ! Ignore SyMBA-specific, not-yet-implemented, or obsolete input parameters
+               case ("NPLMAX", "NTPMAX", "GMTINY", "PARTICLE_FILE", "FRAGMENTATION", "SEED", "YARKOVSKY", "YORP") ! Ignore SyMBA-specific, not-yet-implemented, or obsolete input parameters
                case default
                   write(iomsg,*) "Unknown parameter -> ",param_name
                   iostat = -1
@@ -1211,16 +1211,20 @@ module subroutine io_write_discard(self, param)
 
       associate(tp_discards => self%tp_discards, nsp => self%tp_discards%nbody, pl => self%pl, npl => self%pl%nbody)
          if (nsp == 0) return
-         select case(param%out_stat)
-         case('APPEND')
+         if (lfirst) then
+            select case(param%out_stat)
+            case('APPEND')
+               open(unit = LUN, file = param%discard_out, status = 'OLD', position = 'APPEND', form = 'FORMATTED', iostat = ierr)
+            case('NEW', 'REPLACE', 'UNKNOWN')
+               open(unit = LUN, file = param%discard_out, status = param%out_stat, form = 'FORMATTED', iostat = ierr)
+            case default
+               write(*,*) 'Invalid status code for OUT_STAT: ',trim(adjustl(param%out_stat))
+               call util_exit(FAILURE)
+            end select
+            lfirst = .false.
+         else
             open(unit = LUN, file = param%discard_out, status = 'OLD', position = 'APPEND', form = 'FORMATTED', iostat = ierr)
-         case('NEW', 'REPLACE', 'UNKNOWN')
-            open(unit = LUN, file = param%discard_out, status = param%out_stat, form = 'FORMATTED', iostat = ierr)
-         case default
-            write(*,*) 'Invalid status code for OUT_STAT: ',trim(adjustl(param%out_stat))
-            call util_exit(FAILURE)
-         end select
-         lfirst = .false.
+         end if
          if (param%lgr) call tp_discards%pv2v(param) 
 
          write(LUN, HDRFMT) param%t, nsp, param%lbig_discard

docs/src/setup.f90

@@ -54,8 +54,8 @@ module subroutine setup_construct_system(system, param)
             allocate(symba_pl :: system%pl)
             allocate(symba_tp :: system%tp)
             allocate(symba_tp :: system%tp_discards)
-            allocate(symba_merger :: system%mergeadd_list)
-            allocate(symba_merger :: system%mergesub_list)
+            allocate(symba_merger :: system%pl_adds)
+            allocate(symba_merger :: system%pl_discards)
             allocate(symba_plplenc :: system%plplenc_list)
             allocate(symba_pltpenc :: system%pltpenc_list)
          end select

docs/src/swiftest_classes.f90

@@ -469,11 +469,11 @@ module subroutine fragmentation_initialize(system, param, family, x, v, L_spin,
          real(DP), intent(inout)                 :: Qloss
       end subroutine fragmentation_initialize
 
-      module subroutine fragmentation_regime(Mcb, m1, m2, rad1, rad2, xh1, xh2, vb1, vb2, den1, den2, regime, Mlr, Mslr, mtiny, Qloss)
+      module subroutine fragmentation_regime(Mcb, m1, m2, rad1, rad2, xh1, xh2, vb1, vb2, den1, den2, regime, Mlr, Mslr, Gmtiny, Qloss)
          implicit none
          integer(I4B), intent(out)         :: regime
          real(DP), intent(out)          :: Mlr, Mslr
-         real(DP), intent(in)           :: Mcb, m1, m2, rad1, rad2, den1, den2, mtiny 
+         real(DP), intent(in)           :: Mcb, m1, m2, rad1, rad2, den1, den2, Gmtiny 
          real(DP), dimension(:), intent(in)   :: xh1, xh2, vb1, vb2
          real(DP), intent(out)          :: Qloss !! The residual energy after the collision 
       end subroutine fragmentation_regime

docs/src/symba_classes.f90

@@ -19,7 +19,7 @@ module symba_classes
 
    type, extends(swiftest_parameters) :: symba_parameters
       character(STRMAX)                       :: particle_file  = PARTICLE_OUTFILE !! Name of output particle information file
-      real(DP)                                :: MTINY          = -1.0_DP          !! Smallest mass that is fully gravitating
+      real(DP)                                :: GMTINY          = -1.0_DP          !! Smallest mass that is fully gravitating
       integer(I4B), dimension(:), allocatable :: seed                              !! Random seeds
       logical                                 :: lfragmentation = .false.          !! Do fragmentation modeling instead of simple merger.
    contains
@@ -73,9 +73,9 @@ module symba_classes
    type, extends(helio_pl) :: symba_pl
       logical,                   dimension(:), allocatable :: lcollision !! flag indicating whether body has merged with another this time step
       logical,                   dimension(:), allocatable :: lencounter !! flag indicating whether body is part of an encounter this time step
-      logical,                   dimension(:), allocatable :: lmtiny     !! flag indicating whether this body is below the MTINY cutoff value
-      integer(I4B)                                         :: nplm       !! number of bodies above the MTINY limit
-      integer(I8B)                                         :: nplplm     !! Number of body (all massive)-body (only those above MTINY) comparisons in the flattened upper triangular matrix 
+      logical,                   dimension(:), allocatable :: lmtiny     !! flag indicating whether this body is below the GMTINY cutoff value
+      integer(I4B)                                         :: nplm       !! number of bodies above the GMTINY limit
+      integer(I8B)                                         :: nplplm     !! Number of body (all massive)-body (only those above GMTINY) comparisons in the flattened upper triangular matrix 
       integer(I4B),              dimension(:), allocatable :: nplenc     !! number of encounters with other planets this time step
       integer(I4B),              dimension(:), allocatable :: ntpenc     !! number of encounters with test particles this time step
       integer(I4B),              dimension(:), allocatable :: levelg     !! level at which this body should be moved
@@ -160,8 +160,8 @@ module symba_classes
    !  symba_nbody_system class definitions and method interfaces
    !********************************************************************************************************************************
    type, extends(helio_nbody_system) :: symba_nbody_system
-      class(symba_merger),  allocatable :: mergeadd_list !! List of added bodies in mergers or collisions
-      class(symba_merger),  allocatable :: mergesub_list !! List of subtracted bodies in mergers or collisions
+      class(symba_merger),  allocatable :: pl_adds !! List of added bodies in mergers or collisions
+      class(symba_merger),  allocatable :: pl_discards !! List of subtracted bodies in mergers or collisions
       class(symba_pltpenc), allocatable :: pltpenc_list  !! List of massive body-test particle encounters in a single step 
       class(symba_plplenc), allocatable :: plplenc_list  !! List of massive body-massive body encounters in a single step
       integer(I4B)                      :: irec          !! System recursion level

docs/src/symba_collision.f90

@@ -449,7 +449,7 @@ module subroutine symba_collision_resolve_fragmentations(self, system, param)
                v2_si(:)      = plpl_collisions%v2(:,i) * param%DU2M / param%TU2S    !! The velocity of the parent from inside the step (at collision)
                density_si(:) = mass_si(:) / (4.0_DP / 3._DP * PI * radius_si(:)**3) !! The collective density of the parent and its children
                Mcb_si        = cb%mass * param%MU2KG 
-               mtiny_si      = (param%MTINY / param%GU) * param%MU2KG
+               mtiny_si      = (param%GMTINY / param%GU) * param%MU2KG
 
                mass_res(:) = 0.0_DP
 

docs/src/symba_fragmentation.f90

@@ -34,7 +34,7 @@ module function symba_fragmentation_casedisruption(system, param, family, x, v,
 
       select type(pl => system%pl)
       class is (symba_pl)
-         associate(mergeadd_list => system%mergeadd_list, mergesub_list => system%mergesub_list, cb => system%cb)
+         associate(pl_adds => system%pl_adds, pl_discards => system%pl_discards, cb => system%cb)
             ! Collisional fragments will be uniformly distributed around the pre-impact barycenter
             nfrag = NFRAG_DISRUPT 
             allocate(m_frag(nfrag))
@@ -91,11 +91,11 @@ module function symba_fragmentation_casedisruption(system, param, family, x, v,
                lmask(:) = .false.
                lmask(family(:)) = .true.
                pl%status(family(:)) = MERGED
-               nstart = mergesub_list%nbody + 1
-               nend = mergesub_list%nbody + nfamily
-               call mergesub_list%append(pl, lmask)
+               nstart = pl_discards%nbody + 1
+               nend = pl_discards%nbody + nfamily
+               call pl_discards%append(pl, lmask)
                ! Record how many bodies were subtracted in this event
-               mergesub_list%ncomp(nstart:nend) = nfamily 
+               pl_discards%ncomp(nstart:nend) = nfamily 
 
                allocate(plnew, mold=pl)
                call plnew%setup(nfrag, param)
@@ -133,10 +133,10 @@ module function symba_fragmentation_casedisruption(system, param, family, x, v,
                end if
 
                ! Append the new merged body to the list and record how many we made
-               nstart = mergeadd_list%nbody + 1
-               nend = mergeadd_list%nbody + plnew%nbody
-               call mergeadd_list%append(plnew)
-               mergeadd_list%ncomp(nstart:nend) = plnew%nbody
+               nstart = pl_adds%nbody + 1
+               nend = pl_adds%nbody + plnew%nbody
+               call pl_adds%append(plnew)
+               pl_adds%ncomp(nstart:nend) = plnew%nbody
 
                call plnew%setup(0, param)
                deallocate(plnew)
@@ -179,7 +179,7 @@ module function symba_fragmentation_casehitandrun(system, param, family, x, v, m
 
       select type(pl => system%pl)
       class is (symba_pl)
-         associate(mergeadd_list => system%mergeadd_list, mergesub_list => system%mergesub_list, cb => system%cb)
+         associate(pl_adds => system%pl_adds, pl_discards => system%pl_discards, cb => system%cb)
             mtot = sum(mass(:))
             xcom(:) = (mass(1) * x(:,1) + mass(2) * x(:,2)) / mtot
             vcom(:) = (mass(1) * v(:,1) + mass(2) * v(:,2)) / mtot
@@ -247,11 +247,11 @@ module function symba_fragmentation_casehitandrun(system, param, family, x, v, m
                lmask(:) = .false.
                lmask(family(:)) = .true.
                pl%status(family(:)) = MERGED
-               nstart = mergesub_list%nbody + 1
-               nend = mergesub_list%nbody + nfamily
-               call mergesub_list%append(pl, lmask)
+               nstart = pl_discards%nbody + 1
+               nend = pl_discards%nbody + nfamily
+               call pl_discards%append(pl, lmask)
                ! Record how many bodies were subtracted in this event
-               mergesub_list%ncomp(nstart:nend) = nfamily 
+               pl_discards%ncomp(nstart:nend) = nfamily 
 
                allocate(plnew, mold=pl)
                call plnew%setup(nfrag, param)
@@ -289,10 +289,10 @@ module function symba_fragmentation_casehitandrun(system, param, family, x, v, m
                end if
 
                ! Append the new merged body to the list and record how many we made
-               nstart = mergeadd_list%nbody + 1
-               nend = mergeadd_list%nbody + plnew%nbody
-               call mergeadd_list%append(plnew)
-               mergeadd_list%ncomp(nstart:nend) = plnew%nbody
+               nstart = pl_adds%nbody + 1
+               nend = pl_adds%nbody + plnew%nbody
+               call pl_adds%append(plnew)
+               pl_adds%ncomp(nstart:nend) = plnew%nbody
 
                call plnew%setup(0, param)
                deallocate(plnew)
@@ -334,7 +334,7 @@ module function symba_fragmentation_casemerge(system, param, family, x, v, mass,
 
       select type(pl => system%pl)
       class is (symba_pl)
-         associate(mergeadd_list => system%mergeadd_list, mergesub_list => system%mergesub_list, cb => system%cb)
+         associate(pl_adds => system%pl_adds, pl_discards => system%pl_discards, cb => system%cb)
             status = MERGED
             write(*, '("Merging bodies ",99(I8,",",:))') pl%id(family(:))
             mass_new = sum(mass(:))
@@ -386,11 +386,11 @@ module function symba_fragmentation_casemerge(system, param, family, x, v, mass,
             lmask(:) = .false.
             lmask(family(:)) = .true.
             pl%status(family(:)) = MERGED
-            nstart = mergesub_list%nbody + 1
-            nend = mergesub_list%nbody + nfamily
-            call mergesub_list%append(pl, lmask)
+            nstart = pl_discards%nbody + 1
+            nend = pl_discards%nbody + nfamily
+            call pl_discards%append(pl, lmask)
             ! Record how many bodies were subtracted in this event
-            mergesub_list%ncomp(nstart:nend) = nfamily 
+            pl_discards%ncomp(nstart:nend) = nfamily 
 
             ! Create the new merged body 
             allocate(plnew, mold=pl)
@@ -422,10 +422,10 @@ module function symba_fragmentation_casemerge(system, param, family, x, v, mass,
             end if
 
             ! Append the new merged body to the list and record how many we made
-            nstart = mergeadd_list%nbody + 1
-            nend = mergeadd_list%nbody + plnew%nbody
-            call mergeadd_list%append(plnew)
-            mergeadd_list%ncomp(nstart:nend) = plnew%nbody
+            nstart = pl_adds%nbody + 1
+            nend = pl_adds%nbody + plnew%nbody
+            call pl_adds%append(plnew)
+            pl_adds%ncomp(nstart:nend) = plnew%nbody
 
             call plnew%setup(0, param)
             deallocate(plnew)
@@ -468,7 +468,7 @@ module function symba_fragmentation_casesupercatastrophic(system, param, family,
 
       select type(pl => system%pl)
       class is (symba_pl)
-         associate(mergeadd_list => system%mergeadd_list, mergesub_list => system%mergesub_list, cb => system%cb)
+         associate(pl_adds => system%pl_adds, pl_discards => system%pl_discards, cb => system%cb)
             ! Collisional fragments will be uniformly distributed around the pre-impact barycenter
             nfrag = NFRAG_SUPERCAT
             allocate(m_frag(nfrag))
@@ -521,11 +521,11 @@ module function symba_fragmentation_casesupercatastrophic(system, param, family,
                lmask(:) = .false.
                lmask(family(:)) = .true.
                pl%status(family(:)) = MERGED
-               nstart = mergesub_list%nbody + 1
-               nend = mergesub_list%nbody + nfamily
-               call mergesub_list%append(pl, lmask)
+               nstart = pl_discards%nbody + 1
+               nend = pl_discards%nbody + nfamily
+               call pl_discards%append(pl, lmask)
                ! Record how many bodies were subtracted in this event
-               mergesub_list%ncomp(nstart:nend) = nfamily 
+               pl_discards%ncomp(nstart:nend) = nfamily 
 
                allocate(plnew, mold=pl)
                call plnew%setup(nfrag, param)
@@ -563,10 +563,10 @@ module function symba_fragmentation_casesupercatastrophic(system, param, family,
                end if
 
                ! Append the new merged body to the list and record how many we made
-               nstart = mergeadd_list%nbody + 1
-               nend = mergeadd_list%nbody + plnew%nbody
-               call mergeadd_list%append(plnew)
-               mergeadd_list%ncomp(nstart:nend) = plnew%nbody
+               nstart = pl_adds%nbody + 1
+               nend = pl_adds%nbody + plnew%nbody
+               call pl_adds%append(plnew)
+               pl_adds%ncomp(nstart:nend) = plnew%nbody
 
                call plnew%setup(0, param)
                deallocate(plnew)