Merge branch 'debug'

examples/.gitignore

@@ -9,3 +9,4 @@
 !solar_impact
 !Swifter_Swiftest
 !swiftest_performance
+!rmvs_swifter_comparison

examples/rmvs_swifter_comparison/.gitignore

@@ -0,0 +1,3 @@
+*
+!.gitignore
+!1pl_1tp_encounter

examples/rmvs_swifter_comparison/1pl_1tp_encounter/.gitignore

@@ -0,0 +1,3 @@
+*
+!.gitignore
+!init_cond.py

examples/rmvs_swifter_comparison/1pl_1tp_encounter/init_cond.py

@@ -0,0 +1,38 @@
+#!/usr/bin/env python3
+import numpy as np
+import swiftest 
+
+# 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	= 0.25 * swiftest.JD2S / swiftest.YR2S   # simulation step size
+end_sim = 0.15
+t_print = deltaT  #output interval to print results
+
+iout = int(np.ceil(t_print / deltaT))
+
+radius = np.double(4.25875607065041e-05)
+Gmass = np.double(0.00012002693582795244940133) 
+apl = np.longdouble(1.0)
+atp = np.longdouble(1.01)
+vpl = np.longdouble(2 * np.pi)
+vtp = np.longdouble(2 * np.pi / np.sqrt(atp))
+
+p_pl = np.array([apl, 0.0, 0.0], dtype=np.double)
+v_pl = np.array([0.0, vpl, 0.0], dtype=np.double)
+
+p_tp = np.array([atp, 0.0, 0.0], dtype=np.double)
+v_tp = np.array([0.0, vtp, 0.0], dtype=np.double)
+
+rhill = np.double(apl * 0.0100447248332378922085)
+
+sim = swiftest.Simulation(simdir="swiftest_sim",init_cond_format="XV", general_relativity=False, integrator="RMVS")
+sim.clean()
+sim.add_solar_system_body(["Sun"])
+sim.add_body(name=["Planet"], rh=[p_pl], vh=[v_pl], Gmass=[Gmass], radius=[radius], rhill=[rhill])
+sim.add_body(name=["TestParticle"], rh=[p_tp], vh=[v_tp])
+sim.set_parameter(tstart=t_0, tstop=end_sim, dt=deltaT, istep_out=iout, dump_cadence=0)
+sim.save()
+
+sim.set_parameter(simdir="swifter_sim",codename="Swifter",init_cond_file_type="ASCII",output_format="XV",output_file_type="REAL8")
+sim.save()

python/swiftest/swiftest/init_cond.py

@@ -277,8 +277,8 @@ def vec2xr(param: Dict, **kwargs: Any):
     space_coords = np.array(["x","y","z"])
 
     vector_vars = ["rh","vh","Ip","rot"]
-    scalar_vars = ["name","a","e","inc","capom","omega","capm","Gmass","radius","rhill","J2","J4"]
-    time_vars =  ["rh","vh","Ip","rot","a","e","inc","capom","omega","capm","Gmass","radius","rhill","J2","J4"]
+    scalar_vars = ["name","a","e","inc","capom","omega","capm","Gmass","radius","rhill","j2rp2","j4rp4"]
+    time_vars =  ["rh","vh","Ip","rot","a","e","inc","capom","omega","capm","Gmass","radius","rhill","j2rp2","j4rp4"]
 
     # Check for valid keyword arguments
     kwargs = {k:kwargs[k] for k,v in kwargs.items() if v is not None}

python/swiftest/swiftest/io.py

@@ -16,6 +16,7 @@
 import sys
 import tempfile
 import re
+import os
 
 # This defines features that are new in Swiftest and not in Swifter (for conversion between param.in files)
 newfeaturelist = ("RESTART",
@@ -30,9 +31,14 @@
                   "SEED",
                   "INTERACTION_LOOPS",
                   "ENCOUNTER_CHECK",
+                  "ENCOUNTER_CHECK_PLPL",
+                  "ENCOUNTER_CHECK_PLTP",
                   "TSTART",
                   "DUMP_CADENCE",
                   "ENCOUNTER_SAVE",
+                  "MIN_GMFRAG",
+                  "NFRAG_REDUCTION",
+                  "COLLISION_MODEL",
                   "COARRAY")
 
 # This list defines features that are booleans, so must be converted to/from string when writing/reading from file
