Merge branch 'master' into debug

Makefile.Defines

@@ -45,11 +45,6 @@ COLLRESOLVE_HOME = $(ROOT_DIR)/collresolve/
 
 # DO NOT include in FFLAGS the "-c" option to compile object only
 # this is done explicitly as needed in the Makefile
-ADVIXE_DIR = /apps/cent7/intel/advisor_2019
-ADVIXE_FLAGS = -g -O2 -qopt-report=5 -vecabi=cmdtarget -simd -shared-intel -debug inline-debug-info -DTBB_DEBUG -DTBB_USE_THREADING_TOOLS -xhost -traceback 
-
-VTUNE_FLAGS = -g -O2 -qopt-report=5 -simd -shared-intel -qopenmp -debug inline-debug-info -parallel-source-info=2 -parallel -DTBB_DEBUG -DTBB_USE_THREADING_TOOLS -qopenmp -fp-model no-except -mp1 -xhost -traceback 
-#Be sure to set the environment variable KMP_FORKJOIN_FRAMES=1 for OpenMP debuging in vtune
 
 IDEBUG = -O0 -init=snan,arrays  -nogen-interfaces -no-pie -no-ftz -fpe-all=0 -g -traceback -mp1 -fp-model strict -fpe0 -debug all -align all -pad -ip -prec-div -prec-sqrt -assume protect-parens -CB -no-wrap-margin
 STRICTREAL = -fp-model=precise -prec-div -prec-sqrt -assume protect-parens 
@@ -66,26 +61,27 @@ GMEM =  -fsanitize-address-use-after-scope -fstack-check  -fsanitize=bounds-stri
 GWARNINGS = -Wall -Warray-bounds  -Wimplicit-interface -Wextra -Warray-temporaries
 GPRODUCTION = -O3 -ffree-line-length-none $(GPAR)
 
-
 MKL_ROOT = /apps/spack/bell/apps/intel-parallel-studio/cluster.2019.5-intel-19.0.5-4brgqlf/mkl/lib
 IMKL = -I$(MKLROOT)/include
 LMKL = -L$(MKLROOT)/lib/intel64 -qopt-matmul 
 
-FSTRICTFLAGS  = $(IDEBUG) #$(SIMDVEC) $(PAR)
-FFLAGS  = $(IDEBUG) #$(SIMDVEC) $(PAR)
+FSTRICTFLAGS = $(IPRODUCTION) $(STRICTREAL) $(OPTREPORT) 
+FFLAGS = $(IPRODUCTION) -fp-model=fast $(OPTREPORT) 
 FORTRAN = ifort
 AR     = xiar
 CC    = icc
 
 #FORTRAN = gfortran
 #FFLAGS =  $(GDEBUG)  $(GMEM) #$(GPAR)
-#FSTRICTFLAGS =  $(GDEBUG)  $(GMEM) #$(GPAR)
+#FFLAGS =  $(GPRODUCTION) #-g -fbacktrace -fcheck=all #-Wall 
+#FSTRICTFLAGS=  $(GPRODUCTION) #-g -fbacktrace -fcheck=all #-Wall 
 #AR     = ar 
-#CC     = cc
 
 # DO NOT include in CFLAGS the "-c" option to compile object only
 # this is done explicitly as needed in the Makefile
 
+CC              = icc
+#CC              = cc
 CFLAGS          = -O3 -w -m64 -std=c99
 
 64_BIT_REALS    = -r8

python/swiftest/swiftest/io.py

@@ -58,6 +58,11 @@ def read_swiftest_param(param_file_name, param):
 
         param['ISTEP_OUT'] = int(param['ISTEP_OUT'])
         param['ISTEP_DUMP'] = int(param['ISTEP_DUMP'])
+        param['OUT_TYPE'] = param['OUT_TYPE'].upper()
+        param['OUT_FORM'] = param['OUT_FORM'].upper()
+        param['OUT_STAT'] = param['OUT_STAT'].upper()
+        param['IN_TYPE'] = param['IN_TYPE'].upper()
+        param['IN_FORM'] = param['IN_FORM'].upper()
         param['T0'] = real2float(param['T0'])
         param['TSTOP'] = real2float(param['TSTOP'])
         param['DT'] = real2float(param['DT'])
