diff --git a/examples/Basic_Simulation/initial_conditions.ipynb b/examples/Basic_Simulation/initial_conditions.ipynb index cdcaae0ed..17b85582f 100644 --- a/examples/Basic_Simulation/initial_conditions.ipynb +++ b/examples/Basic_Simulation/initial_conditions.ipynb @@ -2,10 +2,18 @@ "cells": [ { "cell_type": "code", - "execution_count": null, + "execution_count": 1, "id": "2c4f59ea-1251-49f6-af1e-5695d7e25500", "metadata": {}, - "outputs": [], + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "env: OMP_NUM_THREADS=4\n" + ] + } + ], "source": [ "import swiftest\n", "import numpy as np\n", @@ -15,21 +23,468 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 2, "id": "6054c7ab-c748-4b39-9fee-d8b27326f497", "metadata": {}, "outputs": [], "source": [ "# Initialize the simulation object as a variable\n", - "sim = swiftest.Simulation(tstart=0.0, tstop=1.0e6, dt=0.01, tstep_out=1.0e3, fragmentation=True, minimum_fragment_mass = 2.5e-11, mtiny=2.5e-8)" + "sim = swiftest.Simulation(tstart=0.0, tstop=1.0e3, dt=0.01, tstep_out=1.0e0, dump_cadence=2, fragmentation=True, minimum_fragment_mass = 2.5e-11, mtiny=2.5e-8)" ] }, { "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" + ] + }, + { + "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/21)\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",
+       "    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\n",
+       "    nplm           (time) int64 8
" + ], + "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/19)\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.469 0.4169 1.369 0.6314 0.4806\n",
+       "    e              (time, name) float64 0.1092 0.03191 0.03574 0.03611 0.2767\n",
+       "    inc            (time, name) float64 0.2741 70.11 62.39 31.73 47.9\n",
+       "    capom          (time, name) float64 123.3 146.2 205.2 41.36 298.9\n",
+       "    ...             ...\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\n",
+       "    nplm           (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.7527 1.445 0.8756 ... 1.341 0.9409\n",
+       "    e              (time, name) float64 0.267 0.0711 0.04515 ... 0.1502 0.06409\n",
+       "    inc            (time, name) float64 58.34 7.109 33.64 ... 52.18 26.94 7.888\n",
+       "    capom          (time, name) float64 130.7 145.3 68.94 ... 131.8 140.6 81.53\n",
+       "    omega          (time, name) float64 144.5 215.6 104.4 ... 288.9 84.92 180.3\n",
+       "    capm           (time, name) float64 55.73 338.2 71.69 ... 239.2 311.4 187.1\n",
+       "    ntp            int64 10\n",
+       "    npl            int64 0\n",
+       "    nplm           int64 0
" + ], + "text/plain": [ + "\n", + "Dimensions: (name: 10, time: 1)\n", + "Coordinates:\n", + " * name (name) 2\u001b[0m \u001b[43msim\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mrun\u001b[49m\u001b[43m(\u001b[49m\u001b[43m)\u001b[49m\n", + "File \u001b[0;32m~/git_debug/swiftest/python/swiftest/swiftest/simulation_class.py:474\u001b[0m, in \u001b[0;36mSimulation.run\u001b[0;34m(self, **kwargs)\u001b[0m\n\u001b[1;32m 471\u001b[0m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_run_swiftest_driver()\n\u001b[1;32m 473\u001b[0m \u001b[38;5;66;03m# Read in new data\u001b[39;00m\n\u001b[0;32m--> 474\u001b[0m \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mbin2xr\u001b[49m\u001b[43m(\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m 476\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m\n", + "File \u001b[0;32m~/git_debug/swiftest/python/swiftest/swiftest/simulation_class.py:2743\u001b[0m, in \u001b[0;36mSimulation.bin2xr\u001b[0;34m(self)\u001b[0m\n\u001b[1;32m 2741\u001b[0m param_tmp[\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mBIN_OUT\u001b[39m\u001b[38;5;124m'\u001b[39m] \u001b[38;5;241m=\u001b[39m os\u001b[38;5;241m.\u001b[39mpath\u001b[38;5;241m.\u001b[39mjoin(\u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39msim_dir, \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mparam[\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mBIN_OUT\u001b[39m\u001b[38;5;124m'\u001b[39m])\n\u001b[1;32m 2742\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mcodename \u001b[38;5;241m==\u001b[39m \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mSwiftest\u001b[39m\u001b[38;5;124m\"\u001b[39m:\n\u001b[0;32m-> 2743\u001b[0m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mdata \u001b[38;5;241m=\u001b[39m \u001b[43mio\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mswiftest2xr\u001b[49m\u001b[43m(\u001b[49m\u001b[43mparam_tmp\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mverbose\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mverbose\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m 2744\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mverbose: \u001b[38;5;28mprint\u001b[39m(\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mSwiftest simulation data stored as xarray DataSet .data\u001b[39m\u001b[38;5;124m'\u001b[39m)\n\u001b[1;32m 2745\u001b[0m \u001b[38;5;28;01melif\u001b[39;00m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mcodename \u001b[38;5;241m==\u001b[39m \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mSwifter\u001b[39m\u001b[38;5;124m\"\u001b[39m:\n", + "File \u001b[0;32m~/git_debug/swiftest/python/swiftest/swiftest/io.py:854\u001b[0m, in \u001b[0;36mswiftest2xr\u001b[0;34m(param, verbose)\u001b[0m\n\u001b[1;32m 852\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m ((param[\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mOUT_TYPE\u001b[39m\u001b[38;5;124m'\u001b[39m] \u001b[38;5;241m==\u001b[39m \u001b[38;5;124m'\u001b[39m\u001b[38;5;124mNETCDF_DOUBLE\u001b[39m\u001b[38;5;124m'\u001b[39m) \u001b[38;5;129;01mor\u001b[39;00m (param[\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mOUT_TYPE\u001b[39m\u001b[38;5;124m'\u001b[39m] \u001b[38;5;241m==\u001b[39m \u001b[38;5;124m'\u001b[39m\u001b[38;5;124mNETCDF_FLOAT\u001b[39m\u001b[38;5;124m'\u001b[39m)):\n\u001b[1;32m 853\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m verbose: \u001b[38;5;28mprint\u001b[39m(\u001b[38;5;124m'\u001b[39m\u001b[38;5;130;01m\\n\u001b[39;00m\u001b[38;5;124mCreating Dataset from NetCDF file\u001b[39m\u001b[38;5;124m'\u001b[39m)\n\u001b[0;32m--> 854\u001b[0m ds \u001b[38;5;241m=\u001b[39m \u001b[43mxr\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mopen_dataset\u001b[49m\u001b[43m(\u001b[49m\u001b[43mparam\u001b[49m\u001b[43m[\u001b[49m\u001b[38;5;124;43m'\u001b[39;49m\u001b[38;5;124;43mBIN_OUT\u001b[39;49m\u001b[38;5;124;43m'\u001b[39;49m\u001b[43m]\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mmask_and_scale\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;28;43;01mFalse\u001b[39;49;00m\u001b[43m)\u001b[49m\n\u001b[1;32m 855\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m param[\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mOUT_TYPE\u001b[39m\u001b[38;5;124m'\u001b[39m] \u001b[38;5;241m==\u001b[39m \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mNETCDF_DOUBLE\u001b[39m\u001b[38;5;124m\"\u001b[39m:\n\u001b[1;32m 856\u001b[0m ds \u001b[38;5;241m=\u001b[39m fix_types(ds,ftype\u001b[38;5;241m=\u001b[39mnp\u001b[38;5;241m.\u001b[39mfloat64)\n", + "File \u001b[0;32m~/.conda/envs/cent7/2020.11-py38/debug_env/lib/python3.8/site-packages/xarray/backends/api.py:495\u001b[0m, in \u001b[0;36mopen_dataset\u001b[0;34m(filename_or_obj, engine, chunks, cache, decode_cf, mask_and_scale, decode_times, decode_timedelta, use_cftime, concat_characters, decode_coords, drop_variables, backend_kwargs, *args, **kwargs)\u001b[0m\n\u001b[1;32m 483\u001b[0m decoders \u001b[38;5;241m=\u001b[39m _resolve_decoders_kwargs(\n\u001b[1;32m 484\u001b[0m decode_cf,\n\u001b[1;32m 485\u001b[0m open_backend_dataset_parameters\u001b[38;5;241m=\u001b[39mbackend\u001b[38;5;241m.\u001b[39mopen_dataset_parameters,\n\u001b[0;32m (...)\u001b[0m\n\u001b[1;32m 491\u001b[0m decode_coords\u001b[38;5;241m=\u001b[39mdecode_coords,\n\u001b[1;32m 492\u001b[0m )\n\u001b[1;32m 494\u001b[0m overwrite_encoded_chunks \u001b[38;5;241m=\u001b[39m kwargs\u001b[38;5;241m.\u001b[39mpop(\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124moverwrite_encoded_chunks\u001b[39m\u001b[38;5;124m\"\u001b[39m, \u001b[38;5;28;01mNone\u001b[39;00m)\n\u001b[0;32m--> 495\u001b[0m backend_ds \u001b[38;5;241m=\u001b[39m \u001b[43mbackend\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mopen_dataset\u001b[49m\u001b[43m(\u001b[49m\n\u001b[1;32m 496\u001b[0m \u001b[43m \u001b[49m\u001b[43mfilename_or_obj\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m 497\u001b[0m \u001b[43m \u001b[49m\u001b[43mdrop_variables\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mdrop_variables\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m 498\u001b[0m \u001b[43m \u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43mdecoders\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m 499\u001b[0m \u001b[43m \u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43mkwargs\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m 500\u001b[0m \u001b[43m\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m 501\u001b[0m ds \u001b[38;5;241m=\u001b[39m _dataset_from_backend_dataset(\n\u001b[1;32m 502\u001b[0m backend_ds,\n\u001b[1;32m 503\u001b[0m filename_or_obj,\n\u001b[0;32m (...)\u001b[0m\n\u001b[1;32m 510\u001b[0m \u001b[38;5;241m*\u001b[39m\u001b[38;5;241m*\u001b[39mkwargs,\n\u001b[1;32m 511\u001b[0m )\n\u001b[1;32m 512\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m ds\n", + "File \u001b[0;32m~/.conda/envs/cent7/2020.11-py38/debug_env/lib/python3.8/site-packages/xarray/backends/h5netcdf_.py:386\u001b[0m, in \u001b[0;36mH5netcdfBackendEntrypoint.open_dataset\u001b[0;34m(self, filename_or_obj, mask_and_scale, decode_times, concat_characters, decode_coords, drop_variables, use_cftime, decode_timedelta, format, group, lock, invalid_netcdf, phony_dims, decode_vlen_strings)\u001b[0m\n\u001b[1;32m 374\u001b[0m store \u001b[38;5;241m=\u001b[39m H5NetCDFStore\u001b[38;5;241m.\u001b[39mopen(\n\u001b[1;32m 375\u001b[0m filename_or_obj,\n\u001b[1;32m 376\u001b[0m \u001b[38;5;28mformat\u001b[39m\u001b[38;5;241m=\u001b[39m\u001b[38;5;28mformat\u001b[39m,\n\u001b[0;32m (...)\u001b[0m\n\u001b[1;32m 381\u001b[0m decode_vlen_strings\u001b[38;5;241m=\u001b[39mdecode_vlen_strings,\n\u001b[1;32m 382\u001b[0m )\n\u001b[1;32m 384\u001b[0m store_entrypoint \u001b[38;5;241m=\u001b[39m StoreBackendEntrypoint()\n\u001b[0;32m--> 386\u001b[0m ds \u001b[38;5;241m=\u001b[39m \u001b[43mstore_entrypoint\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mopen_dataset\u001b[49m\u001b[43m(\u001b[49m\n\u001b[1;32m 387\u001b[0m \u001b[43m \u001b[49m\u001b[43mstore\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m 388\u001b[0m \u001b[43m \u001b[49m\u001b[43mmask_and_scale\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mmask_and_scale\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m 389\u001b[0m \u001b[43m \u001b[49m\u001b[43mdecode_times\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mdecode_times\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m 390\u001b[0m \u001b[43m \u001b[49m\u001b[43mconcat_characters\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mconcat_characters\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m 391\u001b[0m \u001b[43m \u001b[49m\u001b[43mdecode_coords\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mdecode_coords\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m 392\u001b[0m \u001b[43m \u001b[49m\u001b[43mdrop_variables\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mdrop_variables\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m 393\u001b[0m \u001b[43m \u001b[49m\u001b[43muse_cftime\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43muse_cftime\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m 394\u001b[0m \u001b[43m \u001b[49m\u001b[43mdecode_timedelta\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mdecode_timedelta\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m 395\u001b[0m \u001b[43m\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m 396\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m ds\n", + "File \u001b[0;32m~/.conda/envs/cent7/2020.11-py38/debug_env/lib/python3.8/site-packages/xarray/backends/store.py:24\u001b[0m, in \u001b[0;36mStoreBackendEntrypoint.open_dataset\u001b[0;34m(self, store, mask_and_scale, decode_times, concat_characters, decode_coords, drop_variables, use_cftime, decode_timedelta)\u001b[0m\n\u001b[1;32m 12\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21mopen_dataset\u001b[39m(\n\u001b[1;32m 13\u001b[0m \u001b[38;5;28mself\u001b[39m,\n\u001b[1;32m 14\u001b[0m store,\n\u001b[0;32m (...)\u001b[0m\n\u001b[1;32m 22\u001b[0m decode_timedelta\u001b[38;5;241m=\u001b[39m\u001b[38;5;28;01mNone\u001b[39;00m,\n\u001b[1;32m 23\u001b[0m ):\n\u001b[0;32m---> 24\u001b[0m \u001b[38;5;28mvars\u001b[39m, attrs \u001b[38;5;241m=\u001b[39m \u001b[43mstore\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mload\u001b[49m\u001b[43m(\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m 25\u001b[0m encoding \u001b[38;5;241m=\u001b[39m store\u001b[38;5;241m.\u001b[39mget_encoding()\n\u001b[1;32m 27\u001b[0m \u001b[38;5;28mvars\u001b[39m, attrs, coord_names \u001b[38;5;241m=\u001b[39m conventions\u001b[38;5;241m.\u001b[39mdecode_cf_variables(\n\u001b[1;32m 28\u001b[0m \u001b[38;5;28mvars\u001b[39m,\n\u001b[1;32m 29\u001b[0m attrs,\n\u001b[0;32m (...)\u001b[0m\n\u001b[1;32m 36\u001b[0m decode_timedelta\u001b[38;5;241m=\u001b[39mdecode_timedelta,\n\u001b[1;32m 37\u001b[0m )\n", + "File \u001b[0;32m~/.conda/envs/cent7/2020.11-py38/debug_env/lib/python3.8/site-packages/xarray/backends/common.py:123\u001b[0m, in \u001b[0;36mAbstractDataStore.load\u001b[0;34m(self)\u001b[0m\n\u001b[1;32m 101\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21mload\u001b[39m(\u001b[38;5;28mself\u001b[39m):\n\u001b[1;32m 102\u001b[0m \u001b[38;5;124;03m\"\"\"\u001b[39;00m\n\u001b[1;32m 103\u001b[0m \u001b[38;5;124;03m This loads the variables and attributes simultaneously.\u001b[39;00m\n\u001b[1;32m 104\u001b[0m \u001b[38;5;124;03m A centralized loading function makes it easier to create\u001b[39;00m\n\u001b[0;32m (...)\u001b[0m\n\u001b[1;32m 120\u001b[0m \u001b[38;5;124;03m are requested, so care should be taken to make sure its fast.\u001b[39;00m\n\u001b[1;32m 121\u001b[0m \u001b[38;5;124;03m \"\"\"\u001b[39;00m\n\u001b[1;32m 122\u001b[0m variables \u001b[38;5;241m=\u001b[39m FrozenDict(\n\u001b[0;32m--> 123\u001b[0m (_decode_variable_name(k), v) \u001b[38;5;28;01mfor\u001b[39;00m k, v \u001b[38;5;129;01min\u001b[39;00m \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mget_variables\u001b[49m\u001b[43m(\u001b[49m\u001b[43m)\u001b[49m\u001b[38;5;241m.\u001b[39mitems()\n\u001b[1;32m 124\u001b[0m )\n\u001b[1;32m 125\u001b[0m attributes \u001b[38;5;241m=\u001b[39m FrozenDict(\u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mget_attrs())\n\u001b[1;32m 126\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m variables, attributes\n", + "File \u001b[0;32m~/.conda/envs/cent7/2020.11-py38/debug_env/lib/python3.8/site-packages/xarray/backends/h5netcdf_.py:229\u001b[0m, in \u001b[0;36mH5NetCDFStore.get_variables\u001b[0;34m(self)\u001b[0m\n\u001b[1;32m 228\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21mget_variables\u001b[39m(\u001b[38;5;28mself\u001b[39m):\n\u001b[0;32m--> 229\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[43mFrozenDict\u001b[49m\u001b[43m(\u001b[49m\n\u001b[1;32m 230\u001b[0m \u001b[43m \u001b[49m\u001b[43m(\u001b[49m\u001b[43mk\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mopen_store_variable\u001b[49m\u001b[43m(\u001b[49m\u001b[43mk\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mv\u001b[49m\u001b[43m)\u001b[49m\u001b[43m)\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;28;43;01mfor\u001b[39;49;00m\u001b[43m \u001b[49m\u001b[43mk\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mv\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;129;43;01min\u001b[39;49;00m\u001b[43m \u001b[49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mds\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mvariables\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mitems\u001b[49m\u001b[43m(\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m 231\u001b[0m \u001b[43m \u001b[49m\u001b[43m)\u001b[49m\n", + "File \u001b[0;32m~/.conda/envs/cent7/2020.11-py38/debug_env/lib/python3.8/site-packages/xarray/core/utils.py:476\u001b[0m, in \u001b[0;36mFrozenDict\u001b[0;34m(*args, **kwargs)\u001b[0m\n\u001b[1;32m 475\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21mFrozenDict\u001b[39m(\u001b[38;5;241m*\u001b[39margs, \u001b[38;5;241m*\u001b[39m\u001b[38;5;241m*\u001b[39mkwargs) \u001b[38;5;241m-\u001b[39m\u001b[38;5;241m>\u001b[39m Frozen:\n\u001b[0;32m--> 476\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m Frozen(\u001b[38;5;28;43mdict\u001b[39;49m\u001b[43m(\u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43margs\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43mkwargs\u001b[49m\u001b[43m)\u001b[49m)\n", + "File \u001b[0;32m~/.conda/envs/cent7/2020.11-py38/debug_env/lib/python3.8/site-packages/xarray/backends/h5netcdf_.py:230\u001b[0m, in \u001b[0;36m\u001b[0;34m(.0)\u001b[0m\n\u001b[1;32m 228\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21mget_variables\u001b[39m(\u001b[38;5;28mself\u001b[39m):\n\u001b[1;32m 229\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m FrozenDict(\n\u001b[0;32m--> 230\u001b[0m (k, \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mopen_store_variable\u001b[49m\u001b[43m(\u001b[49m\u001b[43mk\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mv\u001b[49m\u001b[43m)\u001b[49m) \u001b[38;5;28;01mfor\u001b[39;00m k, v \u001b[38;5;129;01min\u001b[39;00m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mds\u001b[38;5;241m.\u001b[39mvariables\u001b[38;5;241m.\u001b[39mitems()\n\u001b[1;32m 231\u001b[0m )\n", + "File \u001b[0;32m~/.conda/envs/cent7/2020.11-py38/debug_env/lib/python3.8/site-packages/xarray/backends/h5netcdf_.py:195\u001b[0m, in \u001b[0;36mH5NetCDFStore.open_store_variable\u001b[0;34m(self, name, var)\u001b[0m\n\u001b[1;32m 192\u001b[0m \u001b[38;5;28;01mimport\u001b[39;00m \u001b[38;5;21;01mh5py\u001b[39;00m\n\u001b[1;32m 194\u001b[0m dimensions \u001b[38;5;241m=\u001b[39m var\u001b[38;5;241m.\u001b[39mdimensions\n\u001b[0;32m--> 195\u001b[0m data \u001b[38;5;241m=\u001b[39m indexing\u001b[38;5;241m.\u001b[39mLazilyIndexedArray(\u001b[43mH5NetCDFArrayWrapper\u001b[49m\u001b[43m(\u001b[49m\u001b[43mname\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[43m)\u001b[49m)\n\u001b[1;32m 196\u001b[0m attrs \u001b[38;5;241m=\u001b[39m _read_attributes(var)\n\u001b[1;32m 198\u001b[0m \u001b[38;5;66;03m# netCDF4 specific encoding\u001b[39;00m\n", + "File \u001b[0;32m~/.conda/envs/cent7/2020.11-py38/debug_env/lib/python3.8/site-packages/xarray/backends/netCDF4_.py:56\u001b[0m, in \u001b[0;36mBaseNetCDF4Array.__init__\u001b[0;34m(self, variable_name, datastore)\u001b[0m\n\u001b[1;32m 53\u001b[0m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mvariable_name \u001b[38;5;241m=\u001b[39m variable_name\n\u001b[1;32m 55\u001b[0m array \u001b[38;5;241m=\u001b[39m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mget_array()\n\u001b[0;32m---> 56\u001b[0m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mshape \u001b[38;5;241m=\u001b[39m \u001b[43marray\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mshape\u001b[49m\n\u001b[1;32m 58\u001b[0m dtype \u001b[38;5;241m=\u001b[39m array\u001b[38;5;241m.\u001b[39mdtype\n\u001b[1;32m 59\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m dtype \u001b[38;5;129;01mis\u001b[39;00m \u001b[38;5;28mstr\u001b[39m:\n\u001b[1;32m 60\u001b[0m \u001b[38;5;66;03m# use object dtype because that's the only way in numpy to\u001b[39;00m\n\u001b[1;32m 61\u001b[0m \u001b[38;5;66;03m# represent variable length strings; it also prevents automatic\u001b[39;00m\n\u001b[1;32m 62\u001b[0m \u001b[38;5;66;03m# string concatenation via conventions.decode_cf_variable\u001b[39;00m\n", + "File \u001b[0;32m~/.conda/envs/cent7/2020.11-py38/debug_env/lib/python3.8/site-packages/h5netcdf/core.py:259\u001b[0m, in \u001b[0;36mBaseVariable.shape\u001b[0;34m(self)\u001b[0m\n\u001b[1;32m 257\u001b[0m \u001b[38;5;124;03m\"\"\"Return current sizes of all variable dimensions.\"\"\"\u001b[39;00m\n\u001b[1;32m 258\u001b[0m \u001b[38;5;66;03m# return actual dimensions sizes, this is in line with netcdf4-python\u001b[39;00m\n\u001b[0;32m--> 259\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[38;5;28mtuple\u001b[39m([\u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_parent\u001b[38;5;241m.\u001b[39m_all_dimensions[d]\u001b[38;5;241m.\u001b[39msize \u001b[38;5;28;01mfor\u001b[39;00m d \u001b[38;5;129;01min\u001b[39;00m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mdimensions])\n", + "File \u001b[0;32m~/.conda/envs/cent7/2020.11-py38/debug_env/lib/python3.8/site-packages/h5netcdf/core.py:259\u001b[0m, in \u001b[0;36m\u001b[0;34m(.0)\u001b[0m\n\u001b[1;32m 257\u001b[0m \u001b[38;5;124;03m\"\"\"Return current sizes of all variable dimensions.\"\"\"\u001b[39;00m\n\u001b[1;32m 258\u001b[0m \u001b[38;5;66;03m# return actual dimensions sizes, this is in line with netcdf4-python\u001b[39;00m\n\u001b[0;32m--> 259\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[38;5;28mtuple\u001b[39m([\u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43m_parent\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43m_all_dimensions\u001b[49m\u001b[43m[\u001b[49m\u001b[43md\u001b[49m\u001b[43m]\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43msize\u001b[49m \u001b[38;5;28;01mfor\u001b[39;00m d \u001b[38;5;129;01min\u001b[39;00m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mdimensions])\n", + "File \u001b[0;32m~/.conda/envs/cent7/2020.11-py38/debug_env/lib/python3.8/site-packages/h5netcdf/dimensions.py:115\u001b[0m, in \u001b[0;36mDimension.size\u001b[0;34m(self)\u001b[0m\n\u001b[1;32m 113\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m reflist \u001b[38;5;129;01mis\u001b[39;00m \u001b[38;5;129;01mnot\u001b[39;00m \u001b[38;5;28;01mNone\u001b[39;00m:\n\u001b[1;32m 114\u001b[0m \u001b[38;5;28;01mfor\u001b[39;00m ref, axis \u001b[38;5;129;01min\u001b[39;00m reflist:\n\u001b[0;32m--> 115\u001b[0m var \u001b[38;5;241m=\u001b[39m \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43m_parent\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43m_h5group\u001b[49m\u001b[43m[\u001b[49m\u001b[38;5;124;43m\"\u001b[39;49m\u001b[38;5;124;43m/\u001b[39;49m\u001b[38;5;124;43m\"\u001b[39;49m\u001b[43m]\u001b[49m\u001b[43m[\u001b[49m\u001b[43mref\u001b[49m\u001b[43m]\u001b[49m\n\u001b[1;32m 116\u001b[0m size \u001b[38;5;241m=\u001b[39m \u001b[38;5;28mmax\u001b[39m(var\u001b[38;5;241m.\u001b[39mshape[axis], size)\n\u001b[1;32m 117\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m size\n", + "File \u001b[0;32mh5py/_objects.pyx:54\u001b[0m, in \u001b[0;36mh5py._objects.with_phil.wrapper\u001b[0;34m()\u001b[0m\n", + "File \u001b[0;32mh5py/_objects.pyx:55\u001b[0m, in \u001b[0;36mh5py._objects.with_phil.wrapper\u001b[0;34m()\u001b[0m\n", + "File \u001b[0;32m~/.conda/envs/cent7/2020.11-py38/debug_env/lib/python3.8/site-packages/h5py/_hl/group.py:337\u001b[0m, in \u001b[0;36mGroup.__getitem__\u001b[0;34m(self, name)\u001b[0m\n\u001b[1;32m 335\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m Group(oid)\n\u001b[1;32m 336\u001b[0m \u001b[38;5;28;01melif\u001b[39;00m otype \u001b[38;5;241m==\u001b[39m h5i\u001b[38;5;241m.\u001b[39mDATASET:\n\u001b[0;32m--> 337\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[43mdataset\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mDataset\u001b[49m\u001b[43m(\u001b[49m\u001b[43moid\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mreadonly\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43m(\u001b[49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mfile\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mmode\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;241;43m==\u001b[39;49m\u001b[43m \u001b[49m\u001b[38;5;124;43m'\u001b[39;49m\u001b[38;5;124;43mr\u001b[39;49m\u001b[38;5;124;43m'\u001b[39;49m\u001b[43m)\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m 338\u001b[0m \u001b[38;5;28;01melif\u001b[39;00m otype \u001b[38;5;241m==\u001b[39m h5i\u001b[38;5;241m.\u001b[39mDATATYPE:\n\u001b[1;32m 339\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m datatype\u001b[38;5;241m.\u001b[39mDatatype(oid)\n", + "File \u001b[0;32mh5py/_objects.pyx:54\u001b[0m, in \u001b[0;36mh5py._objects.with_phil.wrapper\u001b[0;34m()\u001b[0m\n", + "File \u001b[0;32mh5py/_objects.pyx:55\u001b[0m, in \u001b[0;36mh5py._objects.with_phil.wrapper\u001b[0;34m()\u001b[0m\n", + "File \u001b[0;32m~/.conda/envs/cent7/2020.11-py38/debug_env/lib/python3.8/site-packages/h5py/_hl/dataset.py:622\u001b[0m, in \u001b[0;36mDataset.__init__\u001b[0;34m(self, bind, readonly)\u001b[0m\n\u001b[1;32m 619\u001b[0m \u001b[38;5;28;01mraise\u001b[39;00m \u001b[38;5;167;01mValueError\u001b[39;00m(\u001b[38;5;124m\"\u001b[39m\u001b[38;5;132;01m%s\u001b[39;00m\u001b[38;5;124m is not a DatasetID\u001b[39m\u001b[38;5;124m\"\u001b[39m \u001b[38;5;241m%\u001b[39m bind)\n\u001b[1;32m 620\u001b[0m \u001b[38;5;28msuper\u001b[39m()\u001b[38;5;241m.\u001b[39m\u001b[38;5;21m__init__\u001b[39m(bind)\n\u001b[0;32m--> 622\u001b[0m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_dcpl \u001b[38;5;241m=\u001b[39m \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mid\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mget_create_plist\u001b[49m\u001b[43m(\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m 623\u001b[0m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_dxpl \u001b[38;5;241m=\u001b[39m h5p\u001b[38;5;241m.\u001b[39mcreate(h5p\u001b[38;5;241m.\u001b[39mDATASET_XFER)\n\u001b[1;32m 624\u001b[0m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_filters \u001b[38;5;241m=\u001b[39m filters\u001b[38;5;241m.\u001b[39mget_filters(\u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_dcpl)\n", + "\u001b[0;31mKeyboardInterrupt\u001b[0m: " + ] + } + ], "source": [ "# Run the simulation\n", "sim.run()" ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "02a8911d-3b2c-415c-9290-bf1519a3f5c6", + "metadata": {}, + "outputs": [], + "source": [] } ], "metadata": { diff --git a/examples/Basic_Simulation/initial_conditions.py b/examples/Basic_Simulation/initial_conditions.py index 9bb279b51..4d7fe7ae5 100644 --- a/examples/Basic_Simulation/initial_conditions.py +++ b/examples/Basic_Simulation/initial_conditions.py @@ -23,7 +23,7 @@ from numpy.random import default_rng # Initialize the simulation object as a variable -sim = swiftest.Simulation(tstart=0.0, tstop=1.0e3, dt=0.010, tstep_out=1.0e0, fragmentation=True, minimum_fragment_mass = 2.5e-11, mtiny=2.5e-8) +sim = swiftest.Simulation(tstart=0.0, tstop=1.0e3, dt=0.01, tstep_out=1.0e0, dump_cadence=0, fragmentation=True, minimum_fragment_mass = 2.5e-11, mtiny=2.5e-8) # 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/read_old_run.ipynb b/examples/Basic_Simulation/read_old_run.ipynb new file mode 100644 index 000000000..0a7459084 --- /dev/null +++ b/examples/Basic_Simulation/read_old_run.ipynb @@ -0,0 +1,853 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 1, + "id": "86b309f3-d400-4164-88fc-cf56b6cfcf5f", + "metadata": {}, + "outputs": [], + "source": [ + "import swiftest" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "353b2559-83a2-496f-8648-60f6121333a5", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Reading Swiftest file /home/daminton/git_debug/swiftest/examples/Basic_Simulation/simdata/param.in\n", + "\n", + "Creating Dataset from NetCDF file\n", + "Successfully converted 78 output frames.\n", + "Swiftest simulation data stored as xarray DataSet .data\n" + ] + } + ], + "source": [ + "sim = swiftest.Simulation(read_old_output_file=True)" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "id": "c7bceca2-678d-469d-a3a7-06cfa9f29fc0", + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "
<xarray.Dataset>\n",
+       "Dimensions:          (time: 78, name: 25)\n",
+       "Coordinates:\n",
+       "  * time             (time) float64 0.0 2.0 nan 4.0 nan ... nan 76.0 nan 78.0\n",
+       "  * name             (name) <U32 'Sun' 'Mercury' ... 'TestParticle_10'\n",
+       "Data variables: (12/45)\n",
+       "    id               (name) int64 0 1 2 3 4 5 6 7 8 ... 17 18 19 20 21 22 23 24\n",
+       "    npl              (time) int64 14 14 -2147483647 14 ... 14 -2147483647 14\n",
+       "    ntp              (time) int64 10 10 -2147483647 10 ... 10 -2147483647 10\n",
+       "    nplm             (time) int64 13 13 -2147483647 13 ... 13 -2147483647 13\n",
+       "    particle_type    (name) <U32 'Central Body' ... 'Test Particle'\n",
+       "    status           (name) <U32 'ACTIVE' 'ACTIVE' ... 'ACTIVE' 'ACTIVE'\n",
+       "    ...               ...\n",
+       "    Ip3              (time, name) float64 0.07 0.346 0.4 0.3307 ... nan nan nan\n",
+       "    rotx             (time, name) float64 11.21 3.573 0.1765 ... nan nan nan\n",
+       "    roty             (time, name) float64 -38.76 -18.38 -3.661 ... nan nan nan\n",
+       "    rotz             (time, name) float64 82.25 34.36 8.703 ... nan nan nan\n",
+       "    j2rp2            (time) float64 4.754e-12 4.754e-12 nan ... nan 4.754e-12\n",
+       "    j4rp4            (time) float64 -2.247e-18 -2.247e-18 nan ... nan -2.247e-18
" + ], + "text/plain": [ + "\n", + "Dimensions: (time: 78, name: 25)\n", + "Coordinates:\n", + " * time (time) float64 0.0 2.0 nan 4.0 nan ... nan 76.0 nan 78.0\n", + " * name (name)
    • " ], "text/plain": [ "\n", @@ -436,7 +434,7 @@ "Coordinates:\n", " * name (name)
    • " ], "text/plain": [ "\n", @@ -885,7 +881,7 @@ "Coordinates:\n", " * name (name) = 0.0_DP) .and. (rh2 > rmax2)) then tp%status(i) = DISCARDED_RMAX write(idstr, *) tp%id(i) - write(timestr, *) param%t + write(timestr, *) system%t write(*, *) "Particle " // trim(adjustl(tp%info(i)%name)) // " (" // trim(adjustl(idstr)) // ")" // & " too far from the central body at t = " // trim(adjustl(timestr)) tp%ldiscard(i) = .true. tp%lmask(i) = .false. - call tp%info(i)%set_value(status="DISCARDED_RMAX", discard_time=param%t, discard_xh=tp%xh(:,i), & + call tp%info(i)%set_value(status="DISCARDED_RMAX", discard_time=system%t, discard_xh=tp%xh(:,i), & discard_vh=tp%vh(:,i)) else if ((param%rmin >= 0.0_DP) .and. (rh2 < rmin2)) then tp%status(i) = DISCARDED_RMIN write(idstr, *) tp%id(i) - write(timestr, *) param%t + write(timestr, *) system%t write(*, *) "Particle " // trim(adjustl(tp%info(i)%name)) // " (" // trim(adjustl(idstr)) // ")" // & " too close to the central body at t = " // trim(adjustl(timestr)) tp%ldiscard(i) = .true. tp%lmask(i) = .false. - call tp%info(i)%set_value(status="DISCARDED_RMIN", discard_time=param%t, discard_xh=tp%xh(:,i), & + call tp%info(i)%set_value(status="DISCARDED_RMIN", discard_time=system%t, discard_xh=tp%xh(:,i), & discard_vh=tp%vh(:,i), discard_body_id=cb%id) else if (param%rmaxu >= 0.0_DP) then rb2 = dot_product(tp%xb(:, i), tp%xb(:, i)) @@ -159,12 +159,12 @@ subroutine discard_cb_tp(tp, system, param) if ((energy > 0.0_DP) .and. (rb2 > rmaxu2)) then tp%status(i) = DISCARDED_RMAXU write(idstr, *) tp%id(i) - write(timestr, *) param%t + write(timestr, *) system%t write(*, *) "Particle " // trim(adjustl(tp%info(i)%name)) // " (" // trim(adjustl(idstr)) // ")" // & " is unbound and too far from barycenter at t = " // trim(adjustl(timestr)) tp%ldiscard(i) = .true. tp%lmask(i) = .false. - call tp%info(i)%set_value(status="DISCARDED_RMAXU", discard_time=param%t, discard_xh=tp%xh(:,i), & + call tp%info(i)%set_value(status="DISCARDED_RMAXU", discard_time=system%t, discard_xh=tp%xh(:,i), & discard_vh=tp%vh(:,i)) end if end if @@ -194,7 +194,7 @@ subroutine discard_peri_tp(tp, system, param) real(DP), dimension(NDIM) :: dx character(len=STRMAX) :: idstr, timestr - associate(cb => system%cb, ntp => tp%nbody, pl => system%pl, npl => system%pl%nbody, t => param%t) + associate(cb => system%cb, ntp => tp%nbody, pl => system%pl, npl => system%pl%nbody, t => system%t) call tp%get_peri(system, param) do i = 1, ntp if (tp%status(i) == ACTIVE) then @@ -211,11 +211,11 @@ subroutine discard_peri_tp(tp, system, param) (tp%peri(i) <= param%qmin)) then tp%status(i) = DISCARDED_PERI write(idstr, *) tp%id(i) - write(timestr, *) param%t + write(timestr, *) system%t write(*, *) "Particle " // trim(adjustl(tp%info(i)%name)) // " (" // trim(adjustl(idstr)) // ")" // & " perihelion distance too small at t = " // trim(adjustl(timestr)) tp%ldiscard(i) = .true. - call tp%info(i)%set_value(status="DISCARDED_PERI", discard_time=param%t, discard_xh=tp%xh(:,i), & + call tp%info(i)%set_value(status="DISCARDED_PERI", discard_time=system%t, discard_xh=tp%xh(:,i), & discard_vh=tp%vh(:,i), discard_body_id=pl%id(j)) end if end if @@ -246,7 +246,7 @@ subroutine discard_pl_tp(tp, system, param) real(DP), dimension(NDIM) :: dx, dv character(len=STRMAX) :: idstri, idstrj, timestr - associate(ntp => tp%nbody, pl => system%pl, npl => system%pl%nbody, t => param%t, dt => param%dt) + associate(ntp => tp%nbody, pl => system%pl, npl => system%pl%nbody, t => system%t, dt => param%dt) do i = 1, ntp if (tp%status(i) == ACTIVE) then do j = 1, npl @@ -260,12 +260,12 @@ subroutine discard_pl_tp(tp, system, param) pl%ldiscard(j) = .true. write(idstri, *) tp%id(i) write(idstrj, *) pl%id(j) - write(timestr, *) param%t + write(timestr, *) system%t write(*, *) "Test particle " // trim(adjustl(tp%info(i)%name)) // " (" // trim(adjustl(idstri)) // ")" & // " too close to massive body " // trim(adjustl(pl%info(j)%name)) // " (" // trim(adjustl(idstrj)) // ")" & // " at t = " // trim(adjustl(timestr)) tp%ldiscard(i) = .true. - call tp%info(i)%set_value(status="DISCARDED_PLR", discard_time=param%t, discard_xh=tp%xh(:,i), & + call tp%info(i)%set_value(status="DISCARDED_PLR", discard_time=system%t, discard_xh=tp%xh(:,i), & discard_vh=tp%vh(:,i), discard_body_id=pl%id(j)) exit end if diff --git a/src/io/io.f90 b/src/io/io.f90 index de4c547e9..7ee5b9cbd 100644 --- a/src/io/io.f90 +++ b/src/io/io.f90 @@ -35,7 +35,7 @@ module subroutine io_compact_output(self, param, timer) select type(timer) class is (walltimer) - formatted_output = fmt("ILOOP",param%iloop) // fmt("T",param%t) // fmt("NPL",self%pl%nbody) // fmt("NTP",self%tp%nbody) + formatted_output = fmt("ILOOP",param%iloop) // fmt("T",self%t) // fmt("NPL",self%pl%nbody) // fmt("NTP",self%tp%nbody) select type(pl => self%pl) class is (symba_pl) formatted_output = formatted_output // fmt("NPLM",pl%nplm) @@ -250,7 +250,7 @@ module subroutine io_dump_system(self, param) dump_param%in_netcdf = trim(adjustl(DUMP_NC_FILE(idx))) dump_param%nciu%id_chunk = self%pl%nbody + self%tp%nbody dump_param%nciu%time_chunk = 1 - dump_param%T0 = param%t + dump_param%tstart = self%t call dump_param%dump(param_file_name) @@ -270,6 +270,33 @@ module subroutine io_dump_system(self, param) end subroutine io_dump_system + module subroutine io_dump_system_storage(self, param) + !! author: David A. Minton + !! + !! Dumps the time history of the simulation to file. Each time it writes a frame to file, it deallocates the system + !! object from inside. It will only dump frames with systems that are allocated, so this can be called at the end of + !! a simulation for cases when the number of saved frames is not equal to the dump cadence (for instance, if the dump + !! cadence is not divisible by the total number of loops). + implicit none + ! Arguments + class(swiftest_storage(*)), intent(inout) :: self !! Swiftest simulation history storage object + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters + ! Internals + integer(I8B) :: i, iloop_start + + iloop_start = param%iloop - param%istep_out * param%dump_cadence + 1_I8B + do i = 1_I8B, param%dump_cadence + if (allocated(self%frame(i)%system)) then + param%ioutput = int(iloop_start / param%istep_out, kind=I8B) + i + call self%frame(i)%system%write_frame(param) + deallocate(self%frame(i)%system) + end if + end do + + return + end subroutine io_dump_system_storage + + module subroutine io_get_args(integrator, param_file_name, display_style) !! author: David A. Minton !! @@ -461,7 +488,6 @@ module subroutine io_param_reader(self, unit, iotype, v_list, iostat, iomsg) integer, intent(out) :: iostat !! IO status code character(len=*), intent(inout) :: iomsg !! Message to pass if iostat /= 0 ! Internals - logical :: t0_set = .false. !! Is the initial time set in the input file? logical :: tstart_set = .false. !! Is the final time set in the input file? logical :: tstop_set = .false. !! Is the final time set in the input file? logical :: dt_set = .false. !! Is the step size set in the input file? @@ -489,9 +515,8 @@ module subroutine io_param_reader(self, unit, iotype, v_list, iostat, iomsg) select case (param_name) case ("T0") read(param_value, *, err = 667, iomsg = iomsg) param%t0 - t0_set = .true. case ("TSTART") - read(param_value, *, err = 667, iomsg = iomsg) param%t0 + read(param_value, *, err = 667, iomsg = iomsg) param%tstart tstart_set = .true. case ("TSTOP") read(param_value, *, err = 667, iomsg = iomsg) param%tstop @@ -526,8 +551,8 @@ module subroutine io_param_reader(self, unit, iotype, v_list, iostat, iomsg) case ("OUT_STAT") call io_toupper(param_value) param%out_stat = param_value - case ("ISTEP_DUMP") - read(param_value, *, err = 667, iomsg = iomsg) param%istep_dump + case ("DUMP_CADENCE") + read(param_value, *, err = 667, iomsg = iomsg) param%dump_cadence case ("CHK_CLOSE") call io_toupper(param_value) if (param_value == "YES" .or. param_value == 'T') param%lclose = .true. @@ -652,7 +677,7 @@ module subroutine io_param_reader(self, unit, iotype, v_list, iostat, iomsg) iostat = 0 ! Do basic sanity checks on the input values - if ((.not. t0_set) .or. (.not. tstop_set) .or. (.not. dt_set)) then + if ((.not. tstart_set) .or. (.not. tstop_set) .or. (.not. dt_set)) then write(iomsg,*) 'Valid simulation time not set' iostat = -1 return @@ -672,8 +697,13 @@ module subroutine io_param_reader(self, unit, iotype, v_list, iostat, iomsg) iostat = -1 return end if - if ((param%istep_out <= 0) .and. (param%istep_dump <= 0)) then - write(iomsg,*) 'Invalid istep' + if (param%istep_out <= 0) then + write(iomsg,*) 'Invalid ISTEP_OUT. Must be a positive integer' + iostat = -1 + return + end if + if (param%dump_cadence < 0) then + write(iomsg,*) 'Invalid DUMP_CADENCE. Must be a positive integer or 0.' iostat = -1 return end if @@ -858,6 +888,7 @@ module subroutine io_param_writer(self, unit, iotype, v_list, iostat, iomsg) associate(param => self) call io_param_writer_one("T0", param%t0, unit) + call io_param_writer_one("TSTART", param%tstart, unit) call io_param_writer_one("TSTOP", param%tstop, unit) call io_param_writer_one("DT", param%dt, unit) call io_param_writer_one("IN_TYPE", param%in_type, unit) @@ -870,7 +901,7 @@ module subroutine io_param_writer(self, unit, iotype, v_list, iostat, iomsg) end if call io_param_writer_one("IN_FORM", param%in_form, unit) - if (param%istep_dump > 0) call io_param_writer_one("ISTEP_DUMP",param%istep_dump, unit) + if (param%dump_cadence > 0) call io_param_writer_one("DUMP_CADENCE",param%dump_cadence, unit) if (param%istep_out > 0) then call io_param_writer_one("ISTEP_OUT", param%istep_out, unit) call io_param_writer_one("BIN_OUT", param%outfile, unit) @@ -1503,7 +1534,7 @@ module subroutine io_write_frame_system(self, param) logical :: fileExists param%nciu%id_chunk = self%pl%nbody + self%tp%nbody - param%nciu%time_chunk = max(param%istep_dump / param%istep_out, 1) + param%nciu%time_chunk = max(param%dump_cadence / param%istep_out, 1) if (lfirst) then inquire(file=param%outfile, exist=fileExists) diff --git a/src/main/swiftest_driver.f90 b/src/main/swiftest_driver.f90 index 8f02b74f8..a9af0cc71 100644 --- a/src/main/swiftest_driver.f90 +++ b/src/main/swiftest_driver.f90 @@ -23,12 +23,11 @@ program swiftest_driver character(len=:), allocatable :: integrator !! Integrator type code (see swiftest_globals for symbolic names) character(len=:),allocatable :: param_file_name !! Name of the file containing user-defined parameters character(len=:), allocatable :: display_style !! Style of the output display {"STANDARD", "COMPACT", "PROGRESS"}). Default is "STANDARD" - integer(I4B) :: ierr !! I/O error code integer(I8B) :: idump !! Dump cadence counter integer(I8B) :: iout !! Output cadence counter - integer(I8B) :: ioutput_t0 !! The output frame counter at time 0 + integer(I8B) :: istart !! Starting index for loop counter integer(I8B) :: nloops !! Number of steps to take in the simulation - 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 + integer(I8B) :: iframe !! System history frame cindex 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 @@ -40,7 +39,7 @@ program swiftest_driver character(len=64) :: pbarmessage character(*), parameter :: symbacompactfmt = '(";NPLM",ES22.15,$)' - + type(swiftest_storage(nframes=:)), allocatable :: system_history call io_get_args(integrator, param_file_name, display_style) @@ -65,36 +64,42 @@ program swiftest_driver call setup_construct_system(nbody_system, param) call param%read_in(param_file_name) - associate(t => param%t, & + associate(t => nbody_system%t, & t0 => param%t0, & + tstart => param%tstart, & dt => param%dt, & tstop => param%tstop, & iloop => param%iloop, & istep_out => param%istep_out, & - istep_dump => param%istep_dump, & + dump_cadence => param%dump_cadence, & ioutput => param%ioutput, & display_style => param%display_style, & display_unit => param%display_unit) call nbody_system%initialize(param) - t = t0 - iloop = 0 + + ! Set up loop and output cadence variables + t = tstart iout = istep_out - idump = istep_dump nloops = ceiling((tstop - t0) / dt, kind=I8B) - ioutput_t0 = int(t0 / dt / istep_out, kind=I8B) - ioutput = ioutput_t0 - ! Prevent duplicate frames from being written if this is a restarted run + istart = ceiling((tstart - t0) / dt + 1, kind=I8B) + ioutput = int(istart / istep_out, kind=I8B) + + ! Set up system storage for intermittent file dumps + if (dump_cadence == 0) dump_cadence = ceiling(nloops / (1.0_DP * istep_out), kind=I8B) + allocate(swiftest_storage(dump_cadence) :: system_history) + idump = dump_cadence + + ! If this is a new run, compute energy initial conditions (if energy tracking is turned on) and write the initial conditions to file. if (param%lrestart) then - old_t_final = nbody_system%get_old_t_final(param) + if (param%lenergy) call nbody_system%conservation_report(param, lterminal=.true.) else - old_t_final = t0 if (param%lenergy) call nbody_system%conservation_report(param, lterminal=.false.) ! This will save the initial values of energy and momentum - if (istep_out > 0) call nbody_system%write_frame(param) + call nbody_system%write_frame(param) end if write(display_unit, *) " *************** Main Loop *************** " - if (param%lrestart .and. param%lenergy) call nbody_system%conservation_report(param, lterminal=.true.) + if (display_style == "PROGRESS") then call pbar%reset(nloops) write(pbarmessage,fmt=pbarfmt) t0, tstop @@ -103,7 +108,8 @@ program swiftest_driver write(*,*) "SWIFTEST START " // trim(adjustl(param%integrator)) call nbody_system%compact_output(param,integration_timer) end if - do iloop = 1, nloops + + do iloop = istart, nloops !> Step the system forward in time call integration_timer%start() call nbody_system%step(param, t, dt) @@ -119,16 +125,23 @@ program swiftest_driver if (istep_out > 0) then iout = iout - 1 if (iout == 0) then - ioutput = ioutput_t0 + iloop / istep_out - call nbody_system%write_frame(param) + idump = idump - 1 + iframe = dump_cadence - idump + system_history%frame(iframe) = nbody_system + + if (idump == 0) then + call nbody_system%dump(param) + call system_history%dump(param) + idump = dump_cadence + end if - tfrac = (param%t - param%t0) / (param%tstop - param%t0) + tfrac = (t - t0) / (tstop - t0) select type(pl => nbody_system%pl) class is (symba_pl) - write(display_unit, symbastatfmt) param%t, tfrac, pl%nplm, pl%nbody, nbody_system%tp%nbody + write(display_unit, symbastatfmt) t, tfrac, pl%nplm, pl%nbody, nbody_system%tp%nbody class default - write(display_unit, statusfmt) param%t, tfrac, pl%nbody, nbody_system%tp%nbody + write(display_unit, statusfmt) 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:", unit=display_unit, nsubsteps=istep_out) @@ -146,21 +159,11 @@ program swiftest_driver end if end if - !> If the loop counter is at the dump cadence value, dump the state of the system to a file in case it needs to be restarted - if (istep_dump > 0) then - idump = idump - 1 - if (idump == 0) then - call nbody_system%dump(param) - idump = istep_dump - end if - end if end do + ! Dump any remaining history if it exists + call system_history%dump(param) if (display_style == "COMPACT") write(*,*) "SWIFTEST STOP" // trim(adjustl(param%integrator)) end associate - call nbody_system%dealloc() - call util_exit(SUCCESS) - - stop end program swiftest_driver diff --git a/src/modules/swiftest_classes.f90 b/src/modules/swiftest_classes.f90 index 874bfdf31..ab7000e36 100644 --- a/src/modules/swiftest_classes.f90 +++ b/src/modules/swiftest_classes.f90 @@ -32,43 +32,43 @@ module swiftest_classes !> User defined parameters that are read in from the parameters input file. !> Each paramter is initialized to a default values. type :: swiftest_parameters - character(STRMAX) :: integrator = UNKNOWN_INTEGRATOR !! Symbolic name of the nbody integrator used - character(STRMAX) :: param_file_name = "param.in" !! The default name of the parameter input file - integer(I4B) :: maxid = -1 !! The current maximum particle id number - integer(I4B) :: maxid_collision = 0 !! The current maximum collision id number - real(DP) :: t0 = -1.0_DP !! Integration start time - real(DP) :: t = -1.0_DP !! Integration current time - real(DP) :: tstop = -1.0_DP !! Integration stop time - real(DP) :: dt = -1.0_DP !! Time step - integer(I8B) :: iloop = 0_I8B !! Main loop counter - integer(I8B) :: ioutput = 0_I8B !! Output counter - character(STRMAX) :: incbfile = CB_INFILE !! Name of input file for the central body - character(STRMAX) :: inplfile = PL_INFILE !! Name of input file for massive bodies - character(STRMAX) :: intpfile = TP_INFILE !! Name of input file for test particles - character(STRMAX) :: in_netcdf = NC_INFILE !! Name of system input file for NetCDF input - character(STRMAX) :: in_type = "ASCII" !! Data representation type of input data files - character(STRMAX) :: in_form = "XV" !! Format of input data files ("EL" or "XV") - integer(I4B) :: istep_out = -1 !! Number of time steps between binary outputs - character(STRMAX) :: outfile = NETCDF_OUTFILE !! Name of output binary file - character(STRMAX) :: out_type = "NETCDF_DOUBLE" !! Binary format of output file - character(STRMAX) :: out_form = "XVEL" !! Data to write to output file - character(STRMAX) :: out_stat = 'NEW' !! Open status for output binary file - integer(I4B) :: istep_dump = -1 !! Number of time steps between dumps - real(DP) :: rmin = -1.0_DP !! Minimum heliocentric radius for test particle - real(DP) :: rmax = -1.0_DP !! Maximum heliocentric radius for test particle - real(DP) :: rmaxu = -1.0_DP !! Maximum unbound heliocentric radius for test particle - real(DP) :: qmin = -1.0_DP !! Minimum pericenter distance for test particle - character(STRMAX) :: qmin_coord = 'HELIO' !! Coordinate frame to use for qmin - real(DP) :: qmin_alo = -1.0_DP !! Minimum semimajor axis for qmin - real(DP) :: qmin_ahi = -1.0_DP !! Maximum semimajor axis for qmin - real(QP) :: MU2KG = -1.0_QP !! Converts mass units to grams - real(QP) :: TU2S = -1.0_QP !! Converts time units to seconds - real(QP) :: DU2M = -1.0_QP !! Converts distance unit to centimeters - real(DP) :: GU = -1.0_DP !! Universal gravitational constant in the system units - real(DP) :: inv_c2 = -1.0_DP !! Inverse speed of light squared in the system units - character(NAMELEN) :: interaction_loops = "ADAPTIVE" !! Method used to compute interaction loops. Options are "TRIANGULAR", "FLAT", or "ADAPTIVE" - character(NAMELEN) :: encounter_check_plpl = "ADAPTIVE" !! Method used to compute pl-pl encounter checks. Options are "TRIANGULAR", "SORTSWEEP", or "ADAPTIVE" - character(NAMELEN) :: encounter_check_pltp = "ADAPTIVE" !! Method used to compute pl-tp encounter checks. Options are "TRIANGULAR", "SORTSWEEP", or "ADAPTIVE" + character(STRMAX) :: integrator = UNKNOWN_INTEGRATOR !! Symbolic name of the nbody integrator used + character(STRMAX) :: param_file_name = "param.in" !! The default name of the parameter input file + integer(I4B) :: maxid = -1 !! The current maximum particle id number + integer(I4B) :: maxid_collision = 0 !! The current maximum collision id number + real(DP) :: t0 = 0.0_DP !! Integration reference time + real(DP) :: tstart = -1.0_DP !! Integration start time + real(DP) :: tstop = -1.0_DP !! Integration stop time + real(DP) :: dt = -1.0_DP !! Time step + integer(I8B) :: iloop = 0_I8B !! Main loop counter + integer(I8B) :: ioutput = 0_I8B !! Output counter + character(STRMAX) :: incbfile = CB_INFILE !! Name of input file for the central body + character(STRMAX) :: inplfile = PL_INFILE !! Name of input file for massive bodies + character(STRMAX) :: intpfile = TP_INFILE !! Name of input file for test particles + character(STRMAX) :: in_netcdf = NC_INFILE !! Name of system input file for NetCDF input + character(STRMAX) :: in_type = "ASCII" !! Data representation type of input data files + character(STRMAX) :: in_form = "XV" !! Format of input data files ("EL" or "XV") + integer(I4B) :: istep_out = -1 !! Number of time steps between saved outputs + character(STRMAX) :: outfile = NETCDF_OUTFILE !! Name of output binary file + character(STRMAX) :: out_type = "NETCDF_DOUBLE" !! Binary format of output file + character(STRMAX) :: out_form = "XVEL" !! Data to write to output file + character(STRMAX) :: out_stat = 'NEW' !! Open status for output binary file + integer(I4B) :: dump_cadence = 10 !! Number of output steps between dumping simulation data to file + real(DP) :: rmin = -1.0_DP !! Minimum heliocentric radius for test particle + real(DP) :: rmax = -1.0_DP !! Maximum heliocentric radius for test particle + real(DP) :: rmaxu = -1.0_DP !! Maximum unbound heliocentric radius for test particle + real(DP) :: qmin = -1.0_DP !! Minimum pericenter distance for test particle + character(STRMAX) :: qmin_coord = 'HELIO' !! Coordinate frame to use for qmin + real(DP) :: qmin_alo = -1.0_DP !! Minimum semimajor axis for qmin + real(DP) :: qmin_ahi = -1.0_DP !! Maximum semimajor axis for qmin + real(QP) :: MU2KG = -1.0_QP !! Converts mass units to grams + real(QP) :: TU2S = -1.0_QP !! Converts time units to seconds + real(QP) :: DU2M = -1.0_QP !! Converts distance unit to centimeters + real(DP) :: GU = -1.0_DP !! Universal gravitational constant in the system units + real(DP) :: inv_c2 = -1.0_DP !! Inverse speed of light squared in the system units + character(NAMELEN) :: interaction_loops = "ADAPTIVE" !! Method used to compute interaction loops. Options are "TRIANGULAR", "FLAT", or "ADAPTIVE" + character(NAMELEN) :: encounter_check_plpl = "ADAPTIVE" !! Method used to compute pl-pl encounter checks. Options are "TRIANGULAR", "SORTSWEEP", or "ADAPTIVE" + character(NAMELEN) :: encounter_check_pltp = "ADAPTIVE" !! Method used to compute pl-tp encounter checks. Options are "TRIANGULAR", "SORTSWEEP", or "ADAPTIVE" ! The following are used internally, and are not set by the user, but instead are determined by the input value of INTERACTION_LOOPS logical :: lflatten_interactions = .false. !! Use the flattened upper triangular matrix for pl-pl interaction loops @@ -89,18 +89,18 @@ module swiftest_classes logical :: ltides = .false. !! Include tidal dissipation ! Initial values to pass to the energy report subroutine (usually only used in the case of a restart, otherwise these will be updated with initial conditions values) - real(DP) :: Eorbit_orig = 0.0_DP !! Initial orbital energy - real(DP) :: GMtot_orig = 0.0_DP !! Initial system mass - real(DP), dimension(NDIM) :: Ltot_orig = 0.0_DP !! Initial total angular momentum vector - real(DP), dimension(NDIM) :: Lorbit_orig = 0.0_DP !! Initial orbital angular momentum - real(DP), dimension(NDIM) :: Lspin_orig = 0.0_DP !! Initial spin angular momentum vector - real(DP), dimension(NDIM) :: Lescape = 0.0_DP !! Angular momentum of bodies that escaped the system (used for bookeeping) - real(DP) :: GMescape = 0.0_DP !! Mass of bodies that escaped the system (used for bookeeping) - real(DP) :: Ecollisions = 0.0_DP !! Energy lost from system due to collisions - real(DP) :: Euntracked = 0.0_DP !! Energy gained from system due to escaped bodies + real(DP) :: Eorbit_orig = 0.0_DP !! Initial orbital energy + real(DP) :: GMtot_orig = 0.0_DP !! Initial system mass + real(DP), dimension(NDIM) :: Ltot_orig = 0.0_DP !! Initial total angular momentum vector + real(DP), dimension(NDIM) :: Lorbit_orig = 0.0_DP !! Initial orbital angular momentum + real(DP), dimension(NDIM) :: Lspin_orig = 0.0_DP !! Initial spin angular momentum vector + real(DP), dimension(NDIM) :: Lescape = 0.0_DP !! Angular momentum of bodies that escaped the system (used for bookeeping) + real(DP) :: GMescape = 0.0_DP !! Mass of bodies that escaped the system (used for bookeeping) + real(DP) :: Ecollisions = 0.0_DP !! Energy lost from system due to collisions + real(DP) :: Euntracked = 0.0_DP !! Energy gained from system due to escaped bodies logical :: lfirstenergy = .true. !! This is the first time computing energe - logical :: lfirstkick = .true. !! Initiate the first kick in a symplectic step - logical :: lrestart = .false. !! Indicates whether or not this is a restarted run + logical :: lfirstkick = .true. !! Initiate the first kick in a symplectic step + logical :: lrestart = .false. !! Indicates whether or not this is a restarted run character(len=:), allocatable :: display_style !! Style of the output display {"STANDARD", "COMPACT"}). Default is "STANDARD" integer(I4B) :: display_unit !! File unit number for display (either to stdout or to a log file) @@ -159,7 +159,7 @@ module swiftest_classes !******************************************************************************************************************************** ! swiftest_cb class definitions and methods !******************************************************************************************************************************** - !> A concrete lass for the central body in a Swiftest simulation + !> An abstract class for a generic central body in a Swiftest simulation type, abstract, extends(swiftest_base) :: swiftest_cb type(swiftest_particle_info) :: info !! Particle metadata information integer(I4B) :: id = 0 !! External identifier (unique) @@ -346,6 +346,7 @@ module swiftest_classes class(swiftest_tp), allocatable :: tp !! Test particle data structure class(swiftest_tp), allocatable :: tp_discards !! Discarded test particle data structure class(swiftest_pl), allocatable :: pl_discards !! Discarded massive body particle data structure + real(DP) :: t = -1.0_DP !! Integration current time real(DP) :: GMtot = 0.0_DP !! Total system mass - used for barycentric coordinate conversion real(DP) :: ke_orbit = 0.0_DP !! System orbital kinetic energy real(DP) :: ke_spin = 0.0_DP !! System spin kinetic energy @@ -413,12 +414,25 @@ module swiftest_classes procedure :: get_energy_and_momentum => util_get_energy_momentum_system !! Calculates the total system energy and momentum procedure :: rescale => util_rescale_system !! Rescales the system into a new set of units procedure :: validate_ids => util_valid_id_system !! Validate the numerical ids passed to the system and save the maximum value - generic :: write_frame => write_frame_system, write_frame_netcdf !! Generic method call for reading a frame of output data + generic :: write_frame => write_frame_system, write_frame_netcdf !! Generic method call for reading a frame of output data end type swiftest_nbody_system + type storage_frame + class(swiftest_nbody_system), allocatable :: system + contains + procedure :: store => util_copy_store_system !! Stores a snapshot of the nbody system so that later it can be retrieved for saving to file. + generic :: assignment(=) => store + end type + + type, extends(swiftest_base) :: swiftest_storage(nframes) + integer(I8B), len :: nframes + !! A class that that is used to store simulation history data between file output + type(storage_frame), dimension(nframes) :: frame + contains + procedure :: dump => io_dump_system_storage + end type swiftest_storage abstract interface - subroutine abstract_accel(self, system, param, t, lbeg) import swiftest_body, swiftest_nbody_system, swiftest_parameters, DP class(swiftest_body), intent(inout) :: self !! Swiftest body data structure @@ -611,6 +625,13 @@ module subroutine io_dump_system(self, param) class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters end subroutine io_dump_system + + module subroutine io_dump_system_storage(self, param) + implicit none + class(swiftest_storage(*)), intent(inout) :: self !! Swiftest simulation history storage object + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters + end subroutine io_dump_system_storage + module subroutine io_get_args(integrator, param_file_name, display_style) implicit none character(len=:), allocatable, intent(inout) :: integrator !! Symbolic code of the requested integrator @@ -1222,6 +1243,12 @@ module subroutine util_copy_particle_info_arr(source, dest, idx) integer(I4B), dimension(:), intent(in), optional :: idx !! Optional array of indices to draw the source object end subroutine util_copy_particle_info_arr + module subroutine util_copy_store_system(self, system) + implicit none + class(storage_frame), intent(inout) :: self !! Swiftest storage frame object + class(swiftest_nbody_system), intent(in) :: system !! Swiftest n-body system object + end subroutine util_copy_store_system + module subroutine util_dealloc_body(self) implicit none class(swiftest_body), intent(inout) :: self diff --git a/src/netcdf/netcdf.f90 b/src/netcdf/netcdf.f90 index 89824e4ec..abd03f1d6 100644 --- a/src/netcdf/netcdf.f90 +++ b/src/netcdf/netcdf.f90 @@ -793,7 +793,7 @@ module subroutine netcdf_read_hdr_system(self, iu, param) tslot = int(param%ioutput, kind=I4B) + 1 call check( nf90_inquire_dimension(iu%ncid, iu%id_dimid, len=idmax), "netcdf_read_frame_system nf90_inquire_dimension id_dimid" ) - call check( nf90_get_var(iu%ncid, iu%time_varid, param%t, start=[tslot]), "netcdf_read_hdr_system nf90_getvar time_varid" ) + call check( nf90_get_var(iu%ncid, iu%time_varid, self%t, start=[tslot]), "netcdf_read_hdr_system nf90_getvar time_varid" ) allocate(gmtemp(idmax)) allocate(tpmask(idmax)) @@ -1428,7 +1428,7 @@ module subroutine netcdf_write_hdr_system(self, iu, param) tslot = int(param%ioutput, kind=I4B) + 1 - call check( nf90_put_var(iu%ncid, iu%time_varid, param%t, start=[tslot]), "netcdf_write_hdr_system nf90_put_var time_varid" ) + call check( nf90_put_var(iu%ncid, iu%time_varid, self%t, start=[tslot]), "netcdf_write_hdr_system nf90_put_var time_varid" ) call check( nf90_put_var(iu%ncid, iu%npl_varid, self%pl%nbody, start=[tslot]), "netcdf_write_hdr_system nf90_put_var npl_varid" ) call check( nf90_put_var(iu%ncid, iu%ntp_varid, self%tp%nbody, start=[tslot]), "netcdf_write_hdr_system nf90_put_var ntp_varid" ) select type(pl => self%pl) diff --git a/src/rmvs/rmvs_discard.f90 b/src/rmvs/rmvs_discard.f90 index 732cbdea0..60be2f6b0 100644 --- a/src/rmvs/rmvs_discard.f90 +++ b/src/rmvs/rmvs_discard.f90 @@ -29,7 +29,7 @@ module subroutine rmvs_discard_tp(self, system, param) if (self%nbody == 0) return - associate(tp => self, ntp => self%nbody, pl => system%pl, t => param%t) + associate(tp => self, ntp => self%nbody, pl => system%pl, t => system%t) do i = 1, ntp associate(iplperP => tp%plperP(i)) if ((tp%status(i) == ACTIVE) .and. (tp%lperi(i))) then diff --git a/src/setup/setup.f90 b/src/setup/setup.f90 index ae157eaae..655d15b58 100644 --- a/src/setup/setup.f90 +++ b/src/setup/setup.f90 @@ -91,9 +91,7 @@ module subroutine setup_finalize_system(self, param) class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters associate(system => self) - if ((param%out_type == "NETCDF_FLOAT") .or. (param%out_type == "NETCDF_DOUBLE")) then - call param%nciu%close() - end if + call param%nciu%close() end associate return diff --git a/src/setup/symba_collision.f90 b/src/setup/symba_collision.f90 new file mode 100644 index 000000000..e839af1de --- /dev/null +++ b/src/setup/symba_collision.f90 @@ -0,0 +1,1090 @@ +!! Copyright 2022 - David Minton, Carlisle Wishard, Jennifer Pouplin, Jake Elliott, & Dana Singh +!! This file is part of Swiftest. +!! Swiftest is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License +!! as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version. +!! Swiftest is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty +!! of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. +!! You should have received a copy of the GNU General Public License along with Swiftest. +!! If not, see: https://www.gnu.org/licenses. + +submodule (symba_classes) s_symba_collision + use swiftest + +contains + + module function symba_collision_casedisruption(system, param, colliders, frag) result(status) + !! author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton + !! + !! Create the fragments resulting from a non-catastrophic disruption collision + !! + implicit none + ! Arguments + class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object + class(symba_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA additions + class(fraggle_colliders), intent(inout) :: colliders !! Fraggle colliders object + class(fraggle_fragments), intent(inout) :: frag !! Fraggle fragmentation system object + ! Result + integer(I4B) :: status !! Status flag assigned to this outcome + ! Internals + integer(I4B) :: i, nfrag + logical :: lfailure + character(len=STRMAX) :: message + + select case(frag%regime) + case(COLLRESOLVE_REGIME_DISRUPTION) + message = "Disruption between" + case(COLLRESOLVE_REGIME_SUPERCATASTROPHIC) + message = "Supercatastrophic disruption between" + end select + call symba_collision_collider_message(system%pl, colliders%idx, message) + call io_log_one_message(FRAGGLE_LOG_OUT, message) + + ! Collisional fragments will be uniformly distributed around the pre-impact barycenter + call frag%set_mass_dist(colliders, param) + + ! Generate the position and velocity distributions of the fragments + call frag%generate_fragments(colliders, system, param, lfailure) + + if (lfailure) then + call io_log_one_message(FRAGGLE_LOG_OUT, "No fragment solution found, so treat as a pure hit-and-run") + status = ACTIVE + nfrag = 0 + select type(pl => system%pl) + class is (symba_pl) + pl%status(colliders%idx(:)) = status + pl%ldiscard(colliders%idx(:)) = .false. + pl%lcollision(colliders%idx(:)) = .false. + end select + else + ! Populate the list of new bodies + nfrag = frag%nbody + write(message, *) nfrag + call io_log_one_message(FRAGGLE_LOG_OUT, "Generating " // trim(adjustl(message)) // " fragments") + select case(frag%regime) + case(COLLRESOLVE_REGIME_DISRUPTION) + status = DISRUPTION + case(COLLRESOLVE_REGIME_SUPERCATASTROPHIC) + status = SUPERCATASTROPHIC + end select + frag%id(1:nfrag) = [(i, i = param%maxid + 1, param%maxid + nfrag)] + param%maxid = frag%id(nfrag) + call symba_collision_mergeaddsub(system, param, colliders, frag, status) + end if + + return + end function symba_collision_casedisruption + + + module function symba_collision_casehitandrun(system, param, colliders, frag) result(status) + !! author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton + !! + !! Create the fragments resulting from a non-catastrophic hit-and-run collision + !! + implicit none + ! Arguments + class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object + class(symba_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA additions + class(fraggle_colliders), intent(inout) :: colliders !! Fraggle colliders object + class(fraggle_fragments), intent(inout) :: frag !! Fraggle fragmentation system object + ! Result + integer(I4B) :: status !! Status flag assigned to this outcom + ! Internals + integer(I4B) :: i, ibiggest, nfrag, jtarg, jproj + logical :: lpure + character(len=STRMAX) :: message + + message = "Hit and run between" + call symba_collision_collider_message(system%pl, colliders%idx, message) + call io_log_one_message(FRAGGLE_LOG_OUT, trim(adjustl(message))) + + if (colliders%mass(1) > colliders%mass(2)) then + jtarg = 1 + jproj = 2 + else + jtarg = 2 + jproj = 1 + end if + + if (frag%mass_dist(2) > 0.9_DP * colliders%mass(jproj)) then ! Pure hit and run, so we'll just keep the two bodies untouched + call io_log_one_message(FRAGGLE_LOG_OUT, "Pure hit and run. No new fragments generated.") + nfrag = 0 + lpure = .true. + else ! Imperfect hit and run, so we'll keep the largest body and destroy the other + lpure = .false. + call frag%set_mass_dist(colliders, param) + + ! Generate the position and velocity distributions of the fragments + call frag%generate_fragments(colliders, system, param, lpure) + + if (lpure) then + call io_log_one_message(FRAGGLE_LOG_OUT, "Should have been a pure hit and run instead") + nfrag = 0 + else + nfrag = frag%nbody + write(message, *) nfrag + call io_log_one_message(FRAGGLE_LOG_OUT, "Generating " // trim(adjustl(message)) // " fragments") + end if + end if + if (lpure) then ! Reset these bodies back to being active so that nothing further is done to them + status = HIT_AND_RUN_PURE + select type(pl => system%pl) + class is (symba_pl) + pl%status(colliders%idx(:)) = ACTIVE + pl%ldiscard(colliders%idx(:)) = .false. + pl%lcollision(colliders%idx(:)) = .false. + end select + else + ibiggest = colliders%idx(maxloc(system%pl%Gmass(colliders%idx(:)), dim=1)) + frag%id(1) = system%pl%id(ibiggest) + frag%id(2:nfrag) = [(i, i = param%maxid + 1, param%maxid + nfrag - 1)] + param%maxid = frag%id(nfrag) + status = HIT_AND_RUN_DISRUPT + call symba_collision_mergeaddsub(system, param, colliders, frag, status) + end if + + return + end function symba_collision_casehitandrun + + + module function symba_collision_casemerge(system, param, colliders, frag) result(status) + !! author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton + !! + !! Merge massive bodies. + !! + !! Adapted from David E. Kaufmann's Swifter routines symba_merge_pl.f90 and symba_discard_merge_pl.f90 + !! + !! Adapted from Hal Levison's Swift routines symba5_merge.f and discard_mass_merge.f + implicit none + ! Arguments + class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object + class(symba_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA additions + class(fraggle_colliders), intent(inout) :: colliders !! Fraggle colliders object + class(fraggle_fragments), intent(inout) :: frag !! Fraggle fragmentation system object + ! Result + integer(I4B) :: status !! Status flag assigned to this outcome + ! Internals + integer(I4B) :: i, j, k, ibiggest + real(DP) :: pe + real(DP), dimension(NDIM) :: L_spin_new + character(len=STRMAX) :: message + + message = "Merging" + call symba_collision_collider_message(system%pl, colliders%idx, message) + call io_log_one_message(FRAGGLE_LOG_OUT, message) + + select type(pl => system%pl) + class is (symba_pl) + + call frag%set_mass_dist(colliders, param) + ibiggest = colliders%idx(maxloc(pl%Gmass(colliders%idx(:)), dim=1)) + frag%id(1) = pl%id(ibiggest) + frag%xb(:,1) = frag%xbcom(:) + frag%vb(:,1) = frag%vbcom(:) + + if (param%lrotation) then + ! Conserve angular momentum by putting pre-impact orbital momentum into spin of the new body + L_spin_new(:) = colliders%L_orbit(:,1) + colliders%L_orbit(:,2) + colliders%L_spin(:,1) + colliders%L_spin(:,2) + + ! Assume prinicpal axis rotation on 3rd Ip axis + frag%rot(:,1) = L_spin_new(:) / (frag%Ip(3,1) * frag%mass(1) * frag%radius(1)**2) + else ! If spin is not enabled, we will consider the lost pre-collision angular momentum as "escaped" and add it to our bookkeeping variable + param%Lescape(:) = param%Lescape(:) + colliders%L_orbit(:,1) + colliders%L_orbit(:,2) + end if + + ! Keep track of the component of potential energy due to the pre-impact colliders%idx for book-keeping + pe = 0.0_DP + do j = 1, colliders%ncoll + do i = j + 1, colliders%ncoll + pe = pe - pl%Gmass(i) * pl%mass(j) / norm2(pl%xb(:, i) - pl%xb(:, j)) + end do + end do + system%Ecollisions = system%Ecollisions + pe + system%Euntracked = system%Euntracked - pe + + ! Update any encounter lists that have the removed bodies in them so that they instead point to the new + do k = 1, system%plplenc_list%nenc + do j = 1, colliders%ncoll + i = colliders%idx(j) + if (i == ibiggest) cycle + if (system%plplenc_list%id1(k) == pl%id(i)) then + system%plplenc_list%id1(k) = pl%id(ibiggest) + system%plplenc_list%index1(k) = i + end if + if (system%plplenc_list%id2(k) == pl%id(i)) then + system%plplenc_list%id2(k) = pl%id(ibiggest) + system%plplenc_list%index2(k) = i + end if + if (system%plplenc_list%id1(k) == system%plplenc_list%id2(k)) system%plplenc_list%status(k) = INACTIVE + end do + end do + + status = MERGED + + call symba_collision_mergeaddsub(system, param, colliders, frag, status) + + end select + + return + end function symba_collision_casemerge + + + subroutine symba_collision_collider_message(pl, collidx, collider_message) + !! author: David A. Minton + !! + !! Prints a nicely formatted message about which bodies collided, including their names and ids. + !! This subroutine appends the body names and ids to an input message. + implicit none + ! Arguments + class(swiftest_pl), intent(in) :: pl !! Swiftest massive body object + integer(I4B), dimension(:), intent(in) :: collidx !! Index of collisional colliders%idx members + character(*), intent(inout) :: collider_message !! The message to print to the screen. + ! Internals + integer(I4B) :: i, n + character(len=STRMAX) :: idstr + + n = size(collidx) + if (n == 0) return + + do i = 1, n + if (i > 1) collider_message = trim(adjustl(collider_message)) // " and " + collider_message = " " // trim(adjustl(collider_message)) // " " // trim(adjustl(pl%info(collidx(i))%name)) + write(idstr, '(I10)') pl%id(collidx(i)) + collider_message = trim(adjustl(collider_message)) // " (" // trim(adjustl(idstr)) // ") " + end do + + return + end subroutine symba_collision_collider_message + + + module function symba_collision_check_encounter(self, system, param, t, dt, irec) result(lany_collision) + !! author: David A. Minton + !! + !! Check for merger between massive bodies and test particles in SyMBA + !! + !! Adapted from David E. Kaufmann's Swifter routine symba_merge.f90 and symba_merge_tp.f90 + !! + !! Adapted from Hal Levison's Swift routine symba5_merge.f + implicit none + ! Arguments + class(symba_encounter), intent(inout) :: self !! SyMBA pl-tp encounter list object + class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object + class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters + real(DP), intent(in) :: t !! current time + real(DP), intent(in) :: dt !! step size + integer(I4B), intent(in) :: irec !! Current recursion level + ! Result + logical :: lany_collision !! Returns true if cany pair of encounters resulted in a collision + ! Internals + logical, dimension(:), allocatable :: lcollision, lmask + real(DP), dimension(NDIM) :: xr, vr + integer(I4B) :: i, j, k, nenc + real(DP) :: rlim, Gmtot + logical :: isplpl + character(len=STRMAX) :: timestr, idstri, idstrj, message + class(symba_encounter), allocatable :: tmp + + lany_collision = .false. + if (self%nenc == 0) return + + select type(self) + class is (symba_plplenc) + isplpl = .true. + class default + isplpl = .false. + end select + + select type(pl => system%pl) + class is (symba_pl) + select type(tp => system%tp) + class is (symba_tp) + nenc = self%nenc + allocate(lmask(nenc)) + lmask(:) = ((self%status(1:nenc) == ACTIVE) .and. (pl%levelg(self%index1(1:nenc)) >= irec)) + if (isplpl) then + lmask(:) = lmask(:) .and. (pl%levelg(self%index2(1:nenc)) >= irec) + else + lmask(:) = lmask(:) .and. (tp%levelg(self%index2(1:nenc)) >= irec) + end if + if (.not.any(lmask(:))) return + + allocate(lcollision(nenc)) + lcollision(:) = .false. + + if (isplpl) then + do concurrent(k = 1:nenc, lmask(k)) + i = self%index1(k) + j = self%index2(k) + xr(:) = pl%xh(:, i) - pl%xh(:, j) + vr(:) = pl%vb(:, i) - pl%vb(:, j) + rlim = pl%radius(i) + pl%radius(j) + Gmtot = pl%Gmass(i) + pl%Gmass(j) + lcollision(k) = symba_collision_check_one(xr(1), xr(2), xr(3), vr(1), vr(2), vr(3), & + Gmtot, rlim, dt, self%lvdotr(k)) + end do + else + do concurrent(k = 1:nenc, lmask(k)) + i = self%index1(k) + j = self%index2(k) + xr(:) = pl%xh(:, i) - tp%xh(:, j) + vr(:) = pl%vb(:, i) - tp%vb(:, j) + lcollision(k) = symba_collision_check_one(xr(1), xr(2), xr(3), vr(1), vr(2), vr(3), & + pl%Gmass(i), pl%radius(i), dt, self%lvdotr(k)) + end do + end if + + lany_collision = any(lcollision(:)) + if (lany_collision) then + call pl%xh2xb(system%cb) ! Update the central body barycenteric position vector to get us out of DH and into bary + do k = 1, nenc + i = self%index1(k) + j = self%index2(k) + if (lcollision(k)) self%status(k) = COLLISION + self%t(k) = t + self%x1(:,k) = pl%xh(:,i) + system%cb%xb(:) + self%v1(:,k) = pl%vb(:,i) + if (isplpl) then + self%x2(:,k) = pl%xh(:,j) + system%cb%xb(:) + self%v2(:,k) = pl%vb(:,j) + if (lcollision(k)) then + ! Check to see if either of these bodies has been involved with a collision before, and if so, make this a collisional colliders%idx + if (pl%lcollision(i) .or. pl%lcollision(j)) call pl%make_colliders([i,j]) + + ! Set the collision flag for these to bodies to true in case they become involved in another collision later in the step + pl%lcollision([i, j]) = .true. + pl%status([i, j]) = COLLISION + call pl%info(i)%set_value(status="COLLISION", discard_time=t, discard_xh=pl%xh(:,i), discard_vh=pl%vh(:,i)) + call pl%info(j)%set_value(status="COLLISION", discard_time=t, discard_xh=pl%xh(:,j), discard_vh=pl%vh(:,j)) + end if + else + self%x2(:,k) = tp%xh(:,j) + system%cb%xb(:) + self%v2(:,k) = tp%vb(:,j) + if (lcollision(k)) then + tp%status(j) = DISCARDED_PLR + tp%ldiscard(j) = .true. + write(idstri, *) pl%id(i) + write(idstrj, *) tp%id(j) + write(timestr, *) t + call tp%info(j)%set_value(status="DISCARDED_PLR", discard_time=t, discard_xh=tp%xh(:,j), discard_vh=tp%vh(:,j)) + write(message, *) "Particle " // trim(adjustl(tp%info(j)%name)) // " (" // trim(adjustl(idstrj)) // ")" & + // " collided with massive body " // trim(adjustl(pl%info(i)%name)) // " (" // trim(adjustl(idstri)) // ")" & + // " at t = " // trim(adjustl(timestr)) + call io_log_one_message(FRAGGLE_LOG_OUT, message) + end if + end if + end do + end if + end select + end select + + ! Extract the pl-pl or pl-tp encounter list and return the pl-pl or pl-tp collision_list + if (lany_collision) then + select type(self) + class is (symba_plplenc) + call self%extract_collisions(system, param) + class default + allocate(tmp, mold=self) + call self%spill(tmp, lcollision, ldestructive=.true.) ! Remove this encounter pair from the encounter list + end select + end if + + return + end function symba_collision_check_encounter + + + pure elemental function symba_collision_check_one(xr, yr, zr, vxr, vyr, vzr, Gmtot, rlim, dt, lvdotr) result(lcollision) + !! author: David A. Minton + !! + !! Check for a merger between a single pair of particles + !! + !! Adapted from David E. Kaufmann's Swifter routines symba_merge_tp.f90 and symba_merge_pl.f90 + !! + !! Adapted from Hal Levison's Swift routine symba5_merge.f + implicit none + ! Arguments + real(DP), intent(in) :: xr, yr, zr !! Relative position vector components + real(DP), intent(in) :: vxr, vyr, vzr !! Relative velocity vector components + real(DP), intent(in) :: Gmtot !! Sum of G*mass of colliding bodies + real(DP), intent(in) :: rlim !! Collision limit - Typically the sum of the radii of colliding bodies + real(DP), intent(in) :: dt !! Step size + logical, intent(in) :: lvdotr !! Logical flag indicating that these two bodies are approaching in the current substep + ! Result + logical :: lcollision !! Logical flag indicating whether these two bodies will collide or not + ! Internals + real(DP) :: r2, rlim2, a, e, q, vdotr, tcr2, dt2 + + r2 = xr**2 + yr**2 + zr**2 + rlim2 = rlim**2 + + if (r2 <= rlim2) then ! checks if bodies are actively colliding in this time step + lcollision = .true. + else ! if they are not actively colliding in this time step, checks if they are going to collide next time step based on velocities and q + lcollision = .false. + vdotr = xr * vxr + yr * vyr + zr * vzr + if (lvdotr .and. (vdotr > 0.0_DP)) then + tcr2 = r2 / (vxr**2 + vyr**2 + vzr**2) + dt2 = dt**2 + if (tcr2 <= dt2) then + call orbel_xv2aeq(Gmtot, xr, yr, zr, vxr, vyr, vzr, a, e, q) + lcollision = (q < rlim) + end if + end if + end if + + return + end function symba_collision_check_one + + + function symba_collision_consolidate_colliders(pl, cb, param, idx_parent, colliders) result(lflag) + !! author: David A. Minton + !! + !! Loops through the pl-pl collision list and groups families together by index. Outputs the indices of all colliders%idx members, + !! and pairs of quantities (x and v vectors, mass, radius, L_spin, and Ip) that can be used to resolve the collisional outcome. + implicit none + ! Arguments + class(symba_pl), intent(inout) :: pl !! SyMBA massive body object + class(symba_cb), intent(inout) :: cb !! SyMBA central body object + class(symba_parameters), intent(in) :: param !! Current run configuration parameters with SyMBA additions + integer(I4B), dimension(2), intent(inout) :: idx_parent !! Index of the two bodies considered the "parents" of the collision + class(fraggle_colliders), intent(out) :: colliders + ! Result + logical :: lflag !! Logical flag indicating whether a colliders%idx was successfully created or not + ! Internals + type collidx_array + integer(I4B), dimension(:), allocatable :: id + integer(I4B), dimension(:), allocatable :: idx + end type collidx_array + type(collidx_array), dimension(2) :: parent_child_index_array + integer(I4B), dimension(2) :: nchild + integer(I4B) :: i, j, ncolliders, idx_child + real(DP), dimension(2) :: volume, density + real(DP) :: mchild, volchild + real(DP), dimension(NDIM) :: xc, vc, xcom, vcom, xchild, vchild, xcrossv + real(DP), dimension(NDIM,2) :: mxc, vcc + + nchild(:) = pl%kin(idx_parent(:))%nchild + ! If all of these bodies share a parent, but this is still a unique collision, move the last child + ! out of the parent's position and make it the secondary body + if (idx_parent(1) == idx_parent(2)) then + if (nchild(1) == 0) then ! There is only one valid body recorded in this pair (this could happen due to restructuring of the kinship relationships, though it should be rare) + lflag = .false. + call pl%reset_kinship([idx_parent(1)]) + return + end if + idx_parent(2) = pl%kin(idx_parent(1))%child(nchild(1)) + nchild(1) = nchild(1) - 1 + nchild(2) = 0 + pl%kin(idx_parent(:))%nchild = nchild(:) + pl%kin(idx_parent(2))%parent = idx_parent(1) + end if + + colliders%mass(:) = pl%mass(idx_parent(:)) ! Note: This is meant to mass, not G*mass, as the collisional regime determination uses mass values that will be converted to Si + colliders%radius(:) = pl%radius(idx_parent(:)) + volume(:) = (4.0_DP / 3.0_DP) * PI * colliders%radius(:)**3 + + ! Group together the ids and indexes of each collisional parent and its children + do j = 1, 2 + allocate(parent_child_index_array(j)%idx(nchild(j)+ 1)) + allocate(parent_child_index_array(j)%id(nchild(j)+ 1)) + associate(idx_arr => parent_child_index_array(j)%idx, & + id_arr => parent_child_index_array(j)%id, & + ncj => nchild(j), & + pl => pl, & + plkinj => pl%kin(idx_parent(j))) + idx_arr(1) = idx_parent(j) + if (ncj > 0) idx_arr(2:ncj + 1) = plkinj%child(1:ncj) + id_arr(:) = pl%id(idx_arr(:)) + end associate + end do + + ! Consolidate the groups of collsional parents with any children they may have into a single "colliders%idx" index array + ncolliders = 2 + sum(nchild(:)) + allocate(colliders%idx(ncolliders)) + colliders%idx = [parent_child_index_array(1)%idx(:),parent_child_index_array(2)%idx(:)] + + colliders%ncoll = count(pl%lcollision(colliders%idx(:))) + colliders%idx = pack(colliders%idx(:), pl%lcollision(colliders%idx(:))) + colliders%L_spin(:,:) = 0.0_DP + colliders%Ip(:,:) = 0.0_DP + + ! Find the barycenter of each body along with its children, if it has any + do j = 1, 2 + colliders%xb(:, j) = pl%xh(:, idx_parent(j)) + cb%xb(:) + colliders%vb(:, j) = pl%vb(:, idx_parent(j)) + ! Assume principal axis rotation about axis corresponding to highest moment of inertia (3rd Ip) + if (param%lrotation) then + colliders%Ip(:, j) = colliders%mass(j) * pl%Ip(:, idx_parent(j)) + colliders%L_spin(:, j) = colliders%Ip(3, j) * colliders%radius(j)**2 * pl%rot(:, idx_parent(j)) + end if + + if (nchild(j) > 0) then + do i = 1, nchild(j) ! Loop over all children and take the mass weighted mean of the properties + idx_child = parent_child_index_array(j)%idx(i + 1) + if (.not. pl%lcollision(idx_child)) cycle + mchild = pl%mass(idx_child) + xchild(:) = pl%xh(:, idx_child) + cb%xb(:) + vchild(:) = pl%vb(:, idx_child) + volchild = (4.0_DP / 3.0_DP) * PI * pl%radius(idx_child)**3 + volume(j) = volume(j) + volchild + ! Get angular momentum of the child-parent pair and add that to the spin + ! Add the child's spin + if (param%lrotation) then + xcom(:) = (colliders%mass(j) * colliders%xb(:,j) + mchild * xchild(:)) / (colliders%mass(j) + mchild) + vcom(:) = (colliders%mass(j) * colliders%vb(:,j) + mchild * vchild(:)) / (colliders%mass(j) + mchild) + xc(:) = colliders%xb(:, j) - xcom(:) + vc(:) = colliders%vb(:, j) - vcom(:) + xcrossv(:) = xc(:) .cross. vc(:) + colliders%L_spin(:, j) = colliders%L_spin(:, j) + colliders%mass(j) * xcrossv(:) + + xc(:) = xchild(:) - xcom(:) + vc(:) = vchild(:) - vcom(:) + xcrossv(:) = xc(:) .cross. vc(:) + colliders%L_spin(:, j) = colliders%L_spin(:, j) + mchild * xcrossv(:) + + colliders%L_spin(:, j) = colliders%L_spin(:, j) + mchild * pl%Ip(3, idx_child) & + * pl%radius(idx_child)**2 & + * pl%rot(:, idx_child) + colliders%Ip(:, j) = colliders%Ip(:, j) + mchild * pl%Ip(:, idx_child) + end if + + ! Merge the child and parent + colliders%mass(j) = colliders%mass(j) + mchild + colliders%xb(:, j) = xcom(:) + colliders%vb(:, j) = vcom(:) + end do + end if + density(j) = colliders%mass(j) / volume(j) + colliders%radius(j) = (3 * volume(j) / (4 * PI))**(1.0_DP / 3.0_DP) + if (param%lrotation) colliders%Ip(:, j) = colliders%Ip(:, j) / colliders%mass(j) + end do + lflag = .true. + + xcom(:) = (colliders%mass(1) * colliders%xb(:, 1) + colliders%mass(2) * colliders%xb(:, 2)) / sum(colliders%mass(:)) + vcom(:) = (colliders%mass(1) * colliders%vb(:, 1) + colliders%mass(2) * colliders%vb(:, 2)) / sum(colliders%mass(:)) + mxc(:, 1) = colliders%mass(1) * (colliders%xb(:, 1) - xcom(:)) + mxc(:, 2) = colliders%mass(2) * (colliders%xb(:, 2) - xcom(:)) + vcc(:, 1) = colliders%vb(:, 1) - vcom(:) + vcc(:, 2) = colliders%vb(:, 2) - vcom(:) + colliders%L_orbit(:,:) = mxc(:,:) .cross. vcc(:,:) + + ! Destroy the kinship relationships for all members of this colliders%idx + call pl%reset_kinship(colliders%idx(:)) + + return + end function symba_collision_consolidate_colliders + + + module subroutine symba_collision_encounter_extract_collisions(self, system, param) + !! author: David A. Minton + !! + !! Processes the pl-pl encounter list remove only those encounters that led to a collision + !! + implicit none + ! Arguments + class(symba_plplenc), intent(inout) :: self !! SyMBA pl-pl encounter list + class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object + class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters + ! Internals + logical, dimension(:), allocatable :: lplpl_collision + logical, dimension(:), allocatable :: lplpl_unique_parent + integer(I4B), dimension(:), pointer :: plparent + integer(I4B), dimension(:), allocatable :: collision_idx, unique_parent_idx + integer(I4B) :: i, index_coll, ncollisions, nunique_parent, nplplenc + + select type (pl => system%pl) + class is (symba_pl) + associate(plplenc_list => self, idx1 => self%index1, idx2 => self%index2, plparent => pl%kin%parent) + nplplenc = plplenc_list%nenc + allocate(lplpl_collision(nplplenc)) + lplpl_collision(:) = plplenc_list%status(1:nplplenc) == COLLISION + if (.not.any(lplpl_collision)) return + ! Collisions have been detected in this step. So we need to determine which of them are between unique bodies. + + ! Get the subset of pl-pl encounters that lead to a collision + ncollisions = count(lplpl_collision(:)) + allocate(collision_idx(ncollisions)) + collision_idx = pack([(i, i=1, nplplenc)], lplpl_collision) + + ! Get the subset of collisions that involve a unique pair of parents + allocate(lplpl_unique_parent(ncollisions)) + + lplpl_unique_parent(:) = plparent(idx1(collision_idx(:))) /= plparent(idx2(collision_idx(:))) + nunique_parent = count(lplpl_unique_parent(:)) + allocate(unique_parent_idx(nunique_parent)) + unique_parent_idx = pack(collision_idx(:), lplpl_unique_parent(:)) + + ! Scrub all pl-pl collisions involving unique pairs of parents, which will remove all duplicates and leave behind + ! all pairs that have themselves as parents but are not part of the unique parent list. This can hapepn in rare cases + ! due to restructuring of parent/child relationships when there are large numbers of multi-body collisions in a single + ! step + lplpl_unique_parent(:) = .true. + do index_coll = 1, ncollisions + associate(ip1 => plparent(idx1(collision_idx(index_coll))), ip2 => plparent(idx2(collision_idx(index_coll)))) + lplpl_unique_parent(:) = .not. ( any(plparent(idx1(unique_parent_idx(:))) == ip1) .or. & + any(plparent(idx2(unique_parent_idx(:))) == ip1) .or. & + any(plparent(idx1(unique_parent_idx(:))) == ip2) .or. & + any(plparent(idx2(unique_parent_idx(:))) == ip2) ) + end associate + end do + + ! Reassemble collision index list to include only those containing the unique pairs of parents, plus all the non-unique pairs that don't + ! contain a parent body on the unique parent list. + ncollisions = nunique_parent + count(lplpl_unique_parent) + collision_idx = [unique_parent_idx(:), pack(collision_idx(:), lplpl_unique_parent(:))] + + ! Create a mask that contains only the pl-pl encounters that did not result in a collision, and then discard them + lplpl_collision(:) = .false. + lplpl_collision(collision_idx(:)) = .true. + call plplenc_list%spill(system%plplcollision_list, lplpl_collision, ldestructive=.true.) ! Extract any encounters that are not collisions from the list. + end associate + end select + + return + end subroutine symba_collision_encounter_extract_collisions + + + module subroutine symba_collision_make_colliders_pl(self, idx) + !! author: Jennifer L.L. Pouplin, Carlisle A. wishard, and David A. Minton + !! + !! When a single body is involved in more than one collision in a single step, it becomes part of a colliders%idx. + !! The largest body involved in a multi-body collision is the "parent" and all bodies that collide with it are its "children," + !! including those that collide with the children. + !! + !! Adapted from David E. Kaufmann's Swifter routine symba_merge_pl.f90 + !! + !! Adapted from Hal Levison's Swift routine symba5_merge.f + implicit none + ! Arguments + class(symba_pl), intent(inout) :: self !! SyMBA massive body object + integer(I4B), dimension(2), intent(in) :: idx !! Array holding the indices of the two bodies involved in the collision + ! Internals + integer(I4B) :: i, j, index_parent, index_child, p1, p2 + integer(I4B) :: nchild_inherit, nchild_orig, nchild_new + integer(I4B), dimension(:), allocatable :: temp + + associate(pl => self) + p1 = pl%kin(idx(1))%parent + p2 = pl%kin(idx(2))%parent + if (p1 == p2) return ! This is a collision between to children of a shared parent. We will ignore it. + + if (pl%mass(p1) > pl%mass(p2)) then + index_parent = p1 + index_child = p2 + else + index_parent = p2 + index_child = p1 + end if + + ! Expand the child array (or create it if necessary) and copy over the previous lists of children + nchild_orig = pl%kin(index_parent)%nchild + nchild_inherit = pl%kin(index_child)%nchild + nchild_new = nchild_orig + nchild_inherit + 1 + allocate(temp(nchild_new)) + + if (nchild_orig > 0) temp(1:nchild_orig) = pl%kin(index_parent)%child(1:nchild_orig) + ! Find out if the child body has any children of its own. The new parent wil inherit these children + if (nchild_inherit > 0) then + temp(nchild_orig+1:nchild_orig+nchild_inherit) = pl%kin(index_child)%child(1:nchild_inherit) + do i = 1, nchild_inherit + j = pl%kin(index_child)%child(i) + ! Set the childrens' parent to the new parent + pl%kin(j)%parent = index_parent + end do + end if + call pl%reset_kinship([index_child]) + ! Add the new child to its parent + pl%kin(index_child)%parent = index_parent + temp(nchild_new) = index_child + ! Save the new child array to the parent + pl%kin(index_parent)%nchild = nchild_new + call move_alloc(from=temp, to=pl%kin(index_parent)%child) + end associate + + return + end subroutine symba_collision_make_colliders_pl + + + subroutine symba_collision_mergeaddsub(system, param, colliders, frag, status) + !! author: David A. Minton + !! + !! Fills the pl_discards and pl_adds with removed and added bodies + !! + implicit none + ! Arguments + class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object + class(symba_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA additions + class(fraggle_colliders), intent(inout) :: colliders !! Fraggle colliders object + class(fraggle_fragments), intent(inout) :: frag !! Fraggle fragmentation system object + integer(I4B), intent(in) :: status !! Status flag to assign to adds + ! Internals + integer(I4B) :: i, ibiggest, ismallest, iother, nstart, nend, ncolliders, nfrag + logical, dimension(system%pl%nbody) :: lmask + class(symba_pl), allocatable :: plnew, plsub + character(*), parameter :: FRAGFMT = '("Newbody",I0.7)' + character(len=NAMELEN) :: newname + + select type(pl => system%pl) + class is (symba_pl) + select type(pl_discards => system%pl_discards) + class is (symba_merger) + associate(info => pl%info, pl_adds => system%pl_adds, cb => system%cb, npl => pl%nbody) + ! Add the colliders%idx bodies to the subtraction list + ncolliders = colliders%ncoll + nfrag = frag%nbody + + param%maxid_collision = max(param%maxid_collision, maxval(system%pl%info(:)%collision_id)) + param%maxid_collision = param%maxid_collision + 1 + + ! Setup new bodies + allocate(plnew, mold=pl) + call plnew%setup(nfrag, param) + ibiggest = colliders%idx(maxloc(pl%Gmass(colliders%idx(:)), dim=1)) + ismallest = colliders%idx(minloc(pl%Gmass(colliders%idx(:)), dim=1)) + + ! Copy over identification, information, and physical properties of the new bodies from the fragment list + plnew%id(1:nfrag) = frag%id(1:nfrag) + plnew%xb(:, 1:nfrag) = frag%xb(:, 1:nfrag) + plnew%vb(:, 1:nfrag) = frag%vb(:, 1:nfrag) + call pl%vb2vh(cb) + call pl%xh2xb(cb) + do i = 1, nfrag + plnew%xh(:,i) = frag%xb(:, i) - cb%xb(:) + plnew%vh(:,i) = frag%vb(:, i) - cb%vb(:) + end do + plnew%mass(1:nfrag) = frag%mass(1:nfrag) + plnew%Gmass(1:nfrag) = param%GU * frag%mass(1:nfrag) + plnew%radius(1:nfrag) = frag%radius(1:nfrag) + plnew%density(1:nfrag) = frag%mass(1:nfrag) / frag%radius(1:nfrag) + call plnew%set_rhill(cb) + + select case(status) + case(DISRUPTION) + plnew%status(1:nfrag) = NEW_PARTICLE + do i = 1, nfrag + write(newname, FRAGFMT) frag%id(i) + call plnew%info(i)%set_value(origin_type="Disruption", origin_time=system%t, name=newname, & + origin_xh=plnew%xh(:,i), & + origin_vh=plnew%vh(:,i), collision_id=param%maxid_collision) + end do + do i = 1, ncolliders + if (colliders%idx(i) == ibiggest) then + iother = ismallest + else + iother = ibiggest + end if + call pl%info(colliders%idx(i))%set_value(status="Disruption", discard_time=system%t, & + discard_xh=pl%xh(:,i), discard_vh=pl%vh(:,i), discard_body_id=iother) + end do + case(SUPERCATASTROPHIC) + plnew%status(1:nfrag) = NEW_PARTICLE + do i = 1, nfrag + write(newname, FRAGFMT) frag%id(i) + call plnew%info(i)%set_value(origin_type="Supercatastrophic", origin_time=system%t, name=newname, & + origin_xh=plnew%xh(:,i), origin_vh=plnew%vh(:,i), & + collision_id=param%maxid_collision) + end do + do i = 1, ncolliders + if (colliders%idx(i) == ibiggest) then + iother = ismallest + else + iother = ibiggest + end if + call pl%info(colliders%idx(i))%set_value(status="Supercatastrophic", discard_time=system%t, & + discard_xh=pl%xh(:,i), discard_vh=pl%vh(:,i), & + discard_body_id=iother) + end do + case(HIT_AND_RUN_DISRUPT) + call plnew%info(1)%copy(pl%info(ibiggest)) + plnew%status(1) = OLD_PARTICLE + do i = 2, nfrag + write(newname, FRAGFMT) frag%id(i) + call plnew%info(i)%set_value(origin_type="Hit and run fragment", origin_time=system%t, name=newname, & + origin_xh=plnew%xh(:,i), origin_vh=plnew%vh(:,i), & + collision_id=param%maxid_collision) + end do + do i = 1, ncolliders + if (colliders%idx(i) == ibiggest) cycle + iother = ibiggest + call pl%info(colliders%idx(i))%set_value(status="Hit and run fragmention", discard_time=system%t, & + discard_xh=pl%xh(:,i), discard_vh=pl%vh(:,i), & + discard_body_id=iother) + end do + case(MERGED) + call plnew%info(1)%copy(pl%info(ibiggest)) + plnew%status(1) = OLD_PARTICLE + do i = 1, ncolliders + if (colliders%idx(i) == ibiggest) cycle + + iother = ibiggest + call pl%info(colliders%idx(i))%set_value(status="MERGED", discard_time=system%t, discard_xh=pl%xh(:,i), & + discard_vh=pl%vh(:,i), discard_body_id=iother) + end do + end select + + if (param%lrotation) then + plnew%Ip(:, 1:nfrag) = frag%Ip(:, 1:nfrag) + plnew%rot(:, 1:nfrag) = frag%rot(:, 1:nfrag) + end if + + ! if (param%ltides) then + ! plnew%Q = pl%Q(ibiggest) + ! plnew%k2 = pl%k2(ibiggest) + ! plnew%tlag = pl%tlag(ibiggest) + ! end if + + !Copy over or set integration parameters for new bodies + plnew%lcollision(1:nfrag) = .false. + plnew%ldiscard(1:nfrag) = .false. + plnew%levelg(1:nfrag) = pl%levelg(ibiggest) + plnew%levelm(1:nfrag) = pl%levelm(ibiggest) + + ! Log the properties of the new bodies + call fraggle_io_log_pl(plnew, param) + + ! Append the new merged body to the list + nstart = pl_adds%nbody + 1 + nend = pl_adds%nbody + nfrag + call pl_adds%append(plnew, lsource_mask=[(.true., i=1, nfrag)]) + ! Record how many bodies were added in this event + pl_adds%ncomp(nstart:nend) = plnew%nbody + + ! Add the discarded bodies to the discard list + pl%status(colliders%idx(:)) = MERGED + pl%ldiscard(colliders%idx(:)) = .true. + pl%lcollision(colliders%idx(:)) = .true. + lmask(:) = .false. + lmask(colliders%idx(:)) = .true. + + call plnew%setup(0, param) + deallocate(plnew) + + allocate(plsub, mold=pl) + call pl%spill(plsub, lmask, ldestructive=.false.) + + nstart = pl_discards%nbody + 1 + nend = pl_discards%nbody + ncolliders + call pl_discards%append(plsub, lsource_mask=[(.true., i = 1, ncolliders)]) + + ! Record how many bodies were subtracted in this event + pl_discards%ncomp(nstart:nend) = ncolliders + + call plsub%setup(0, param) + deallocate(plsub) + end associate + end select + end select + + return + end subroutine symba_collision_mergeaddsub + + + module subroutine symba_collision_resolve_fragmentations(self, system, param) + !! author: David A. Minton + !! + !! Process list of collisions, determine the collisional regime, and then create fragments. + !! + implicit none + ! Arguments + class(symba_plplenc), intent(inout) :: self !! SyMBA pl-pl encounter list + class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object + class(symba_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA additions + ! Internals + ! Internals + integer(I4B), dimension(2) :: idx_parent !! Index of the two bodies considered the "parents" of the collision + logical :: lgoodcollision + integer(I4B) :: i + type(fraggle_colliders) :: colliders !! Fraggle colliders object + type(fraggle_fragments) :: frag !! Fraggle fragmentation system object + + associate(plplcollision_list => self, ncollisions => self%nenc, idx1 => self%index1, idx2 => self%index2) + select type(pl => system%pl) + class is (symba_pl) + select type (cb => system%cb) + class is (symba_cb) + do i = 1, ncollisions + idx_parent(1) = pl%kin(idx1(i))%parent + idx_parent(2) = pl%kin(idx2(i))%parent + lgoodcollision = symba_collision_consolidate_colliders(pl, cb, param, idx_parent, colliders) + if ((.not. lgoodcollision) .or. any(pl%status(idx_parent(:)) /= COLLISION)) cycle + + call colliders%regime(frag, system, param) + + select case (frag%regime) + case (COLLRESOLVE_REGIME_DISRUPTION, COLLRESOLVE_REGIME_SUPERCATASTROPHIC) + plplcollision_list%status(i) = symba_collision_casedisruption(system, param, colliders, frag) + case (COLLRESOLVE_REGIME_HIT_AND_RUN) + plplcollision_list%status(i) = symba_collision_casehitandrun(system, param, colliders, frag) + case (COLLRESOLVE_REGIME_MERGE, COLLRESOLVE_REGIME_GRAZE_AND_MERGE) + plplcollision_list%status(i) = symba_collision_casemerge(system, param, colliders, frag) + case default + write(*,*) "Error in symba_collision, unrecognized collision regime" + call util_exit(FAILURE) + end select + end do + end select + end select + end associate + + return + end subroutine symba_collision_resolve_fragmentations + + + module subroutine symba_collision_resolve_mergers(self, system, param) + !! author: David A. Minton + !! + !! Process list of collisions and merge colliding bodies together. + !! + implicit none + ! Arguments + class(symba_plplenc), intent(inout) :: self !! SyMBA pl-pl encounter list + class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object + class(symba_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA additions + ! Internals + integer(I4B), dimension(2) :: idx_parent !! Index of the two bodies considered the "parents" of the collision + logical :: lgoodcollision + integer(I4B) :: i + type(fraggle_colliders) :: colliders !! Fraggle colliders object + type(fraggle_fragments) :: frag !! Fraggle fragmentation system object + + associate(plplcollision_list => self, ncollisions => self%nenc, idx1 => self%index1, idx2 => self%index2) + select type(pl => system%pl) + class is (symba_pl) + select type(cb => system%cb) + class is (symba_cb) + do i = 1, ncollisions + idx_parent(1) = pl%kin(idx1(i))%parent + idx_parent(2) = pl%kin(idx2(i))%parent + lgoodcollision = symba_collision_consolidate_colliders(pl, cb, param, idx_parent, colliders) + if (.not. lgoodcollision) cycle + if (any(pl%status(idx_parent(:)) /= COLLISION)) cycle ! One of these two bodies has already been resolved + + frag%regime = COLLRESOLVE_REGIME_MERGE + frag%mtot = sum(colliders%mass(:)) + frag%mass_dist(1) = frag%mtot + frag%mass_dist(2) = 0.0_DP + frag%mass_dist(3) = 0.0_DP + frag%xbcom(:) = (colliders%mass(1) * colliders%xb(:,1) + colliders%mass(2) * colliders%xb(:,2)) / frag%mtot + frag%vbcom(:) = (colliders%mass(1) * colliders%vb(:,1) + colliders%mass(2) * colliders%vb(:,2)) / frag%mtot + plplcollision_list%status(i) = symba_collision_casemerge(system, param, colliders, frag) + end do + end select + end select + end associate + + return + end subroutine symba_collision_resolve_mergers + + + module subroutine symba_collision_resolve_plplenc(self, system, param, t, dt, irec) + !! author: David A. Minton + !! + !! Process the pl-pl collision list, then modifiy the massive bodies based on the outcome of the collision + !! + implicit none + ! Arguments + class(symba_plplenc), intent(inout) :: self !! SyMBA pl-pl encounter list + class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA additions + real(DP), intent(in) :: t !! Current simulation time + real(DP), intent(in) :: dt !! Current simulation step size + integer(I4B), intent(in) :: irec !! Current recursion level + ! Internals + real(DP) :: Eorbit_before, Eorbit_after + logical :: lplpl_collision + character(len=STRMAX) :: timestr + class(symba_parameters), allocatable :: tmp_param + + associate(plplenc_list => self, plplcollision_list => system%plplcollision_list) + select type(pl => system%pl) + class is (symba_pl) + select type(param) + class is (symba_parameters) + if (plplcollision_list%nenc == 0) return ! No collisions to resolve + ! Make sure that the heliocentric and barycentric coordinates are consistent with each other + call pl%vb2vh(system%cb) + call pl%xh2xb(system%cb) + + ! Get the energy before the collision is resolved + if (param%lenergy) then + call system%get_energy_and_momentum(param) + Eorbit_before = system%te + end if + + do + write(timestr,*) t + call io_log_one_message(FRAGGLE_LOG_OUT, "") + call io_log_one_message(FRAGGLE_LOG_OUT, "***********************************************************" // & + "***********************************************************") + call io_log_one_message(FRAGGLE_LOG_OUT, "Collision between massive bodies detected at time t = " // & + trim(adjustl(timestr))) + call io_log_one_message(FRAGGLE_LOG_OUT, "***********************************************************" // & + "***********************************************************") + allocate(tmp_param, source=param) + if (param%lfragmentation) then + call plplcollision_list%resolve_fragmentations(system, param) + else + call plplcollision_list%resolve_mergers(system, param) + end if + + ! Destroy the collision list now that the collisions are resolved + call plplcollision_list%setup(0_I8B) + + if ((system%pl_adds%nbody == 0) .and. (system%pl_discards%nbody == 0)) exit + + ! Save the add/discard information to file + call system%write_discard(tmp_param) + + ! Rearrange the arrays: Remove discarded bodies, add any new bodies, resort, and recompute all indices and encounter lists + call pl%rearray(system, tmp_param) + + ! Destroy the add/discard list so that we don't append the same body multiple times if another collision is detected + call system%pl_discards%setup(0, param) + call system%pl_adds%setup(0, param) + deallocate(tmp_param) + + ! Check whether or not any of the particles that were just added are themselves in a collision state. This will generate a new plplcollision_list + lplpl_collision = plplenc_list%collision_check(system, param, t, dt, irec) + + if (.not.lplpl_collision) exit + end do + + if (param%lenergy) then + call system%get_energy_and_momentum(param) + Eorbit_after = system%te + system%Ecollisions = system%Ecollisions + (Eorbit_after - Eorbit_before) + end if + + end select + end select + end associate + + return + end subroutine symba_collision_resolve_plplenc + + + module subroutine symba_collision_resolve_pltpenc(self, system, param, t, dt, irec) + !! author: David A. Minton + !! + !! Process the pl-tp collision list, then modifiy the massive bodies based on the outcome of the collision + !! + implicit none + ! Arguments + class(symba_pltpenc), intent(inout) :: self !! SyMBA pl-pl encounter list + class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA additions + real(DP), intent(in) :: t !! Current simulation tim + real(DP), intent(in) :: dt !! Current simulation step size + integer(I4B), intent(in) :: irec !! Current recursion level + + ! Make sure coordinate systems are all synced up due to being inside the recursion at this point + call system%pl%vb2vh(system%cb) + call system%tp%vb2vh(system%cb%vb) + call system%pl%b2h(system%cb) + call system%tp%b2h(system%cb) + + ! Discard the collider + call system%tp%discard(system, param) + + return + end subroutine symba_collision_resolve_pltpenc + +end submodule s_symba_collision \ No newline at end of file diff --git a/src/symba/symba_collision.f90 b/src/symba/symba_collision.f90 index 1423adc7e..e839af1de 100644 --- a/src/symba/symba_collision.f90 +++ b/src/symba/symba_collision.f90 @@ -761,7 +761,7 @@ subroutine symba_collision_mergeaddsub(system, param, colliders, frag, status) plnew%status(1:nfrag) = NEW_PARTICLE do i = 1, nfrag write(newname, FRAGFMT) frag%id(i) - call plnew%info(i)%set_value(origin_type="Disruption", origin_time=param%t, name=newname, & + call plnew%info(i)%set_value(origin_type="Disruption", origin_time=system%t, name=newname, & origin_xh=plnew%xh(:,i), & origin_vh=plnew%vh(:,i), collision_id=param%maxid_collision) end do @@ -771,14 +771,14 @@ subroutine symba_collision_mergeaddsub(system, param, colliders, frag, status) else iother = ibiggest end if - call pl%info(colliders%idx(i))%set_value(status="Disruption", discard_time=param%t, & + call pl%info(colliders%idx(i))%set_value(status="Disruption", discard_time=system%t, & discard_xh=pl%xh(:,i), discard_vh=pl%vh(:,i), discard_body_id=iother) end do case(SUPERCATASTROPHIC) plnew%status(1:nfrag) = NEW_PARTICLE do i = 1, nfrag write(newname, FRAGFMT) frag%id(i) - call plnew%info(i)%set_value(origin_type="Supercatastrophic", origin_time=param%t, name=newname, & + call plnew%info(i)%set_value(origin_type="Supercatastrophic", origin_time=system%t, name=newname, & origin_xh=plnew%xh(:,i), origin_vh=plnew%vh(:,i), & collision_id=param%maxid_collision) end do @@ -788,7 +788,7 @@ subroutine symba_collision_mergeaddsub(system, param, colliders, frag, status) else iother = ibiggest end if - call pl%info(colliders%idx(i))%set_value(status="Supercatastrophic", discard_time=param%t, & + call pl%info(colliders%idx(i))%set_value(status="Supercatastrophic", discard_time=system%t, & discard_xh=pl%xh(:,i), discard_vh=pl%vh(:,i), & discard_body_id=iother) end do @@ -797,14 +797,14 @@ subroutine symba_collision_mergeaddsub(system, param, colliders, frag, status) plnew%status(1) = OLD_PARTICLE do i = 2, nfrag write(newname, FRAGFMT) frag%id(i) - call plnew%info(i)%set_value(origin_type="Hit and run fragment", origin_time=param%t, name=newname, & + call plnew%info(i)%set_value(origin_type="Hit and run fragment", origin_time=system%t, name=newname, & origin_xh=plnew%xh(:,i), origin_vh=plnew%vh(:,i), & collision_id=param%maxid_collision) end do do i = 1, ncolliders if (colliders%idx(i) == ibiggest) cycle iother = ibiggest - call pl%info(colliders%idx(i))%set_value(status="Hit and run fragmention", discard_time=param%t, & + call pl%info(colliders%idx(i))%set_value(status="Hit and run fragmention", discard_time=system%t, & discard_xh=pl%xh(:,i), discard_vh=pl%vh(:,i), & discard_body_id=iother) end do @@ -815,7 +815,7 @@ subroutine symba_collision_mergeaddsub(system, param, colliders, frag, status) if (colliders%idx(i) == ibiggest) cycle iother = ibiggest - call pl%info(colliders%idx(i))%set_value(status="MERGED", discard_time=param%t, discard_xh=pl%xh(:,i), & + call pl%info(colliders%idx(i))%set_value(status="MERGED", discard_time=system%t, discard_xh=pl%xh(:,i), & discard_vh=pl%vh(:,i), discard_body_id=iother) end do end select @@ -1019,7 +1019,6 @@ module subroutine symba_collision_resolve_plplenc(self, system, param, t, dt, ir call io_log_one_message(FRAGGLE_LOG_OUT, "***********************************************************" // & "***********************************************************") allocate(tmp_param, source=param) - tmp_param%t = t if (param%lfragmentation) then call plplcollision_list%resolve_fragmentations(system, param) else diff --git a/src/symba/symba_discard.f90 b/src/symba/symba_discard.f90 index f00b222cb..76bcbaf41 100644 --- a/src/symba/symba_discard.f90 +++ b/src/symba/symba_discard.f90 @@ -44,7 +44,7 @@ subroutine symba_discard_cb_pl(pl, system, param) pl%lcollision(i) = .false. pl%status(i) = DISCARDED_RMAX write(idstr, *) pl%id(i) - write(timestr, *) param%t + write(timestr, *) system%t write(message, *) trim(adjustl(pl%info(i)%name)) // " (" // trim(adjustl(idstr)) // ")" // & " too far from the central body at t = " // trim(adjustl(timestr)) call io_log_one_message(FRAGGLE_LOG_OUT, "") @@ -54,14 +54,14 @@ subroutine symba_discard_cb_pl(pl, system, param) call io_log_one_message(FRAGGLE_LOG_OUT, "***********************************************************" // & "***********************************************************") call io_log_one_message(FRAGGLE_LOG_OUT, "") - call pl%info(i)%set_value(status="DISCARDED_RMAX", discard_time=param%t, discard_xh=pl%xh(:,i), & + call pl%info(i)%set_value(status="DISCARDED_RMAX", discard_time=system%t, discard_xh=pl%xh(:,i), & discard_vh=pl%vh(:,i)) else if ((param%rmin >= 0.0_DP) .and. (rh2 < rmin2)) then pl%ldiscard(i) = .true. pl%lcollision(i) = .false. pl%status(i) = DISCARDED_RMIN write(idstr, *) pl%id(i) - write(timestr, *) param%t + write(timestr, *) system%t write(message, *) trim(adjustl(pl%info(i)%name)) // " (" // trim(adjustl(idstr)) // ")" // & " too close to the central body at t = " // trim(adjustl(timestr)) call io_log_one_message(FRAGGLE_LOG_OUT, "") @@ -71,7 +71,7 @@ subroutine symba_discard_cb_pl(pl, system, param) call io_log_one_message(FRAGGLE_LOG_OUT, "************************************************************" // & "************************************************************") call io_log_one_message(FRAGGLE_LOG_OUT, "") - call pl%info(i)%set_value(status="DISCARDED_RMIN", discard_time=param%t, discard_xh=pl%xh(:,i), & + call pl%info(i)%set_value(status="DISCARDED_RMIN", discard_time=system%t, discard_xh=pl%xh(:,i), & discard_vh=pl%vh(:,i), discard_body_id=cb%id) else if (param%rmaxu >= 0.0_DP) then rb2 = dot_product(pl%xb(:,i), pl%xb(:,i)) @@ -82,7 +82,7 @@ subroutine symba_discard_cb_pl(pl, system, param) pl%lcollision(i) = .false. pl%status(i) = DISCARDED_RMAXU write(idstr, *) pl%id(i) - write(timestr, *) param%t + write(timestr, *) system%t write(message, *) trim(adjustl(pl%info(i)%name)) // " (" // trim(adjustl(idstr)) // ")" // & " is unbound and too far from barycenter at t = " // trim(adjustl(timestr)) call io_log_one_message(FRAGGLE_LOG_OUT, "") @@ -92,7 +92,7 @@ subroutine symba_discard_cb_pl(pl, system, param) call io_log_one_message(FRAGGLE_LOG_OUT, "************************************************************" // & "************************************************************") call io_log_one_message(FRAGGLE_LOG_OUT, "") - call pl%info(i)%set_value(status="DISCARDED_RMAXU", discard_time=param%t, discard_xh=pl%xh(:,i), & + call pl%info(i)%set_value(status="DISCARDED_RMAXU", discard_time=system%t, discard_xh=pl%xh(:,i), & discard_vh=pl%vh(:,i)) end if end if @@ -325,11 +325,11 @@ subroutine symba_discard_peri_pl(pl, system, param) pl%ldiscard(i) = .true. pl%lcollision(i) = .false. pl%status(i) = DISCARDED_PERI - write(timestr, *) param%t + write(timestr, *) system%t write(idstr, *) pl%id(i) write(*, *) trim(adjustl(pl%info(i)%name)) // " (" // trim(adjustl(idstr)) // & ") perihelion distance too small at t = " // trim(adjustl(timestr)) - call pl%info(i)%set_value(status="DISCARDED_PERI", discard_time=param%t, & + call pl%info(i)%set_value(status="DISCARDED_PERI", discard_time=system%t, & discard_xh=pl%xh(:,i), discard_vh=pl%vh(:,i), discard_body_id=system%cb%id) end if end if diff --git a/src/symba/symba_util.f90 b/src/symba/symba_util.f90 index 0087c24e4..f0f300bc1 100644 --- a/src/symba/symba_util.f90 +++ b/src/symba/symba_util.f90 @@ -482,6 +482,7 @@ module subroutine symba_util_final_pl(self) return end subroutine symba_util_final_pl + module subroutine symba_util_final_system(self) !! author: David A. Minton !! @@ -495,6 +496,7 @@ module subroutine symba_util_final_system(self) return end subroutine symba_util_final_system + module subroutine symba_util_final_tp(self) !! author: David A. Minton !! diff --git a/src/util/util_copy.f90 b/src/util/util_copy.f90 index 2266396fb..51e210787 100644 --- a/src/util/util_copy.f90 +++ b/src/util/util_copy.f90 @@ -78,5 +78,18 @@ module subroutine util_copy_particle_info_arr(source, dest, idx) end subroutine util_copy_particle_info_arr + module subroutine util_copy_store_system(self, system) + !! author: David A. Minton + !! + !! Stores a snapshot of the nbody system so that later it can be retrieved for saving to file. + implicit none + class(storage_frame), intent(inout) :: self !! Swiftest storage frame object + class(swiftest_nbody_system), intent(in) :: system !! Swiftest n-body system object + + if (allocated(self%system)) deallocate(self%system) + allocate(self%system, source=system) + return + + end subroutine util_copy_store_system end submodule s_util_copy diff --git a/src/walltime/walltime.f90 b/src/walltime/walltime.f90 index 104b63d95..491a2e478 100644 --- a/src/walltime/walltime.f90 +++ b/src/walltime/walltime.f90 @@ -314,7 +314,6 @@ module subroutine walltime_interaction_time_this_loop(self, param, ninteractions end select self%is_on = .true. - write(tstr,*) param%t select case(self%stage) case(1) self%stage1_ninteractions = ninteractions