diff --git a/examples/.gitignore b/examples/.gitignore new file mode 100644 index 000000000..643407776 --- /dev/null +++ b/examples/.gitignore @@ -0,0 +1 @@ +*/param.* diff --git a/examples/Basic_Simulation/initial_conditions.ipynb b/examples/Basic_Simulation/initial_conditions.ipynb new file mode 100644 index 000000000..cef068fa0 --- /dev/null +++ b/examples/Basic_Simulation/initial_conditions.ipynb @@ -0,0 +1,155 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": null, + "id": "2c4f59ea-1251-49f6-af1e-5695d7e25500", + "metadata": {}, + "outputs": [], + "source": [ + "import swiftest\n", + "import numpy as np\n", + "from numpy.random import default_rng" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6054c7ab-c748-4b39-9fee-d8b27326f497", + "metadata": {}, + "outputs": [], + "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", + "sim.get_parameter()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "1c122676-bacb-447c-bc37-5ef8019be0d0", + "metadata": {}, + "outputs": [], + "source": [ + "# Add the modern planets and the Sun using the JPL Horizons Database\n", + "sim.add_solar_system_body([\"Sun\",\"Mercury\",\"Venus\",\"Earth\",\"Mars\",\"Jupiter\",\"Saturn\",\"Uranus\",\"Neptune\",\"Pluto\"])" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "97fe9f16-bc2e-443c-856b-7dacb1267f2d", + "metadata": {}, + "outputs": [], + "source": [ + "# Add 5 user-defined massive bodies\n", + "npl = 5\n", + "density_pl = 3000.0 / (sim.param['MU2KG'] / sim.param['DU2M'] ** 3)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "566a742e-3935-484d-9c7a-87310b6aaa3a", + "metadata": {}, + "outputs": [], + "source": [ + "name_pl = [\"MassiveBody_01\", \"MassiveBody_02\", \"MassiveBody_03\", \"MassiveBody_04\", \"MassiveBody_05\"]\n", + "a_pl = default_rng().uniform(0.3, 1.5, npl)\n", + "e_pl = default_rng().uniform(0.0, 0.3, npl)\n", + "inc_pl = default_rng().uniform(0.0, 90, npl)\n", + "capom_pl = default_rng().uniform(0.0, 360.0, npl)\n", + "omega_pl = default_rng().uniform(0.0, 360.0, npl)\n", + "capm_pl = default_rng().uniform(0.0, 360.0, npl)\n", + "GM_pl = (np.array([6e23, 8e23, 1e24, 3e24, 5e24]) / sim.param['MU2KG']) * sim.GU\n", + "R_pl = np.full(npl, (3 * (GM_pl / sim.GU) / (4 * np.pi * density_pl)) ** (1.0 / 3.0))\n", + "Rh_pl = a_pl * ((GM_pl) / (3 * sim.GU)) ** (1.0 / 3.0)\n", + "Ip1_pl = [0.4, 0.4, 0.4, 0.4, 0.4]\n", + "Ip2_pl = [0.4, 0.4, 0.4, 0.4, 0.4]\n", + "Ip3_pl = [0.4, 0.4, 0.4, 0.4, 0.4]\n", + "rotx_pl = [0.0, 0.0, 0.0, 0.0, 0.0]\n", + "roty_pl = [0.0, 0.0, 0.0, 0.0, 0.0]\n", + "rotz_pl = [0.0, 0.0, 0.0, 0.0, 0.0]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d2c14121-e5b4-45bd-99a7-acde0ee77955", + "metadata": {}, + "outputs": [], + "source": [ + "sim.add_body(name_pl, a_pl, e_pl, inc_pl, capom_pl, omega_pl, capm_pl, GMpl=GM_pl, Rpl=R_pl, rhill=Rh_pl, Ip1=Ip1_pl, Ip2=Ip2_pl, Ip3=Ip3_pl, rotx=rotx_pl, roty=roty_pl, rotz=rotz_pl)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b3f979b3-5238-492f-8589-0cf5d9a3c2bc", + "metadata": {}, + "outputs": [], + "source": [ + "# Add 10 user-defined test particles\n", + "ntp = 10" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5777c1bb-8c91-452a-8869-ce6f951b33a3", + "metadata": {}, + "outputs": [], + "source": [ + "name_tp = [\"TestParticle_01\", \"TestParticle_02\", \"TestParticle_03\", \"TestParticle_04\", \"TestParticle_05\", \"TestParticle_06\", \"TestParticle_07\", \"TestParticle_08\", \"TestParticle_09\", \"TestParticle_10\"]\n", + "a_tp = default_rng().uniform(0.3, 1.5, ntp)\n", + "e_tp = default_rng().uniform(0.0, 0.3, ntp)\n", + "inc_tp = default_rng().uniform(0.0, 90, ntp)\n", + "capom_tp = default_rng().uniform(0.0, 360.0, ntp)\n", + "omega_tp = default_rng().uniform(0.0, 360.0, ntp)\n", + "capm_tp = default_rng().uniform(0.0, 360.0, ntp)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8df4bcb7-46a6-402a-b10f-f37c55a0fdee", + "metadata": {}, + "outputs": [], + "source": [ + "sim.add_body(name_tp, a_tp, e_tp, inc_tp, capom_tp, omega_tp, capm_tp)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "fd30638b-6868-4162-bc05-1011bb255968", + "metadata": {}, + "outputs": [], + "source": [ + "# Run the simulation\n", + "sim.run()" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "swiftest", + "language": "python", + "name": "swiftest" + }, + "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.8.5" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/examples/Basic_Simulation/param.in b/examples/Basic_Simulation/param.in deleted file mode 100644 index 0ee870562..000000000 --- a/examples/Basic_Simulation/param.in +++ /dev/null @@ -1,37 +0,0 @@ -! 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/python/swiftest/swiftest/simulation_class.py b/python/swiftest/swiftest/simulation_class.py index c6fe5e8da..c58d68db8 100644 --- a/python/swiftest/swiftest/simulation_class.py +++ b/python/swiftest/swiftest/simulation_class.py @@ -2325,12 +2325,14 @@ def _combine_and_fix_dsnew(self,dsnew): 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.") + self.data = xr.combine_by_coords([self.data, dsnew]) + if self.param['OUT_TYPE'] == "NETCDF_DOUBLE": dsnew = io.fix_types(dsnew, ftype=np.float64) + self.data = io.fix_types(self.data, ftype=np.float64) elif self.param['OUT_TYPE'] == "NETCDF_FLOAT": dsnew = io.fix_types(dsnew, ftype=np.float32) - - self.data = xr.combine_by_coords([self.data, dsnew]) + self.data = io.fix_types(self.data, ftype=np.float32) def get_nvals(ds): if "name" in ds.dims: @@ -2348,6 +2350,8 @@ def get_nvals(ds): dsnew = get_nvals(dsnew) self.data = get_nvals(self.data) + self.data = self.data.sortby("id") + return dsnew def read_param(self,