diff --git a/examples/Basic_Simulation/initial_conditions.py b/examples/Basic_Simulation/initial_conditions.py index 1edc13879..b7c6a323c 100644 --- a/examples/Basic_Simulation/initial_conditions.py +++ b/examples/Basic_Simulation/initial_conditions.py @@ -30,11 +30,10 @@ from numpy.random import default_rng # Initialize the simulation object as a variable -sim = swiftest.Simulation() -sim.set_simulation_time(tstart=0.0, tstop=10.0, dt=0.005, tstep_out=1.0) +sim = swiftest.Simulation(tstart=0.0, tstop=10.0, dt=0.005, tstep_out=1.0, fragmentation=True, minimum_fragment_gmass = 1e-9) +sim.get_parameter() # Add parameter attributes to the simulation object sim.param['GMTINY'] = 1e-6 -sim.param['MIN_GMFRAG'] = 1e-9 # Add the modern planets and the Sun using the JPL Horizons Database sim.add_solar_system_body(["Sun","Mercury","Venus","Earth","Mars","Jupiter","Saturn","Uranus","Neptune","Pluto"]) diff --git a/examples/Basic_Simulation/test_io.ipynb b/examples/Basic_Simulation/test_io.ipynb new file mode 100644 index 000000000..ce92a8229 --- /dev/null +++ b/examples/Basic_Simulation/test_io.ipynb @@ -0,0 +1,143 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 1, + "id": "86c845ce-1801-46ca-8a8a-1cabb266e6a6", + "metadata": {}, + "outputs": [], + "source": [ + "import swiftest\n", + "import xarray as xr\n", + "import numpy as np\n", + "import os" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "d716c371-8eb4-4fc1-82af-8b5c444c831e", + "metadata": {}, + "outputs": [], + "source": [ + "sim = swiftest.Simulation()" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "id": "83cebbc1-387b-4ef5-b96e-76856b6672e5", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "integrator symba\n", + "codename Swiftest\n", + "driver_executable /Users/daminton/git/swiftest/build/swiftest_driver\n", + "t0 0.0 y\n", + "tstart 0.0 y\n", + "tstop NOT SET\n", + "dt NOT SET\n", + "istep_out NOT SET\n", + "istep_dump NOT SET\n", + "init_cond_file_type NETCDF_DOUBLE\n", + "init_cond_format EL\n", + "init_cond_file_name init_cond.nc\n", + "output_file_type NETCDF_DOUBLE\n", + "output_file_name bin.nc\n", + "output_format XVEL\n", + "rmin 0.004650467260962157 AU\n", + "rmax 10000.0 AU\n", + "qmin_coord HELIO\n", + "MU: MSun 1.988409870698051e+30 kg / MSun\n", + "DU: AU 149597870700.0 m / AU\n", + "TU: y 31557600.0 s / y\n", + "close_encounter_check True\n", + "general_relativity True\n", + "fragmentation True\n", + "rotation True\n", + "compute_conservation_values False\n", + "extra_force False\n", + "big_discard False\n", + "rhill_present False\n", + "restart False\n", + "interaction_loops TRIANGULAR\n", + "encounter_check_loops TRIANGULAR\n", + "ephemeris_date 2027-04-30\n" + ] + }, + { + "data": { + "text/plain": [ + "{'! VERSION': 'Swiftest parameter input',\n", + " 'T0': 0.0,\n", + " 'TSTART': 0.0,\n", + " 'IN_TYPE': 'NETCDF_DOUBLE',\n", + " 'IN_FORM': 'EL',\n", + " 'NC_IN': 'init_cond.nc',\n", + " 'OUT_TYPE': 'NETCDF_DOUBLE',\n", + " 'BIN_OUT': 'bin.nc',\n", + " 'OUT_FORM': 'XVEL',\n", + " 'CHK_RMIN': 0.004650467260962157,\n", + " 'CHK_RMAX': 10000.0,\n", + " 'CHK_QMIN_COORD': 'HELIO',\n", + " 'CHK_QMIN': 0.004650467260962157,\n", + " 'CHK_QMIN_RANGE': '0.004650467260962157 10000.0',\n", + " 'MU2KG': 1.988409870698051e+30,\n", + " 'DU2M': 149597870700.0,\n", + " 'TU2S': 31557600.0,\n", + " 'CHK_CLOSE': True,\n", + " 'GR': True,\n", + " 'FRAGMENTATION': True,\n", + " 'ROTATION': True,\n", + " 'ENERGY': False,\n", + " 'EXTRA_FORCE': False,\n", + " 'BIG_DISCARD': False,\n", + " 'RHILL_PRESENT': False,\n", + " 'RESTART': False,\n", + " 'INTERACTION_LOOPS': 'TRIANGULAR',\n", + " 'ENCOUNTER_CHECK': 'TRIANGULAR'}" + ] + }, + "execution_count": 3, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "sim.get_parameter()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ec7452d6-4c9b-4df3-acc0-b11c32264b91", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.9.13" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/python/swiftest/swiftest/simulation_class.py b/python/swiftest/swiftest/simulation_class.py index a731603ce..a359f8caa 100644 --- a/python/swiftest/swiftest/simulation_class.py +++ b/python/swiftest/swiftest/simulation_class.py @@ -14,6 +14,7 @@ from swiftest import init_cond from swiftest import tool from swiftest import constants +from swiftest import __file__ as _pyfile import os import datetime import xarray as xr @@ -66,12 +67,12 @@ def __init__(self, rmax: float = 10000.0, gmtiny: float | None = None, mtiny: float | None = None, - min_fragment_mass: float | None = None, - min_fragment_gmass: float | None = None, qmin_coord: Literal["HELIO","BARY"] = "HELIO", close_encounter_check: bool = True, general_relativity: bool = True, - fragmentation: bool = True, + fragmentation: bool = False, + minimum_fragment_mass: float | None = None, + minimum_fragment_gmass: float | None = None, rotation: bool = True, compute_conservation_values: bool = False, extra_force: bool = False, @@ -201,6 +202,12 @@ def __init__(self, fragmentation : bool, default True If set to True, this turns on the Fraggle fragment generation code and `rotation` must also be True. This argument only applies to Swiftest-SyMBA simulations. It will be ignored otherwise. + minimum_fragment_gmass : float, optional + If fragmentation is turned on, this sets the mimimum G*mass of a collisional fragment that can be generated. + *Note.* Only set one of minimum_fragment_gmass or minimum_fragment_mass + minimum_fragment_mass : float, optional + If fragmentation is turned on, this sets the mimimum mass of a collisional fragment that can be generated. + *Note.* Only set one of minimum_fragment_gmass or minimum_fragment_mass rotation : bool, default True If set to True, this turns on rotation tracking and radius, rotation vector, and moments of inertia values must be included in the initial conditions. @@ -284,6 +291,8 @@ def __init__(self, close_encounter_check=close_encounter_check, general_relativity=general_relativity, fragmentation=fragmentation, + minimum_fragment_gmass=minimum_fragment_gmass, + minimum_fragment_mass=minimum_fragment_mass, rotation=rotation, compute_conservation_values=compute_conservation_values, extra_force=extra_force, @@ -605,9 +614,11 @@ def set_integrator(self, Returns ------- - None + integrator_dict: dict + A dictionary containing the subset of the parameter dictonary that was updated by this setter """ + # TODO: Improve how it finds the executable binary if integrator is None and codename is None: return @@ -627,6 +638,13 @@ def set_integrator(self, self.param['! VERSION'] = f"{self.codename} parameter input" update_list.append("codename") + if self.codename == "Swiftest": + self.binary_path = os.path.realpath(os.path.join(os.path.dirname(os.path.realpath(_pyfile)),os.pardir,os.pardir,os.pardir,"build")) + self.driver_executable = os.path.join(self.binary_path,"swiftest_driver") + else: + self.binary_path = "NOT SET" + self.driver_executable = "NOT SET" + update_list.append("driver_executable") if integrator is not None: valid_integrator = ["symba","rmvs","whm","helio"] @@ -670,7 +688,8 @@ def get_integrator(self,arg_list: str | List[str] | None = None, verbose: bool | valid_var = {"codename": "! VERSION"} valid_instance_vars = {"integrator": self.integrator, - "codename": self.codename} + "codename": self.codename, + "driver_executable": self.driver_executable} try: self.integrator @@ -684,29 +703,25 @@ def get_integrator(self,arg_list: str | List[str] | None = None, verbose: bool | print(f"codename is not set") return {} + if verbose is None: + verbose = self.verbose if not bool(kwargs) and arg_list is None: arg_list = list(valid_instance_vars.keys()) + integrator = self._get_instance_var(arg_list, valid_instance_vars, verbose, **kwargs) valid_arg, integrator_dict = self._get_valid_arg_list(arg_list, valid_var) - if verbose is None: - verbose = self.verbose - if verbose: - for arg in arg_list: - if arg in valid_arg: - key = valid_var[arg] - print(f"{arg:<{self._getter_column_width}} {integrator_dict[key]}") - elif arg in valid_instance_vars: - print(f"{arg:<{self._getter_column_width}} {valid_instance_vars[arg]}") return integrator_dict def set_feature(self, close_encounter_check: bool | None = None, general_relativity: bool | None = None, fragmentation: bool | None = None, + minimum_fragment_gmass: float | None = None, + minimum_fragment_mass: float | None = None, rotation: bool | None = None, compute_conservation_values: bool | None = None, extra_force: bool | None = None, @@ -732,6 +747,12 @@ def set_feature(self, fragmentation : bool, optional If set to True, this turns on the Fraggle fragment generation code and `rotation` must also be True. This argument only applies to Swiftest-SyMBA simulations. It will be ignored otherwise. + minimum_fragment_gmass : float, optional + If fragmentation is turned on, this sets the mimimum G*mass of a collisional fragment that can be generated. + *Note.* Only set one of minimum_fragment_gmass or minimum_fragment_mass + minimum_fragment_mass : float, optional + If fragmentation is turned on, this sets the mimimum mass of a collisional fragment that can be generated. + *Note.* Only set one of minimum_fragment_gmass or minimum_fragment_mass rotation : bool, optional If set to True, this turns on rotation tracking and radius, rotation vector, and moments of inertia values must be included in the initial conditions. @@ -800,8 +821,27 @@ def set_feature(self, update_list.append("general_relativity") if fragmentation is not None: - self.param['FRAGMENTATION'] = fragmentation - update_list.append("fragmentation") + if self.codename != "Swiftest" and self.integrator != "symba" and fragmentation: + print("Fragmentation is only available on Swiftest SyMBA.") + self.param['FRAGMENTATION'] = False + else: + self.param['FRAGMENTATION'] = fragmentation + update_list.append("fragmentation") + if fragmentation: + if "MIN_GMFRAG" not in self.param and minimum_fragment_mass is None and minimum_fragment_gmass is None: + print("Minimum fragment mass is not set. Set it using minimum_fragment_gmass or minimum_fragment_mass") + else: + update_list.append("minimum_fragment_gmass") + update_list.append("minimum_fragment_mass") + if minimum_fragment_gmass is not None and minimum_fragment_mass is not None: + print("Warning! Only set either minimum_fragment_mass or minimum_fragment_gmass, but not both!") + + if minimum_fragment_gmass is not None: + self.param["MIN_GMFRAG"] = minimum_fragment_gmass + update_list.append("minimum_fragment_gmass") + elif minimum_fragment_mass is not None: + self.param["MIN_GMFRAG"] = minimum_fragment_mass * self.GU + update_list.append("minimum_fragment_gmass") if rotation is not None: self.param['ROTATION'] = rotation @@ -891,7 +931,8 @@ def get_feature(self, arg_list: str | List[str] | None = None, verbose: bool | N "rhill_present": "RHILL_PRESENT", "restart": "RESTART", "interaction_loops": "INTERACTION_LOOPS", - "encounter_check_loops": "ENCOUNTER_CHECK" + "encounter_check_loops": "ENCOUNTER_CHECK", + "minimum_fragment_gmass" : "MIN_GMFRAG" } valid_arg, feature_dict = self._get_valid_arg_list(arg_list, valid_var) @@ -902,7 +943,16 @@ def get_feature(self, arg_list: str | List[str] | None = None, verbose: bool | N if verbose: for arg in valid_arg: key = valid_var[arg] - print(f"{arg:<{self._getter_column_width}} {feature_dict[key]}") + if key in feature_dict: + if arg == "minimum_fragment_gmass": + if self.param['FRAGMENTATION']: + print(f"{arg:<{self._getter_column_width}} {feature_dict[key]} {self.DU_name}^3 / {self.TU_name}^2 ") + print(f"{'minimum_fragment_mass':<{self._getter_column_width}} {feature_dict[key] / self.GU} {self.MU_name}") + else: + print(f"{arg:<{self._getter_column_width}} {feature_dict[key]}") + else: + print(f"{arg:<{self._getter_column_width}} NOT SET") + return feature_dict