Added an example involving multiple test particles encountering multiple planets

Makefile.Defines

@@ -51,29 +51,32 @@ ADVIXE_FLAGS = -g -O2 -qopt-report=5 -vec -vecabi=cmdtarget -simd -shared-intel
 VTUNE_FLAGS = -g -O2 -vec -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 -warn all -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 
+IDEBUG = -O0 -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 =  -mp1 -prec-div -prec-sqrt -ftz -assume protect-parens -fp-speculation safe
-SIMDVEC = -vec -simd -xhost -fma -align all -assume contiguous_assumed_shape -vecabi=cmdtarget -qopt-zmm-usage=high -fp-model no-except
+STRICTREAL = -mp1 -fp-model strict -prec-div -prec-sqrt -assume protect-parens 
+SIMDVEC = -simd -xhost -align all -assume contiguous_assumed_shape -vecabi=cmdtarget -prec-div -prec-sqrt -assume protect-parens
 PAR = -parallel -qopenmp 
 HEAPARR = -heap-arrays 1048576 
-OPTIMIZE = -qopt-report=5
+OPTIMIZE = -qopt-report=5 
 
-#FFLAGS  = -init=snan,arrays -traceback -no-wrap-margin -O3 -g -CB -nogen-interfaces -no-pie -fp-speculation=safe $(SIMDVEC) $(PAR) #$(HEAPARR)
-FORTRAN = ifort
-FFLAGS  = $(IDEBUG) 
-#FFLAGS  = -init=snan,arrays -traceback -no-wrap-margin -O3 -g -shared-intel -debug inline-debug-info -qopt-report=5 $(SIMDVEC) $(STRICTREAL) $(PAR)
 
+#debug flags
 
-GDEBUG = -ggdb -g3 -Og -fbacktrace -fbounds-check 
+GDEBUG = -g -O0 -fbacktrace -fbounds-check 
 GPRODUCTION = -O3
 GPAR = -fopenmp -ftree-parallelize-loops=4
 GMEM = -fsanitize=undefined -fsanitize=address -fsanitize=leak 
 GWARNINGS = -Wall -Warray-bounds  -Wimplicit-interface -Wextra -Warray-temporaries
 
-#FORTRAN = gfortran
-#FFLAGS = $(GDEBUG)
+#FFLAGS  = -init=snan,arrays -traceback -no-wrap-margin -O3 -g -CB -nogen-interfaces -no-pie -fp-speculation=safe $(SIMDVEC) $(PAR) #$(HEAPARR)
+#FFLAGS  = $(IDEBUG) $(HEAPARR)
+FFLAGS  = -init=snan,arrays -no-wrap-margin -fp-model strict -fp-model no-except -traceback -g $(OPTIMIZE) -O3  $(SIMDVEC) $(PAR) $(HEAPARR) 
+FORTRAN = ifort
+#AR     = xiar
 
+#FORTRAN = gfortran
+#FFLAGS =  -ffree-line-length-none $(GDEBUG) #$(GMEM)
+AR     = ar
 
 # DO NOT include in CFLAGS the "-c" option to compile object only
 # this is done explicitly as needed in the Makefile

examples/rmvs_swifter_comparison/9pl_18tp_encounters/cb.swiftest.in

@@ -0,0 +1,4 @@
+39.47692640889762629
+0.0046504672609621575315
+4.7535806948127356533e-12
+-2.2473967953572827815e-18

examples/rmvs_swifter_comparison/9pl_18tp_encounters/config.swiftest.in

