From d4e17769bb218788853392a4050ec5e49b6225c7 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Mon, 30 Aug 2021 14:17:37 -0400 Subject: [PATCH] Convert angles to degrees before putting them into NetCDF. Also write out final state of variables on the mass loss exception. --- examples/symba_mars_disk/testnetcdf.ipynb | 6957 ++++++++++++++++----- src/io/io.f90 | 4 +- src/netcdf/netcdf.f90 | 8 +- 3 files changed, 5368 insertions(+), 1601 deletions(-) diff --git a/examples/symba_mars_disk/testnetcdf.ipynb b/examples/symba_mars_disk/testnetcdf.ipynb index 50ac91b3b..be4bf6a11 100644 --- a/examples/symba_mars_disk/testnetcdf.ipynb +++ b/examples/symba_mars_disk/testnetcdf.ipynb @@ -2,7 +2,7 @@ "cells": [ { "cell_type": "code", - "execution_count": 10, + "execution_count": 1, "metadata": {}, "outputs": [ { @@ -11,7 +11,7 @@ "'/home/daminton/git/swiftest/examples/symba_mars_disk'" ] }, - "execution_count": 10, + "execution_count": 1, "metadata": {}, "output_type": "execute_result" } @@ -26,7 +26,7 @@ }, { "cell_type": "code", - "execution_count": 11, + "execution_count": 2, "metadata": {}, "outputs": [ { @@ -36,7 +36,7 @@ "Reading Swiftest file param.in\n", "\n", "Creating Dataset\n", - "Successfully converted 2 output frames.\n", + "Successfully converted 6 output frames.\n", "Swiftest simulation data stored as xarray DataSet .ds\n" ] } @@ -48,7 +48,7 @@ }, { "cell_type": "code", - "execution_count": 12, + "execution_count": 3, "metadata": {}, "outputs": [ { @@ -405,25 +405,28 @@ " stroke: currentColor;\n", " fill: currentColor;\n", "}\n", - "
<xarray.DataArray 'npl' (time: 2)>\n",
-       "array([1500, 1518], dtype=int32)\n",
+       "
<xarray.DataArray 'status' (id: 1521)>\n",
+       "array(['ACTIVE', 'ACTIVE', 'ACTIVE', ..., 'ACTIVE', 'ACTIVE', 'ACTIVE'],\n",
+       "      dtype='<U17')\n",
        "Coordinates:\n",
-       "  * time     (time) float64 0.0 6e+03
"
+       "  * id       (id) int32 0 1 2 3 4 5 6 7 ... 1514 1515 1516 1517 1518 1519 1520
" ], "text/plain": [ - "\n", - "array([1500, 1518], dtype=int32)\n", + "\n", + "array(['ACTIVE', 'ACTIVE', 'ACTIVE', ..., 'ACTIVE', 'ACTIVE', 'ACTIVE'],\n", + " dtype='
<xarray.DataArray 'xhx' (time: 2)>\n",
-       "array([      0.      , 2742814.426863])\n",
+       "
<xarray.DataArray 'status' ()>\n",
+       "array('Supercatastrophic', dtype='<U17')\n",
        "Coordinates:\n",
-       "  * time     (time) float64 0.0 6e+03\n",
-       "    id       int32 1510
"
+       "    id       int32 726
" ], "text/plain": [ - "\n", - "array([ 0. , 2742814.426863])\n", + "\n", + "array('Supercatastrophic', dtype='\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "
<xarray.DataArray 'discard_time' (id: 1521)>\n",
+       "array([ 0.000000e+000,  6.013470e-154,  6.013470e-154, ..., -1.797693e+308,\n",
+       "       -1.797693e+308, -1.797693e+308])\n",
+       "Coordinates:\n",
+       "  * id       (id) int32 0 1 2 3 4 5 6 7 ... 1514 1515 1516 1517 1518 1519 1520
" + ], + "text/plain": [ + "\n", + "array([ 0.000000e+000, 6.013470e-154, 6.013470e-154, ..., -1.797693e+308,\n", + " -1.797693e+308, -1.797693e+308])\n", + "Coordinates:\n", + " * id (id) int32 0 1 2 3 4 5 6 7 ... 1514 1515 1516 1517 1518 1519 1520" + ] + }, + "execution_count": 6, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "sim.ds['discard_time']" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [], + "source": [ + "dslost = sim.ds.where(sim.ds['status'] != \"ACTIVE\", drop=True)" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "
<xarray.DataArray 'discard_time' (id: 2)>\n",
+       "array([2400., 2400.])\n",
+       "Coordinates:\n",
+       "  * id       (id) int32 230 726
" + ], + "text/plain": [ + "\n", + "array([2400., 2400.])\n", + "Coordinates:\n", + " * id (id) int32 230 726" + ] + }, + "execution_count": 8, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "dslost['discard_time']" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": {}, + "outputs": [], + "source": [ + "lastloss = dslost.where(dslost['discard_time'] == 98400.0, drop=True)" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": {}, + "outputs": [], + "source": [ + "dsactive = sim.ds.where(sim.ds['status'] == \"ACTIVE\", drop=True)" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "
<xarray.Dataset>\n",
+       "Dimensions:          (id: 1519, time: 6)\n",
+       "Coordinates:\n",
+       "  * time             (time) float64 0.0 600.0 1.2e+03 1.8e+03 2.4e+03 3e+03\n",
+       "  * id               (id) int32 0 1 2 3 4 5 6 ... 1515 1516 1517 1518 1519 1520\n",
+       "Data variables: (12/57)\n",
+       "    npl              (time, id) float64 1.5e+03 1.5e+03 ... 1.518e+03 1.518e+03\n",
+       "    ntp              (time, id) float64 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0\n",
+       "    name             (id) object 'Mars' 'Body2' ... 'Newbody0001520'\n",
+       "    particle_type    (id) object 'Central Body' ... 'Massive Body'\n",
+       "    xhx              (time, id) float64 0.0 -2.358e+06 ... -3.489e+06 -3.476e+06\n",
+       "    xhy              (time, id) float64 0.0 8.604e+06 ... -8.532e+06 -8.571e+06\n",
+       "    ...               ...\n",
+       "    discard_xhy      (id) float64 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0\n",
+       "    discard_xhz      (id) float64 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0\n",
+       "    discard_vhx      (id) float64 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0\n",
+       "    discard_vhy      (id) float64 0.0 2.122e-314 2.122e-314 ... 0.0 0.0 0.0\n",
+       "    discard_vhz      (id) float64 0.0 3.024e-153 3.024e-153 ... 0.0 0.0 0.0\n",
+       "    discard_body_id  (id) float64 -2.147e+09 -2.147e+09 ... -2.147e+09
" + ], + "text/plain": [ + "\n", + "Dimensions: (id: 1519, time: 6)\n", + "Coordinates:\n", + " * time (time) float64 0.0 600.0 1.2e+03 1.8e+03 2.4e+03 3e+03\n", + " * id (id) int32 0 1 2 3 4 5 6 ... 1515 1516 1517 1518 1519 1520\n", + "Data variables: (12/57)\n", + " npl (time, id) float64 1.5e+03 1.5e+03 ... 1.518e+03 1.518e+03\n", + " ntp (time, id) float64 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0\n", + " name (id) object 'Mars' 'Body2' ... 'Newbody0001520'\n", + " particle_type (id) object 'Central Body' ... 'Massive Body'\n", + " xhx (time, id) float64 0.0 -2.358e+06 ... -3.489e+06 -3.476e+06\n", + " xhy (time, id) float64 0.0 8.604e+06 ... -8.532e+06 -8.571e+06\n", + " ... ...\n", + " discard_xhy (id) float64 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0\n", + " discard_xhz (id) float64 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0\n", + " discard_vhx (id) float64 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0\n", + " discard_vhy (id) float64 0.0 2.122e-314 2.122e-314 ... 0.0 0.0 0.0\n", + " discard_vhz (id) float64 0.0 3.024e-153 3.024e-153 ... 0.0 0.0 0.0\n", + " discard_body_id (id) float64 -2.147e+09 -2.147e+09 ... -2.147e+09" + ] + }, + "execution_count": 11, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "dsactive" + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "metadata": {}, + "outputs": [], + "source": [ + "lastadd = sim.ds.where(sim.ds['origin_time'] == 2400.0, drop=True)" + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "
<xarray.Dataset>\n",
+       "Dimensions:          (id: 20, time: 6)\n",
+       "Coordinates:\n",
+       "  * time             (time) float64 0.0 600.0 1.2e+03 1.8e+03 2.4e+03 3e+03\n",
+       "  * id               (id) int32 1501 1502 1503 1504 1505 ... 1517 1518 1519 1520\n",
+       "Data variables: (12/57)\n",
+       "    npl              (time, id) float64 1.5e+03 1.5e+03 ... 1.518e+03 1.518e+03\n",
+       "    ntp              (time, id) float64 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0\n",
+       "    name             (id) object 'Newbody0001501' ... 'Newbody0001520'\n",
+       "    particle_type    (id) object 'Massive Body' ... 'Massive Body'\n",
+       "    xhx              (time, id) float64 0.0 0.0 0.0 ... -3.489e+06 -3.476e+06\n",
+       "    xhy              (time, id) float64 0.0 0.0 0.0 ... -8.532e+06 -8.571e+06\n",
+       "    ...               ...\n",
+       "    discard_xhy      (id) float64 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0\n",
+       "    discard_xhz      (id) float64 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0\n",
+       "    discard_vhx      (id) float64 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0\n",
+       "    discard_vhy      (id) float64 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0\n",
+       "    discard_vhz      (id) float64 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0\n",
+       "    discard_body_id  (id) float64 -2.147e+09 -2.147e+09 ... -2.147e+09
" + ], + "text/plain": [ + "\n", + "Dimensions: (id: 20, time: 6)\n", + "Coordinates:\n", + " * time (time) float64 0.0 600.0 1.2e+03 1.8e+03 2.4e+03 3e+03\n", + " * id (id) int32 1501 1502 1503 1504 1505 ... 1517 1518 1519 1520\n", + "Data variables: (12/57)\n", + " npl (time, id) float64 1.5e+03 1.5e+03 ... 1.518e+03 1.518e+03\n", + " ntp (time, id) float64 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0\n", + " name (id) object 'Newbody0001501' ... 'Newbody0001520'\n", + " particle_type (id) object 'Massive Body' ... 'Massive Body'\n", + " xhx (time, id) float64 0.0 0.0 0.0 ... -3.489e+06 -3.476e+06\n", + " xhy (time, id) float64 0.0 0.0 0.0 ... -8.532e+06 -8.571e+06\n", + " ... ...\n", + " discard_xhy (id) float64 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0\n", + " discard_xhz (id) float64 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0\n", + " discard_vhx (id) float64 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0\n", + " discard_vhy (id) float64 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0\n", + " discard_vhz (id) float64 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0\n", + " discard_body_id (id) float64 -2.147e+09 -2.147e+09 ... -2.147e+09" + ] + }, + "execution_count": 14, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "lastadd" + ] + }, + { + "cell_type": "code", + "execution_count": 15, + "metadata": {}, + "outputs": [], + "source": [ + "lastframe = dsactive.isel(time=-1, drop=True)\n", + "firstframe = sim.ds.isel(time=0, drop=True)\n", + "firstframe = firstframe.where(firstframe['a'] < 1e20, drop=True)\n", + "midframe = sim.ds.isel(time=-2, drop=True)\n", + "midframe = midframe.where(midframe['a'] < 1e20, drop=True)" + ] + }, + { + "cell_type": "code", + "execution_count": 16, + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "
<xarray.DataArray 'id' (id: 1519)>\n",
+       "array([   0,    1,    2, ..., 1518, 1519, 1520], dtype=int32)\n",
+       "Coordinates:\n",
+       "  * id       (id) int32 0 1 2 3 4 5 6 7 ... 1514 1515 1516 1517 1518 1519 1520
" + ], + "text/plain": [ + "\n", + "array([ 0, 1, 2, ..., 1518, 1519, 1520], dtype=int32)\n", + "Coordinates:\n", + " * id (id) int32 0 1 2 3 4 5 6 7 ... 1514 1515 1516 1517 1518 1519 1520" + ] + }, + "execution_count": 16, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "lastframe.id" + ] + }, + { + "cell_type": "code", + "execution_count": 17, + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "
<xarray.DataArray 'id' (id: 0)>\n",
+       "array([], dtype=int32)\n",
+       "Coordinates:\n",
+       "  * id       (id) int32
" + ], + "text/plain": [ + "\n", + "array([], dtype=int32)\n", + "Coordinates:\n", + " * id (id) int32 " + ] + }, + "execution_count": 17, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "lastloss.id" + ] + }, + { + "cell_type": "code", + "execution_count": 18, + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "
<xarray.DataArray 'Gmass' ()>\n",
+       "array(4.28396188e+13)
" + ], + "text/plain": [ + "\n", + "array(4.28396188e+13)" + ] + }, + "execution_count": 18, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "lastframe['Gmass'].sum()" + ] + }, + { + "cell_type": "code", + "execution_count": 19, + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "
<xarray.DataArray 'Gmass' ()>\n",
+       "array(4.28396188e+13)
" + ], + "text/plain": [ + "\n", + "array(4.28396188e+13)" + ] + }, + "execution_count": 19, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "firstframe['Gmass'].sum()" + ] + }, + { + "cell_type": "code", + "execution_count": 20, + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "
<xarray.DataArray 'Gmass' ()>\n",
+       "array(4.28396188e+13)
" + ], + "text/plain": [ + "\n", + "array(4.28396188e+13)" + ] + }, + "execution_count": 20, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "midframe['Gmass'].sum()" + ] + }, + { + "cell_type": "code", + "execution_count": 22, + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "
<xarray.Dataset>\n",
+       "Dimensions:          ()\n",
+       "Coordinates:\n",
+       "    id               int32 2\n",
+       "Data variables: (12/57)\n",
+       "    npl              float64 1.518e+03\n",
+       "    ntp              float64 0.0\n",
+       "    name             object 'Body3'\n",
+       "    particle_type    object 'Massive Body'\n",
+       "    xhx              float64 -3.025e+06\n",
+       "    xhy              float64 -1.02e+07\n",
+       "    ...               ...\n",
+       "    discard_xhy      float64 0.0\n",
+       "    discard_xhz      float64 0.0\n",
+       "    discard_vhx      float64 0.0\n",
+       "    discard_vhy      float64 2.122e-314\n",
+       "    discard_vhz      float64 3.024e-153\n",
+       "    discard_body_id  float64 -2.147e+09
" + ], + "text/plain": [ + "\n", + "Dimensions: ()\n", + "Coordinates:\n", + " id int32 2\n", + "Data variables: (12/57)\n", + " npl float64 1.518e+03\n", + " ntp float64 0.0\n", + " name object 'Body3'\n", + " particle_type object 'Massive Body'\n", + " xhx float64 -3.025e+06\n", + " xhy float64 -1.02e+07\n", + " ... ...\n", + " discard_xhy float64 0.0\n", + " discard_xhz float64 0.0\n", + " discard_vhx float64 0.0\n", + " discard_vhy float64 2.122e-314\n", + " discard_vhz float64 3.024e-153\n", + " discard_body_id float64 -2.147e+09" + ] + }, + "execution_count": 22, + "metadata": {}, + "output_type": "execute_result" } ], "source": [ - "for n in sim.ds['radius']:\n", - " print(n.values)" + "lastframe.sel(id=2)" ] }, { diff --git a/src/io/io.f90 b/src/io/io.f90 index a7e1d5d9f..637d1705f 100644 --- a/src/io/io.f90 +++ b/src/io/io.f90 @@ -69,6 +69,9 @@ module subroutine io_conservation_report(self, param, lterminal) Merror = (GMtot_now - param%GMtot_orig) / param%GMtot_orig if (Merror < -10 * epsilon(Merror)) then write(*,*) 'Mass loss! Halting!' + call param%nciu%open(param) + call pl%write_frame(param%nciu, param) + call param%nciu%close(param) call util_exit(FAILURE) end if write(*, EGYTERMFMT) Lerror, Ecoll_error, Etotal_error, Merror @@ -1919,7 +1922,6 @@ module subroutine io_write_frame_system(self, param) call cb%write_frame(iu, param) call pl%write_frame(iu, param) call tp%write_frame(iu, param) - close(iu, err = 667, iomsg = errmsg) else if ((param%out_type == NETCDF_FLOAT_TYPE) .or. (param%out_type == NETCDF_DOUBLE_TYPE)) then call cb%write_frame(param%nciu, param) diff --git a/src/netcdf/netcdf.f90 b/src/netcdf/netcdf.f90 index 6672ee7d0..64b30a44a 100644 --- a/src/netcdf/netcdf.f90 +++ b/src/netcdf/netcdf.f90 @@ -282,10 +282,10 @@ module subroutine netcdf_write_frame_base(self, iu, param) if ((param%out_form == EL) .or. (param%out_form == XVEL)) then call check( nf90_put_var(iu%ncid, iu%a_varid, self%a(j), start=[idslot, tslot]) ) call check( nf90_put_var(iu%ncid, iu%e_varid, self%e(j), start=[idslot, tslot]) ) - call check( nf90_put_var(iu%ncid, iu%inc_varid, self%inc(j), start=[idslot, tslot]) ) - call check( nf90_put_var(iu%ncid, iu%capom_varid, self%capom(j), start=[idslot, tslot]) ) - call check( nf90_put_var(iu%ncid, iu%omega_varid, self%omega(j), start=[idslot, tslot]) ) - call check( nf90_put_var(iu%ncid, iu%capm_varid, self%capm(j), start=[idslot, tslot]) ) + call check( nf90_put_var(iu%ncid, iu%inc_varid, self%inc(j) * RAD2DEG, start=[idslot, tslot]) ) + call check( nf90_put_var(iu%ncid, iu%capom_varid, self%capom(j) * RAD2DEG, start=[idslot, tslot]) ) + call check( nf90_put_var(iu%ncid, iu%omega_varid, self%omega(j) * RAD2DEG, start=[idslot, tslot]) ) + call check( nf90_put_var(iu%ncid, iu%capm_varid, self%capm(j) * RAD2DEG, start=[idslot, tslot]) ) end if select type(self)