@@ -74,11 +79,14 @@ def read_swiftest_param(param_file_name, param):
         param['RHILL_PRESENT'] = param['RHILL_PRESENT'].upper()
         param['FRAGMENTATION'] = param['FRAGMENTATION'].upper()
         if param['FRAGMENTATION'] == 'YES' and param['PARTICLE_OUT'] == '':
-            param['PARTICLE_OUT'] = 'particle.dat'
+            if param['OUT_TYPE'] == 'REAL8' or param['OUT_TYPE'] == 'REAL4':
+                param['PARTICLE_OUT'] = 'particle.dat'
         param['ROTATION'] = param['ROTATION'].upper()
         param['TIDES'] = param['TIDES'].upper()
         param['ENERGY'] = param['ENERGY'].upper()
         param['GR'] = param['GR'].upper()
+        param['INTERACTION_LOOPS'] = param['INTERACTION_LOOPS'].upper()
+        param['ENCOUNTER_CHECK'] = param['ENCOUNTER_CHECK'].upper()
         if 'GMTINY' in param:
             param['GMTINY'] = real2float(param['GMTINY'])
     except IOError:
@@ -1517,6 +1525,11 @@ def swiftest2swifter_param(swiftest_param, J2=0.0, J4=0.0):
        swifter_param['OUT_TYPE'] = 'REAL4'
     if swifter_param['OUT_FORM'] == 'XVEL':
        swifter_param['OUT_FORM'] = 'XV'
+    IN_FORM = swifter_param.pop("IN_FORM", None)
+    INTERACTION_LOOPS = swifter_param.pop("INTERACTION_LOOPS", None)
+    ENCOUNTER_CHECK = swifter_param.pop("ENCOUNTER_CHECK", None)
+    ENCOUNTER_CHECK_PLPL = swifter_param.pop("ENCOUNTER_CHECK_PLPL", None)
+    ENCOUNTER_CHECK_PLTP = swifter_param.pop("ENCOUNTER_CHECK_PLTP", None)
     swifter_param['! VERSION'] = "Swifter parameter file converted from Swiftest"
 
     return swifter_param

python/swiftest/swiftest/simulation_class.py

@@ -34,12 +34,12 @@ def __init__(self, codename="Swiftest", param_file=""):
             'CHK_QMIN': "-1.0",
             'CHK_QMIN_COORD': "HELIO",
             'CHK_QMIN_RANGE': "-1.0 -1.0",
-            'ENC_OUT': "enc.dat",
+            'ENC_OUT': "",
             'MU2KG': constants.MSun,
             'TU2S': constants.JD2S,
             'DU2M': constants.AU2M,
             'EXTRA_FORCE': "NO",
-            'DISCARD_OUT': "discard.out",
+            'DISCARD_OUT': "",
             'PARTICLE_OUT' : "",
             'BIG_DISCARD': "NO",
             'CHK_CLOSE': "YES",
@@ -48,7 +48,9 @@ def __init__(self, codename="Swiftest", param_file=""):
             'ROTATION': "NO",
             'TIDES': "NO",
             'ENERGY': "NO",
-            'GR': "NO",
+            'GR': "YES",
+            'INTERACTION_LOOPS': "ADAPTIVE",
+            'ENCOUNTER_CHECK': "ADAPTIVE"
         }
         self.codename = codename
         if param_file != "" :

src/io/io.f90

@@ -721,6 +721,8 @@ module subroutine io_param_reader(self, unit, iotype, v_list, iostat, iomsg)
                   read(param_value, *, err = 667, iomsg = iomsg) param%Euntracked
                case ("MAXID")
                   read(param_value, *, err = 667, iomsg = iomsg) param%maxid 