@@ -0,0 +1,29 @@
+! Swiftest input file generated using init_cond.py
+T0             0 
+TSTOP          1.0
+DT             0.0027378507871321013
+CB_IN          cb.swiftest.in
+PL_IN          pl.swiftest.in
+TP_IN          tp.in
+IN_TYPE        ASCII
+ISTEP_OUT      1
+ISTEP_DUMP     1
+BIN_OUT        bin.swiftest.dat
+OUT_TYPE       REAL8
+OUT_FORM       XV
+OUT_STAT       REPLACE
+CHK_CLOSE      yes
+CHK_RMIN       0.004650467260962157
+CHK_RMAX       1000.0
+CHK_EJECT      1000.0
+CHK_QMIN       0.004650467260962157
+CHK_QMIN_COORD HELIO
+CHK_QMIN_RANGE 0.004650467260962157 1000.0
+ENC_OUT        enc.swiftest.dat
+EXTRA_FORCE    no
+BIG_DISCARD    no
+ROTATION       no
+GR             no
+MU2KG          1.988409870698051e+30
+DU2M           149597870700.0
+TU2S           31557600.0

examples/rmvs_swifter_comparison/9pl_18tp_encounters/init_cond.py

@@ -0,0 +1,245 @@
+"""
+For testing RMVS, the code generates clones of test particles based on one that is fated to impact Mercury.
+To use the script, modify the variables just after the  "if __name__ == '__main__':" line
+"""
+import numpy as np
+from astroquery.jplhorizons import Horizons
+import astropy.constants as const 
+import swiftestio as swio
+from scipy.io import FortranFile
+import sys
+
+#Values from JPL Horizons
+AU2M = np.longdouble(const.au.value)
+GMSunSI = np.longdouble(const.GM_sun.value)
+Rsun = np.longdouble(const.R_sun.value)
+GC = np.longdouble(const.G.value)
+JD = 86400
+year = np.longdouble(365.25 * JD)
+c = np.longdouble(299792458.0)
+
+MU2KG    = np.longdouble(GMSunSI / GC) #Conversion from mass unit to kg
+DU2M     = np.longdouble(AU2M)         #Conversion from radius unit to centimeters
+TU2S     = np.longdouble(year)         #Conversion from time unit to seconds
+GU       = np.longdouble(GC / (DU2M**3 / (MU2KG * TU2S**2)))
+GMSun = np.longdouble(GMSunSI / (DU2M**3 / TU2S**2))
+
+# Solar oblatenes values: From Mecheri et al. (2004), using Corbard (b) 2002 values (Table II)
+J2 = np.longdouble(2.198e-7) * (Rsun / DU2M)**2
+J4 = np.longdouble(-4.805e-9) * (Rsun / DU2M)**4
+
+npl = 9
+ntp = 2 * npl
+
+# Planet ids
+plid = {
+   'mercury'     : '1',
+   'venus'       : '2',
+   'earthmoon'   : '3',
+   'mars'        : '4',
+   'jupiter'     : '5',
+   'saturn'      : '6',
+   'uranus'      : '7',
+   'neptune'     : '8',
+   'plutocharon' : '9'
+}
+
+#Planet Msun/M ratio
+MSun_over_Mpl = {
+   'mercury'     : np.longdouble(6023600.0),
+   'venus'       : np.longdouble(408523.71),
+   'earthmoon'   : np.longdouble(328900.56),
+   'mars'        : np.longdouble(3098708.),
+   'jupiter'     : np.longdouble(1047.3486),
+   'saturn'      : np.longdouble(3497.898),
+   'uranus'      : np.longdouble(22902.98),
+   'neptune'     : np.longdouble(19412.24),
+   'plutocharon' : np.longdouble(1.35e8)
+}
+
+#Planet radii in meters
+Rpl = {
+   'mercury'     : np.longdouble(2439.4e3),
+   'venus'       : np.longdouble(6051.8e3),
+   'earthmoon'   : np.longdouble(6371.0084e3), # Earth only for radius
+   'mars'        : np.longdouble(3389.50e3),
+   'jupiter'     : np.longdouble(69911e3),
+   'saturn'      : np.longdouble(58232.0e3),
+   'uranus'      : np.longdouble(25362.e3),
+   'neptune'     : np.longdouble(24622.e3),
+   'plutocharon' : np.longdouble(1188.3e3)
+}
+
+THIRDLONG = np.longdouble(1.0) / np.longdouble(3.0)
+
+# Both codes use the same tp input file
+tpin = "tp.in"
+
+swifter_input  = "param.swifter.in"
+swifter_pl     = "pl.swifter.in"
+swifter_bin    = "bin.swifter.dat"
+swifter_enc    = "enc.swifter.dat"
+
+swiftest_input = "config.swiftest.in"
+swiftest_pl    = "pl.swiftest.in"
+swiftest_cb    = "cb.swiftest.in"
+swiftest_bin   = "bin.swiftest.dat"
+swiftest_enc   = "enc.swiftest.dat"
+
+# Simple initial conditions of a circular planet with one test particle in a close encounter state
+# Simulation start, stop, and output cadence times
+t_0	  = 0 # simulation start time
+deltaT	= 1.00 * JD / TU2S   # simulation step size
+end_sim = year / TU2S # simulation end time
+t_print = deltaT  #output interval to print results
+
+iout = int(np.ceil(t_print / deltaT))
+rmin = Rsun / DU2M
+rmax = 1000.0
+
+# Dates to fetch planet ephemerides from JPL Horizons
+tstart = '2021-06-15'
+tend = '2021-06-16'
+tstep = '1d'
+
+#######################################################
+# Start generating initial conditions for the planets #
+#######################################################
+pldata = {}
+tpdata = {}
+plvec = {}
+tpvec = {}
+Rhill = {}
+
+for key,val in plid.items():
+   pldata[key] = Horizons(id=val, id_type='majorbody',location='@sun',
+            epochs={'start': tstart, 'stop': tend,
+            'step': tstep})
+   plvec[key] = np.array([
+         pldata[key].vectors()['x'][0],
+         pldata[key].vectors()['y'][0],
+         pldata[key].vectors()['z'][0],
+         pldata[key].vectors()['vx'][0],
+         pldata[key].vectors()['vy'][0],
+         pldata[key].vectors()['vz'][0]
+         ])
+   Rhill[key] = pldata[key].elements()['a'][0] * (3 * MSun_over_Mpl[key])**(-1.0 / 3.0)
+
+# For each planet, we will initialize a pair of test particles. One on its way in, and one on its way out. We will also initialize two additional particles that don't encounter anything
+tpnames = range(101, 101 + ntp)
+tpxv1 = np.empty((6))
+tpxv2 = np.empty((6))
+for idx,key in enumerate(plid):
+    rstart = 2 * Rpl[key] / DU2M # Start the test particles at a multiple of the planet radius away
+    vstart = 1.5 * np.sqrt(2 * GMSun / MSun_over_Mpl[key] / rstart)  # Start the test particle velocities at a multiple of the escape speed
+    xvstart = np.array([rstart / np.sqrt(2.0), rstart / np.sqrt(2.0), 0.0, vstart, 0.0, 0.0])
+    # The positions and velocities of each pair of test particles will be in reference to a planet
+    tpxv1 = plvec[key] + xvstart
+    tpxv2 = plvec[key] - xvstart
+    tpvec[tpnames[idx]] = tpxv1
+    tpvec[tpnames[idx + npl]] = tpxv2
+
+# TP file
+tpfile = open(tpin, 'w')
+print(ntp, file=tpfile)
+for tpid, tp in tpvec.items():
+    print(tpid, file=tpfile)
+    print(f'{tp[0]} {tp[1]} {tp[2]}', file=tpfile)
+    print(f'{tp[3]} {tp[4]} {tp[5]}', file=tpfile)
+tpfile.close()
+
+# Swifter PL file
+plfile = open(swifter_pl, 'w')
+print(npl + 1, file=plfile)
+print(1,GMSun,file=plfile)
+print('0.0 0.0 0.0',file=plfile)
+print('0.0 0.0 0.0',file=plfile)
+for key, pl in plvec.items():
+    print(f'{int(plid[key])+1} {GMSun / MSun_over_Mpl[key]} {Rhill[key]} ! {key}', file=plfile)
+    print(f'{Rpl[key] / DU2M}', file=plfile)
+    print(f'{pl[0]} {pl[1]} {pl[2]}', file=plfile)
+    print(f'{pl[3]} {pl[4]} {pl[5]}', file=plfile)
+plfile.close()
+
+# Swiftest Central body and pl file
+cbfile = open(swiftest_cb, 'w')
+print(GMSun, file=cbfile)
+print(rmin, file=cbfile)
+print(J2, file=cbfile)
+print(J4, file=cbfile)
+cbfile.close()
+plfile = open(swiftest_pl, 'w')
+print(npl, file=plfile)
+for key, pl in plvec.items():
+    print(f'{int(plid[key]) + 1} {GMSun / MSun_over_Mpl[key]} ! {key}', file=plfile)
+    print(f'{Rpl[key] / DU2M}', file=plfile)
+    print(f'{pl[0]} {pl[1]} {pl[2]}', file=plfile)
+    print(f'{pl[3]} {pl[4]} {pl[5]}', file=plfile)
+plfile.close()
+
+# Swifter parameter file
+sys.stdout = open(swifter_input, "w")
+print(f'! Swifter input file generated using init_cond.py')
+print(f'T0            {t_0} ')
+print(f'TSTOP         {end_sim}')
+print(f'DT            {deltaT}')
+print(f'PL_IN         {swifter_pl}')
+print(f'TP_IN         {tpin}')
+print(f'IN_TYPE       ASCII')
+print(f'ISTEP_OUT     {iout:d}')
+print(f'ISTEP_DUMP    {iout:d}')
+print(f'BIN_OUT       {swifter_bin}')
+print(f'OUT_TYPE      REAL8')
+print(f'OUT_FORM      XV')
+print(f'OUT_STAT      NEW')
+print(f'J2            {J2}')
+print(f'J4            {J4}')
+print(f'CHK_CLOSE     yes')
+print(f'CHK_RMIN      {rmin}')
+print(f'CHK_RMAX      {rmax}')
+print(f'CHK_EJECT     {rmax}')
+print(f'CHK_QMIN      {rmin}')
+print(f'CHK_QMIN_COORD HELIO')
+print(f'CHK_QMIN_RANGE {rmin} {rmax}')
+print(f'ENC_OUT        {swifter_enc}')
+print(f'EXTRA_FORCE    no')
+print(f'BIG_DISCARD    no')
+print(f'RHILL_PRESENT  yes')
+sys.stdout = sys.__stdout__
+
+# Swiftest configuration file
+sys.stdout = open(swiftest_input, "w")
+print(f'! Swiftest input file generated using init_cond.py')
+print(f'T0             {t_0} ')
+print(f'TSTOP          {end_sim}')
+print(f'DT             {deltaT}')
+print(f'CB_IN          {swiftest_cb}')
+print(f'PL_IN          {swiftest_pl}')
+print(f'TP_IN          {tpin}')
+print(f'IN_TYPE        ASCII')
+print(f'ISTEP_OUT      {iout:d}')
+print(f'ISTEP_DUMP     {iout:d}')
+print(f'BIN_OUT        {swiftest_bin}')
+print(f'OUT_TYPE       REAL8')
+print(f'OUT_FORM       XV')
+print(f'OUT_STAT       REPLACE')
+print(f'CHK_CLOSE      yes')
+print(f'CHK_RMIN       {rmin}')
+print(f'CHK_RMAX       {rmax}')
+print(f'CHK_EJECT      {rmax}')
+print(f'CHK_QMIN       {rmin}')
+print(f'CHK_QMIN_COORD HELIO')
+print(f'CHK_QMIN_RANGE {rmin} {rmax}')
+print(f'ENC_OUT        {swiftest_enc}')
+print(f'EXTRA_FORCE    no')
+print(f'BIG_DISCARD    no')
+print(f'ROTATION       no')
+print(f'GR             no')
+print(f'MU2KG          {MU2KG}')
+print(f'DU2M           {DU2M}')
+print(f'TU2S           {TU2S}')
+
+
+
+
+

