diff --git a/examples/rmvs_swifter_comparison/1pl_1tp_encounter/init_cond.py b/examples/rmvs_swifter_comparison/1pl_1tp_encounter/init_cond.py index 34072613e..ca1be07bb 100755 --- a/examples/rmvs_swifter_comparison/1pl_1tp_encounter/init_cond.py +++ b/examples/rmvs_swifter_comparison/1pl_1tp_encounter/init_cond.py @@ -85,7 +85,7 @@ print(f'TP_IN {swifter_tp}') print(f'IN_TYPE ASCII') print(f'ISTEP_OUT {iout:d}') -print(f'ISTEP_DUMP {iout:d}') +print(f'ISTEP_DUMP {10*iout:d}') print(f'BIN_OUT {swifter_bin}') print(f'OUT_TYPE REAL8') print(f'OUT_FORM XV') diff --git a/examples/rmvs_swifter_comparison/1pl_1tp_encounter/param.swiftest.in b/examples/rmvs_swifter_comparison/1pl_1tp_encounter/param.swiftest.in index d81a1d316..5dd9ef3d8 100644 --- a/examples/rmvs_swifter_comparison/1pl_1tp_encounter/param.swiftest.in +++ b/examples/rmvs_swifter_comparison/1pl_1tp_encounter/param.swiftest.in @@ -7,7 +7,7 @@ PL_IN pl.swiftest.in TP_IN tp.swiftest.in IN_TYPE REAL8 ISTEP_OUT 1 -ISTEP_DUMP 1 +ISTEP_DUMP 10 BIN_OUT bin.swiftest.nc OUT_TYPE NETCDF_DOUBLE OUT_FORM XVEL diff --git a/examples/rmvs_swifter_comparison/1pl_1tp_encounter/swiftest_vs_swifter.ipynb b/examples/rmvs_swifter_comparison/1pl_1tp_encounter/swiftest_vs_swifter.ipynb index aa7e8e5e3..43c398d4a 100644 --- a/examples/rmvs_swifter_comparison/1pl_1tp_encounter/swiftest_vs_swifter.ipynb +++ b/examples/rmvs_swifter_comparison/1pl_1tp_encounter/swiftest_vs_swifter.ipynb @@ -43,8 +43,8 @@ "output_type": "stream", "text": [ "Reading Swiftest file param.swiftest.in\n", - "Reading in time 1.506e-01\n", - "Creating Dataset\n", + "\n", + "Creating Dataset from NetCDF file\n", "Successfully converted 221 output frames.\n", "Swiftest simulation data stored as xarray DataSet .ds\n" ] @@ -75,23 +75,23 @@ }, { "cell_type": "code", - "execution_count": 6, + "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "text/plain": [ - "[,\n", - " ]" + "[,\n", + " ]" ] }, - "execution_count": 6, + "execution_count": 7, "metadata": {}, "output_type": "execute_result" }, { "data": { - "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYYAAAERCAYAAAB/4wAeAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAZLklEQVR4nO3de5ScdZng8e/T3SERCehykcQGEw1gEmAQekHQ5Q4jOAMDDBxYVmTBYfG6s6w7yw5nR1lXxUVGnJWZPUFhgUWyCyNkRS6KkIUDghPuAWSCEqFDlBBFyCC35Nk/qrqp6vSluuvyVld9P+f06a56bw9F6nnqeX/v763ITCRJGtJTdACSpPZiYZAkVbEwSJKqWBgkSVUsDJKkKhYGSVKVaVcYIuKyiHg+IlY2aH+3RMSLEXHjiOfvioiHyj/PRcQNjTieJLW7aVcYgP8JfKSB+7sQ+NjIJzPzX2TmXpm5F/AT4HsNPKYkta1pVxgy807gN5XPRcT7yp/87y9/0n//JPb3Y+DlsZZHxGzgUOCGKYYsSdNKX9EBNMgS4OzMXBUR+wF/SymZN8JxwI8z86UG7U+S2tq0LwwRsRVwAHBtRAw9PbO87Hjgv4yy2ZrM/MMaD3EK8O1645Sk6WLaFwZKp8NeLI8FVMnM71HH2EBEbAvsS6lrkKSuMO3GGEYqn+J5OiJOBIiSP2jQ7k8EbszMVxu0P0lqe9OuMETENZSuEtotIgYj4kzgVODMiHgYeAw4dhL7uwu4FjisvL/KU0wnA9c0LnpJan/hbbclSZWmXccgSWquaTX4vN122+W8efOKDkOSppX777//hczcvtb1p1VhmDdvHitWrCg6DEmaViLil5NZ31NJkqQqFgZJUhULgySpioVBklTFwiBJqlJoYYiIj0TEkxHxVEScW2QskqSSwgpDRPQClwBHAYuAUyJiUVHxSJJKipzHsC/wVGb+AiAillK6x9HjjT7Qd69awha/fnD48XN9O3H32w6Z8v7mvvksH/r9csDbiUhqrnW9O7J8yyNZNHdrvvDHi1tyzCILw7uBZyseDwL7jVwpIs4CzgLYeeedp3SgBS/dy8CG0t23e0g20cM9sw4iY2oN01H/dANHvvIDNhETryxJdXh8iz1ZvuWRLT1mkYVhtKy62UfwzFxC6RvaGBgYmNJH9H0/fRlwWenBXX9Nz4/PZ+mZ+8CMWVPZHSy7Gn7+bnrOaXhzI0lVdgf+d4uPWeTg8yCwU8XjfuC5ph+1b2bp95t1fMVCJqPXNUma/oosDP8A7BIR8yNiC0rfffB/m37U4cLw2tT3kZtgiqehJKndFXYqKTPfjIjPALcCvcBlmflY0w/cVz59VFfHsAnCjkFSZyr07qqZeRNwU0sPOlwY7BgkaTTdl90aMsZgYZDUubovu9kxSNK4ui+72TFI0ri6L7sNdQwb7RgkaTTdl928XFWSxtV92a0hl6umhUFSx+q+7NawjsF5DJI6UxcWhkZNcOu+l05Sd+i+7OYYgySNq/uymx2DJI2r+7Jbrx2DJI2n+7JbTw/0bmHHIElj6M7s1jfLjkGSxtCd2a1vpvMYJGkM3ZndGtIxOI9BUmfq0sJQb8fgqSRJnas7s5tjDJI0pu7MbnYMkjSm7sxudgySNKbuzG52DJI0pu7Mbn2zLAySNIbuzG59M+s8leQ8BkmdqzuzW0M6BucxSOpMXVoY6u0YPJUkqXN1Z3ZzjEGSxtSd2c2OQZLG1J3ZbahjyJza9hYGSR2skOwWESdGxGMRsSkiBloewNDXe258fWrbWxgkdbCisttK4HjgzkKOPvwtblMcZ7AwSOpgfUUcNDOfAIiiLvnsq/PrPZ3HIKmDtX12i4izImJFRKxYt25dY3baN6v0245BkjbTtI4hIm4Ddhxl0XmZuazW/WTmEmAJwMDAwBRHi0cYLgxT7Ric4CapczWtMGTm4c3ad936HGOQpLF0Z3ZrSMfQnS+dpM5X1OWqx0XEILA/8IOIuLWlAdQ9+GxhkNS5iroq6Xrg+iKODTj4LEnj6M7sZscgSWPqzuxWd8fgPAZJnas7s5sdgySNqTuzW0PGGJzHIKkzdWlhsGOQpLF0Z3bzqiRJGlN3ZrehjuGV9fD6K5Pf3sIgqYN1Z3br6YUZW8I9fwN//X54Y5Kdg4VBUgfr3ux2yjWw+5/Cq7+D116a3LYWBkkdrHuz23sPhvcdUvp7MmMNmYDzGCR1ru7OblO5md7Q90RbGCR1qO7OblO5/XZuKv22MEjqUN2d3abUMQwVBie4SepMXV4Y7BgkaaTuzm5TmujmGIOkztbd2W0qt8awY5DU4bo7u02lY7AwSOpw3Z3d7BgkaTPdnd3sGCRpM92d3abUMTj4LKmzdXd2q6tjcB6DpM7U3YWh1zEGSRqpu7NbTw/0buEYgyRVMLv1zbJjkKQKZre+mXYMklTB7GbHIElVCsluEXFhRPwsIh6JiOsj4h1FxAHYMUjSCEVltx8Bu2fmnsA/Av+poDjsGCRphEKyW2b+MDPfLD+8F+gvIg5gCh2DE9wkdbZ2yG5nADePtTAizoqIFRGxYt26dY0/+pQ7Bie4SepMTSsMEXFbRKwc5efYinXOA94Erh5rP5m5JDMHMnNg++23b3ygjjFIUpW+Zu04Mw8fb3lEfBz4I+CwzKHzMwXonQlvTqITsTBI6nBNKwzjiYiPAP8ROCgzXykihmF9Mx18lqQKRWW3bwGzgR9FxEMR8T8KiqM8xuCpJEkaUkjHkJkLijjuqOwYJKmK2c2OQZKqmN3sGCSpitltaB5DrRdGOcFNUoczu/XNBBI2vlHb+k5wk9ThLAyT/XpPTyVJ6nBmt75Jfr2nhUFShzO72TFIUhWz23BhsGOQJLAwVJxKsmOQJLAw2DFI0ghmNzsGSapidpv04LMT3CR1NrPblC9XdYKbpM5kYfByVUmqYnZzgpskVTG72TFIUhWzmx2DJFUxu9kxSFIVs5sdgyRVmTC7RcQOozy3W3PCKUBPL/TMcB6DJJXVkt3uioiThh5ExL8Hrm9eSAUY+ha3WjiPQVKH66thnYOBJRFxIvAu4Alg32YG1XJ9M+HJH8Cbv4ejLoTecV4WTyVJ6nATZrfMXAvcAuwPzAOuzMwNTY6rtRb/Cbz5Oqy4DH67evx1LQySOlwtYww/AvYDdgeOBr4REV9vdmAt9dGL4A//a+nv3Dj+uhYGSR2ulux2M/CXmfliZq4EDgB+19ywCtBTPn206c3x17MwSOpwtWS32cCtEXFXRHwa2DYzv9TkuFrPwiBJQG1jDOdn5mLg08Bc4P9FxG1Nj6zVLAySBExugtvzwK+A9cBmcxsmIyK+FBGPRMRDEfHDiJhbz/4aoqe39HuTYwySulstg8+fjIjlwI+B7YA/y8w96zzuhZm5Z2buBdwI/FWd+6tfzR2DE9wkdbZa5jG8B/jzzHyoUQfNzJcqHr4dyEbte8omfSrJCW6SOtOEhSEzz23GgSPiy8BplK5wOmSc9c4CzgLYeeedmxFKiWMMkgQ08SZ6EXFbRKwc5edYgMw8LzN3Aq4GPjPWfjJzSWYOZObA9ttv36xwKwpDrWMMdgySOlMtp5KmJDMPr3HV7wI/AL7QrFhqMjz4bMcgqbsVkt0iYpeKh8cAPysijiqeSpIkoIkdwwQuKN+6exPwS+DsguJ4i4VBkoCCCkNmnlDEccc16TEGC4OkzmR2G+IYgyQBFoa3OMFNkgALw1scY5AkwMLwFmc+SxJgYXjLZAaf7RYkdTAz3JDJDD5bGCR1MDPckMmcSrIwSOpgZrghFgZJAiwMb3GMQZIAC8NbhpJ9LfMYLAySOpgZbkhEqWvwVJKkLmeGq1RzYXAOg6TOZWGo1NPnGIOkrmeGq9TT66kkSV3PDFfJMQZJsjBUsTBIkoWhioVBkiwMVXp6HXyW1PXMcJVq6hic4Caps5nhKjmPQZIsDFUcY5AkC0MVxxgkycJQxY5BkiwMVSwMkmRhqGJhkCQLQxVvoidJFoYqNd1Ez3kMkjpboRkuIj4fERkR2xUZxzDnMUhScYUhInYCjgCeKSqGzTjGIEmFdgzfAP4CyAJjqOYYgyQVUxgi4hhgTWY+XMO6Z0XEiohYsW7duuYG5hf1SBJ9zdpxRNwG7DjKovOAvwSOrGU/mbkEWAIwMDDQ3O7CU0mS1LzCkJmHj/Z8ROwBzAcejtIgbj/wQETsm5m/alY8NbEwSFLzCsNYMvNRYIehxxGxGhjIzBdaHctmHGOQJOcxVIke5zFI6not7xhGysx5RccwzFNJkmTHUMUJbpJkYajiGIMkWRiq+EU9kmRhqOIYgyRZGKpYGCTJwlClpw9yY+mS1LFYGCR1uMIvV20rPeWXY9NG6B3jpbEwSG3rjTfeYHBwkFdffbXoUAoxa9Ys+vv7mTFjRl37sTBU6ukt/d705jiFwQluUrsaHBxk9uzZzJs3j+iyy8ozk/Xr1zM4OMj8+fPr2pcZrtJwxzDOOIPzGKS29eqrr7Ltttt2XVEAiAi23XbbhnRLFoZKNRcGXzapXXVjURjSqP92M1ylyjGGsVgYJHU4M1ylyjGGsVgYJFU44IADRn3+9NNP57rrrmtxNI1hhqvkqSRJk3TPPfcUHULDeVVSpZoKg1clSXrLVlttxYYNG8hMPvvZz3L77bczf/58crz5UG3ODFfJjkHSFF1//fU8+eSTPProo1x66aXTupMww1UaHmMYb/DZjkHS5u68805OOeUUent7mTt3LoceemjRIU2ZGa6S8xgk1aFTLpW1MFTyVJKkKTrwwANZunQpGzduZO3atdxxxx1FhzRlDj5XsjBImqLjjjuO22+/nT322INdd92Vgw46qOiQpszCUMkJbpImacOGDUDpNNK3vvWtgqNpDDNcJSe4SZKFoYqnkiTJwlDFwiBJFoYqNY0xOI9BUmczw1VyjEGSLAxVnOAmScUUhoj4YkSsiYiHyj9HFxHHZhxjkFSnZ599lkMOOYSFCxeyePFivvnNb262Tmbyuc99jgULFrDnnnvywAMPFBDp2Iqcx/CNzPx6gcffnIVBUp36+vq46KKL2HvvvXn55ZfZZ599OOKII1i0aNHwOjfffDOrVq1i1apV3HfffXzyk5/kvvvuKzDqak5wq1TTTfQsDNJ0cP73H+Px515q6D4Xzd2aL/zx4nHXmTNnDnPmzAFg9uzZLFy4kDVr1lQVhmXLlnHaaacREXzwgx/kxRdfZO3atcPbFa3IDPeZiHgkIi6LiHeOtVJEnBURKyJixbp165obkR2DpAZavXo1Dz74IPvtt1/V82vWrGGnnXYaftzf38+aNWtaHd6YmtYxRMRtwI6jLDoP+DvgS0CWf18EnDHafjJzCbAEYGBgoLnffGFhkDrGRJ/sm23Dhg2ccMIJXHzxxWy99dZVy0b7Ep92ujNr0wpDZh5ey3oRcSlwY7PimJSJCkMm4DwGSeN74403OOGEEzj11FM5/vjjN1ve39/Ps88+O/x4cHCQuXPntjLEcRV1VVLlibTjgJVFxLGZiSa4DVV5C4OkMWQmZ555JgsXLuScc84ZdZ1jjjmGK6+8kszk3nvvZZtttmmb8QUobvD5v0XEXpROJa0G/k1BcVSbaIJbbir9bqOWT1J7ufvuu7nqqqvYY4892GuvvQD4yle+wjPPPAPA2WefzdFHH81NN93EggUL2HLLLbn88ssLjHhzhRSGzPxYEced0ISnkoYKgx2DpNF9+MMfHnUMoVJEcMkll7Qooskzw1WyMEiShaHKhGMMFgZJnc8MV6nmMQZfNkmdywxXKQKi18IgqauZ4Ubq6bMwSOpqZriRLAySupwZbqSePie4SZqyM844gx122IHdd999+Lnf/OY3HHHEEeyyyy4cccQR/Pa3vx1e9tWvfpUFCxaw2267ceutt466z/G2bwYz3Eg9tYwxOMFN0uhOP/10brnllqrnLrjgAg477DBWrVrFYYcdxgUXXADA448/ztKlS3nssce45ZZb+NSnPsXGjZt/MB1r+2bxttsjeSpJ6gw3nwu/erSx+9xxDzhq/KR84IEHsnr16qrnli1bxvLlywH4+Mc/zsEHH8zXvvY1li1bxsknn8zMmTOZP38+CxYs4Kc//Sn7779/Tds3ixluJAuDpAb79a9/PXwvpDlz5vD8888Dtd9+e6ztm8WOYaRxxxgsDNK0McEn+3bQrrffNsONVNMYgy+bpNq9613vYu3atQCsXbuWHXbYAaj99ttjbd8sZriRPJUkqcGOOeYYrrjiCgCuuOIKjj322OHnly5dymuvvcbTTz/NqlWr2HfffWvevlk8lTRSTx+s+hFcst/myza+Xv6j+FZPUns65ZRTWL58OS+88AL9/f2cf/75nHvuuZx00kl85zvfYeedd+baa68FYPHixZx00kksWrSIvr4+LrnkEnp7S7fm+cQnPsHZZ5/NwMDAmNs3S0x0e9h2MjAwkCtWrGjuQR64Ep66bezlvTPh8C/ANv3NjUPSpD3xxBMsXLiw6DAKNdprEBH3Z+ZArfuwYxhp79NKP5LUpTxZLkmqYmGQ1FGm0+nxRmvUf7uFQVLHmDVrFuvXr+/K4pCZrF+/nlmzZtW9L8cYJHWM/v5+BgcHWbduXdGhFGLWrFn099d/YYyFQVLHmDFjBvPnzy86jGnPU0mSpCoWBklSFQuDJKnKtJr5HBHrgF9OcfPtgBcaGE4rGHNrGHNrGHNrjBbzezJz+1p3MK0KQz0iYsVkpoS3A2NuDWNuDWNujUbE7KkkSVIVC4MkqUo3FYYlRQcwBcbcGsbcGsbcGnXH3DVjDJKk2nRTxyBJqoGFQZJUpSMKQ0R8JCKejIinIuLcUZZHRPxNefkjEbF3rdu2W8wRsVNE3BERT0TEYxHxb9s53orlvRHxYETc2Ip46405It4REddFxM/Kr/X+0yDmf1f+N7EyIq6JiPpvs9mYmN8fET+JiNci4vOT2bbdYi7q/VdPzBXLa38PZua0/gF6gZ8D7wW2AB4GFo1Y52jgZkpf1vxB4L5at23DmOcAe5f/ng38Y7NjrifeiuXnAN8Fbmz3fxflZVcAnyj/vQXwjnaOGXg38DTwtvLj/wOc3iYx7wD8c+DLwOcns20bxtzy91+9MVcsr/k92Akdw77AU5n5i8x8HVgKHDtinWOBK7PkXuAdETGnxm3bKubMXJuZDwBk5svAE5SSQlvGCxAR/cBHgW83Oc6GxBwRWwMHAt8ByMzXM/PFdo65vKwPeFtE9AFbAs+1Q8yZ+Xxm/gPwxmS3bbeYC3r/1RUzTP492AmF4d3AsxWPB9n8f9RY69SybTPUE/OwiJgHfAC4r/EhTi6WCda5GPgLYFOT4htNPTG/F1gHXF5uvb8dEW9vZrATxDPhOpm5Bvg68AywFvhdZv6wibGOG08Ltq1HQ47bwvcf1B/zxUziPdgJhSFGeW7kNbhjrVPLts1QT8ylhRFbAX8P/HlmvtTA2EYz5Xgj4o+A5zPz/saHNa56XuM+YG/g7zLzA8A/Aa04/13P6/xOSp8g5wNzgbdHxL9qcHyjqec91M7vv/F30Nr3H9QR81Teg51QGAaBnSoe97N5Cz3WOrVs2wz1xExEzKD0j/LqzPxeE+OcMJYa1vkQcExErKbU/h4aEf+reaFOGE8t6wwCg5k59EnwOkqFotnqiflw4OnMXJeZbwDfAw5oYqwTxdPsbetR13ELeP9BfTFP/j3Y7EGTZv9Q+nT3C0qflIYGZRaPWOejVA/Y/bTWbdsw5gCuBC6eDq/xiHUOpnWDz3XFDNwF7Fb++4vAhe0cM7Af8BilsYWgNHj+2XaIuWLdL1I9kNu2779xYm75+6/emEcsq+k92LL/sCa/aEdTujrg58B55efOBs6u+J95SXn5o8DAeNu2c8zAhym1kI8AD5V/jm7XeKfyj7IdYgb2AlaUX+cbgHdOg5jPB34GrASuAma2Scw7UvrE+xLwYvnvrcfatp1jLur9V+/rXLGPmt6D3hJDklSlE8YYJEkNZGGQJFWxMEiSqlgYJElVLAySpCoWBnWt8h1UP1XxeG5EXNekY/1JRPzVBOt8PSIObcbxpcnwclV1rfK9bm7MzN1bcKx7gGMy84Vx1nkPcGlmHtnseKTx2DGom10AvC8iHoqICyNiXkSsBIiI0yPihoj4fkQ8HRGfiYhzyjfVuzci/ll5vfdFxC0RcX9E3BUR7x95kIjYFXgtM1+IiNnl/c0oL9s6IlZHxIzM/CWwbUTs2MLXQNqMhUHd7Fzg55m5V2b+h1GW7w78S0q3PP4y8EqWbqr3E+C08jpLKN16Yh/g88DfjrKfDwGVt2peTunWFgAnA3+fpfsbUV7vQ3X+d0l16Ss6AKmN3VFO5C9HxO+A75effxTYs3yHzQOAayOGb345c5T9zKF0G+8h36Z0C+QbgH8N/FnFsucp3R1VKoyFQRrbaxV/b6p4vInSe6cHeDEz95pgP78Hthl6kJl3l09bHQT0ZubKinVnldeXCuOpJHWzlyl9PeOUZOk+/E9HxIkw/H3MfzDKqk8AC0Y8dyVwDXD5iOd3pXQTPKkwFgZ1rcxcD9wdESsj4sIp7uZU4MyIeJjSba9H+2rKO4EPRMX5JuBq4J2UigMwfJ//BZTu6ioVxstVpRaIiG8C38/M28qP/xQ4NjM/VrHOcZS+aP4/FxSmBDjGILXKVyh9mQ4R8d+BoyjdX79SH3BRi+OSNmPHIEmq4hiDJKmKhUGSVMXCIEmqYmGQJFWxMEiSqvx/wtUxbYbYMtUAAAAASUVORK5CYII=\n", + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYYAAAERCAYAAAB/4wAeAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAZd0lEQVR4nO3df5RkZXng8e/TP4YmMgSE4TA4jIwoyi8doQOCLkQEo5jAInpWZDWsrBPX1ZMsx8264SRqsrruIhvcaJIzitmFGHBFUEElBiHiQUwcFGUQUZFRB1CGQYUJgkP3s3/c6raqp7q7uuvHrbn3+zmnT/26de9DMfU+932f+74VmYkkSTNGyg5AkjRcTAySpBYmBklSCxODJKmFiUGS1MLEIElqsdslhoj4SEQ8GBGbe7S/6yPiZxFx3ZznvxQRtzf+7o+IT/bieJI07Ha7xAD8H+BlPdzfRcDr5j6Zmf8qM9dn5nrgVuDqHh5TkobWbpcYMvNm4OHm5yLi0MaZ/22NM/3nLGF/XwAene/1iFgJnAJ8cpkhS9JuZazsAHpkI/CmzPxuRBwP/CVFY94LZwFfyMxHerQ/SRpqu31iiIi9gBOBj0fEzNN7NF57JfCnbd52X2b+VoeHOAf4cLdxStLuYrdPDBTDYT9r1AJaZObVdFEbiIj9gOMoeg2SVAu7XY1hrsYQz70R8WqAKDyvR7t/NXBdZj7eo/1J0tArNTEs59LTiLiC4iqhZ0fE1og4HzgXOD8ivgHcCZy5hP19Cfg48JLG/pqHmF4DXNHpviSpCqLMZbcj4iRgB3BZZh5VWiCSpFml9hjaXXoqSSrXblV83n///fOQQw4pOwxJ2q3cdtttD2Xmqk63H/rEEBEbgA0Aa9euZdOmTSVHJEm7l4j4wVK2H/qrkjJzY2ZOZubkqlUdJzxJ0jINfWKQJA1W2Zertrv0VJJUolJrDJl5TpnHlyTtyqEkSVILE4MkqYWJQZLUon6J4ZH74dufKTsKSRpa9UsMt34QPvY6KHGNKEkaZvVLDA9/H3IKpp8sOxJJGkr1TAwAU78sNw5JGlL1SgzT0/DwvcV9E4MktVWvxPDo/TD1RHF/ame5sUjSkKpXYpgZRgJ7DJI0DxODJKlFzRLDvb+671CSJLVVs8Rgj0GSFlOzxHAvjDQWlDUxSFJb9UkMmUWPYb9nFo8dSpKktuqTGHY8CDv/BVY9u3hsj0GS2qpPYpipL6x6TnFrYpCktuqTGH7auCJp/8OKW4eSJKmt+iSGh78PMQr7HVo8tscgSW3VKzHsczCs2Kt4bI9BktqqV2J46jNgdLx4bI9BktqqR2LIhO2NxDBiYpCkhdQjMfzip/DEz2HfdTC6onjOoSRJaqvUxBARL4uIuyPiexHx9r4daOZSVYeSJGlRpSWGiBgFPgi8HDgCOCcijujLwWYWz3vqM5p6DCYGSWpnrMRjHwd8LzO/DxARVwJnAt/q9YFuuvUrnEzw+qt/whQPcQXwsa/cw9V33rqs/R28815O+sUXgOxpnJI0109GD+KGp7yCIw7am3f8zpEDOWaZieFpwI+aHm8Fjp+7UURsADYArF27dlkH2nPqUR4aPYCdsQIymSYYZWpZ+wL4rceu5bTHPsvjscey9yFJnbhrxdHc8JRXDPSYZSaGaPPcLqfgmbkR2AgwOTm5rFP0F7x5I0xP8bGR0eKJP1vBq553AK867YTl7A4+9bdwz9OYuKDnnRtJavF84GMDPmaZxeetwMFNj9cA9/ftaDNJAYo6QzdXJWVC1OOCLkn1U2br9lXgWRGxLiJWAK8BPj2QI4+Od1d8zmmIdh0eSdr9lTaUlJlPRsRbgL8HRoGPZOadAzn46IoeJAZ7DJKqqcwaA5n5WeCzAz9w10NJ08WCfJJUQfU87e12KGl6yh6DpMqqZ+vmUJIkzauerdvoeA+Gkur50Umqvnq2br3oMYxYY5BUTSaG5fByVUkVVtPE4FCSJM2nnq2bxWdJmlc9W7eezGOo50cnqfrq2br1ZB6DxWdJ1VTTxOBQkiTNp56tm0NJkjSverZuXa+u6rLbkqqrnq1b10NJUzBSz49OUvXVs3VzHoMkzauerZvFZ0maVz1bt9EVML2zqBUsh4lBUoXVs3UbHS9ulzuc5DwGSRVW08Sworhd7nCSPQZJFVbP1q3rxODlqpKqq56tW7dDSS67LanCapoYuu0xTPlDPZIqy8SwHNYYJFVYPVu3ngwl1fOjk1R9pbRuEfHqiLgzIqYjYnLgAdhjkKR5ldW6bQZeCdxcytFnE4PzGCRprrEyDpqZdwFEWVf2zA4lebmqJM019K1bRGyIiE0RsWnbtm292alDSZI0r771GCLiBuDANi9dmJmf6nQ/mbkR2AgwOTm5zMWN5uh2KMl5DJIqrG+JITNP7de+uzYzlDTdRWJwHoOkiqrneEgvJrg5lCSposq6XPWsiNgKnAB8JiL+fqAB9GQoycQgqZrKuirpGuCaMo4N9OCqJBODpOqqZ+vmVUmSNK96tm5dT3CbdoKbpMqqaWLoxVCSl6tKqqaaJgaHkiRpPvVs3VxdVZLmVc/WbWS0qBF869Nw038v1j5aCn+oR1KF1TMxADzrpfDzH8IX3wuPPby099pjkFRh9W3dXnslvPS/Ffef/MXS3mtikFRh9W7dxvYsbnc+3vl7ZoadTAySKqrerdvYHsXtUnoM01PFrfMYJFVUvRPD+HJ6DNPFrfMYJFVUvRPD2ERxu5Qew2xiqPdHJ6m66t26ddVjqPdHJ6m66t26LavH0KgxOI9BUkXVOzHYY5CkXdS7dbPGIEm7qHfrtqweg/MYJFVbvVu35fQYZucx1Pujk1Rd9W7dZhKDNQZJmlXv1m1kBEb3sMYgSU1s3cYn7DFIUhNbt7E9l9djcB6DpIoqJTFExEUR8e2I+GZEXBMR+5QRB7CMHoPFZ0nVVlbr9g/AUZn5XOA7wH8tKY7l9xhMDJIqqpTWLTM/n5lPNh5+BVhTRhzAMnoMzmOQVG3D0Lq9AfhcaUcf2xOetPgsSTPG+rXjiLgBOLDNSxdm5qca21wIPAl8dIH9bAA2AKxdu7b3gY5PwOOPdL69E9wkVVzfEkNmnrrQ6xHxu8BvAy/JnBmfabufjcBGgMnJyXm3W7axPeHJBzvf3h6DpIrrW2JYSES8DPgvwMmZ+VgZMcwan4CdFp8laUZZrdsHgJXAP0TE7RHx1yXFYY1BkuYopceQmc8s47htjU8sMTH4Qz2Sqs3T3jGXxJCkZrZuYxPFBLf569+tnMcgqeIWbd0i4vw5j0cj4h39C2nAxieKXsDUzs62n+0xRP9ikqQSdXLa+5KI+GxErI6IoyhmKq/sc1yDM9b4FbdOl8WYncdgjUFSNS1afM7M10bEvwHuAB4DzsnMW/oe2aCMN/1Yz8SvL769NQZJFdfJUNKzgN8HPgFsAV4XEb/W57gGZ6k9BhODpIrrpHW7FvjjzPw94GTgu8BX+xrVII0v8ec9TQySKq6TeQzHZeYjAI2lKy6OiE/3N6wBWnKPwXkMkqqtkxrDIxFxInDInO2/26+gBsoegyS1WDQxRMTlwKHA7UDjdJkELutfWAO05B6D8xgkVVsnQ0mTwBELrYC6W1t2j8F5DJKqqZPT3s20/12FanAegyS1mLfHEBHXUgwZrQS+FRH/DDwx83pmntH/8AbAGoMktVhoKOl9jdtjgUuBHwLVGz9xHoMktZg3MWTmFwEi4jeB84GHgSuBqzLzJ4MIbiDsMUhSi0Vbt8x8V2YeCfxH4CDgi43fc66G5fYYnMcgqaKWctr7IPBjYDtwQH/CKcHoeHH2b49BkoDO1kr6DxHxj8AXgP2BN2bmc/sd2MBELO3nPU0Mkiquk3kMTwf+IDNv73Ms5RmfgLuuhScegVf8r6IXMR/nMUiquE5qDG+vdFIAOPrVRYP/tcvgpz9YeNvZxGCNQVI1OR4C8PL/Aaf9aXF/+smFt52d4OZHJ6mabN1mzAwfTS/yE5/WGCRVnK3bjJFGuWWxHoOJQVLFldK6RcSfRcQ3I+L2iPh8RBxURhwtRho9hikTg6R6K6t1uygzn5uZ64HrgD8pKY5fmZmw1mmPwQlukiqqlMQw84twDU+hWKyvXNYYJAnobB5DX0TEu4HXAz8HXrzAdhuADQBr167tX0BLrjE4j0FSNfXttDciboiIzW3+zgTIzAsz82Dgo8Bb5ttPZm7MzMnMnFy1alW/wrXGIEkNfesxZOapHW76d8BngHf0K5aOdFpj8Id6JFVcWVclPavp4RnAt8uIo4U1BkkCyqsxvDcing1MAz8A3lRSHL/iPAZJAkpKDJl5dhnHXdBsYphaeDsTg6SKs3WbMZMYpjocSnIeg6SKMjHMmK0xOJQkqd5s3WbMDiV1Wnx2HoOkajIxzFhKjcHegqQKs4Wb0WmNYXrKOQySKs3EMGMpl6vaY5BUYbZwM5Yywc3EIKnCbOFmWGOQJMDE8CsRRe1g0XkM6RwGSZVmYmg2MtZBjWHKS1UlVZqJodnouMVnSbVnC9dsZNTEIKn2bOGajYx3tlaS8xgkVZiJoVknNYbpKXsMkirNFq6ZNQZJMjG06KjGkCYGSZVmC9es4xqDH5uk6rKFa9bpPIYRPzZJ1WUL12y0k8Rgj0FStdnCNeuox2BikFRttnDNRrwqSZJs4ZqNjPlDPZJqz8TQbHTMZbcl1V6pLVxEvC0iMiL2LzOOWSNjHfxQj/MYJFVbaS1cRBwMnAb8sKwYdtFxjcFltyVVV5mnvn8O/CGQJcbQamQMpjqZx2CNQVJ1lZIYIuIM4L7M/EYH226IiE0RsWnbtm39DcxltyWJsX7tOCJuAA5s89KFwB8BL+1kP5m5EdgIMDk52d/exeh4BzUGE4OkautbYsjMU9s9HxFHA+uAb0QxVr8G+FpEHJeZP+5XPB1xgpsk9S8xzCcz7wAOmHkcEVuAycx8aNCx7GJkfPEag/MYJFWcp77NXHZbkgbfY5grMw8pO4ZZndYYvCpJUoWVnhiGSsc1hvHBxCNpKOzcuZOtW7fy+OOPlx3KgiYmJlizZg3j4921USaGZs5jkNTG1q1bWblyJYcccggxpBNcM5Pt27ezdetW1q1b19W+HCxv5lVJktp4/PHH2W+//YY2KQBEBPvtt19PejW2cM2cxyBpHsOcFGb0KkZbuGYjY0XDPz09/zYmBkkVZwvXbKRRclloOCmnnccgaVlOPPHEts+fd955XHXVVQOOZn4mhmadJIZpV1eVtDxf/vKXyw6hI16V1Gy0cYnXQnUGh5IkLdNee+3Fjh07yEze+ta3cuONN7Ju3Toyh2eRabDH0Gq2x7DAr7iZGCR16ZprruHuu+/mjjvu4EMf+tDQ9SRs4Zp1XGPwY5O0fDfffDPnnHMOo6OjHHTQQZxyyillh9TCFq7ZTGKYWmgoyQlukro3zJe/mhiazdYY7DFI6p+TTjqJK6+8kqmpKR544AFuuummskNqYfG5mUNJkgbgrLPO4sYbb+Too4/msMMO4+STTy47pBYmhmYdJQaX3Za0PDt27ACKYaQPfOADJUczP1u4Zp3UGKanTAySKs0WrplDSZJkYmhh8VmSTAwtZi5DNTFIqjFbuGYjjR6D8xgk1ZiJoZk1BkkyMbToqMbg5aqSyvGGN7yBAw44gKOOOqqvx7GFa2aNQdIQO++887j++uv7fpxSJrhFxDuBNwLbGk/9UWZ+toxYWnRSY3Aeg1Rr77r2Tr51/yM93ecRB+3NO37nyEW3O+mkk9iyZUtPj91OmTOf/zwz31fi8XdljUGSXBKjhfMYJC2ikzP73V2ZLdxbIuKbEfGRiNh3vo0iYkNEbIqITdu2bZtvs96wxiBJ/UsMEXFDRGxu83cm8FfAocB64AHg4vn2k5kbM3MyMydXrVrVr3ALzmOQpP4lhsw8NTOPavP3qcz8SWZOZeY08CHguH7FsSSL1RhmfpfVHoOkEpxzzjmccMIJ3H333axZs4ZLL720L8cp66qk1Zn5QOPhWcDmMuLYxWI1hpwubk0MkkpwxRVXDOQ4ZRWf/2dErAcS2AL8XklxtFqsxjCbGIb3J/kkqVulJIbMfF0Zx13UYjWG2cRgjUFSdTkm0myxGsP0VHHrUJKkCrOFazabGKbav26NQVIN2MI1GxkpGv3pxYaS/NgkVZct3Fwj44sXn53HIKnCTAxzjYx1UHz2Y5M0WD/60Y948YtfzOGHH86RRx7J+9///r4dy7WS5hoZs8YgaeiMjY1x8cUXc8wxx/Doo49y7LHHctppp3HEEUf0/lg93+PubnSsgxqD8xik2vrc2+HHd/R2nwceDS9/74KbrF69mtWrVwOwcuVKDj/8cO67776+JAZPfecaGXPms6ShtmXLFr7+9a9z/PHH92X/9hjmGhmHqcXmMVh8lmprkTP7ftuxYwdnn302l1xyCXvvvXdfjuGp71wjo/YYJA2lnTt3cvbZZ3Puuefyyle+sm/HsYWba3TceQyShk5mcv7553P44YdzwQUX9PVYDiXNNTIG3/k8fLDN2N3UL4tbE4OkAbvlllu4/PLLOfroo1m/fj0A73nPezj99NN7fiwTw1wveDPc84X5X1/zG7DupMHFI0nAi170InLmN2H6zMQw17G/W/xJUk05JiJJamFikKQODGoYpxu9itHEIEmLmJiYYPv27UOdHDKT7du3MzEx0fW+rDFI0iLWrFnD1q1b2bZtW9mhLGhiYoI1a9Z0vR8TgyQtYnx8nHXr1pUdxsA4lCRJamFikCS1MDFIklrEMFfZ54qIbcAPlvn2/YGHehjOIBjzYBjzYBjzYLSL+emZuarTHexWiaEbEbEpMyfLjmMpjHkwjHkwjHkwehGzQ0mSpBYmBklSizolho1lB7AMxjwYxjwYxjwYXcdcmxqDJKkzdeoxSJI6YGKQJLWoRGKIiJdFxN0R8b2IeHub1yMi/nfj9W9GxDGdvnfYYo6IgyPipoi4KyLujIjfH+Z4m14fjYivR8R1g4i325gjYp+IuCoivt34rE/YDWL+T41/E5sj4oqI6H6Zzd7E/JyIuDUinoiIty3lvcMWc1nfv25ibnq98+9gZu7Wf8AocA/wDGAF8A3giDnbnA58DgjgBcA/dfreIYx5NXBM4/5K4Dv9jrmbeJtevwD4O+C6Yf930Xjt/wL/vnF/BbDPMMcMPA24F9iz8fj/AecNScwHAL8BvBt421LeO4QxD/z7123MTa93/B2sQo/hOOB7mfn9zPwlcCVw5pxtzgQuy8JXgH0iYnWH7x2qmDPzgcz8GkBmPgrcRdEoDGW8ABGxBngF8OE+x9mTmCNib+Ak4FKAzPxlZv5smGNuvDYG7BkRY8CvAfcPQ8yZ+WBmfhXYudT3DlvMJX3/uooZlv4drEJieBrwo6bHW9n1f9R823Ty3n7oJuZZEXEI8Hzgn3of4tJiWWSbS4A/BKb7FF873cT8DGAb8DeNrveHI+Ip/Qx2kXgW3SYz7wPeB/wQeAD4eWZ+vo+xLhjPAN7bjZ4cd4DfP+g+5ktYwnewCokh2jw39xrc+bbp5L390E3MxYsRewGfAP4gMx/pYWztLDveiPht4MHMvK33YS2om894DDgG+KvMfD7wL8Agxr+7+Zz3pTiDXAccBDwlIv5tj+Nrp5vv0DB//xbewWC/f9BFzMv5DlYhMWwFDm56vIZdu9DzbdPJe/uhm5iJiHGKf5Qfzcyr+xjnorF0sM0LgTMiYgtF9/eUiPjb/oW6aDydbLMV2JqZM2eCV1Ekin7rJuZTgXszc1tm7gSuBk7sY6yLxdPv93ajq+OW8P2D7mJe+new30WTfv9RnN19n+JMaaYoc+ScbV5Ba8Hunzt97xDGHMBlwCW7w2c8Z5vfZHDF565iBr4EPLtx/53ARcMcM3A8cCdFbSEoiudvHYaYm7Z9J62F3KH9/i0Q88C/f93GPOe1jr6DA/sP6/OHdjrF1QH3ABc2nnsT8Kam/5kfbLx+BzC50HuHOWbgRRRdyG8Ctzf+Th/WeJfzj3IYYgbWA5san/MngX13g5jfBXwb2AxcDuwxJDEfSHHG+wjws8b9ved77zDHXNb3r9vPuWkfHX0HXRJDktSiCjUGSVIPmRgkSS1MDJKkFiYGSVILE4MkqYWJQbXVWEH1zU2PD4qIq/p0rH8dEX+yyDbvi4hT+nF8aSm8XFW11Vjr5rrMPGoAx/oycEZmPrTANk8HPpSZL+13PNJC7DGozt4LHBoRt0fERRFxSERsBoiI8yLikxFxbUTcGxFviYgLGovqfSUintrY7tCIuD4ibouIL0XEc+YeJCIOA57IzIciYmVjf+ON1/aOiC0RMZ6ZPwD2i4gDB/gZSLswMajO3g7ck5nrM/M/t3n9KOC1FEsevxt4LItF9W4FXt/YZiPF0hPHAm8D/rLNfl4INC/V/I8US1sAvAb4RBbrG9HY7oVd/ndJXRkrOwBpiN3UaMgfjYifA9c2nr8DeG5jhc0TgY9HzC5+uUeb/aymWMZ7xocplkD+JPDvgDc2vfYgxeqoUmlMDNL8nmi6P930eJriuzMC/Cwz1y+yn18Avz7zIDNvaQxbnQyMZubmpm0nGttLpXEoSXX2KMXPMy5LFuvw3xsRr4bZ32N+XptN7wKeOee5y4ArgL+Z8/xhFIvgSaUxMai2MnM7cEtEbI6Ii5a5m3OB8yPiGxTLXrf7acqbgedH03gT8FFgX4rkAMyu8/9MilVdpdJ4uao0ABHxfuDazLyh8fhVwJmZ+bqmbc6i+KH5Py4pTAmwxiANynsofkyHiPgL4OUU6+s3GwMuHnBc0i7sMUiSWlhjkCS1MDFIklqYGCRJLUwMkqQWJgZJUov/D7bkwCZGuX/yAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] @@ -103,7 +103,7 @@ } ], "source": [ - "swiftdiff['vx'].plot.line(x=\"time (y)\")" + "swiftdiff['vhx'].plot.line(x=\"time (y)\")" ] }, { diff --git a/src/io/io.f90 b/src/io/io.f90 index 226811d61..29562763b 100644 --- a/src/io/io.f90 +++ b/src/io/io.f90 @@ -27,7 +27,7 @@ module subroutine io_conservation_report(self, param, lterminal) "; DM/M0 = ", ES12.5)' associate(system => self, pl => self%pl, cb => self%cb, npl => self%pl%nbody) - if ((param%out_type == REAL4_TYPE) .or. (param%out_type == REAL8_TYPE) .and. (param%energy_out /= "")) then + if (((param%out_type == REAL4_TYPE) .or. (param%out_type == REAL8_TYPE)) .and. (param%energy_out /= "")) then if (param%lfirstenergy .and. (param%out_stat /= "OLD")) then open(unit = EGYIU, file = param%energy_out, form = "formatted", status = "replace", action = "write", err = 667, iomsg = errmsg) write(EGYIU,EGYHEADER, err = 667, iomsg = errmsg) @@ -35,8 +35,10 @@ module subroutine io_conservation_report(self, param, lterminal) open(unit = EGYIU, file = param%energy_out, form = "formatted", status = "old", action = "write", position = "append", err = 667, iomsg = errmsg) end if end if + call pl%vb2vh(cb) call pl%xh2xb(cb) + call system%get_energy_and_momentum(param) ke_orbit_now = system%ke_orbit ke_spin_now = system%ke_spin @@ -56,7 +58,7 @@ module subroutine io_conservation_report(self, param, lterminal) param%lfirstenergy = .false. end if - if ((param%out_type == REAL4_TYPE) .or. (param%out_type == REAL8_TYPE) .and. (param%energy_out /= "")) then + if (((param%out_type == REAL4_TYPE) .or. (param%out_type == REAL8_TYPE)) .and. (param%energy_out /= "")) then write(EGYIU,EGYFMT, err = 667, iomsg = errmsg) param%t, Eorbit_now, param%Ecollisions, Ltot_now, GMtot_now close(EGYIU, err = 667, iomsg = errmsg) end if @@ -71,11 +73,11 @@ module subroutine io_conservation_report(self, param, lterminal) if (abs(Merror) > 100 * epsilon(Merror)) then write(*,*) "Severe error! Mass not conserved! Halting!" call pl%xv2el(cb) - call param%nciu%open(param) + !call param%nciu%open(param) call self%write_hdr(param%nciu, param) call cb%write_frame(param%nciu, param) call pl%write_frame(param%nciu, param) - call param%nciu%close(param) + call param%nciu%close() call util_exit(FAILURE) end if end if @@ -200,9 +202,7 @@ module subroutine io_dump_particle_info_base(self, param, idx) close(unit = LUN, err = 667, iomsg = errmsg) !else if ((param%out_type == NETCDF_FLOAT_TYPE) .or. (param%out_type == NETCDF_DOUBLE_TYPE)) then if ((param%out_type == NETCDF_FLOAT_TYPE) .or. (param%out_type == NETCDF_DOUBLE_TYPE)) then - call param%nciu%open(param) call self%write_particle_info(param%nciu) - call param%nciu%close(param) end if return @@ -275,41 +275,52 @@ module subroutine io_dump_system(self, param) real(DP) :: tfrac character(*), parameter :: statusfmt = '("Time = ", ES12.5, "; fraction done = ", F6.3, "; Number of active pl, tp = ", I5, ", ", I5)' character(*), parameter :: symbastatfmt = '("Time = ", ES12.5, "; fraction done = ", F6.3, "; Number of active plm, pl, tp = ", I5, ", ", I5, ", ", I5)' - logical, save :: lfirst = .true. - - if (lfirst) then - lfirst = .false. - if (param%lenergy) call self%conservation_report(param, lterminal=.false.) - else - allocate(dump_param, source=param) - param_file_name = trim(adjustl(DUMP_PARAM_FILE(idx))) + + allocate(dump_param, source=param) + param_file_name = trim(adjustl(DUMP_PARAM_FILE(idx))) + dump_param%in_form = XV + dump_param%out_form = XV + dump_param%out_stat = 'APPEND' + if ((param%out_type == REAL8_TYPE) .or. (param%out_type == REAL4_TYPE)) then dump_param%incbfile = trim(adjustl(DUMP_CB_FILE(idx))) dump_param%inplfile = trim(adjustl(DUMP_PL_FILE(idx))) dump_param%intpfile = trim(adjustl(DUMP_TP_FILE(idx))) - dump_param%in_form = XV - dump_param%out_stat = 'APPEND' dump_param%in_type = REAL8_TYPE - dump_param%T0 = param%t + else if ((param%out_type == NETCDF_FLOAT_TYPE) .or. (param%out_type == NETCDF_DOUBLE_TYPE)) then + dump_param%outfile = trim(adjustl(DUMP_NC_FILE(idx))) + dump_param%in_type = NETCDF_DOUBLE_TYPE + end if + dump_param%T0 = param%t + dump_param%ioutput = 0 - call dump_param%dump(param_file_name) + call dump_param%dump(param_file_name) - dump_param%out_form = XV + if ((param%out_type == REAL8_TYPE) .or. (param%out_type == REAL4_TYPE)) then call self%cb%dump(dump_param) call self%pl%dump(dump_param) call self%tp%dump(dump_param) + else if ((param%out_type == NETCDF_FLOAT_TYPE) .or. (param%out_type == NETCDF_DOUBLE_TYPE)) then + call dump_param%nciu%initialize(dump_param) + call self%write_hdr(dump_param%nciu, dump_param) + call self%cb%write_frame(dump_param%nciu, dump_param) + call self%pl%write_frame(dump_param%nciu, dump_param) + call self%tp%write_frame(dump_param%nciu, dump_param) + call dump_param%nciu%close() + ! Syncrhonize the disk and memory buffer of the NetCDF file (e.g. commit the frame files stored in memory to disk) + call param%nciu%sync() + end if - idx = idx + 1 - if (idx > NDUMPFILES) idx = 1 + idx = idx + 1 + if (idx > NDUMPFILES) idx = 1 - tfrac = (param%t - param%t0) / (param%tstop - param%t0) - - select type(pl => self%pl) - class is (symba_pl) - write(*, symbastatfmt) param%t, tfrac, pl%nplm, pl%nbody, self%tp%nbody - class default - write(*, statusfmt) param%t, tfrac, pl%nbody, self%tp%nbody - end select - end if + tfrac = (param%t - param%t0) / (param%tstop - param%t0) + + select type(pl => self%pl) + class is (symba_pl) + write(*, symbastatfmt) param%t, tfrac, pl%nplm, pl%nbody, self%tp%nbody + class default + write(*, statusfmt) param%t, tfrac, pl%nbody, self%tp%nbody + end select return end subroutine io_dump_system @@ -393,7 +404,7 @@ module function io_get_old_t_final_system(self, param) result(old_t_final) class(swiftest_nbody_system), allocatable :: tmpsys class(swiftest_parameters), allocatable :: tmpparam integer(I4B) :: ierr, iu = LUN - character(len=STRMAX) :: errmsg + character(len=STRMAX) :: errmsg old_t_final = 0.0_DP allocate(tmpsys, source=self) @@ -1143,7 +1154,42 @@ module subroutine io_param_writer_one_QP(param_name, param_value, unit) end subroutine io_param_writer_one_QP - module subroutine io_read_in_body(self, param) + module subroutine io_read_in_base(self,param) + !! author: Carlisle A. Wishard and David A. Minton + !! + !! Reads in either a central body, test particle, or massive body object. For the swiftest_body types (non-central body), it allocates array space for them + implicit none + class(swiftest_base), intent(inout) :: self !! Swiftest base object + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters + ! Internals + integer(I4B) :: ierr !! Error code: returns 0 if the read is successful + + select case(param%in_type) + case(NETCDF_DOUBLE_TYPE, NETCDF_FLOAT_TYPE) + select type(self) + class is (swiftest_body) + if (self%nbody == 0) return + call self%setup(self%nbody, param) + end select + + ierr = self%read_frame(param%nciu, param) + if (ierr == 0) return + case default + select type(self) + class is (swiftest_body) + call io_read_in_body(self, param) + class is (swiftest_cb) + call io_read_in_cb(self, param) + end select + return + end select + + 667 continue + write(*,*) "Error reading body in io_read_in_base" + end subroutine io_read_in_base + + + subroutine io_read_in_body(self, param) !! author: The Purdue Swiftest Team - David A. Minton, Carlisle A. Wishard, Jennifer L.L. Pouplin, and Jacob R. Elliott !! !! Read in either test particle or massive body data @@ -1157,20 +1203,19 @@ module subroutine io_read_in_body(self, param) ! Internals integer(I4B) :: iu = LUN integer(I4B) :: i, nbody - logical :: is_ascii, is_pl + logical :: is_ascii character(len=:), allocatable :: infile real(DP) :: t character(STRMAX) :: errmsg integer(I4B) :: ierr ! Select the appropriate polymorphic class (test particle or massive body) + select type(self) class is (swiftest_pl) infile = param%inplfile - is_pl = .true. class is (swiftest_tp) infile = param%intpfile - is_pl = .false. end select is_ascii = (param%in_type == 'ASCII') @@ -1207,7 +1252,7 @@ module subroutine io_read_in_body(self, param) end subroutine io_read_in_body - module subroutine io_read_in_cb(self, param) + subroutine io_read_in_cb(self, param) !! author: David A. Minton !! !! Reads in central body data @@ -1268,6 +1313,30 @@ module subroutine io_read_in_cb(self, param) end subroutine io_read_in_cb + module subroutine io_read_in_system(self, param) + !! author: David A. Minton and Carlisle A. Wishard + !! + !! Reads in the system from input files + implicit none + ! Arguments + class(swiftest_nbody_system), intent(inout) :: self + class(swiftest_parameters), intent(inout) :: param + ! Internals + integer(I4B) :: ierr + + if ((param%in_type == NETCDF_DOUBLE_TYPE) .or. (param%in_type == NETCDF_FLOAT_TYPE)) then + ierr = self%read_frame(param%nciu, param) + if (ierr /=0) call util_exit(FAILURE) + else + call self%cb%read_in(param) + call self%pl%read_in(param) + call self%tp%read_in(param) + end if + + return + end subroutine io_read_in_system + + function io_read_encounter(t, id1, id2, Gmass1, Gmass2, radius1, radius2, & xh1, xh2, vh1, vh2, enc_out, out_type) result(ierr) !! author: David A. Minton @@ -1787,9 +1856,7 @@ module subroutine io_write_discard(self, param) ! Record the discarded body metadata information to file if ((param%out_type == NETCDF_FLOAT_TYPE) .or. (param%out_type == NETCDF_DOUBLE_TYPE)) then - call param%nciu%open(param) call tp_discards%write_particle_info(param%nciu) - call param%nciu%close(param) end if if (param%discard_out == "") return @@ -2094,29 +2161,20 @@ module subroutine io_write_frame_system(self, param) errmsg = param%outfile // " not found! You must specify OUT_STAT = NEW, REPLACE, or UNKNOWN" goto 667 end if + call param%nciu%open(param) case('NEW') if (fileExists) then - errmsg = param%outfile // " Alread Exists! You must specify OUT_STAT = OLD, REPLACE, or UNKNOWN" + errmsg = param%outfile // " Alread Exists! You must specify OUT_STAT = APPEND, REPLACE, or UNKNOWN" goto 667 end if + call param%nciu%initialize(param) case('REPLACE', 'UNKNOWN') - if (fileExists) then - open(file=param%outfile, unit=iu, status='OLD') - close (unit=BINUNIT, status="delete") - end if - end select - - select case(param%out_stat) - case('APPEND') - call param%nciu%open(param) - case('NEW', 'REPLACE', 'UNKNOWN') call param%nciu%initialize(param) - call param%nciu%close(param) - call param%nciu%open(param) end select + lfirst = .false. else - call param%nciu%open(param) + !call param%nciu%open(param) end if end if @@ -2142,7 +2200,6 @@ module subroutine io_write_frame_system(self, param) call cb%write_frame(param%nciu, param) call pl%write_frame(param%nciu, param) call tp%write_frame(param%nciu, param) - call param%nciu%close(param) end if return diff --git a/src/main/swiftest_driver.f90 b/src/main/swiftest_driver.f90 index bd4d50855..755468868 100644 --- a/src/main/swiftest_driver.f90 +++ b/src/main/swiftest_driver.f90 @@ -57,11 +57,13 @@ program swiftest_driver ioutput = ioutput_t0 ! Prevent duplicate frames from being written if this is a restarted run if ((param%lrestart) .and. ((param%out_type == REAL8_TYPE) .or. param%out_type == REAL4_TYPE)) then - old_t_final = nbody_system%get_old_t_final(param) + old_t_final = nbody_system%get_old_t_final_bin(param) + else if ((param%lrestart) .and. ((param%out_type == NETCDF_DOUBLE_TYPE) .or. param%out_type == NETCDF_FLOAT_TYPE)) then + old_t_final = nbody_system%get_old_t_final_netcdf(param) else old_t_final = t0 if (istep_out > 0) call nbody_system%write_frame(param) - call nbody_system%dump(param) + if (param%lenergy) call nbody_system%conservation_report(param, lterminal=.false.) ! This will save the initial values of energy and momentum end if !> Define the maximum number of threads @@ -104,6 +106,8 @@ program swiftest_driver end do end associate + call nbody_system%finalize(param) + call util_exit(SUCCESS) stop diff --git a/src/modules/swiftest_classes.f90 b/src/modules/swiftest_classes.f90 index 4776dc0ba..0f341a392 100644 --- a/src/modules/swiftest_classes.f90 +++ b/src/modules/swiftest_classes.f90 @@ -79,6 +79,7 @@ module swiftest_classes procedure :: close => netcdf_close !! Closes an open NetCDF file procedure :: initialize => netcdf_initialize_output !! Initialize a set of parameters used to identify a NetCDF output object procedure :: open => netcdf_open !! Opens a NetCDF file + procedure :: sync => netcdf_sync !! Syncrhonize the disk and memory buffer of the NetCDF file (e.g. commit the frame files stored in memory to disk) end type netcdf_parameters !******************************************************************************************************************************** @@ -200,9 +201,12 @@ module swiftest_classes !! The minimal methods that all systems must have procedure :: dump => io_dump_base !! Dump contents to file procedure :: dump_particle_info => io_dump_particle_info_base !! Dump contents of particle information metadata to file + procedure :: read_in => io_read_in_base !! Read in body initial conditions from a file procedure :: write_frame_netcdf => netcdf_write_frame_base !! I/O routine for writing out a single frame of time-series data for all bodies in the system in NetCDF format + procedure :: read_frame_netcdf => netcdf_read_frame_base !! I/O routine for writing out a single frame of time-series data for all bodies in the system in NetCDF format procedure :: write_particle_info_netcdf => netcdf_write_particle_info_base !! Writes out the particle information metadata to NetCDF file generic :: write_frame => write_frame_netcdf !! Set up generic procedure that will switch between NetCDF or Fortran binary depending on arguments + generic :: read_frame => read_frame_netcdf !! Set up generic procedure that will switch between NetCDF or Fortran binary depending on arguments generic :: write_particle_info => write_particle_info_netcdf end type swiftest_base @@ -236,10 +240,10 @@ module swiftest_classes real(DP), dimension(NDIM) :: L0 = 0.0_DP !! Initial angular momentum of the central body real(DP), dimension(NDIM) :: dL = 0.0_DP !! Change in angular momentum of the central body contains - procedure :: read_in => io_read_in_cb !! I/O routine for reading in central body data - procedure :: read_frame => io_read_frame_cb !! I/O routine for reading out a single frame of time-series data for the central body + procedure :: read_frame_bin => io_read_frame_cb !! I/O routine for reading out a single frame of time-series data for the central body procedure :: write_frame_bin => io_write_frame_cb !! I/O routine for writing out a single frame of time-series data for the central body generic :: write_frame => write_frame_bin !! Write a frame (either binary or NetCDF, using generic procedures) + generic :: read_frame => read_frame_bin !! Write a frame (either binary or NetCDF, using generic procedures) end type swiftest_cb !******************************************************************************************************************************** @@ -283,8 +287,7 @@ module swiftest_classes procedure :: drift => drift_body !! Loop through bodies and call Danby drift routine on heliocentric variables procedure :: v2pv => gr_vh2pv_body !! Converts from velocity to psudeovelocity for GR calculations using symplectic integrators procedure :: pv2v => gr_pv2vh_body !! Converts from psudeovelocity to velocity for GR calculations using symplectic integrators - procedure :: read_in => io_read_in_body !! Read in body initial conditions from a file - procedure :: read_frame => io_read_frame_body !! I/O routine for writing out a single frame of time-series data for the central body + procedure :: read_frame_bin => io_read_frame_body !! I/O routine for writing out a single frame of time-series data for the central body procedure :: write_frame_bin => io_write_frame_body !! I/O routine for writing out a single frame of time-series data for the central body procedure :: accel_obl => obl_acc_body !! Compute the barycentric accelerations of bodies due to the oblateness of the central body procedure :: el2xv => orbel_el2xv_vec !! Convert orbital elements to position and velocity vectors @@ -299,6 +302,7 @@ module swiftest_classes procedure :: rearrange => util_sort_rearrange_body !! Rearranges the order of array elements of body based on an input index array. Used in sorting methods procedure :: spill => util_spill_body !! "Spills" bodies from one object to another depending on the results of a mask (uses the PACK intrinsic) generic :: write_frame => write_frame_bin !! Add the generic write frame for Fortran binary files + generic :: read_frame => read_frame_bin !! Add the generic read frame for Fortran binary files end type swiftest_body !******************************************************************************************************************************** @@ -419,14 +423,21 @@ module swiftest_classes procedure :: discard => discard_system !! Perform a discard step on the system procedure :: conservation_report => io_conservation_report !! Compute energy and momentum and print out the change with time procedure :: dump => io_dump_system !! Dump the state of the system to a file - procedure :: get_old_t_final => io_get_old_t_final_system !! Validates the dump file to check whether the dump file initial conditions duplicate the last frame of the binary output. - procedure :: read_frame => io_read_frame_system !! Read in a frame of input data from file - procedure :: read_particle_info => io_read_particle_info_system !! Read in particle metadata from file - procedure :: write_discard => io_write_discard !! Write out information about discarded test particles - procedure :: write_frame => io_write_frame_system !! Append a frame of output data to file + procedure :: get_old_t_final_bin => io_get_old_t_final_system !! Validates the dump file to check whether the dump file initial conditions duplicate the last frame of the binary output. + procedure :: get_old_t_final_netcdf => netcdf_get_old_t_final_system !! Validates the dump file to check whether the dump file initial conditions duplicate the last frame of the netcdf output. + procedure :: read_frame_bin => io_read_frame_system !! Read in a frame of input data from file + procedure :: write_frame_bin => io_write_frame_system !! Append a frame of output data to file + procedure :: read_frame_netcdf => netcdf_read_frame_system !! Read in a frame of input data from file + procedure :: write_frame_netcdf => netcdf_write_frame_system !! Write a frame of input data from file procedure :: write_hdr_bin => io_write_hdr_system !! Write a header for an output frame in Fortran binary format + !procedure :: read_hdr_bin => io_read_hdr !! Read a header for an output frame in Fortran binary format + procedure :: read_hdr_netcdf => netcdf_read_hdr_system !! Read a header for an output frame in NetCDF format procedure :: write_hdr_netcdf => netcdf_write_hdr_system !! Write a header for an output frame in NetCDF format + procedure :: read_in => io_read_in_system !! Reads the initial conditions for an nbody system + procedure :: read_particle_info => io_read_particle_info_system !! Read in particle metadata from file + procedure :: write_discard => io_write_discard !! Write out information about discarded test particles procedure :: obl_pot => obl_pot_system !! Compute the contribution to the total gravitational potential due solely to the oblateness of the central body + procedure :: finalize => setup_finalize_system !! Runs any finalization subroutines when ending the simulation. procedure :: initialize => setup_initialize_system !! Initialize the system from input files procedure :: init_particle_info => setup_initialize_particle_info_system !! Initialize the system from input files procedure :: step_spin => tides_step_spin_system !! Steps the spins of the massive & central bodies due to tides. @@ -435,6 +446,9 @@ module swiftest_classes 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_hdr => write_hdr_bin, write_hdr_netcdf !! Generic method call for writing headers + generic :: read_hdr => read_hdr_netcdf !! Generic method call for reading headers + generic :: read_frame => read_frame_bin, read_frame_netcdf !! Generic method call for reading a frame of output data + generic :: write_frame => write_frame_bin, write_frame_netcdf !! Generic method call for writing a frame of output data end type swiftest_nbody_system type :: swiftest_encounter @@ -834,17 +848,12 @@ end subroutine io_param_writer_one_QP end interface io_param_writer_one interface - module subroutine io_read_in_body(self, param) - implicit none - class(swiftest_body), intent(inout) :: self !! Swiftest body object - class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters - end subroutine io_read_in_body - module subroutine io_read_in_cb(self, param) + module subroutine io_read_in_base(self,param) implicit none - class(swiftest_cb), intent(inout) :: self !! Swiftest central body object - class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters - end subroutine io_read_in_cb + class(swiftest_base), intent(inout) :: self !! Swiftest base object + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters + end subroutine io_read_in_base module subroutine io_read_in_param(self, param_file_name) implicit none @@ -858,6 +867,12 @@ module subroutine io_read_in_particle_info(self, iu) integer(I4B), intent(in) :: iu !! Open file unit number end subroutine io_read_in_particle_info + module subroutine io_read_in_system(self, param) + implicit none + class(swiftest_nbody_system), intent(inout) :: self + class(swiftest_parameters), intent(inout) :: param + end subroutine io_read_in_system + module function io_read_frame_body(self, iu, param) result(ierr) implicit none class(swiftest_body), intent(inout) :: self !! Swiftest body object @@ -1012,12 +1027,18 @@ module pure subroutine kick_getacch_int_one_tp(rji2, xr, yr, zr, Gmpl, ax, ay, a real(DP), intent(inout) :: ax, ay, az !! Acceleration vector components of test particle end subroutine kick_getacch_int_one_tp - module subroutine netcdf_close(self, param) + module subroutine netcdf_close(self) implicit none class(netcdf_parameters), intent(inout) :: self !! Parameters used to identify a particular NetCDF dataset - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters end subroutine netcdf_close + module function netcdf_get_old_t_final_system(self, param) result(old_t_final) + implicit none + class(swiftest_nbody_system), intent(in) :: self + class(swiftest_parameters), intent(inout) :: param + real(DP) :: old_t_final + end function netcdf_get_old_t_final_system + module subroutine netcdf_initialize_output(self, param) implicit none class(netcdf_parameters), intent(inout) :: self !! Parameters used to for writing a NetCDF dataset to file @@ -1030,6 +1051,34 @@ module subroutine netcdf_open(self, param) class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters end subroutine netcdf_open + module subroutine netcdf_sync(self) + implicit none + class(netcdf_parameters), intent(inout) :: self !! Parameters used to identify a particular NetCDF dataset + end subroutine netcdf_sync + + module function netcdf_read_frame_base(self, iu, param) result(ierr) + implicit none + class(swiftest_base), intent(inout) :: self !! Swiftest base object + class(netcdf_parameters), intent(inout) :: iu !! Parameters used to for writing a NetCDF dataset to file + class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters + integer(I4B) :: ierr !! Error code: returns 0 if the read is successful + end function netcdf_read_frame_base + + module function netcdf_read_frame_system(self, iu, param) result(ierr) + implicit none + class(swiftest_nbody_system), intent(inout) :: self !! Swiftest system object + class(netcdf_parameters), intent(inout) :: iu !! Parameters used to for reading a NetCDF dataset to file + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters + integer(I4B) :: ierr !! Error code: returns 0 if the read is successful + end function netcdf_read_frame_system + + module subroutine netcdf_read_hdr_system(self, iu, param) + implicit none + class(swiftest_nbody_system), intent(inout) :: self !! Swiftest nbody system object + class(netcdf_parameters), intent(inout) :: iu !! Parameters used to for reading a NetCDF dataset to file + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters + end subroutine netcdf_read_hdr_system + module subroutine netcdf_write_frame_base(self, iu, param) implicit none class(swiftest_base), intent(in) :: self !! Swiftest base object @@ -1039,8 +1088,8 @@ end subroutine netcdf_write_frame_base module subroutine netcdf_write_frame_system(self, iu, param) implicit none - class(swiftest_nbody_system), intent(in) :: self !! Swiftest system object - integer(I4B), intent(inout) :: iu !! Parameters used to for writing a NetCDF dataset to file + class(swiftest_nbody_system), intent(inout) :: self !! Swiftest system object + class(netcdf_parameters), intent(inout) :: iu !! Parameters used to for writing a NetCDF dataset to file class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters end subroutine netcdf_write_frame_system @@ -1141,6 +1190,12 @@ module subroutine setup_encounter(self, n) integer(I4B), intent(in) :: n !! Number of encounters to allocate space for end subroutine setup_encounter + module subroutine setup_finalize_system(self, param) + implicit none + class(swiftest_nbody_system), intent(inout) :: self !! Swiftest system object + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters + end subroutine setup_finalize_system + module subroutine setup_initialize_particle_info_system(self, param) implicit none class(swiftest_nbody_system), intent(inout) :: self !! Swiftest nbody system object diff --git a/src/modules/swiftest_globals.f90 b/src/modules/swiftest_globals.f90 index 169cdcd2f..37a1b3a2b 100644 --- a/src/modules/swiftest_globals.f90 +++ b/src/modules/swiftest_globals.f90 @@ -114,9 +114,10 @@ module swiftest_globals !> Standard file names integer(I4B), parameter :: NDUMPFILES = 2 - character(*), dimension(2), parameter :: DUMP_CB_FILE = ['dump_cb1.bin', 'dump_cb2.bin' ] - character(*), dimension(2), parameter :: DUMP_PL_FILE = ['dump_pl1.bin', 'dump_pl2.bin' ] - character(*), dimension(2), parameter :: DUMP_TP_FILE = ['dump_tp1.bin', 'dump_tp2.bin' ] + character(*), dimension(2), parameter :: DUMP_CB_FILE = ['dump_cb1.bin', 'dump_cb2.bin' ] + character(*), dimension(2), parameter :: DUMP_PL_FILE = ['dump_pl1.bin', 'dump_pl2.bin' ] + character(*), dimension(2), parameter :: DUMP_TP_FILE = ['dump_tp1.bin', 'dump_tp2.bin' ] + character(*), dimension(2), parameter :: DUMP_NC_FILE = ['dump_bin1.nc', 'dump_bin2.nc' ] character(*), dimension(2), parameter :: DUMP_PARAM_FILE = ['dump_param1.in', 'dump_param2.in'] !> Default file names that can be changed by the user in the parameters file diff --git a/src/netcdf/netcdf.f90 b/src/netcdf/netcdf.f90 index b26a1ba63..cb3c43198 100644 --- a/src/netcdf/netcdf.f90 +++ b/src/netcdf/netcdf.f90 @@ -19,20 +19,50 @@ subroutine check(status) return end subroutine check - module subroutine netcdf_close(self, param) + module subroutine netcdf_close(self) !! author: Carlisle A. Wishard, Dana Singh, and David A. Minton !! !! Closes a NetCDF file implicit none ! Arguments class(netcdf_parameters), intent(inout) :: self !! Parameters used to identify a particular NetCDF dataset - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters call check( nf90_close(self%ncid) ) return end subroutine netcdf_close + module function netcdf_get_old_t_final_system(self, param) result(old_t_final) + !! author: David A. Minton + !! + !! Validates the dump file to check whether the dump file initial conditions duplicate the last frame of the netcdf output. + !! + implicit none + ! Arguments + class(swiftest_nbody_system), intent(in) :: self + class(swiftest_parameters), intent(inout) :: param + ! Result + real(DP) :: old_t_final + ! Internals + class(swiftest_nbody_system), allocatable :: tmpsys + class(swiftest_parameters), allocatable :: tmpparam + Integer(I4B) :: ierr + + old_t_final = 0.0_DP + allocate(tmpsys, source=self) + allocate(tmpparam, source=param) + ierr = 0 + do + ierr = tmpsys%read_frame(param%nciu, tmpparam) + end do + + if (is_iostat_end(ierr)) then + old_t_final = tmpparam%t + return + end if + + end function netcdf_get_old_t_final_system + module subroutine netcdf_initialize_output(self, param) !! author: Carlisle A. Wishard, Dana Singh, and David A. Minton @@ -44,15 +74,22 @@ module subroutine netcdf_initialize_output(self, param) class(netcdf_parameters), intent(inout) :: self !! Parameters used to identify a particular NetCDF dataset class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters ! Internals - logical :: fileExists - integer(I4B) :: old_mode, nvar, varid, vartype + integer(I4B) :: old_mode, nvar, varid, vartype, old_unit real(DP) :: dfill real(SP) :: sfill + logical :: fileExists + character(len=STRMAX) :: errmsg dfill = ieee_value(dfill, IEEE_QUIET_NAN) sfill = ieee_value(sfill, IEEE_QUIET_NAN) - !! Create the new output file, deleting any previously existing output file of the same name + ! Check if the file exists, and if it does, delete it + inquire(file=param%outfile, exist=fileExists) + if (fileExists) then + open(unit=LUN, file=param%outfile, status="old", err=667, iomsg=errmsg) + close(unit=LUN, status="delete") + end if + call check( nf90_create(param%outfile, NF90_NETCDF4, self%ncid) ) ! Define the NetCDF dimensions with particle name as the record dimension @@ -159,7 +196,14 @@ module subroutine netcdf_initialize_output(self, param) end select end do + ! Take the file out of define mode + call check( nf90_enddef(self%ncid) ) + return + + 667 continue + write(*,*) "Error creating NetCDF output file. " // trim(adjustl(errmsg)) + call util_exit(FAILURE) end subroutine netcdf_initialize_output @@ -256,6 +300,220 @@ module subroutine netcdf_open(self, param) end subroutine netcdf_open + module function netcdf_read_frame_base(self, iu, param) result(ierr) + !! author: Carlisle A. Wishard, Dana Singh, and David A. Minton + !! + !! Read a frame of output of either test particle or massive body data from the binary output file + !! Note: If outputting to orbital elements, but sure that the conversion is done prior to calling this method + implicit none + ! Arguments + class(swiftest_base), intent(inout) :: self !! Swiftest base object + class(netcdf_parameters), intent(inout) :: iu !! Parameters used to for writing a NetCDF dataset to file + class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters + ! Internals + integer(I4B) :: i, j, tslot, strlen, idslot, idmax, n + integer(I4B), dimension(:), allocatable :: ind + character(len=:), allocatable :: charstring + integer(I4B) :: ierr !! Error code: returns 0 if the read is successful + real(DP), dimension(:), allocatable :: real_temp + logical, dimension(:), allocatable :: validmask, tpmask + + tslot = int(param%ioutput, kind=I4B) + 1 + + select type(self) + class is (swiftest_body) + if (self%nbody == 0) return + call check( nf90_inquire_dimension(iu%ncid, iu%id_dimid, len=idmax) ) + allocate(ind(idmax)) + allocate(real_temp(idmax)) + allocate(validmask(idmax)) + allocate(tpmask(idmax)) + ind(:) = [(i, i = 1, idmax)] + + ! First filter out only the id slots that contain valid bodies + if (param%in_form == XV) then + call check( nf90_get_var(iu%ncid, iu%xhx_varid, real_temp(:)) ) + else + call check( nf90_get_var(iu%ncid, iu%a_varid, real_temp(:)) ) + end if + validmask(:) = real_temp(:) == real_temp(:) + + ! Filter out the central body, which is always in id dimension array position 1 + validmask(1) = .false. + + ! Next, filter only bodies that don't have mass (test particles) + call check( nf90_get_var(iu%ncid, iu%Gmass_varid, real_temp(:)) ) + tpmask(:) = real_temp(:) /= real_temp(:) + + select type(self) + class is (swiftest_pl) + ind(:) = pack(ind(:), (.not.tpmask(:) .and. validmask(:))) + n = count(.not.tpmask(:)) + class is (swiftest_tp) + ind(:) = pack(ind(:), (tpmask(:) .and. validmask(:))) + n = count(tpmask(:)) + end select + + do i = j, n + idslot = ind(j) + + if (param%in_form == XV) then + call check( nf90_get_var(iu%ncid, iu%xhx_varid, self%xh(1, j), start=[idslot, tslot]) ) + call check( nf90_get_var(iu%ncid, iu%xhy_varid, self%xh(2, j), start=[idslot, tslot]) ) + call check( nf90_get_var(iu%ncid, iu%xhz_varid, self%xh(3, j), start=[idslot, tslot]) ) + call check( nf90_get_var(iu%ncid, iu%vhx_varid, self%vh(1, j), start=[idslot, tslot]) ) + call check( nf90_get_var(iu%ncid, iu%vhy_varid, self%vh(2, j), start=[idslot, tslot]) ) + call check( nf90_get_var(iu%ncid, iu%vhz_varid, self%vh(3, j), start=[idslot, tslot]) ) + else if (param%in_form == EL) then + call check( nf90_get_var(iu%ncid, iu%a_varid, self%a(j), start=[idslot, tslot]) ) + call check( nf90_get_var(iu%ncid, iu%e_varid, self%e(j), start=[idslot, tslot]) ) + call check( nf90_get_var(iu%ncid, iu%inc_varid, self%inc(j), start=[idslot, tslot]) ) + call check( nf90_get_var(iu%ncid, iu%capom_varid, self%capom(j), start=[idslot, tslot]) ) + call check( nf90_get_var(iu%ncid, iu%omega_varid, self%omega(j), start=[idslot, tslot]) ) + call check( nf90_get_var(iu%ncid, iu%capm_varid, self%capm(j), start=[idslot, tslot]) ) + end if + + select type(self) + class is (swiftest_pl) ! Additional output if the passed polymorphic object is a massive body + call check( nf90_get_var(iu%ncid, iu%Gmass_varid, self%Gmass(j), start=[idslot, tslot]) ) + if (param%lrhill_present) then + call check( nf90_get_var(iu%ncid, iu%rhill_varid, self%rhill(j), start=[idslot, tslot]) ) + end if + if (param%lclose) then + call check( nf90_get_var(iu%ncid, iu%radius_varid, self%radius(j), start=[idslot, tslot]) ) + end if + if (param%lrotation) then + call check( nf90_get_var(iu%ncid, iu%Ip1_varid, self%Ip(1, j), start=[idslot, tslot]) ) + call check( nf90_get_var(iu%ncid, iu%Ip2_varid, self%Ip(2, j), start=[idslot, tslot]) ) + call check( nf90_get_var(iu%ncid, iu%Ip3_varid, self%Ip(3, j), start=[idslot, tslot]) ) + call check( nf90_get_var(iu%ncid, iu%rotx_varid, self%rot(1, j), start=[idslot, tslot]) ) + call check( nf90_get_var(iu%ncid, iu%roty_varid, self%rot(2, j), start=[idslot, tslot]) ) + call check( nf90_get_var(iu%ncid, iu%rotz_varid, self%rot(3, j), start=[idslot, tslot]) ) + end if + if (param%ltides) then + call check( nf90_get_var(iu%ncid, iu%k2_varid, self%k2(j), start=[idslot, tslot]) ) + call check( nf90_get_var(iu%ncid, iu%Q_varid, self%Q(j), start=[idslot, tslot]) ) + end if + + end select + end do + + class is (swiftest_cb) + + idslot = 1 + call check( nf90_get_var(iu%ncid, iu%id_varid, self%id, start=[idslot]) ) + + call check( nf90_get_var(iu%ncid, iu%Gmass_varid, self%Gmass, start=[idslot, tslot]) ) + call check( nf90_get_var(iu%ncid, iu%radius_varid, self%radius, start=[idslot, tslot]) ) + if (param%lrotation) then + call check( nf90_get_var(iu%ncid, iu%Ip1_varid, self%Ip(1), start=[idslot, tslot]) ) + call check( nf90_get_var(iu%ncid, iu%Ip2_varid, self%Ip(2), start=[idslot, tslot]) ) + call check( nf90_get_var(iu%ncid, iu%Ip3_varid, self%Ip(3), start=[idslot, tslot]) ) + call check( nf90_get_var(iu%ncid, iu%rotx_varid, self%rot(1), start=[idslot, tslot]) ) + call check( nf90_get_var(iu%ncid, iu%roty_varid, self%rot(2), start=[idslot, tslot]) ) + call check( nf90_get_var(iu%ncid, iu%rotz_varid, self%rot(3), start=[idslot, tslot]) ) + end if + if (param%ltides) then + call check( nf90_get_var(iu%ncid, iu%k2_varid, self%k2, start=[idslot, tslot]) ) + call check( nf90_get_var(iu%ncid, iu%Q_varid, self%Q, start=[idslot, tslot]) ) + end if + + + end select + + !call self%read_particle_info(iu) ! THIS NEEDS TO BE IMPLEMENTED + + return + + end function netcdf_read_frame_base + + + module function netcdf_read_frame_system(self, iu, param) result(ierr) + !! author: The Purdue Swiftest Team - David A. Minton, Carlisle A. Wishard, Jennifer L.L. Pouplin, and Jacob R. Elliott + !! + !! Read a frame (header plus records for each massive body and active test particle) from a output binary file + implicit none + ! Arguments + class(swiftest_nbody_system), intent(inout) :: self !! Swiftest system object + class(netcdf_parameters), intent(inout) :: iu !! Parameters used to identify a particular NetCDF dataset + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters + integer(I4B) :: ierr !! Error code: returns 0 if the read is successful + + call iu%open(param) + + call self%read_hdr(iu, param) + call self%cb%read_in(param) + call self%pl%read_in(param) + call self%tp%read_in(param) + call iu%close() + + ierr = 0 + return + + 667 continue + write(*,*) "Error reading system frame in netcdf_read_frame_system" + + end function netcdf_read_frame_system + + module subroutine netcdf_read_hdr_system(self, iu, param) + !! author: David A. Minton + !! + !! Reads header information (variables that change with time, but not particle id). + !! This subroutine significantly improves the output over the original binary file, allowing us to track energy, momentum, and other quantities that + !! previously were handled as separate output files. + implicit none + ! Arguments + class(swiftest_nbody_system), intent(inout) :: self !! Swiftest nbody system object + class(netcdf_parameters), intent(inout) :: iu !! Parameters used to for writing a NetCDF dataset to file + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters + ! Internals + integer(I4B) :: tslot, old_mode + + tslot = int(param%ioutput, kind=I4B) + 1 + + call check( nf90_open(param%outfile, nf90_nowrite, iu%ncid) ) + call check( nf90_set_fill(iu%ncid, nf90_nofill, old_mode) ) + + call check( nf90_get_var(iu%ncid, iu%time_varid, param%t, start=[tslot]) ) + call check( nf90_get_var(iu%ncid, iu%npl_varid, self%pl%nbody, start=[tslot]) ) + call check( nf90_get_var(iu%ncid, iu%ntp_varid, self%tp%nbody, start=[tslot]) ) + + if (param%lenergy) then + call check( nf90_get_var(iu%ncid, iu%KE_orb_varid, self%ke_orbit, start=[tslot]) ) + call check( nf90_get_var(iu%ncid, iu%KE_spin_varid, self%ke_spin, start=[tslot]) ) + call check( nf90_get_var(iu%ncid, iu%PE_varid, self%pe, start=[tslot]) ) + call check( nf90_get_var(iu%ncid, iu%L_orbx_varid, self%Lorbit(1), start=[tslot]) ) + call check( nf90_get_var(iu%ncid, iu%L_orby_varid, self%Lorbit(2), start=[tslot]) ) + call check( nf90_get_var(iu%ncid, iu%L_orbz_varid, self%Lorbit(3), start=[tslot]) ) + call check( nf90_get_var(iu%ncid, iu%L_spinx_varid, self%Lspin(1), start=[tslot]) ) + call check( nf90_get_var(iu%ncid, iu%L_spiny_varid, self%Lspin(2), start=[tslot]) ) + call check( nf90_get_var(iu%ncid, iu%L_spinz_varid, self%Lspin(3), start=[tslot]) ) + call check( nf90_get_var(iu%ncid, iu%L_escapex_varid, param%Lescape(1), start=[tslot]) ) + call check( nf90_get_var(iu%ncid, iu%L_escapey_varid, param%Lescape(2), start=[tslot]) ) + call check( nf90_get_var(iu%ncid, iu%L_escapez_varid, param%Lescape(3), start=[tslot]) ) + call check( nf90_get_var(iu%ncid, iu%Ecollisions_varid, param%Ecollisions, start=[tslot]) ) + call check( nf90_get_var(iu%ncid, iu%Euntracked_varid, param%Euntracked, start=[tslot]) ) + call check( nf90_get_var(iu%ncid, iu%GMescape_varid, param%GMescape, start=[tslot]) ) + end if + + return + end subroutine netcdf_read_hdr_system + + + module subroutine netcdf_sync(self) + !! author: David A. Minton + !! + !! Syncrhonize the disk and memory buffer of the NetCDF file (e.g. commit the frame files stored in memory to disk) + !! + implicit none + ! Arguments + class(netcdf_parameters), intent(inout) :: self !! Parameters used to identify a particular NetCDF dataset + + call check( nf90_sync(self%ncid) ) + + return + end subroutine netcdf_sync + module subroutine netcdf_write_frame_base(self, iu, param) !! author: Carlisle A. Wishard, Dana Singh, and David A. Minton !! @@ -356,7 +614,24 @@ module subroutine netcdf_write_frame_base(self, iu, param) return end subroutine netcdf_write_frame_base - + module subroutine netcdf_write_frame_system(self, iu, param) + !! author: The Purdue Swiftest Team - David A. Minton, Carlisle A. Wishard, Jennifer L.L. Pouplin, and Jacob R. Elliott + !! + !! Write a frame (header plus records for each massive body and active test particle) to a output binary file + implicit none + ! Arguments + class(swiftest_nbody_system), intent(inout) :: self !! Swiftest system object + class(netcdf_parameters), intent(inout) :: iu !! Parameters used to identify a particular NetCDF dataset + class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters + + call self%write_hdr(iu, param) + call self%cb%write_frame(iu, param) + call self%pl%write_frame(iu, param) + call self%tp%write_frame(iu, param) + + return + end subroutine netcdf_write_frame_system + module subroutine netcdf_write_particle_info_base(self, iu) !! author: Carlisle A. Wishard, Dana Singh, and David A. Minton !! diff --git a/src/setup/setup.f90 b/src/setup/setup.f90 index 7be3eb2c4..bd8930ae3 100644 --- a/src/setup/setup.f90 +++ b/src/setup/setup.f90 @@ -126,6 +126,27 @@ module subroutine setup_encounter(self, n) end subroutine setup_encounter + module subroutine setup_finalize_system(self, param) + !! author: David A. Minton + !! + !! Runs any finalization subroutines when ending the simulation. + !! + implicit none + ! Arguments + class(swiftest_nbody_system), intent(inout) :: self !! Swiftest system object + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters + integer(I4B) :: ierr + + associate(system => self) + if ((param%out_type == NETCDF_FLOAT_TYPE) .or. (param%out_type == NETCDF_DOUBLE_TYPE)) then + call param%nciu%close() + end if + end associate + + return + end subroutine setup_finalize_system + + module subroutine setup_initialize_particle_info_system(self, param) !! author: David A. Minton !! @@ -163,28 +184,31 @@ module subroutine setup_initialize_system(self, param) ! Arguments class(swiftest_nbody_system), intent(inout) :: self !! Swiftest system object class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters - - call self%cb%read_in(param) - call self%pl%read_in(param) - call self%tp%read_in(param) - call self%validate_ids(param) - call self%set_msys() - call self%pl%set_mu(self%cb) - call self%tp%set_mu(self%cb) - if (param%in_form == EL) then - call self%pl%el2xv(self%cb) - call self%tp%el2xv(self%cb) - end if - call self%pl%flatten(param) - if (.not.param%lrhill_present) call self%pl%set_rhill(self%cb) - self%pl%lfirst = param%lfirstkick - self%tp%lfirst = param%lfirstkick - - if (param%lrestart) then - call self%read_particle_info(param) - else - call self%init_particle_info(param) - end if + integer(I4B) :: ierr + + associate(system => self, cb => self%cb, pl => self%pl, tp => self%tp) + + call system%read_in(param) + call system%validate_ids(param) + call system%set_msys() + call pl%set_mu(cb) + call tp%set_mu(cb) + if (param%in_form == EL) then + call pl%el2xv(cb) + call tp%el2xv(cb) + end if + call pl%flatten(param) + if (.not.param%lrhill_present) call pl%set_rhill(cb) + pl%lfirst = param%lfirstkick + tp%lfirst = param%lfirstkick + + if (param%lrestart) then + call system%read_particle_info(param) + else + call system%init_particle_info(param) + end if + end associate + return end subroutine setup_initialize_system diff --git a/src/symba/symba_io.f90 b/src/symba/symba_io.f90 index bea4e9de3..632d4e64b 100644 --- a/src/symba/symba_io.f90 +++ b/src/symba/symba_io.f90 @@ -176,9 +176,7 @@ module subroutine symba_io_write_discard(self, param) ! Record the discarded body metadata information to file if ((param%out_type == NETCDF_FLOAT_TYPE) .or. (param%out_type == NETCDF_DOUBLE_TYPE)) then - call param%nciu%open(param) call pl_discards%write_particle_info(param%nciu) - call param%nciu%close(param) end if if (param%discard_out == "") return