diff --git a/examples/Basic_Simulation/.swiftest/param.in b/examples/Basic_Simulation/.swiftest/param.in new file mode 100644 index 000000000..0ee870562 --- /dev/null +++ b/examples/Basic_Simulation/.swiftest/param.in @@ -0,0 +1,37 @@ +! VERSION Swiftest input file +T0 0.0 +TSTART 0.0 +TSTOP 10.0 +DT 0.005 +ISTEP_OUT 200 +ISTEP_DUMP 200 +NC_IN init_cond.nc +IN_TYPE NETCDF_DOUBLE +IN_FORM EL +BIN_OUT bin.nc +OUT_FORM XVEL +OUT_TYPE NETCDF_DOUBLE +OUT_STAT REPLACE +CHK_QMIN 0.004650467260962157 +CHK_RMIN 0.004650467260962157 +CHK_RMAX 10000.0 +CHK_EJECT 10000.0 +CHK_QMIN_COORD HELIO +CHK_QMIN_RANGE 0.004650467260962157 10000.0 +MU2KG 1.988409870698051e+30 +TU2S 31557600.0 +DU2M 149597870700.0 +GMTINY 9.869231602224408e-07 +RESTART NO +CHK_CLOSE YES +GR YES +FRAGMENTATION YES +ROTATION YES +ENERGY NO +EXTRA_FORCE NO +BIG_DISCARD NO +RHILL_PRESENT NO +INTERACTION_LOOPS TRIANGULAR +ENCOUNTER_CHECK TRIANGULAR +TIDES NO +MIN_GMFRAG 9.869231602224408e-10 diff --git a/examples/Basic_Simulation/initial_conditions.ipynb b/examples/Basic_Simulation/initial_conditions.ipynb index cef068fa0..384c55c1d 100644 --- a/examples/Basic_Simulation/initial_conditions.ipynb +++ b/examples/Basic_Simulation/initial_conditions.ipynb @@ -2,7 +2,7 @@ "cells": [ { "cell_type": "code", - "execution_count": null, + "execution_count": 1, "id": "2c4f59ea-1251-49f6-af1e-5695d7e25500", "metadata": {}, "outputs": [], @@ -14,10 +14,100 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 2, "id": "6054c7ab-c748-4b39-9fee-d8b27326f497", "metadata": {}, - "outputs": [], + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "codename Swiftest\n", + "integrator symba\n", + "param_file /home/daminton/git_debug/swiftest/examples/Basic_Simulation/.swiftest/param.in\n", + "driver_executable /home/daminton/git_debug/swiftest/bin/swiftest_driver\n", + "gmtiny 9.869231602224408e-07 AU^3 / y^2 \n", + "mtiny 2.5e-08 MSun\n", + "t0 0.0 y\n", + "tstart 0.0 y\n", + "tstop 10.0 y\n", + "dt 0.005 y\n", + "istep_out 200 \n", + "istep_dump 200 \n", + "tstep_out 1.0 y\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", + "restart REPLACE\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", + "fragmentation True\n", + "minimum_fragment_gmass 9.869231602224408e-10 AU^3 / y^2 \n", + "minimum_fragment_mass 2.5e-11 MSun\n", + "rotation True\n", + "general_relativity True\n", + "compute_conservation_values False\n", + "rhill_present False\n", + "extra_force False\n", + "big_discard False\n", + "interaction_loops TRIANGULAR\n", + "encounter_check_loops TRIANGULAR\n", + "restart False\n", + "ephemeris_date 2027-04-30\n" + ] + }, + { + "data": { + "text/plain": [ + "{'GMTINY': 9.869231602224408e-07,\n", + " 'T0': 0.0,\n", + " 'TSTART': 0.0,\n", + " 'TSTOP': 10.0,\n", + " 'DT': 0.005,\n", + " 'ISTEP_OUT': 200,\n", + " 'ISTEP_DUMP': 200,\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", + " 'OUT_STAT': 'REPLACE',\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", + " 'FRAGMENTATION': True,\n", + " 'MIN_GMFRAG': 9.869231602224408e-10,\n", + " 'ROTATION': True,\n", + " 'GR': True,\n", + " 'ENERGY': False,\n", + " 'RHILL_PRESENT': False,\n", + " 'EXTRA_FORCE': False,\n", + " 'BIG_DISCARD': False,\n", + " 'INTERACTION_LOOPS': 'TRIANGULAR',\n", + " 'ENCOUNTER_CHECK': 'TRIANGULAR',\n", + " 'RESTART': False}" + ] + }, + "execution_count": 2, + "metadata": {}, + "output_type": "execute_result" + } + ], "source": [ "# Initialize the simulation object as a variable\n", "sim = swiftest.Simulation(tstart=0.0, tstop=10.0, dt=0.005, tstep_out=1.0, fragmentation=True, minimum_fragment_mass = 2.5e-11, mtiny=2.5e-8)\n", @@ -26,10 +116,459 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 3, "id": "1c122676-bacb-447c-bc37-5ef8019be0d0", "metadata": {}, - "outputs": [], + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Creating the Sun as a central body\n", + "Fetching ephemerides data for Mercury from JPL/Horizons\n", + "Fetching ephemerides data for Venus from JPL/Horizons\n", + "Fetching ephemerides data for Earth from JPL/Horizons\n", + "Fetching ephemerides data for Mars from JPL/Horizons\n", + "Fetching ephemerides data for Jupiter from JPL/Horizons\n", + "Fetching ephemerides data for Saturn from JPL/Horizons\n", + "Fetching ephemerides data for Uranus from JPL/Horizons\n", + "Fetching ephemerides data for Neptune from JPL/Horizons\n", + "Fetching ephemerides data for Pluto from JPL/Horizons\n", + "Writing initial conditions to file /home/daminton/git_debug/swiftest/examples/Basic_Simulation/.swiftest/init_cond.nc\n", + "Writing parameter inputs to file /home/daminton/git_debug/swiftest/examples/Basic_Simulation/.swiftest/param.in\n" + ] + }, + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "
<xarray.Dataset>\n",
+       "Dimensions:        (name: 10, time: 1)\n",
+       "Coordinates:\n",
+       "  * name           (name) <U32 'Sun' 'Mercury' 'Venus' ... 'Neptune' 'Pluto'\n",
+       "  * time           (time) float64 0.0\n",
+       "Data variables: (12/20)\n",
+       "    particle_type  (name) <U32 'Central Body' 'Massive Body' ... 'Massive Body'\n",
+       "    id             (name) int64 0 1 2 3 4 5 6 7 8 9\n",
+       "    a              (time, name) float64 nan 0.3871 0.7233 ... 19.24 30.04 39.37\n",
+       "    e              (time, name) float64 nan 0.2056 0.006718 ... 0.008956 0.2487\n",
+       "    inc            (time, name) float64 nan 7.003 3.394 ... 0.773 1.771 17.17\n",
+       "    capom          (time, name) float64 nan 48.3 76.6 ... 74.01 131.8 110.3\n",
+       "    ...             ...\n",
+       "    roty           (time, name) float64 -38.76 -18.38 ... -2.177e+03 261.3\n",
+       "    rotz           (time, name) float64 82.25 34.36 8.703 ... 2.33e+03 -38.57\n",
+       "    j2rp2          (time, name) float64 4.754e-12 nan nan nan ... nan nan nan\n",
+       "    j4rp4          (time, name) float64 -2.247e-18 nan nan nan ... nan nan nan\n",
+       "    ntp            (time) int64 0\n",
+       "    npl            (time) int64 9
" + ], + "text/plain": [ + "\n", + "Dimensions: (name: 10, time: 1)\n", + "Coordinates:\n", + " * name (name) \n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "
<xarray.Dataset>\n",
+       "Dimensions:        (name: 5, time: 1)\n",
+       "Coordinates:\n",
+       "  * name           (name) <U14 'MassiveBody_01' ... 'MassiveBody_05'\n",
+       "  * time           (time) float64 0.0\n",
+       "Data variables: (12/18)\n",
+       "    particle_type  (name) <U14 'Massive Body' 'Massive Body' ... 'Massive Body'\n",
+       "    id             (name) int64 10 11 12 13 14\n",
+       "    a              (time, name) float64 1.311 0.4269 0.999 0.8757 0.6402\n",
+       "    e              (time, name) float64 0.2241 0.1367 0.1002 0.1051 0.07477\n",
+       "    inc            (time, name) float64 59.21 8.384 79.9 38.54 18.64\n",
+       "    capom          (time, name) float64 248.3 69.55 197.0 300.4 103.7\n",
+       "    ...             ...\n",
+       "    Ip3            (time, name) float64 0.4 0.4 0.4 0.4 0.4\n",
+       "    rotx           (time, name) float64 0.0 0.0 0.0 0.0 0.0\n",
+       "    roty           (time, name) float64 0.0 0.0 0.0 0.0 0.0\n",
+       "    rotz           (time, name) float64 0.0 0.0 0.0 0.0 0.0\n",
+       "    ntp            (time) int64 0\n",
+       "    npl            (time) int64 4
" + ], + "text/plain": [ + "\n", + "Dimensions: (name: 5, time: 1)\n", + "Coordinates:\n", + " * name (name) \n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "
<xarray.Dataset>\n",
+       "Dimensions:        (name: 10, time: 1)\n",
+       "Coordinates:\n",
+       "  * name           (name) <U15 'TestParticle_01' ... 'TestParticle_10'\n",
+       "  * time           (time) float64 0.0\n",
+       "Data variables:\n",
+       "    particle_type  (name) <U15 'Test Particle' ... 'Test Particle'\n",
+       "    id             (name) int64 15 16 17 18 19 20 21 22 23 24\n",
+       "    a              (time, name) float64 0.5963 0.7501 1.242 ... 0.6646 0.4335\n",
+       "    e              (time, name) float64 0.2692 0.1215 0.144 ... 0.08782 0.1085\n",
+       "    inc            (time, name) float64 31.01 21.79 71.4 ... 49.32 11.98 82.34\n",
+       "    capom          (time, name) float64 14.83 71.52 313.1 ... 55.01 264.6 27.41\n",
+       "    omega          (time, name) float64 263.9 53.24 270.9 ... 129.9 55.52 322.0\n",
+       "    capm           (time, name) float64 203.5 276.4 200.6 ... 296.5 199.3 11.76\n",
+       "    ntp            int64 10\n",
+       "    npl            int64 0
" + ], + "text/plain": [ + "\n", + "Dimensions: (name: 10, time: 1)\n", + "Coordinates:\n", + " * name (name) (tstop - tstart): - warnings.warn("dt must be smaller than tstop-tstart") - warnings.warn(f"Setting dt = {tstop - tstart} instead of {dt}") + msg = "dt must be smaller than tstop-tstart" + msg +=f"\nSetting dt = {tstop - tstart} instead of {dt}" + warnings.warn(msg,stacklevel=2) dt = tstop - tstart if dt is not None: @@ -535,7 +557,7 @@ def set_simulation_time(self, if istep_out is None and tstep_out is None: istep_out = self.param.pop("ISTEP_OUT", None) elif istep_out is not None and tstep_out is not None: - warnings.warn("istep_out and tstep_out cannot both be set") + warnings.warn("istep_out and tstep_out cannot both be set",stacklevel=2) return {} else: update_list.append("istep_out") @@ -643,7 +665,7 @@ def set_parameter(self, verbose: bool = True, **kwargs): default_arguments = { "codename" : "Swiftest", "integrator": "symba", - "param_file": "param.in", + "param_file": Path.cwd() / ".swiftest" / "param.in", "t0": 0.0, "tstart": 0.0, "tstop": None, @@ -696,9 +718,17 @@ def set_parameter(self, verbose: bool = True, **kwargs): param_file = kwargs.pop("param_file",None) + # Extract the simulation directory and create it if it doesn't exist if param_file is not None: - self.param_file = os.path.realpath(param_file) - self.sim_dir = os.path.dirname(self.param_file) + self.param_file = Path.cwd() / param_file + self.sim_dir = self.param_file.parent + if self.sim_dir.exists(): + if not self.sim_dir.is_dir(): + msg = f"Cannot create the {self.sim_dir} directory: File exists." + msg += "\nDelete the file or change the location of param_file" + warnings.warn(msg,stacklevel=2) + else: + self.sim_dir.mkdir(parents=True, exist_ok=False) # Setters returning parameter dictionary values @@ -783,7 +813,7 @@ def set_integrator(self, if codename is not None: valid_codename = ["Swiftest", "Swifter", "Swift"] if codename.title() not in valid_codename: - warnings.warn(f"{codename} is not a valid codename. Valid options are ",",".join(valid_codename)) + warnings.warn(f"{codename} is not a valid codename. Valid options are ",",".join(valid_codename),stacklevel=2) try: self.codename except: @@ -794,10 +824,10 @@ def set_integrator(self, self.param['! VERSION'] = f"{self.codename} input file" 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,"bin")) - self.driver_executable = os.path.join(self.binary_path,"swiftest_driver") - if not os.path.exists(self.driver_executable): - warnings.warn(f"Cannot find the Swiftest driver in {self.binary_path}") + self.binary_path = Path(_pyfile).parent.parent.parent.parent / "bin" + self.driver_executable = self.binary_path / "swiftest_driver" + if not self.driver_executable.exists(): + warnings.warn(f"Cannot find the Swiftest driver in {str(self.binary_path)}",stacklevel=2) self.driver_executable = None else: self.binary_path = "NOT IMPLEMENTED FOR THIS CODE" @@ -807,7 +837,7 @@ def set_integrator(self, if integrator is not None: valid_integrator = ["symba","rmvs","whm","helio"] if integrator.lower() not in valid_integrator: - warnings.warn(f"{integrator} is not a valid integrator. Valid options are ",",".join(valid_integrator)) + warnings.warn(f"{integrator} is not a valid integrator. Valid options are ",",".join(valid_integrator),stacklevel=2) try: self.integrator except: @@ -818,9 +848,9 @@ def set_integrator(self, if mtiny is not None or gmtiny is not None: if self.integrator != "symba": - warnings.warn("mtiny and gmtiny are only used by SyMBA.") + warnings.warn("mtiny and gmtiny are only used by SyMBA.",stacklevel=2) if mtiny is not None and gmtiny is not None: - warnings.warn("Only set mtiny or gmtiny, not both!") + warnings.warn("Only set mtiny or gmtiny, not both.",stacklevel=2) elif gmtiny is not None: self.param['GMTINY'] = gmtiny update_list.append("gmtiny") @@ -859,19 +889,19 @@ def get_integrator(self,arg_list: str | List[str] | None = None, verbose: bool | valid_instance_vars = {"codename": self.codename, "integrator": self.integrator, - "param_file": self.param_file, - "driver_executable": self.driver_executable} + "param_file": str(self.param_file), + "driver_executable": str(self.driver_executable)} try: self.integrator except: - warnings.warn(f"integrator is not set") + warnings.warn(f"integrator is not set",stacklevel=2) return {} try: self.codename except: - warnings.warn(f"codename is not set") + warnings.warn(f"codename is not set",stacklevel=2) return {} if verbose is None: @@ -1008,19 +1038,19 @@ def set_feature(self, if fragmentation is not None: if self.codename != "Swiftest" and self.integrator != "symba" and fragmentation: - warnings.warn("Fragmentation is only available on Swiftest SyMBA.") + warnings.warn("Fragmentation is only available on Swiftest SyMBA.",stacklevel=2) 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: - warnings.warn("Minimum fragment mass is not set. Set it using minimum_fragment_gmass or minimum_fragment_mass") + warnings.warn("Minimum fragment mass is not set. Set it using minimum_fragment_gmass or minimum_fragment_mass",stacklevel=2) else: update_list.append("minimum_fragment_gmass") if minimum_fragment_gmass is not None and minimum_fragment_mass is not None: - warnings.warn("Only set either minimum_fragment_mass or minimum_fragment_gmass, but not both!") + warnings.warn("Only set either minimum_fragment_mass or minimum_fragment_gmass, but not both!",stacklevel=2) if minimum_fragment_gmass is not None: self.param["MIN_GMFRAG"] = minimum_fragment_gmass @@ -1065,8 +1095,9 @@ def set_feature(self, if interaction_loops is not None: valid_vals = ["TRIANGULAR", "FLAT", "ADAPTIVE"] if interaction_loops not in valid_vals: - warnings.warn(f"{interaction_loops} is not a valid option for interaction loops.") - warnings.warn(f"Must be one of {valid_vals}") + msg = f"{interaction_loops} is not a valid option for interaction loops." + msg += f"\nMust be one of {valid_vals}" + warnings.warn(msg,stacklevel=2) if "INTERACTION_LOOPS" not in self.param: self.param["INTERACTION_LOOPS"] = valid_vals[0] else: @@ -1076,8 +1107,9 @@ def set_feature(self, if encounter_check_loops is not None: valid_vals = ["TRIANGULAR", "SORTSWEEP", "ADAPTIVE"] if encounter_check_loops not in valid_vals: - warnings.warn(f"{encounter_check_loops} is not a valid option for interaction loops.") - warnings.warn(f"Must be one of {valid_vals}") + msg = f"{encounter_check_loops} is not a valid option for interaction loops." + msg += f"\nMust be one of {valid_vals}" + warnings.warn(msg,stacklevel=2) if "ENCOUNTER_CHECK" not in self.param: self.param["ENCOUNTER_CHECK"] = valid_vals[0] else: @@ -1208,13 +1240,14 @@ def set_init_cond_files(self, return {} def ascii_file_input_error_msg(codename): - warnings.warn(f"in set_init_cond_files: init_cond_file_name must be a dictionary of the form: ") - warnings.warn('{') + msg = f"in set_init_cond_files: init_cond_file_name must be a dictionary of the form: " + msg += "\n {" if codename == "Swiftest": - warnings.warn('"CB" : *path to central body initial conditions file*,') - warnings.warn('"PL" : *path to massive body initial conditions file*,') - warnings.warn('"TP" : *path to test particle initial conditions file*') - warnings.warn('}') + msg += '\n"CB" : *path to central body initial conditions file*,' + msg += '\n"PL" : *path to massive body initial conditions file*,' + msg += '\n"TP" : *path to test particle initial conditions file*' + msg += '\n}' + warnings.warn(msg,stacklevel=2) return {} if init_cond_format is None: @@ -1234,21 +1267,21 @@ def ascii_file_input_error_msg(codename): else: init_cond_keys = ["PL", "TP"] if init_cond_file_type != "ASCII": - warnings.warn(f"{init_cond_file_type} is not supported by {self.codename}. Using ASCII instead") + warnings.warn(f"{init_cond_file_type} is not supported by {self.codename}. Using ASCII instead",stacklevel=2) init_cond_file_type = "ASCII" if init_cond_format != "XV": - warnings.warn(f"{init_cond_format} is not supported by {self.codename}. Using XV instead") + warnings.warn(f"{init_cond_format} is not supported by {self.codename}. Using XV instead",stacklevel=2) init_cond_format = "XV" valid_formats = {"EL", "XV"} if init_cond_format not in valid_formats: - warnings.warn(f"{init_cond_format} is not a valid input format") + warnings.warn(f"{init_cond_format} is not a valid input format",stacklevel=2) else: self.param['IN_FORM'] = init_cond_format valid_types = {"NETCDF_DOUBLE", "NETCDF_FLOAT", "ASCII"} if init_cond_file_type not in valid_types: - warnings.warn(f"{init_cond_file_type} is not a valid input type") + warnings.warn(f"{init_cond_file_type} is not a valid input type",stackevel=2) else: self.param['IN_TYPE'] = init_cond_file_type @@ -1275,7 +1308,7 @@ def ascii_file_input_error_msg(codename): elif type(init_cond_file_name) is dict: # Oops, accidentally passed a dictionary instead of the expected single string or path-like for NetCDF # input type. - warnings.warn(f"Only a single input file is used for NetCDF files") + warnings.warn(f"Only a single input file is used for NetCDF files",stacklevel=2) else: self.param["NC_IN"] = init_cond_file_name @@ -1416,7 +1449,7 @@ def set_output_files(self, if output_file_type is None: output_file_type = "NETCDF_DOUBLE" elif output_file_type not in ["NETCDF_DOUBLE", "NETCDF_FLOAT"]: - warnings.warn(f"{output_file_type} is not compatible with Swiftest. Setting to NETCDF_DOUBLE") + warnings.warn(f"{output_file_type} is not compatible with Swiftest. Setting to NETCDF_DOUBLE",stacklevel=2) output_file_type = "NETCDF_DOUBLE" elif self.codename == "Swifter": if output_file_type is None: @@ -1424,7 +1457,7 @@ def set_output_files(self, if output_file_type is None: output_file_type = "REAL8" elif output_file_type not in ["REAL4", "REAL8", "XDR4", "XDR8"]: - warnings.warn(f"{output_file_type} is not compatible with Swifter. Setting to REAL8") + warnings.warn(f"{output_file_type} is not compatible with Swifter. Setting to REAL8",stacklevel=2) output_file_type = "REAL8" elif self.codename == "Swift": if output_file_type is None: @@ -1432,7 +1465,7 @@ def set_output_files(self, if output_file_type is None: output_file_type = "REAL4" if output_file_type not in ["REAL4"]: - warnings.warn(f"{output_file_type} is not compatible with Swift. Setting to REAL4") + warnings.warn(f"{output_file_type} is not compatible with Swift. Setting to REAL4",stacklevel=2) output_file_type = "REAL4" self.param['OUT_TYPE'] = output_file_type @@ -1445,7 +1478,7 @@ def set_output_files(self, self.param['BIN_OUT'] = output_file_name if output_format != "XV" and self.codename != "Swiftest": - warnings.warn(f"{output_format} is not compatible with {self.codename}. Setting to XV") + warnings.warn(f"{output_format} is not compatible with {self.codename}. Setting to XV",stacklevel=2) output_format = "XV" self.param["OUT_FORM"] = output_format @@ -1620,7 +1653,7 @@ def set_unit_system(self, self.param['MU2KG'] = 1000.0 self.MU_name = "g" else: - warnings.warn(f"{MU} not a recognized unit system. Using MSun as a default.") + warnings.warn(f"{MU} not a recognized unit system. Using MSun as a default.",stacklevel=2) self.param['MU2KG'] = constants.MSun self.MU_name = "MSun" @@ -1643,7 +1676,7 @@ def set_unit_system(self, self.param['DU2M'] = 100.0 self.DU_name = "cm" else: - warnings.warn(f"{DU} not a recognized unit system. Using AU as a default.") + warnings.warn(f"{DU} not a recognized unit system. Using AU as a default.",stacklevel=2) self.param['DU2M'] = constants.AU2M self.DU_name = "AU" @@ -1663,7 +1696,7 @@ def set_unit_system(self, self.param['TU2S'] = 1.0 self.TU_name = "s" else: - warnings.warn(f"{TU} not a recognized unit system. Using YR as a default.") + warnings.warn(f"{TU} not a recognized unit system. Using YR as a default.",stacklevel=2) self.param['TU2S'] = constants.YR2S self.TU_name = "y" @@ -1846,7 +1879,7 @@ def set_distance_range(self, if qmin_coord is not None: valid_qmin_coord = ["HELIO","BARY"] if qmin_coord.upper() not in valid_qmin_coord: - warnings.warn(f"qmin_coord = {qmin_coord} is not a valid option. Must be one of",','.join(valid_qmin_coord)) + warnings.warn(f"qmin_coord = {qmin_coord} is not a valid option. Must be one of",','.join(valid_qmin_coord),stacklevel=2) self.param['CHK_QMIN_COORD'] = valid_qmin_coord[0] else: self.param['CHK_QMIN_COORD'] = qmin_coord.upper() @@ -1969,7 +2002,7 @@ def add_solar_system_body(self, if type(ephemeris_id) is int: ephemeris_id = [ephemeris_id] if len(ephemeris_id) != len(name): - warnings.warn(f"The length of ephemeris_id ({len(ephemeris_id)}) does not match the length of name ({len(name)})") + warnings.warn(f"The length of ephemeris_id ({len(ephemeris_id)}) does not match the length of name ({len(name)})",stacklevel=2) return None else: ephemeris_id = [None] * len(name) @@ -1982,11 +2015,11 @@ def add_solar_system_body(self, try: datetime.datetime.fromisoformat(date) except: - warnings.warn(f"{date} is not a valid date format. Must be 'YYYY-MM-DD'. Setting to {self.ephemeris_date}") + warnings.warn(f"{date} is not a valid date format. Must be 'YYYY-MM-DD'. Setting to {self.ephemeris_date}",stacklevel=2) date = self.ephemeris_date if source.upper() != "HORIZONS": - warnings.warn("Currently only the JPL Horizons ephemeris service is supported") + warnings.warn("Currently only the JPL Horizons ephemeris service is supported",stacklevel=2) body_list = [] for i,n in enumerate(name): @@ -2091,8 +2124,9 @@ def set_ephemeris_date(self, datetime.datetime.fromisoformat(ephemeris_date) except: valid_date_args = ['"MBCL"', '"TODAY"', '"YYYY-MM-DD"'] - warnings.warn(f"{ephemeris_date} is not a valid format. Valid options include:", ', '.join(valid_date_args)) - warnings.warn("Using MBCL for date.") + msg = f"{ephemeris_date} is not a valid format. Valid options include:", ', '.join(valid_date_args) + msg += "\nUsing MBCL for date." + warnings.warn(msg,stacklevel=2) ephemeris_date = minton_bcl self.ephemeris_date = ephemeris_date @@ -2127,7 +2161,7 @@ def get_ephemeris_date(self, arg_list: str | List[str] | None = None, verbose: b try: self.ephemeris_date except: - warnings.warn(f"ephemeris_date is not set") + warnings.warn(f"ephemeris_date is not set",stacklevel=2) return valid_arg = {"ephemeris_date": self.ephemeris_date} @@ -2322,8 +2356,9 @@ def _combine_and_fix_dsnew(self,dsnew): dsnew = dsnew.swap_dims({"id" : "name"}) dsnew = dsnew.reset_coords("id") else: - warnings.warn("Non-unique names detected for bodies. The Dataset will be dimensioned by integer id instead of name.") - warnings.warn("Consider using unique names instead.") + msg = "Non-unique names detected for bodies. The Dataset will be dimensioned by integer id instead of name." + msg +="\nConsider using unique names instead." + print(msg) self.data = xr.combine_by_coords([self.data, dsnew]) @@ -2393,7 +2428,7 @@ def read_param(self, elif codename == "Swift": self.param = io.read_swift_param(param_file, verbose=verbose) else: - warnings.warn(f'{codename} is not a recognized code name. Valid options are "Swiftest", "Swifter", or "Swift".') + warnings.warn(f'{codename} is not a recognized code name. Valid options are "Swiftest", "Swifter", or "Swift".',stacklevel=2) return False return True @@ -2444,7 +2479,7 @@ def write_param(self, elif codename == "Swift": io.write_swift_param(param, param_file) else: - warnings.warn('Cannot process unknown code type. Call the read_param method with a valid code name. Valid options are "Swiftest", "Swifter", or "Swift".') + warnings.warn('Cannot process unknown code type. Call the read_param method with a valid code name. Valid options are "Swiftest", "Swifter", or "Swift".',stacklevel=2) return def convert(self, param_file, newcodename="Swiftest", plname="pl.swiftest.in", tpname="tp.swiftest.in", @@ -2474,10 +2509,10 @@ def convert(self, param_file, newcodename="Swiftest", plname="pl.swiftest.in", t """ oldparam = self.param if self.codename == newcodename: - warnings.warn(f"This parameter configuration is already in {newcodename} format") + warnings.warn(f"This parameter configuration is already in {newcodename} format",stacklevel=2) return oldparam if newcodename != "Swift" and newcodename != "Swifter" and newcodename != "Swiftest": - warnings.warn(f'{newcodename} is an invalid code type. Valid options are "Swiftest", "Swifter", or "Swift".') + warnings.warn(f'{newcodename} is an invalid code type. Valid options are "Swiftest", "Swifter", or "Swift".',stacklevel=2) return oldparam goodconversion = True if self.codename == "Swifter": @@ -2498,7 +2533,7 @@ def convert(self, param_file, newcodename="Swiftest", plname="pl.swiftest.in", t if goodconversion: self.write_param(param_file) else: - warnings.warn(f"Conversion from {self.codename} to {newcodename} is not supported.") + warnings.warn(f"Conversion from {self.codename} to {newcodename} is not supported.",stacklevel=2) return oldparam def bin2xr(self): @@ -2525,9 +2560,9 @@ def bin2xr(self): self.data = io.swifter2xr(param_tmp, verbose=self.verbose) if self.verbose: print('Swifter simulation data stored as xarray DataSet .data') elif self.codename == "Swift": - warnings.warn("Reading Swift simulation data is not implemented yet") + warnings.warn("Reading Swift simulation data is not implemented yet",stacklevel=2) else: - warnings.warn('Cannot process unknown code type. Call the read_param method with a valid code name. Valid options are "Swiftest", "Swifter", or "Swift".') + warnings.warn('Cannot process unknown code type. Call the read_param method with a valid code name. Valid options are "Swiftest", "Swifter", or "Swift".',stacklevel=2) return def follow(self, codestyle="Swifter"): @@ -2558,7 +2593,7 @@ def follow(self, codestyle="Swifter"): i_list = [i for i in line.split(" ") if i.strip()] nskp = int(i_list[0]) except IOError: - warnings.warn('No follow.in file found') + warnings.warn('No follow.in file found',stacklevel=2) ifol = None nskp = None fol = tool.follow_swift(self.data, ifol=ifol, nskp=nskp) @@ -2603,7 +2638,8 @@ def save(self, param = self.param if codename == "Swiftest": - io.swiftest_xr2infile(ds=self.data, param=param, in_type=self.param['IN_TYPE'], framenum=framenum) + infile_name = Path(self.sim_dir) / param['NC_IN'] + io.swiftest_xr2infile(ds=self.data, param=param, in_type=self.param['IN_TYPE'], infile_name=infile_name, framenum=framenum) self.write_param(param_file=param_file) elif codename == "Swifter": if codename == "Swiftest": @@ -2613,7 +2649,7 @@ def save(self, io.swifter_xr2infile(self.data, swifter_param, framenum) self.write_param(param_file, param=swifter_param) else: - warnings.warn(f'Saving to {codename} not supported') + warnings.warn(f'Saving to {codename} not supported',stacklevel=2) return @@ -2658,7 +2694,7 @@ def initial_conditions_from_bin(self, framenum=-1, new_param=None, new_param_fil elif self.param['OUT_TYPE'] == 'NETCDF_FLOAT': new_param['IN_TYPE'] = 'NETCDF_FLOAT' else: - warnings.warn(f"{self.param['OUT_TYPE']} is an invalid OUT_TYPE file") + warnings.warn(f"{self.param['OUT_TYPE']} is an invalid OUT_TYPE file",stacklevel=2) return if self.param['BIN_OUT'] != new_param['BIN_OUT'] and restart: diff --git a/src/io/io.f90 b/src/io/io.f90 index 86cf55728..08ab7f868 100644 --- a/src/io/io.f90 +++ b/src/io/io.f90 @@ -9,8 +9,80 @@ submodule (swiftest_classes) s_io use swiftest + contains + module subroutine io_pbar_reset(self, nloops) + !! author: David A. Minton + !! + !! Resets the progress bar to the beginning + implicit none + ! Arguments + class(progress_bar),intent(inout) :: self + integer(I8B), intent(in) :: nloops + ! Internals + character(len=2) :: numchar + character(len=5) :: barfmt + + if (.not.allocated(self%barstr)) then + allocate(character(self%PBARSIZE+4) :: self%barstr) + end if + write(numchar,'(I2)') self%PBARSIZE + self%fmt = '(A1,"[",A' // numchar // ',"]",$)' + barfmt = '(A' // numchar // ')' + write(self%barstr,barfmt) " " + self%nloops = nloops + self%spinner = 0 + + call self%update(0) + + return + end subroutine io_pbar_reset + + module subroutine io_pbar_update(self,i) + !! author: David A. Minton + !! + !! Updates the progress bar with new values and causes the "spinner" to flip. + implicit none + ! Arguments + class(progress_bar), intent(inout) :: self + integer(I8B), intent(in) :: i + ! Internals + integer(I4B) :: k + real(DP) :: frac + integer(I4B) :: pos !! The current integer position of the progress bar + + ! Compute the current position + frac = real(i,kind=DP) / real(self%nloops,kind=DP) + pos = 1 + min(int(ceiling(frac * self%PBARSIZE),kind=I4B),self%PBARSIZE) + + ! Fill in the bar character up to the current position + do k = 1, pos + self%barstr(k:k) = self%barchar + end do + + self%spinner = self%spinner + 1 + if (self%spinner > 4) self%spinner = 1 + select case(self%spinner) + case(1) + self%barstr(pos:pos) = "/" + case(2) + self%barstr(pos:pos) = "-" + case(3) + self%barstr(pos:pos) = "\" + case(4) + self%barstr(pos:pos) = "|" + end select + + write(*,fmt=self%fmt) char(13),self%barstr + + + return + end subroutine io_pbar_update + + + + module subroutine io_conservation_report(self, param, lterminal) !! author: The Purdue Swiftest Team - David A. Minton, Carlisle A. Wishard, Jennifer L.L. Pouplin, and Jacob R. Elliott !! diff --git a/src/main/swiftest_driver.f90 b/src/main/swiftest_driver.f90 index 467403269..0eca02509 100644 --- a/src/main/swiftest_driver.f90 +++ b/src/main/swiftest_driver.f90 @@ -31,6 +31,7 @@ program swiftest_driver real(DP) :: old_t_final = 0.0_DP !! Output time at which writing should start, in order to prevent duplicate lines being written for restarts type(walltimer) :: integration_timer !! Object used for computing elapsed wall time real(DP) :: tfrac + type(progress_bar) :: pbar !! Object used to print out a progress bar character(*), parameter :: statusfmt = '("Time = ", ES12.5, "; fraction done = ", F6.3, ' // & '"; Number of active pl, tp = ", I5, ", ", I5)' character(*), parameter :: symbastatfmt = '("Time = ", ES12.5, "; fraction done = ", F6.3, ' // & @@ -88,6 +89,7 @@ program swiftest_driver !$ write(*,'(a,i3,/)') ' Number of threads = ', nthreads write(*, *) " *************** Main Loop *************** " if (param%lrestart .and. param%lenergy) call nbody_system%conservation_report(param, lterminal=.true.) + call pbar%reset(nloops) do iloop = 1, nloops !> Step the system forward in time call integration_timer%start() @@ -100,6 +102,7 @@ program swiftest_driver call nbody_system%discard(param) !> If the loop counter is at the output cadence value, append the data file with a single frame + call pbar%update(iloop) if (istep_out > 0) then iout = iout - 1 if (iout == 0) then @@ -108,15 +111,15 @@ program swiftest_driver tfrac = (param%t - param%t0) / (param%tstop - param%t0) - select type(pl => nbody_system%pl) - class is (symba_pl) - write(*, symbastatfmt) param%t, tfrac, pl%nplm, pl%nbody, nbody_system%tp%nbody - class default - write(*, statusfmt) param%t, tfrac, pl%nbody, nbody_system%tp%nbody - end select + ! select type(pl => nbody_system%pl) + ! class is (symba_pl) + ! write(*, symbastatfmt) param%t, tfrac, pl%nplm, pl%nbody, nbody_system%tp%nbody + ! class default + ! write(*, statusfmt) param%t, tfrac, pl%nbody, nbody_system%tp%nbody + ! end select if (param%lenergy) call nbody_system%conservation_report(param, lterminal=.true.) - call integration_timer%report(message="Integration steps:", nsubsteps=istep_out) - call integration_timer%reset() + !call integration_timer%report(message="Integration steps:", nsubsteps=istep_out) + !call integration_timer%reset() iout = istep_out end if diff --git a/src/modules/swiftest_classes.f90 b/src/modules/swiftest_classes.f90 index 6c0ac2be3..49735a424 100644 --- a/src/modules/swiftest_classes.f90 +++ b/src/modules/swiftest_classes.f90 @@ -414,6 +414,21 @@ module swiftest_classes generic :: read_particle_info => read_particle_info_bin, read_particle_info_netcdf !! Genereric method call for reading in the particle information metadata end type swiftest_nbody_system + type :: progress_bar + !! author: David A. Minton + !! + !! Implements a class for a simple progress bar that can print on the screen. + integer(I4B) :: PBARSIZE = 80 !! Number of characters acros for a whole progress bar + integer(I8B) :: nloops !! The total number of loops that the progrees bar is executing + character(len=:), allocatable :: barstr !! The string that prints out as the progress bar + integer(I4B) :: spinner !! Position of the "spinner" that indicates that progress is being made + character(len=1) :: barchar = "=" !! The progress bar character + character(len=32) :: fmt !! The format string that is used to define the progress bar itself + contains + procedure :: reset => io_pbar_reset !! Resets the progress bar to the beginning + procedure :: update => io_pbar_update !! Updates the progress bar with new values and causes the "spinner" to flip. + end type progress_bar + abstract interface subroutine abstract_accel(self, system, param, t, lbeg) @@ -582,6 +597,18 @@ pure module subroutine gr_vh2pv_body(self, param) class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters end subroutine gr_vh2pv_body + module subroutine io_pbar_reset(self, nloops) + implicit none + class(progress_bar),intent(inout) :: self + integer(I8B), intent(in) :: nloops + end subroutine io_pbar_reset + + module subroutine io_pbar_update(self,i) + implicit none + class(progress_bar), intent(inout) :: self + integer(I8B), intent(in) :: i + end subroutine io_pbar_update + module subroutine io_conservation_report(self, param, lterminal) implicit none class(swiftest_nbody_system), intent(inout) :: self !! Swiftest nbody system object