diff --git a/examples/.gitignore b/examples/.gitignore index 90ecccf87..dbdf17577 100644 --- a/examples/.gitignore +++ b/examples/.gitignore @@ -9,3 +9,4 @@ !solar_impact !Swifter_Swiftest !swiftest_performance +!rmvs_swifter_comparison diff --git a/examples/rmvs_swifter_comparison/.gitignore b/examples/rmvs_swifter_comparison/.gitignore new file mode 100644 index 000000000..ecf4c57c8 --- /dev/null +++ b/examples/rmvs_swifter_comparison/.gitignore @@ -0,0 +1,3 @@ +* +!.gitignore +!1pl_1tp_encounter diff --git a/examples/rmvs_swifter_comparison/1pl_1tp_encounter/.gitignore b/examples/rmvs_swifter_comparison/1pl_1tp_encounter/.gitignore new file mode 100644 index 000000000..80c7e7bd6 --- /dev/null +++ b/examples/rmvs_swifter_comparison/1pl_1tp_encounter/.gitignore @@ -0,0 +1,3 @@ +* +!.gitignore +!init_cond.py diff --git a/examples/rmvs_swifter_comparison/1pl_1tp_encounter/init_cond.py b/examples/rmvs_swifter_comparison/1pl_1tp_encounter/init_cond.py new file mode 100755 index 000000000..9d997f720 --- /dev/null +++ b/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() \ No newline at end of file diff --git a/python/swiftest/swiftest/init_cond.py b/python/swiftest/swiftest/init_cond.py index 30fb46c6d..a4c47548c 100644 --- a/python/swiftest/swiftest/init_cond.py +++ b/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} diff --git a/python/swiftest/swiftest/io.py b/python/swiftest/swiftest/io.py index 0c0007488..f7517e5af 100644 --- a/python/swiftest/swiftest/io.py +++ b/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 diff --git a/python/swiftest/swiftest/simulation_class.py b/python/swiftest/swiftest/simulation_class.py index 32171e72d..59f4e6669 100644 --- a/python/swiftest/swiftest/simulation_class.py +++ b/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) diff --git a/src/swiftest/swiftest_io.f90 b/src/swiftest/swiftest_io.f90 index d63a1a95d..8cb720ec0 100644 --- a/src/swiftest/swiftest_io.f90 +++ b/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