@@ -1227,7 +1233,7 @@ def swiftest_xr2infile(ds, param, in_type="NETCDF_DOUBLE", infile_name=None,fram
         print(f"{in_type} is an unknown file type")
 
 
-def swifter_xr2infile(ds, param, framenum=-1):
+def swifter_xr2infile(ds, param, simdir=os.getcwd, framenum=-1):
     """
     Writes a set of Swifter input files from a single frame of a Swiftest xarray dataset
 
@@ -1266,7 +1272,7 @@ def swifter_xr2infile(ds, param, framenum=-1):
 
     if param['IN_TYPE'] == 'ASCII':
         # Swiftest Central body file
-        plfile = open(param['PL_IN'], 'w')
+        plfile = open(os.path.join(simdir,param['PL_IN']), 'w')
         print(pl.id.count().values + 1, file=plfile)
         print(cb.id.values[0], GMSun, file=plfile)
         print('0.0 0.0 0.0', file=plfile)
@@ -1279,21 +1285,20 @@ def swifter_xr2infile(ds, param, framenum=-1):
                 print(i.values, pli['Gmass'].values, file=plfile)
             if param['CHK_CLOSE']:
                 print(pli['radius'].values, file=plfile)
-            print(pli['rh'].values[0,0], pli['ry'].values[0,1], pli['rh'].values[0,2], file=plfile)
-            print(pli['vh'].values[0,0], pli['vh'].values[0,1], pli['vh'].values[0,2], file=plfile)
+            print(pli['rh'].values[0], pli['rh'].values[1], pli['rh'].values[2], file=plfile)
+            print(pli['vh'].values[0], pli['vh'].values[1], pli['vh'].values[2], file=plfile)
         plfile.close()
 
         # TP file
-        tpfile = open(param['TP_IN'], 'w')
+        tpfile = open(os.path.join(simdir,param['TP_IN']), 'w')
         print(tp.id.count().values, file=tpfile)
         for i in tp.id:
             tpi = tp.sel(id=i)
             print(i.values, file=tpfile)
-            print(tpi['rh'].values[0,0], tpi['ry'].values[0,1], tpi['rh'].values[0,2], file=tpfile)
-            print(tpi['vh'].values[0,0], tpi['vh'].values[0,1], tpi['vh'].values[0,2], file=tpfile)
+            print(tpi['rh'].values[0], tpi['rh'].values[1], tpi['rh'].values[2], file=tpfile)
+            print(tpi['vh'].values[0], tpi['vh'].values[1], tpi['vh'].values[2], file=tpfile)
         tpfile.close()
     else:
-        # Now make Swiftest files
         print(f"{param['IN_TYPE']} is an unknown input file type")
 
     return
@@ -1848,6 +1853,8 @@ def swiftest2swifter_param(swiftest_param, J2=0.0, J4=0.0):
            swifter_param['C'] =  swiftest.einsteinC * np.longdouble(TU2S) / np.longdouble(DU2M)
     for key in newfeaturelist:
        tmp = swifter_param.pop(key, None)
+    if "ISTEP_DUMP" not in swifter_param:
+        swifter_param["ISTEP_DUMP"] = swifter_param["ISTEP_OUT"]
     swifter_param['J2'] = J2
     swifter_param['J4'] = J4
     if swifter_param['OUT_STAT'] == "REPLACE":
@@ -1858,11 +1865,6 @@ 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

@@ -2223,9 +2223,9 @@ def add_solar_system_body(self,
             values = list(np.hsplit(np.array(body_list[0],dtype=np.dtype(object)),17))
         else:
             values = list(np.squeeze(np.hsplit(np.array(body_list,np.dtype(object)),17)))
-        keys = ["id","name","a","e","inc","capom","omega","capm","rh","vh","Gmass","radius","rhill","Ip","rot","J2","J4"]
+        keys = ["id","name","a","e","inc","capom","omega","capm","rh","vh","Gmass","radius","rhill","Ip","rot","j2rp2","j4rp4"]
         kwargs = dict(zip(keys,values))
-        scalar_floats = ["a","e","inc","capom","omega","capm","Gmass","radius","rhill","J2","J4"]
+        scalar_floats = ["a","e","inc","capom","omega","capm","Gmass","radius","rhill","j2rp2","j4rp4"]
         vector_floats = ["rh","vh","Ip","rot"]
         scalar_ints = ["id"]
 
@@ -2434,7 +2434,7 @@ def add_body(self,
         Hill's radius values if these are massive bodies
         rot: (3) or (n,3) array-like of float, optional
         Rotation rate vectors if these are massive bodies with rotation enabled.
-        Ip: (3) or (n,3) array-like of flaot, optional
+        Ip: (3) or (n,3) array-like of float, optional
         Principal axes moments of inertia vectors if these are massive bodies with rotation enabled.
 
         Returns
@@ -2451,7 +2451,7 @@ def input_to_array(val,t,n=None):
             elif t == "i":
                 t = np.int64
             elif t == "s":
-                t = np.str
+                t = str
 
             if val is None:
                 return None, n
@@ -2551,7 +2551,7 @@ def input_to_array_3d(val,n=None):
                 Gmass = self.GU * mass
 
         dsnew = init_cond.vec2xr(self.param, name=name, a=a, e=e, inc=inc, capom=capom, omega=omega, capm=capm, id=id,
-                                 Gmass=Gmass, radius=radius, rhill=rhill, Ip=Ip, rh=rh, vh=vh,rot=rot, J2=J2, J4=J4, time=time)
+                                 Gmass=Gmass, radius=radius, rhill=rhill, Ip=Ip, rh=rh, vh=vh,rot=rot, j2rp2=J2, j4rp4=J4, time=time)
 
         dsnew = self._combine_and_fix_dsnew(dsnew)
         self.save(verbose=False)
@@ -2984,12 +2984,11 @@ def save(self,
             else:
                 shutil.copy2(self.binary_source, self.driver_executable)  
         elif codename == "Swifter":
-            if codename == "Swiftest":
-                swifter_param = io.swiftest2swifter_param(param)
-            else:
-                swifter_param = param
-            io.swifter_xr2infile(self.data, swifter_param, framenum)
-            self.write_param(param_file, param=swifter_param,**kwargs)
+            swifter_param = io.swiftest2swifter_param(param)
+            if "rhill" in self.data:
+                swifter_param['RHILL_PRESENT'] = 'YES'
+            io.swifter_xr2infile(self.data, swifter_param, self.simdir, framenum)
+            self.write_param(codename=codename,param_file=param_file,param=swifter_param,**kwargs)
         else:
             warnings.warn(f'Saving to {codename} not supported',stacklevel=2)
 

src/swiftest/swiftest_io.f90

@@ -2930,6 +2930,15 @@ module subroutine swiftest_io_read_in_system(self, nc, param)
       if (.not.param%loblatecb) then
          if (allocated(self%pl%aobl)) deallocate(self%pl%aobl)
          if (allocated(self%tp%aobl)) deallocate(self%tp%aobl)
+      else
+         if (self%pl%nbody > 0) then
+            if (.not. allocated(self%pl%aobl)) allocate(self%pl%aobl(NDIM,self%pl%nbody))
+            self%pl%aobl(:,:) = 0.0_DP
+         end if
+         if (self%tp%nbody > 0) then
+            if (.not. allocated(self%tp%aobl)) allocate(self%tp%aobl(NDIM,self%tp%nbody))
+            self%tp%aobl(:,:) = 0.0_DP
+         end if
       end if
 
       return