+               case ("MAXID_COLLISION")
+                  read(param_value, *, err = 667, iomsg = iomsg) param%maxid_collision
                case ("PARTICLE_OUT")
                   param%particle_out = param_value
                case ("NPLMAX", "NTPMAX", "GMTINY", "MIN_GMFRAG", "FRAGMENTATION", "SEED", "YARKOVSKY", "YORP") ! Ignore SyMBA-specific, not-yet-implemented, or obsolete input parameters
@@ -1020,6 +1022,7 @@ module subroutine io_param_writer(self, unit, iotype, v_list, iostat, iomsg)
          end if
          call io_param_writer_one("FIRSTKICK",param%lfirstkick, unit)
          call io_param_writer_one("MAXID",param%maxid, unit)
+         call io_param_writer_one("MAXID_COLLISION",param%maxid_collision, unit)
 
          iostat = 0
          iomsg = "UDIO not implemented"

src/modules/swiftest_classes.f90

@@ -96,6 +96,7 @@ module swiftest_classes
       integer(I4B)         :: integrator     = UNKNOWN_INTEGRATOR !! Symbolic name of the nbody integrator  used
       character(STRMAX)    :: param_file_name = "param.in"        !! The default name of the parameter input file
       integer(I4B)         :: maxid          = -1                 !! The current maximum particle id number 
+      integer(I4B)         :: maxid_collision = 0                !! The current maximum collision 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

src/modules/symba_classes.f90

@@ -173,7 +173,6 @@ module symba_classes
       class(symba_plplenc), allocatable :: plplenc_list       !! List of massive body-massive body encounters in a single step
       class(symba_plplenc), allocatable :: plplcollision_list !! List of massive body-massive body collisions in a single step
       integer(I4B)                      :: irec               !! System recursion level
-      integer(I4B)                      :: collision_counter  !! Counter for number of collisions over total simulation 
    contains
       procedure :: write_discard    => symba_io_write_discard             !! Write out information about discarded and merged planets and test particles in SyMBA
       procedure :: initialize       => symba_setup_initialize_system      !! Performs SyMBA-specific initilization steps

src/netcdf/netcdf.f90

@@ -1057,6 +1057,7 @@ module subroutine netcdf_write_particle_info_base(self, iu)
          call check( nf90_put_var(iu%ncid, iu%origin_type_varid, emptystr, start=[1, idslot], count=[NAMELEN, 1]) )
          call check( nf90_put_var(iu%ncid, iu%origin_type_varid, charstring, start=[1, idslot], count=[strlen, 1]) )
 
+         call check( nf90_put_var(iu%ncid, iu%collision_id_varid, self%info%collision_id, start=[idslot]) )
          call check( nf90_put_var(iu%ncid, iu%origin_time_varid, self%info%origin_time, start=[idslot]) )
          call check( nf90_put_var(iu%ncid, iu%origin_xhx_varid, self%info%origin_xh(1), start=[idslot]) )
          call check( nf90_put_var(iu%ncid, iu%origin_xhy_varid, self%info%origin_xh(2), start=[idslot]) )

src/setup/setup.f90

@@ -106,12 +106,12 @@ module subroutine setup_initialize_particle_info_system(self, param)
 
       associate(cb => self%cb, pl => self%pl, npl => self%pl%nbody, tp => self%tp, ntp => self%tp%nbody)
 
-         call cb%info%set_value(particle_type=CB_TYPE_NAME, status="ACTIVE", origin_type="Initial conditions", origin_time=param%t0, origin_xh=[0.0_DP, 0.0_DP, 0.0_DP], origin_vh=[0.0_DP, 0.0_DP, 0.0_DP], collision_id=0)
+         call cb%info%set_value(particle_type=CB_TYPE_NAME, status="ACTIVE", origin_type="Initial conditions", origin_time=param%t0, origin_xh=[0.0_DP, 0.0_DP, 0.0_DP], origin_vh=[0.0_DP, 0.0_DP, 0.0_DP])
          do i = 1, self%pl%nbody
