From 2e06f9391157b2d87cf2798b9cbe0027621898e2 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Sat, 31 Jul 2021 06:57:05 -0400 Subject: [PATCH] Started adding collisional family code for pl-pl mergers before they shut down the clusters for maintenence --- .../swiftest_vs_swifter.ipynb | 483 +++++++++++++++++- src/symba/symba_collision.f90 | 50 +- 2 files changed, 505 insertions(+), 28 deletions(-) diff --git a/examples/symba_swifter_comparison/1pl_1pl_encounter/swiftest_vs_swifter.ipynb b/examples/symba_swifter_comparison/1pl_1pl_encounter/swiftest_vs_swifter.ipynb index 57dd1934a..34c978f58 100644 --- a/examples/symba_swifter_comparison/1pl_1pl_encounter/swiftest_vs_swifter.ipynb +++ b/examples/symba_swifter_comparison/1pl_1pl_encounter/swiftest_vs_swifter.ipynb @@ -21,9 +21,9 @@ "output_type": "stream", "text": [ "Reading Swifter file param.swifter.in\n", - "Reading in time 6.845e-04\n", + "Reading in time 1.506e-01\n", "Creating Dataset\n", - "Successfully converted 2 output frames.\n", + "Successfully converted 221 output frames.\n", "Swifter simulation data stored as xarray DataSet .ds\n" ] } @@ -44,25 +44,9 @@ "text": [ "Reading Swiftest file param.swiftest.in\n", "Reading in time 1.506e-01\n", - "Creating Dataset\n" - ] - }, - { - "ename": "MergeError", - "evalue": "conflicting values for variable 'Mass' on objects to be combined. You can skip this check by specifying compat='override'.", - "output_type": "error", - "traceback": [ - "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", - "\u001b[0;31mMergeError\u001b[0m Traceback (most recent call last)", - "\u001b[0;32m\u001b[0m in \u001b[0;36m\u001b[0;34m\u001b[0m\n\u001b[1;32m 1\u001b[0m \u001b[0mswiftestsim\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mswiftest\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mSimulation\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mparam_file\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m\"param.swiftest.in\"\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 2\u001b[0;31m \u001b[0mswiftestsim\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mbin2xr\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m", - "\u001b[0;32m~/git/swiftest/python/swiftest/swiftest/simulation_class.py\u001b[0m in \u001b[0;36mbin2xr\u001b[0;34m(self)\u001b[0m\n\u001b[1;32m 135\u001b[0m \u001b[0;32mdef\u001b[0m \u001b[0mbin2xr\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 136\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcodename\u001b[0m \u001b[0;34m==\u001b[0m \u001b[0;34m\"Swiftest\"\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 137\u001b[0;31m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mds\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mio\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mswiftest2xr\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mparam\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 138\u001b[0m \u001b[0mprint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'Swiftest simulation data stored as xarray DataSet .ds'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 139\u001b[0m \u001b[0;32melif\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcodename\u001b[0m \u001b[0;34m==\u001b[0m \u001b[0;34m\"Swifter\"\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", - "\u001b[0;32m~/git/swiftest/python/swiftest/swiftest/io.py\u001b[0m in \u001b[0;36mswiftest2xr\u001b[0;34m(param)\u001b[0m\n\u001b[1;32m 632\u001b[0m \u001b[0mtpds\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mtpda\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mto_dataset\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mdim\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m'vec'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 633\u001b[0m \u001b[0mprint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'\\nCreating Dataset'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 634\u001b[0;31m \u001b[0mds\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mxr\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcombine_by_coords\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mcbds\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mplds\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mtpds\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 635\u001b[0m \u001b[0mprint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34mf\"Successfully converted {ds.sizes['time']} output frames.\"\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 636\u001b[0m \u001b[0;32mreturn\u001b[0m \u001b[0mds\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", - "\u001b[0;32m~/.conda/envs/cent7/2020.02-py37/swiftestOOF/lib/python3.7/site-packages/xarray/core/combine.py\u001b[0m in \u001b[0;36mcombine_by_coords\u001b[0;34m(datasets, compat, data_vars, coords, fill_value, join, combine_attrs)\u001b[0m\n\u001b[1;32m 816\u001b[0m \u001b[0mfill_value\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mfill_value\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 817\u001b[0m \u001b[0mjoin\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mjoin\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 818\u001b[0;31m \u001b[0mcombine_attrs\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mcombine_attrs\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 819\u001b[0m )\n", - "\u001b[0;32m~/.conda/envs/cent7/2020.02-py37/swiftestOOF/lib/python3.7/site-packages/xarray/core/merge.py\u001b[0m in \u001b[0;36mmerge\u001b[0;34m(objects, compat, join, fill_value, combine_attrs)\u001b[0m\n\u001b[1;32m 893\u001b[0m \u001b[0mjoin\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 894\u001b[0m \u001b[0mcombine_attrs\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mcombine_attrs\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 895\u001b[0;31m \u001b[0mfill_value\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mfill_value\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 896\u001b[0m )\n\u001b[1;32m 897\u001b[0m \u001b[0mmerged\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mDataset\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_construct_direct\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m**\u001b[0m\u001b[0mmerge_result\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_asdict\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", - "\u001b[0;32m~/.conda/envs/cent7/2020.02-py37/swiftestOOF/lib/python3.7/site-packages/xarray/core/merge.py\u001b[0m in \u001b[0;36mmerge_core\u001b[0;34m(objects, compat, join, combine_attrs, priority_arg, explicit_coords, indexes, fill_value)\u001b[0m\n\u001b[1;32m 625\u001b[0m \u001b[0mprioritized\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0m_get_priority_vars_and_indexes\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0maligned\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mpriority_arg\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcompat\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mcompat\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 626\u001b[0m variables, out_indexes = merge_collected(\n\u001b[0;32m--> 627\u001b[0;31m \u001b[0mcollected\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mprioritized\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcompat\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mcompat\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcombine_attrs\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mcombine_attrs\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 628\u001b[0m )\n\u001b[1;32m 629\u001b[0m \u001b[0massert_unique_multiindex_level_names\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mvariables\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", - "\u001b[0;32m~/.conda/envs/cent7/2020.02-py37/swiftestOOF/lib/python3.7/site-packages/xarray/core/merge.py\u001b[0m in \u001b[0;36mmerge_collected\u001b[0;34m(grouped, prioritized, compat, combine_attrs)\u001b[0m\n\u001b[1;32m 232\u001b[0m \u001b[0mvariables\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m[\u001b[0m\u001b[0mvariable\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0mvariable\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0m_\u001b[0m \u001b[0;32min\u001b[0m \u001b[0melements_list\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 233\u001b[0m \u001b[0;32mtry\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 234\u001b[0;31m \u001b[0mmerged_vars\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mname\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0munique_variable\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mname\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mvariables\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcompat\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 235\u001b[0m \u001b[0;32mexcept\u001b[0m \u001b[0mMergeError\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 236\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0mcompat\u001b[0m \u001b[0;34m!=\u001b[0m \u001b[0;34m\"minimal\"\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", - "\u001b[0;32m~/.conda/envs/cent7/2020.02-py37/swiftestOOF/lib/python3.7/site-packages/xarray/core/merge.py\u001b[0m in \u001b[0;36munique_variable\u001b[0;34m(name, variables, compat, equals)\u001b[0m\n\u001b[1;32m 140\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0mequals\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 141\u001b[0m raise MergeError(\n\u001b[0;32m--> 142\u001b[0;31m \u001b[0;34mf\"conflicting values for variable {name!r} on objects to be combined. \"\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 143\u001b[0m \u001b[0;34m\"You can skip this check by specifying compat='override'.\"\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 144\u001b[0m )\n", - "\u001b[0;31mMergeError\u001b[0m: conflicting values for variable 'Mass' on objects to be combined. You can skip this check by specifying compat='override'." + "Creating Dataset\n", + "Successfully converted 221 output frames.\n", + "Swiftest simulation data stored as xarray DataSet .ds\n" ] } ], @@ -73,7 +57,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 4, "metadata": {}, "outputs": [], "source": [ @@ -82,7 +66,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 5, "metadata": {}, "outputs": [], "source": [ @@ -91,18 +75,463 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 6, "metadata": {}, - "outputs": [], + "outputs": [ + { + "data": { + "text/plain": [ + "[,\n", + " ]" + ] + }, + "execution_count": 6, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZUAAAEGCAYAAACtqQjWAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAg50lEQVR4nO3df5QdZZ3n8fcnnV+ixPCj0Q4NpjXBSQfYEHshOgoow5rEWXrU1Ul0DaBOJkp2dpZ1d+J6dv2xB4czDjvKGslEiRJXyTAySvQEYnRUPGqARhgkQEwTGNOhJTEMCBMJSe53/6jqcHP7duf27ap7O9Wf1zn39K16nqfqWze5/e2qp+p5FBGYmZllYUKzAzAzs+JwUjEzs8w4qZiZWWacVMzMLDNOKmZmlpmJzQ6gmU499dSYOXNms8MwMzuu3Hvvvb+JiNZqZeM6qcycOZOenp5mh2FmdlyR9M9Dlfnyl5mZZcZJxczMMuOkYmZmmRnXfSrVHDx4kL6+Pp5//vlmh9IUU6dOpb29nUmTJjU7FDM7DjmpVOjr6+PEE09k5syZSGp2OA0VEezbt4++vj46OjqaHY6ZHYd8+avC888/zymnnDLuEgqAJE455ZRxe5ZmZqPnpFLFeEwoA8bzsZvZ6OWaVCQtlLRdUq+kVVXKJen6tPwBSfPLytZJ2iPpwYo2fyfp/vT1uKT70/UzJf2urGxNnsdmZlavR379W37S+5tmh5GL3JKKpBZgNbAI6ASWSuqsqLYImJ2+lgM3lJV9BVhYud2I+OOImBcR84BbgX8oK350oCwiVmR1LHl7wxveUHX9FVdcwTe+8Y0GR2NmeVuydivv/dJd/Pb5g80OJXN5nqmcD/RGxM6IeAHYAHRX1OkG1kdiKzBdUhtARNwJPDXUxpVcp3k3cHMu0TfQT3/602aHYGYN9LIpyT1Sf9/T1+RIspdnUjkd2FW23JeuG2mdobwJeDIidpSt65B0n6QfSXpTtUaSlkvqkdSzd+/eGneVr5e97GVAcvfVypUr6ezs5G1vext79uxpcmRmlofOtmkAfPVnj1MqFWv23TyTSrUe38pPr5Y6Q1nK0Wcp/cCZEXEecDXwdUnTBm08Ym1EdEVEV2tr1fHQmuab3/wm27dv5xe/+AVf/OIXfQZjVlADv+Qe37eflTf/nGf2F+cyWJ7PqfQBZ5QttwNP1FFnEEkTgXcArxtYFxEHgAPp+3slPQqcBRw3I0beeeedLF26lJaWFmbMmMFb3vKWZodkZjmICDrbpnHZvBl8ZvN2vv/w97jorFbOPv3l/N4rT+TMU07gpBMmM/2ESUyZ2NLscEckz6RyDzBbUgewG1gCvKeizkZgpaQNwAXAMxHRX8O2/wB4JCKOXJCU1Ao8FRGHJb2apPN/ZwbH0VC+pdes+EoBLRPEiotew0VntXLz3b/iR7/cy3cfenJQ3ZdMauHEqROZ1DKBSS1iYssEJk4QkycmP1smiIjk7Cci0p/p2dBRy0GpBKUIDpeCi1/bysfeVnnv1OjlllQi4pCklcBmoAVYFxHbJK1Iy9cAm4DFQC+wH7hyoL2km4GLgVMl9QEfj4gb0+IlDO6gvxD4lKRDwGFgRUQM2dE/Fl144YX87d/+LcuWLWPPnj384Ac/4D3vqczDZna8K0UwIf37cU7bND7VfTYA/3rgENuffJb+p5/n6d+9wNP7D/L0/hd47sAhDh4ODh4ucehw8MLhEocOlzh4OEkQEskLMfB3qSTEwPpkeYJggpJE9MqXvySXY8t1mJaI2ESSOMrXrSl7H8BVQ7RdOsx2r6iy7laSW4yPW29/+9v5x3/8R8455xzOOussLrroomaHZGY5KAVQ5arES6dMZP6ZJ8GZjY8pKx77awx47rnngOQvic9//vNNjsbM8hZlZypF42FazMwaLCK5DFVETipmZg1W8pmKmZllpRRR2Ds9nVTMzBqsFNWf/C4CJxUzs0Zzn4qZmWWlFMGEgv72LehhHd927drFm9/8ZubMmcPcuXP53Oc+N6hORPBnf/ZnzJo1i3PPPZef//znTYjUzOqRdNQX80zFz6mMQRMnTuS6665j/vz5PPvss7zuda/j0ksvpbPzxSEVbr/9dnbs2MGOHTu46667+NCHPsRdd93VxKjNrFYFG5j4KD5TGYPa2tqYPz+ZBPPEE09kzpw57N69+6g6t912G8uWLUMSCxYs4Omnn6a/v5Zh08ys2cJnKuPTJ7+9jYee+G2m2+ycMY2P//u5Ndd//PHHue+++7jggguOWr97927OOOPFAZ7b29vZvXs3bW1tmcVqZvkI8HMq1njPPfcc73znO/nsZz/LtGlHTw2TDJt2tKLe925WNO5TGadGckaRtYMHD/LOd76T9773vbzjHe8YVN7e3s6uXS9OmtnX18eMGTMaGaKZ1alUKu4fgT5TGYMigg984APMmTOHq6++umqdyy67jPXr1xMRbN26lZe//OW+9GV2nEieqG92FPnwmcoY9JOf/ISvfvWrnHPOOcybNw+AT3/60/zqV78CYMWKFSxevJhNmzYxa9YsTjjhBL785S83MWIzG4lkQMlmR5EPJ5Ux6I1vfGPVPpNykli9enWDIjKzLAXF7VPx5S8zswYreZgWMzPLSpH7VJxUzMwaLMJ3f9VF0kJJ2yX1SlpVpVySrk/LH5A0v6xsnaQ9kh6saPMJSbsl3Z++FpeVfTTd1nZJb83z2MzM6uXphOsgqQVYDSwCOoGlkjorqi0CZqev5cANZWVfARYOsfm/iYh56WtTur9OYAkwN233hTQGM7MxxX0q9Tkf6I2InRHxArAB6K6o0w2sj8RWYLqkNoCIuBN4agT76wY2RMSBiHgM6E1jMDMbU9ynUp/TgV1ly33pupHWqWZlerlsnaSTRrItScsl9Ujq2bt3bw27arz3v//9nHbaaZx99tlH1j311FNceumlzJ49m0svvZR/+Zd/OVL2l3/5l8yaNYvXvva1bN68ueo2h2tvZo0VPlOpS7VPrPLhi1rqVLoBeA0wD+gHrhvJtiJibUR0RURXa2vrMXbVHFdccQV33HHHUeuuvfZaLrnkEnbs2MEll1zCtddeC8BDDz3Ehg0b2LZtG3fccQcf/vCHOXz48KBtDtXezBqvFOHphOvQB5xRttwOPFFHnaNExJMRcTgiSsAXefES14i3NVZdeOGFnHzyyUetu+2227j88ssBuPzyy/nWt751ZP2SJUuYMmUKHR0dzJo1i7vvvnvQNodqb2aNV+QzlTyfqL8HmC2pA9hN0on+noo6G0kuZW0ALgCeiYhhJwWR1FZW5+3AwN1hG4GvS/o/wAySzv/Bv11H4vZV8OtfjGoTg7zyHFg08rOEJ5988sjYXm1tbezZswdIhsBfsGDBkXoDQ+DX2t7MGq/I0wnnllQi4pCklcBmoAVYFxHbJK1Iy9cAm4DFJJ3q+4ErB9pLuhm4GDhVUh/w8Yi4EfgrSfNILm09Dvxpur1tkm4BHgIOAVdFxODrQAXjIfDNjj+lAj+nkuvYX+ntvpsq1q0pex/AVUO0XTrE+vcNs79rgGvqCraaOs4o8vKKV7yC/v5+2tra6O/v57TTTgNqHwJ/qPZm1nh+TsWa7rLLLuOmm24C4KabbqK7u/vI+g0bNnDgwAEee+wxduzYwfnnD76Teqj2ZtZ4SUd9MbOKk8oYtHTpUl7/+tezfft22tvbufHGG1m1ahVbtmxh9uzZbNmyhVWrkgEK5s6dy7vf/W46OztZuHAhq1evpqUleebzgx/8ID09PQBDtjezxivydMI61hDrRdbV1RUDv3QHPPzww8yZM6dJEY0N/gzM8nXuJzbzjvntfOKy5s0uOxqS7o2IrmplPlMxM2uwIt9S7KRiZtZgHqZlnBnPlwTH87GbNUqpwNMJO6lUmDp1Kvv27RuXv1wjgn379jF16tRmh2JWaEWeTthz1Fdob2+nr6+PsTrYZN6mTp1Ke3t7s8MwKzQ//DiOTJo0iY6OjmaHYWYF5ocfzcwsM8mZSrOjyIeTiplZg5WiuH0qTipmZg0WBe5TcVIxM2uggTtL3adiZmajVkqfVvDlLzMzG7VSeqZSzJTipGJm1lADz1VPKOj1LycVM7MGOnKmUsyc4qRiZtZI4T4VMzPLSsl3f9VP0kJJ2yX1Sho01aAS16flD0iaX1a2TtIeSQ9WtPmMpEfS+t+UND1dP1PS7yTdn77W5HlsZmb1eLGjvphZJbekIqkFWA0sAjqBpZI6K6otAmanr+XADWVlXwEWVtn0FuDsiDgX+CXw0bKyRyNiXvpakcmBmJllaGD884Je/cr1TOV8oDcidkbEC8AGoLuiTjewPhJbgemS2gAi4k7gqcqNRsR3I+JQurgV8JC6ZnbciFLy030qI3c6sKtsuS9dN9I6w3k/cHvZcoek+yT9SNKbqjWQtFxSj6Se8Tq8vZk1j/tU6lftI6uc+aqWOtU3Ln0MOAR8LV3VD5wZEecBVwNflzRt0MYj1kZEV0R0tba21rIrM7PMHEkqBc0qeSaVPuCMsuV24Ik66gwi6XLgD4H3RjqQTkQciIh96ft7gUeBs+qO3swsBwPDtBQzpeSbVO4BZkvqkDQZWAJsrKizEViW3gW2AHgmIvqH26ikhcBfAJdFxP6y9a3pzQFIejVJ5//O7A7HzGz0goGHH4uZVnKb+TEiDklaCWwGWoB1EbFN0oq0fA2wCVgM9AL7gSsH2ku6GbgYOFVSH/DxiLgR+DwwBdiS/qNsTe/0uhD4lKRDwGFgRUQM6ug3M2umoj/8mOt0whGxiSRxlK9bU/Y+gKuGaLt0iPWzhlh/K3Br3cGamTWAO+rNzCwzR/pUnFTMzGy0SqVi96k4qZiZNUFR+1ScVMzMGsh9KmZmlhlPJ2xmZpnxJF1mZpaZCHfUm5lZRl58+LG5ceTFScXMrIHcp2JmZpnx3V9mZpaZgaRS1HGKnVTMzBrIfSpmZpaZoo9S7KRiZtZAL8782ORAclLQwzIzG5tKfk7FzMyy4umEzcwsQwO3FBczrTipmJk1kB9+NDOzzAxM0uVbiusgaaGk7ZJ6Ja2qUi5J16flD0iaX1a2TtIeSQ9WtDlZ0hZJO9KfJ5WVfTTd1nZJb83z2MzM6lEq9rOP+SUVSS3AamAR0AksldRZUW0RMDt9LQduKCv7CrCwyqZXAd+PiNnA99Nl0m0vAeam7b6QxmBmNmZEuE+lXucDvRGxMyJeADYA3RV1uoH1kdgKTJfUBhARdwJPVdluN3BT+v4m4I/K1m+IiAMR8RjQm8ZgZjZmDJyoOKmM3OnArrLlvnTdSOtUekVE9AOkP08bybYkLZfUI6ln7969xzwIM7MseUDJ+lX7yKKOOlnuj4hYGxFdEdHV2tpa567MzOpz5DkVn6mMWB9wRtlyO/BEHXUqPTlwiSz9uWcU2zIzayhPJ1y/e4DZkjokTSbpRN9YUWcjsCy9C2wB8MzApa1hbAQuT99fDtxWtn6JpCmSOkg6/+/O4kDMzLJS9I76iXltOCIOSVoJbAZagHURsU3SirR8DbAJWEzSqb4fuHKgvaSbgYuBUyX1AR+PiBuBa4FbJH0A+BXwrnR72yTdAjwEHAKuiojDeR2fmVk9ij70fW5JBSAiNpEkjvJ1a8reB3DVEG2XDrF+H3DJEGXXANfUG6+ZWd78RL2ZmWXGfSpmZpaZgT4VFfSReicVM7MGOnL5q6C/fQt6WGZmY9O4n05Y0mlV1r02n3DMzIrNT9TDjyW9e2BB0n8FvplfSGZmxVX06YRruaX4YmCtpHcBrwAexgM1mpnVJcb7dMLpE+53AK8HZpKMKvxcznGZmRVSabw/US9pC9APnE0yntY6SXdGxEfyDs7MrGjGfUc9cDvwPyLi6Yh4EHgD8Ey+YZmZFZMffoQTgc2SfizpKuCUiPjfOcdlZlZIR85UCnr7Vy19Kp+MiLkkY3TNAH4k6Xu5R2ZmVkBHzlSaHEdeRvLw4x7g18A+Xpxt0czMRmDcTycs6UOSfgh8HzgV+JOIODfvwMzMiqjoDz/W8pzKq4A/j4j7c47FzKzwij6d8DGTSkSsakQgZmbjQfjuLzMzy0qpVOyHH51UzMwa6MWO+qaGkRsnFTOzBip6n0quSUXSQknbJfVKGtQ3o8T1afkDkuYfq62kv5N0f/p6XNL96fqZkn5XVrYmz2MzM6tH+O6v+khqAVYDlwJ9wD2SNkbEQ2XVFgGz09cFwA3ABcO1jYg/LtvHdRw9ZMyjETEvr2MyMxutog99n+eZyvlAb0TsjIgXgA1Ad0WdbpJRjyMitgLTJbXV0lbJv8i7gZtzPAYzs0wdmU64mDkl16RyOrCrbLkvXVdLnVravgl4MiJ2lK3rkHSfpB9JelO1oCQtl9QjqWfv3r21H42ZWQY8SnH9qn1iUWOdWtou5eizlH7gzIg4D7ga+LqkaYM2ErE2Iroioqu1tXXI4M3M8lD0UYpz61MhObs4o2y5HXiixjqTh2sraSLwDuB1A+si4gBwIH1/r6RHgbOAntEeiJlZVqLgk3TleaZyDzBbUoekycASYGNFnY3AsvQusAXAM+lMk8dq+wfAIxHRN7BCUmvawY+kV5N0/u/M6+DMzOpRKvh0wrmdqUTEIUkrgc1AC7AuIrZJWpGWrwE2AYuBXmA/cOVwbcs2v4TBHfQXAp+SdAg4DKyIiKfyOj4zs3qM++mERyMiNpEkjvJ1a8reB8k8LTW1LSu7osq6W4FbRxGumVnu4sjDj82NIy9+ot7MrIEiAsnPqZiZWQZKUdxLX+CkYmbWUKWIwnbSg5OKmVlD+UzFzMwyE0RhO+nBScXMrKHCZypmZpaVUslnKmZmlhH3qZiZWWZK4TMVMzPLkM9UzMwsE6WIwk7QBU4qZmYNlVz+Km5WcVIxM2ugpKO+2VHkx0nFzKyBIoo7mCQ4qZiZNVS4T8XMzLKSdNQXN6s4qZiZNVApijuVMDipmJk1lO/+MjOz7ARMKPBv3lwPTdJCSdsl9UpaVaVckq5Pyx+QNP9YbSV9QtJuSfenr8VlZR9N62+X9NY8j83MrB5F71OZmNeGJbUAq4FLgT7gHkkbI+KhsmqLgNnp6wLgBuCCGtr+TUT8dcX+OoElwFxgBvA9SWdFxOG8jtHMbKQ8oGT9zgd6I2JnRLwAbAC6K+p0A+sjsRWYLqmtxraVuoENEXEgIh4DetPtmJmNGZ5OuH6nA7vKlvvSdbXUOVbblenlsnWSThrB/pC0XFKPpJ69e/eO5HjMzEYtefix2VHkJ8+kUu1jixrrDNf2BuA1wDygH7huBPsjItZGRFdEdLW2tlZpYmaWn8B9KvXqA84oW24HnqixzuSh2kbEkwMrJX0R+M4I9mdm1lSlkvtU6nUPMFtSh6TJJJ3oGyvqbASWpXeBLQCeiYj+4dqmfS4D3g48WLatJZKmSOog6fy/O6+DMzOrR9En6crtTCUiDklaCWwGWoB1EbFN0oq0fA2wCVhM0qm+H7hyuLbppv9K0jySS1uPA3+attkm6RbgIeAQcJXv/DKzsaZU8AEl87z8RURsIkkc5evWlL0P4Kpa26br3zfM/q4Brqk3XjOzvHlASTMzy0zgPhUzM8uIpxM2M7PMlIJCP6jipGJm1kDuUzEzs8wUfUBJJxUzswaKwGcqZmaWDU/SZWZmmfF0wmZmlplwn4qZmWUlPJ2wmZllxXd/mZlZZoo+oKSTiplZA4WnEzYzs6yU/JyKmZllpejTCTupmJk1UKnkPhUzM8uIh743M7PMRBR65HsnFTOzRvJzKqMgaaGk7ZJ6Ja2qUi5J16flD0iaf6y2kj4j6ZG0/jclTU/Xz5T0O0n3p681eR6bmVk9PJ1wnSS1AKuBRUAnsFRSZ0W1RcDs9LUcuKGGtluAsyPiXOCXwEfLtvdoRMxLXyvyOTIzs/oloxQ3O4r85Hmmcj7QGxE7I+IFYAPQXVGnG1gfia3AdEltw7WNiO9GxKG0/VagPcdjMDPLVPiJ+rqdDuwqW+5L19VSp5a2AO8Hbi9b7pB0n6QfSXpTtaAkLZfUI6ln7969tR2JmVlGfPdX/ap9bFFjnWO2lfQx4BDwtXRVP3BmRJwHXA18XdK0QRuJWBsRXRHR1draeoxDMDPLVtE76ifmuO0+4Iyy5XbgiRrrTB6uraTLgT8ELomIAIiIA8CB9P29kh4FzgJ6sjgYM7Ms+Jbi+t0DzJbUIWkysATYWFFnI7AsvQtsAfBMRPQP11bSQuAvgMsiYv/AhiS1ph38SHo1Sef/zhyPz8xsxJI56oubVXI7U4mIQ5JWApuBFmBdRGyTtCItXwNsAhYDvcB+4Mrh2qab/jwwBdiSdnZtTe/0uhD4lKRDwGFgRUQ8ldfxmZnVo1TwUYrzvPxFRGwiSRzl69aUvQ/gqlrbputnDVH/VuDW0cRrZpa3ovep+Il6M7MGKnk6YTMzy4qfUzEzs8yEn1MxM7OsJB31xc0qTipmZg3k6YTNzCwzyYCSxc0qTipmZo1U8IcfnVTMzBrIA0qamVlmSh77y8zMsuIn6s3MLDOBH340M7OM+OFHMzPLTMl3f5mZWVaS51SaHUV+nFTMzBrIA0qamVkm0tnP3adiZmajV0pyivtUzMxs9ErpmUpxU4qTiplZwwwklQkFvv6Va1KRtFDSdkm9klZVKZek69PyByTNP1ZbSSdL2iJpR/rzpLKyj6b1t0t6a57HZmY2UmlO8d1f9ZDUAqwGFgGdwFJJnRXVFgGz09dy4IYa2q4Cvh8Rs4Hvp8uk5UuAucBC4AvpdszMxoQYB30qE3Pc9vlAb0TsBJC0AegGHiqr0w2sj+SWiK2SpktqA2YO07YbuDhtfxPwQ+Av0vUbIuIA8Jik3jSGn2V9YI/8+rc8tO7DvObwY1lv2swKLIANkw/z0h3z4KI1zQ4nF3kmldOBXWXLfcAFNdQ5/RhtXxER/QAR0S/ptLJtba2yraNIWk5yVsSZZ545gsN50dSJLUw/YRIvOeATITMbmZdObqH15BOaHUZu8kwq1c7vosY6tbStZ39ExFpgLUBXV9extlnVzFNfysw//3I9Tc3MCi3Pjvo+4Iyy5XbgiRrrDNf2yfQSGenPPSPYn5mZ5SjPpHIPMFtSh6TJJJ3oGyvqbASWpXeBLQCeSS9tDdd2I3B5+v5y4Lay9UskTZHUQdL5f3deB2dmZoPldvkrIg5JWglsBlqAdRGxTdKKtHwNsAlYDPQC+4Erh2ubbvpa4BZJHwB+BbwrbbNN0i0knfmHgKsi4nBex2dmZoNpYCya8airqyt6enqaHYaZ2XFF0r0R0VWtzE/Um5lZZpxUzMwsM04qZmaWGScVMzPLzLjuqJe0F/jnUWziVOA3GYWTB8c3emM9xrEeH4z9GMd6fDD2YnxVRLRWKxjXSWW0JPUMdQfEWOD4Rm+sxzjW44OxH+NYjw+OjxgH+PKXmZllxknFzMwy46QyOmubHcAxOL7RG+sxjvX4YOzHONbjg+MjRsB9KmZmliGfqZiZWWacVMzMLDNOKlVIWihpu6ReSauqlEvS9Wn5A5Lm19q2mfFJOkPSDyQ9LGmbpP+cR3yjibGsvEXSfZK+M9biS6e9/oakR9LP8vVjMMb/kv4bPyjpZklTmxDf70n6maQDkj4ykrbNjrFR35XRfIZpea7fk7pEhF9lL5Kh9h8FXg1MBv4J6Kyosxi4nWS2yQXAXbW2bXJ8bcD89P2JwC+zjm+0MZaVXw18HfjOWIsPuAn4YPp+MjB9LMVIMo32Y8BL0uVbgCuaEN9pwL8FrgE+MpK2YyDG3L8ro4mvEd+Tel8+UxnsfKA3InZGxAvABqC7ok43sD4SW4HpSmahrKVt0+KLiP6I+DlARDwLPEzyCyhro/kMkdQOvA34Ug6xjSo+SdOAC4EbASLihYh4eizFmJZNBF4iaSJwAtnPgnrM+CJiT0TcAxwcadtmx9ig78poPsNGfE/q4qQy2OnArrLlPgb/ZxqqTi1tmxnfEZJmAucBd2UcX037P0adzwL/HSjlENto43s1sBf4cnrZ4UuSXjqWYoyI3cBfk0xi108yo+p3mxBfHm1HIpP95PhdGW18nyXf70ldnFQGU5V1lfddD1WnlrajNZr4kkLpZcCtwJ9HxG8zjK2m/Q9XR9IfAnsi4t7swxp+3zXWmQjMB26IiPOAfwXy6BMYzWd4EslfvB3ADOClkv5jE+LLo+1IjHo/OX9X6o6vQd+TujipDNYHnFG23M7gSwdD1amlbTPjQ9Ikki/J1yLiHzKOLYsYfx+4TNLjJJcD3iLp/42h+PqAvogY+Kv1GyRJJmujifEPgMciYm9EHAT+AXhDE+LLo+1IjGo/DfiujCa+RnxP6tPsTp2x9iL5S3QnyV95A51ncyvqvI2jO0jvrrVtk+MTsB747Fj9DCvqXEw+HfWjig/4MfDa9P0ngM+MpRiBC4BtJH0pIrmx4D81Or6yup/g6E7w3L8nGcSY+3dlNPFVlOXyPan7uJodwFh8kdxV80uSOzM+lq5bAaxI3wtYnZb/Augaru1YiQ94I8np9QPA/elr8ViKsWIbuX1ZRvlvPA/oST/HbwEnjcEYPwk8AjwIfBWY0oT4Xkny1/hvgafT99Ma9T0ZTYyN+q6M5jNsxPeknpeHaTEzs8y4T8XMzDLjpGJmZplxUjEzs8w4qZiZWWacVMzMLDNOKmYZSUcv/nDZ8gxJ38hpX38k6X8do85fS3pLHvs3G4pvKTbLSDpG1Hci4uwG7OunwGUR8Zth6rwK+GJE/Lu84zEb4DMVs+xcC7xG0v2SPiNppqQHASRdIelbkr4t6TFJKyVdnQ5KuVXSyWm910i6Q9K9kn4s6fcqdyLpLOBARPxG0onp9ialZdMkPS5pUkT8M3CKpFc28DOwcc5JxSw7q4BHI2JeRPy3KuVnA+8hGfL8GmB/JINS/gxYltZZSzKkyuuAjwBfqLKd3wfKh2X/IcmQLQBLgFsjGfOLtN7vj/K4zGo2sdkBmI0jP0iTwLOSngG+na7/BXBuOiLuG4C/l44MYDulynbaSIbfH/AlkiHQvwVcCfxJWdkekpGKzRrCScWscQ6UvS+VLZdIvosTgKcjYt4xtvM74OUDCxHxk/RS20VAS0Q8WFZ3alrfrCF8+cssO8+STD1bl0jm63hM0rvgyBz0/6ZK1YeBWRXr1gM3A1+uWH8WyaCSZg3hpGKWkYjYB/xE0oOSPlPnZt4LfEDSP5EMX19tmt07gfNUdo0M+BpwEkliAY7MBzKLZERls4bwLcVmxyFJnwO+HRHfS5f/A9AdEe8rq/N2YH5E/M8mhWnjkPtUzI5PnyaZjAtJ/xdYRDI3R7mJwHUNjsvGOZ+pmJlZZtynYmZmmXFSMTOzzDipmJlZZpxUzMwsM04qZmaWmf8P4LgpItRxF68AAAAASUVORK5CYII=\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], "source": [ "swiftdiff['vx'].plot.line(x=\"time (y)\")" ] }, { "cell_type": "code", - "execution_count": null, + "execution_count": 7, "metadata": {}, - "outputs": [], + "outputs": [ + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "
<xarray.DataArray 'vx' (time (y): 221)>\n",
+       "array([ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,\n",
+       "        0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,\n",
+       "        0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,\n",
+       "        0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,\n",
+       "        0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,\n",
+       "        0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,\n",
+       "        0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,\n",
+       "        0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,\n",
+       "        0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,\n",
+       "        0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,\n",
+       "        0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,\n",
+       "        0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,\n",
+       "        0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,\n",
+       "        0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,\n",
+       "        0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,\n",
+       "        0.,  0., nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan,\n",
+       "       nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan])\n",
+       "Coordinates:\n",
+       "    id        float64 100.0\n",
+       "  * time (y)  (time (y)) float64 0.0 0.0006845 0.001369 ... 0.1492 0.1499 0.1506
" + ], + "text/plain": [ + "\n", + "array([ 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,\n", + " 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,\n", + " 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,\n", + " 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,\n", + " 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,\n", + " 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,\n", + " 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,\n", + " 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,\n", + " 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,\n", + " 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,\n", + " 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,\n", + " 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,\n", + " 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,\n", + " 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,\n", + " 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,\n", + " 0., 0., nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan,\n", + " nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan])\n", + "Coordinates:\n", + " id float64 100.0\n", + " * time (y) (time (y)) float64 0.0 0.0006845 0.001369 ... 0.1492 0.1499 0.1506" + ] + }, + "execution_count": 7, + "metadata": {}, + "output_type": "execute_result" + } + ], "source": [ "swiftdiff['vx'].sel(id=100)" ] diff --git a/src/symba/symba_collision.f90 b/src/symba/symba_collision.f90 index 2eaa521ee..bdf83a595 100644 --- a/src/symba/symba_collision.f90 +++ b/src/symba/symba_collision.f90 @@ -19,6 +19,54 @@ module subroutine symba_collision_check_plplenc(self, system, param, t, dt, irec real(DP), intent(in) :: t !! current time real(DP), intent(in) :: dt !! step size integer(I4B), intent(in) :: irec !! Current recursion level + ! Internals + logical, dimension(:), allocatable :: lcollision, lmask + real(DP), dimension(NDIM) :: xr, vr + integer(I4B) :: k + real(DP) :: rlim, mtot + + if (self%nenc == 0) return + + select type(pl => system%pl) + class is (symba_pl) + select type(tp => system%tp) + class is (symba_tp) + associate(plplenc_list => self, nplplenc => self%nenc, ind1 => self%index1, ind2 => self%index2) + allocate(lmask(nplplenc)) + lmask(:) = ((plplenc_list%status(1:nplplenc) == ACTIVE) & + .and. (pl%levelg(ind1(1:nplplenc)) >= irec) & + .and. (pl%levelg(ind2(1:nplplenc)) >= irec)) + if (.not.any(lmask(:))) return + + allocate(lcollision(nplplenc)) + lcollision(:) = .false. + + do concurrent(k = 1:nplplenc, lmask(k)) + xr(:) = pl%xh(:, ind1(k)) - pl%xh(:, ind2(k)) + vr(:) = pl%vb(:, ind1(k)) - pl%vb(:, ind2(k)) + rlim = pl%radius(ind1(k)) + pl%radius(ind2(k)) + mtot = pl%Gmass(ind1(k)) + pl%Gmass(ind2(k)) + lcollision(k) = symba_collision_check_one(xr(1), xr(2), xr(3), vr(1), vr(2), vr(3), mtot, rlim, dt, plplenc_list%lvdotr(k)) + end do + + if (any(lcollision(:))) then + do k = 1, nplplenc + if (plplenc_list%status(k) /= COLLISION) cycle + plplenc_list%status(k) = COLLISION + plplenc_list%xh1(:,k) = pl%xh(:,ind1(k)) + plplenc_list%vb1(:,k) = pl%vb(:,ind1(k)) + plplenc_list%xh2(:,k) = pl%xh(:,ind2(k)) + plplenc_list%vb2(:,k) = pl%vb(:,ind2(k)) + if (pl%lcollision(ind1(k)) .or. pl%lcollision(ind2(k))) call plplenc_list%make_family(system) + pl%lcollision(ind1(k)) = .true. + pl%lcollision(ind2(k)) = .true. + end do + end if + end associate + end select + end select + + return return end subroutine symba_collision_check_plplenc @@ -51,7 +99,7 @@ module subroutine symba_collision_check_pltpenc(self, system, param, t, dt, irec class is (symba_pl) select type(tp => system%tp) class is (symba_tp) - associate(pltpenc_list => self, npltpenc => self%nenc, plind => self%index1, tpind => self%index2 ) + associate(pltpenc_list => self, npltpenc => self%nenc, plind => self%index1, tpind => self%index2) allocate(lmask(npltpenc)) lmask(:) = ((pltpenc_list%status(1:npltpenc) == ACTIVE) & .and. (pl%levelg(plind(1:npltpenc)) >= irec) &