examples/rmvs_swifter_comparison/9pl_18tp_encounters/param.swifter.in

@@ -0,0 +1,26 @@
+! Swifter input file generated using init_cond.py
+T0            0 
+TSTOP         1.0
+DT            0.0027378507871321013
+PL_IN         pl.swifter.in
+TP_IN         tp.in
+IN_TYPE       ASCII
+ISTEP_OUT     1
+ISTEP_DUMP    1
+BIN_OUT       bin.swifter.dat
+OUT_TYPE      REAL8
+OUT_FORM      XV
+OUT_STAT      NEW
+J2            4.7535806948127355e-12
+J4            -2.2473967953572827e-18
+CHK_CLOSE     yes
+CHK_RMIN      0.004650467260962157
+CHK_RMAX      1000.0
+CHK_EJECT     1000.0
+CHK_QMIN      0.004650467260962157
+CHK_QMIN_COORD HELIO
+CHK_QMIN_RANGE 0.004650467260962157 1000.0
+ENC_OUT        enc.swifter.dat
+EXTRA_FORCE    no
+BIG_DISCARD    no
+RHILL_PRESENT  yes

examples/rmvs_swifter_comparison/9pl_18tp_encounters/pl.swifter.in

@@ -0,0 +1,40 @@
+10
+1 39.47692640889762629
+0.0 0.0 0.0
+0.0 0.0 0.0
+2 6.553709809565314e-06 0.0014751258227142052 ! mercury
+1.6306381826061646e-05
+0.008059842448018334 -0.4616051037329109 -0.03846017738329229
+0.02248719132054853 0.001934639213990692 -0.001904656977422976
+3 9.663313399581539e-05 0.006759134232034941 ! venus
+4.0453784346544176e-05
+-0.5115875215389065 0.5030818749037324 0.03642547299277956
+-0.01425515725454357 -0.01452868630179309 0.0006232072038298823
+4 0.00012002693582795246 0.010044625087011915 ! earthmoon
+4.25875607065041e-05
+-0.1090020607540907 -1.009893805009766 4.823302918632528e-05
+0.01682491922568941 -0.001910549762056979 3.992660742687128e-08
+5 1.2739802010675942e-05 0.0072467897902424765 ! mars
+2.2657408050928896e-05
+-1.342897929331636 0.9778655112682739 0.05343398538723887
+-0.007712315645393206 -0.01011917844182223 -2.287744801261131e-05
+6 0.037692251088985676 0.3552720805286442 ! jupiter
+0.0004673261703049093
+3.923184193414315 -3.168419770483168 -0.0746147877972047
+0.004655552638985802 0.006232623300954468 -0.0001300429201057457
+7 0.011285899820091273 0.4376460836930155 ! saturn
+0.00038925687730393614
+6.185794462795267 -7.804174837804826 -0.110498432926239
+0.004066833203985018 0.003458637040736611 -0.0002219310939327014
+8 0.001723658947826773 0.46946272948265794 ! uranus
+0.00016953449859497232
+14.9290976575471 12.92949673572929 -0.1454099139559955
+-0.002599557960646664 0.002795888198858545 4.391864857782088e-05
+9 0.0020336100526728304 0.78119478483336 ! neptune
+0.00016458790412449367
+29.54416169025338 -4.716921603714237 -0.5838030174427992
+0.0004792636209523189 0.00312573757291745 -7.53264045199501e-05
+10 2.924216771029454e-07 0.05379680851617536 ! plutocharon
+7.943294877391593e-06
+14.54448346259197 -31.05223519593471 -0.8828000265625595
+0.002923077617691739 0.0006625916902153526 -0.0009142553677224461