-            call pl%info(i)%set_value(particle_type=PL_TYPE_NAME, status="ACTIVE", origin_type="Initial conditions", origin_time=param%t0, origin_xh=self%pl%xh(:,i), origin_vh=self%pl%vh(:,i), collision_id=0)
+            call pl%info(i)%set_value(particle_type=PL_TYPE_NAME, status="ACTIVE", origin_type="Initial conditions", origin_time=param%t0, origin_xh=self%pl%xh(:,i), origin_vh=self%pl%vh(:,i))
          end do
          do i = 1, self%tp%nbody
-            call tp%info(i)%set_value(particle_type=TP_TYPE_NAME, status="ACTIVE", origin_type="Initial conditions", origin_time=param%t0, origin_xh=self%tp%xh(:,i), origin_vh=self%tp%vh(:,i), collision_id=0)
+            call tp%info(i)%set_value(particle_type=TP_TYPE_NAME, status="ACTIVE", origin_type="Initial conditions", origin_time=param%t0, origin_xh=self%tp%xh(:,i), origin_vh=self%tp%vh(:,i))
          end do
 
       end associate

src/symba/symba_collision.f90

@@ -718,7 +718,9 @@ subroutine symba_collision_mergeaddsub(system, param, colliders, frag, status)
                ncolliders = colliders%ncoll
                nfrag = frag%nbody
 
-               system%collision_counter = system%collision_counter + 1
+               param%maxid_collision = max(param%maxid_collision, maxval(system%pl%info(:)%collision_id))
+               param%maxid_collision = param%maxid_collision + 1
+
                ! Setup new bodies
                allocate(plnew, mold=pl)
                call plnew%setup(nfrag, param)
@@ -746,7 +748,7 @@ subroutine symba_collision_mergeaddsub(system, param, colliders, frag, status)
                   plnew%status(1:nfrag) = NEW_PARTICLE
                   do i = 1, nfrag
                      write(newname, FRAGFMT) frag%id(i)
-                     call plnew%info(i)%set_value(origin_type="Disruption", origin_time=param%t, name=newname, origin_xh=plnew%xh(:,i), origin_vh=plnew%vh(:,i), collision_id=system%collision_counter)
+                     call plnew%info(i)%set_value(origin_type="Disruption", origin_time=param%t, name=newname, origin_xh=plnew%xh(:,i), origin_vh=plnew%vh(:,i), collision_id=param%maxid_collision)
                   end do
                   do i = 1, ncolliders
                      if (colliders%idx(i) == ibiggest) then
@@ -760,7 +762,7 @@ subroutine symba_collision_mergeaddsub(system, param, colliders, frag, status)
                   plnew%status(1:nfrag) = NEW_PARTICLE
                   do i = 1, nfrag
                      write(newname, FRAGFMT) frag%id(i)
-                     call plnew%info(i)%set_value(origin_type="Supercatastrophic", origin_time=param%t, name=newname, origin_xh=plnew%xh(:,i), origin_vh=plnew%vh(:,i), collision_id=system%collision_counter)
+                     call plnew%info(i)%set_value(origin_type="Supercatastrophic", origin_time=param%t, name=newname, origin_xh=plnew%xh(:,i), origin_vh=plnew%vh(:,i), collision_id=param%maxid_collision)
                   end do
                   do i = 1, ncolliders
                      if (colliders%idx(i) == ibiggest) then
@@ -775,7 +777,7 @@ subroutine symba_collision_mergeaddsub(system, param, colliders, frag, status)
                   plnew%status(1) = OLD_PARTICLE
                   do i = 2, nfrag
                      write(newname, FRAGFMT) frag%id(i)
-                     call plnew%info(i)%set_value(origin_type="Hit and run fragment", origin_time=param%t, name=newname, origin_xh=plnew%xh(:,i), origin_vh=plnew%vh(:,i), collision_id=system%collision_counter)
+                     call plnew%info(i)%set_value(origin_type="Hit and run fragment", origin_time=param%t, name=newname, origin_xh=plnew%xh(:,i), origin_vh=plnew%vh(:,i), collision_id=param%maxid_collision)
                   end do
                   do i = 1, ncolliders
                      if (colliders%idx(i) == ibiggest) cycle