From b92eb1f20f5381795efec9751ea120272772e9bd Mon Sep 17 00:00:00 2001 From: David A Minton Date: Sat, 28 Aug 2021 02:14:14 -0400 Subject: [PATCH] Fixed bugs in info tracking. Code seems to work now. --- examples/symba_energy_momentum/Untitled.ipynb | 470 ++ .../param.disruption_headon.in | 6 +- .../param.disruption_off_axis.in | 8 +- .../param.supercatastrophic_headon.in | 6 +- examples/symba_mars_disk/mars.in | 6016 ++++++++--------- examples/symba_mars_disk/testnetcdf.ipynb | 433 +- src/netcdf/netcdf.f90 | 2 +- src/symba/symba_util.f90 | 10 - src/util/util_append.f90 | 36 +- src/util/util_fill.f90 | 28 +- src/util/util_sort.f90 | 12 +- 11 files changed, 3969 insertions(+), 3058 deletions(-) create mode 100644 examples/symba_energy_momentum/Untitled.ipynb diff --git a/examples/symba_energy_momentum/Untitled.ipynb b/examples/symba_energy_momentum/Untitled.ipynb new file mode 100644 index 000000000..3b730328e --- /dev/null +++ b/examples/symba_energy_momentum/Untitled.ipynb @@ -0,0 +1,470 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "'/home/daminton/git/swiftest/examples/symba_energy_momentum'" + ] + }, + "execution_count": 1, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "import swiftest\n", + "import os\n", + "import xarray as xr\n", + "import numpy as np\n", + "os.getcwd()" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Reading Swiftest file param.disruption_off_axis.in\n", + "\n", + "Creating Dataset\n", + "Successfully converted 102 output frames.\n", + "Swiftest simulation data stored as xarray DataSet .ds\n" + ] + } + ], + "source": [ + "sim = swiftest.Simulation(param_file=\"param.disruption_off_axis.in\")\n", + "sim.bin2xr()" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "
<xarray.DataArray 'origin_type' (id: 15)>\n",
+       "array(['Initial conditions', 'Initial conditions', 'Initial conditions',\n",
+       "       'Disruption', 'Disruption', 'Disruption', 'Disruption',\n",
+       "       'Disruption', 'Disruption', 'Disruption', 'Disruption',\n",
+       "       'Disruption', 'Disruption', 'Disruption', 'Disruption'],\n",
+       "      dtype='<U18')\n",
+       "Coordinates:\n",
+       "  * id       (id) int32 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14
" + ], + "text/plain": [ + "\n", + "array(['Initial conditions', 'Initial conditions', 'Initial conditions',\n", + " 'Disruption', 'Disruption', 'Disruption', 'Disruption',\n", + " 'Disruption', 'Disruption', 'Disruption', 'Disruption',\n", + " 'Disruption', 'Disruption', 'Disruption', 'Disruption'],\n", + " dtype='
<xarray.DataArray 'particle_type' (id: 1521)>\n",
-       "array([b'Central Body', b'Massive Body', b'Massive Body', ..., b'Massive Body',\n",
-       "       b'Massive Body', b'Massive Body'], dtype='|S32')\n",
+       "
<xarray.DataArray 'origin_type' (id: 1521)>\n",
+       "array(['Initial conditions', 'Body                    Initial',\n",
+       "       'Body                    Initial', ..., '', '', ''], dtype='<U32')\n",
        "Coordinates:\n",
-       "  * id       (id) int32 0 1 2 3 4 5 6 7 ... 1514 1515 1516 1517 1518 1519 1520
"
+       "  * id       (id) int32 0 1 2 3 4 5 6 7 ... 1514 1515 1516 1517 1518 1519 1520
" ], "text/plain": [ - "\n", - "array([b'Central Body', b'Massive Body', b'Massive Body', ..., b'Massive Body',\n", - " b'Massive Body', b'Massive Body'], dtype='|S32')\n", + "\n", + "array(['Initial conditions', 'Body Initial',\n", + " 'Body Initial', ..., '', '', ''], dtype='\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "
<xarray.Dataset>\n",
+       "Dimensions:        ()\n",
+       "Coordinates:\n",
+       "    time           float64 0.0\n",
+       "    id             int32 2\n",
+       "Data variables: (12/48)\n",
+       "    npl            int32 1500\n",
+       "    ntp            int32 0\n",
+       "    name           <U32 'Body3'\n",
+       "    particle_type  <U12 'Massive Body'\n",
+       "    xhx            float64 -8.012e+06\n",
+       "    xhy            float64 -6.936e+06\n",
+       "    ...             ...\n",
+       "    origin_xhx     float64 6.013e-154\n",
+       "    origin_xhy     float64 6.013e-154\n",
+       "    origin_xhz     float64 0.0\n",
+       "    origin_vhx     float64 7.029e+06\n",
+       "    origin_vhy     float64 6.052e+06\n",
+       "    origin_vhz     float64 3.879e+03
" + ], + "text/plain": [ + "\n", + "Dimensions: ()\n", + "Coordinates:\n", + " time float64 0.0\n", + " id int32 2\n", + "Data variables: (12/48)\n", + " npl int32 ...\n", + " ntp int32 ...\n", + " name self%nbody) if (n == 0) return allocate(ind(n)) - call util_sort(self%id(1:n), ind(1:n)) + call util_sort(self%id(1:n), ind) do i = 1, n j = ind(i) diff --git a/src/symba/symba_util.f90 b/src/symba/symba_util.f90 index 93d4155e2..b8428dd99 100644 --- a/src/symba/symba_util.f90 +++ b/src/symba/symba_util.f90 @@ -79,7 +79,6 @@ module subroutine symba_util_append_pl(self, source, lsource_mask) call util_append(self%peri, source%peri, nold, nsrc, lsource_mask) call util_append(self%atp, source%atp, nold, nsrc, lsource_mask) call util_append(self%kin, source%kin, nold, nsrc, lsource_mask) - call util_append(self%info, source%info, nold, nsrc, lsource_mask) call util_append_pl(self, source, lsource_mask) ! Note: helio_pl does not have its own append method, so we skip back to the base class end associate @@ -148,7 +147,6 @@ module subroutine symba_util_append_tp(self, source, lsource_mask) call util_append(self%nplenc, source%nplenc, nold, nsrc, lsource_mask) call util_append(self%levelg, source%levelg, nold, nsrc, lsource_mask) call util_append(self%levelm, source%levelm, nold, nsrc, lsource_mask) - call util_append(self%info, source%info, nold, nsrc, lsource_mask) call util_append_tp(self, source, lsource_mask) ! Note: helio_tp does not have its own append method, so we skip back to the base class end associate @@ -229,7 +227,6 @@ module subroutine symba_util_fill_pl(self, inserts, lfill_list) call util_fill(keeps%peri, inserts%peri, lfill_list) call util_fill(keeps%atp, inserts%atp, lfill_list) call util_fill(keeps%kin, inserts%kin, lfill_list) - call util_fill(keeps%info, inserts%info, lfill_list) call util_fill_pl(keeps, inserts, lfill_list) ! Note: helio_pl does not have its own fill method, so we skip back to the base class class default @@ -260,7 +257,6 @@ module subroutine symba_util_fill_tp(self, inserts, lfill_list) call util_fill(keeps%nplenc, inserts%nplenc, lfill_list) call util_fill(keeps%levelg, inserts%levelg, lfill_list) call util_fill(keeps%levelm, inserts%levelm, lfill_list) - call util_fill(keeps%info, inserts%info, lfill_list) call util_fill_tp(keeps, inserts, lfill_list) ! Note: helio_tp does not have its own fill method, so we skip back to the base class class default @@ -607,7 +603,6 @@ module subroutine symba_util_resize_pl(self, nnew) call util_resize(self%peri, nnew) call util_resize(self%atp, nnew) call util_resize(self%kin, nnew) - call util_resize(self%info, nnew) call util_resize_pl(self, nnew) @@ -627,7 +622,6 @@ module subroutine symba_util_resize_tp(self, nnew) call util_resize(self%nplenc, nnew) call util_resize(self%levelg, nnew) call util_resize(self%levelm, nnew) - call util_resize(self%info, nnew) call util_resize_tp(self, nnew) @@ -779,7 +773,6 @@ module subroutine symba_util_sort_rearrange_pl(self, ind) call util_sort_rearrange(pl%isperi, ind, npl) call util_sort_rearrange(pl%peri, ind, npl) call util_sort_rearrange(pl%atp, ind, npl) - call util_sort_rearrange(pl%info, ind, npl) call util_sort_rearrange(pl%kin, ind, npl) call util_sort_rearrange_pl(pl,ind) @@ -803,7 +796,6 @@ module subroutine symba_util_sort_rearrange_tp(self, ind) call util_sort_rearrange(tp%nplenc, ind, ntp) call util_sort_rearrange(tp%levelg, ind, ntp) call util_sort_rearrange(tp%levelm, ind, ntp) - call util_sort_rearrange(tp%info, ind, ntp) call util_sort_rearrange_tp(tp,ind) end associate @@ -880,7 +872,6 @@ module subroutine symba_util_spill_pl(self, discards, lspill_list, ldestructive) call util_spill(keeps%isperi, discards%isperi, lspill_list, ldestructive) call util_spill(keeps%peri, discards%peri, lspill_list, ldestructive) call util_spill(keeps%atp, discards%atp, lspill_list, ldestructive) - call util_spill(keeps%info, discards%info, lspill_list, ldestructive) call util_spill(keeps%kin, discards%kin, lspill_list, ldestructive) call util_spill_pl(keeps, discards, lspill_list, ldestructive) @@ -945,7 +936,6 @@ module subroutine symba_util_spill_tp(self, discards, lspill_list, ldestructive) call util_spill(keeps%nplenc, discards%nplenc, lspill_list, ldestructive) call util_spill(keeps%levelg, discards%levelg, lspill_list, ldestructive) call util_spill(keeps%levelm, discards%levelm, lspill_list, ldestructive) - call util_spill(keeps%info, discards%info, lspill_list, ldestructive) call util_spill_tp(keeps, discards, lspill_list, ldestructive) class default diff --git a/src/util/util_append.f90 b/src/util/util_append.f90 index d6b18d4e8..98a98ac0a 100644 --- a/src/util/util_append.f90 +++ b/src/util/util_append.f90 @@ -128,23 +128,35 @@ module subroutine util_append_arr_info(arr, source, nold, nsrc, lsource_mask) logical, dimension(:), intent(in) :: lsource_mask !! Logical mask indicating which elements to append to ! Internals integer(I4B) :: nnew - class(swiftest_particle_info), allocatable :: tmp + class(swiftest_particle_info), dimension(:), allocatable :: arr_tmp, source_tmp if (.not. allocated(source)) return nnew = count(lsource_mask(1:nsrc)) - if (.not.allocated(arr)) then - select type(source) - class is (symba_particle_info) - allocate(symba_particle_info :: arr(nold+nnew)) - class default - allocate(swiftest_particle_info :: arr(nold+nnew)) - end select - else - call util_resize(arr, nold + nnew) - end if - arr(nold + 1:nold + nnew) = pack(source(1:nsrc), lsource_mask(1:nsrc)) + select type(source) + class is (symba_particle_info) + allocate(symba_particle_info :: arr_tmp(nold+nnew)) + if (nold > 0) then + arr_tmp(1:nold) = arr(1:nold) + deallocate(arr) + end if + class is (swiftest_particle_info) + allocate(swiftest_particle_info :: arr_tmp(nold+nnew)) + if (nold > 0) then + arr_tmp(1:nold) = arr(1:nold) + deallocate(arr) + end if + end select + + select type(source) + class is (symba_particle_info) + arr_tmp(nold + 1:nold + nnew) = pack(source(1:nsrc), lsource_mask(1:nsrc)) + class is (swiftest_particle_info) + arr_tmp(nold + 1:nold + nnew) = pack(source(1:nsrc), lsource_mask(1:nsrc)) + end select + + call move_alloc(arr_tmp, arr) return end subroutine util_append_arr_info diff --git a/src/util/util_fill.f90 b/src/util/util_fill.f90 index b1bf951d0..f3f6a3a95 100644 --- a/src/util/util_fill.f90 +++ b/src/util/util_fill.f90 @@ -92,12 +92,34 @@ module subroutine util_fill_arr_info(keeps, inserts, lfill_list) ! Arguments class(swiftest_particle_info), dimension(:), allocatable, intent(inout) :: keeps !! Array of values to keep class(swiftest_particle_info), dimension(:), allocatable, intent(in) :: inserts !! Array of values to insert into keep - logical, dimension(:), intent(in) :: lfill_list !! Logical array of bodies to merge into the keeps + logical, dimension(:), intent(in) :: lfill_list !! Logical array of bodies to merge into the keeps + ! Internals + class(swiftest_particle_info), dimension(:), allocatable :: ktmp, itmp + integer(I4B) :: nk, ni if (.not.allocated(keeps) .or. .not.allocated(inserts)) return - keeps(:) = unpack(keeps(:), .not.lfill_list(:), keeps(:)) - keeps(:) = unpack(inserts(:), lfill_list(:), keeps(:)) + nk = size(keeps) + ni = size(inserts) + + select type(keeps) + class is (symba_particle_info) + allocate(symba_particle_info :: ktmp(nk)) + class is (swiftest_particle_info) + allocate(swiftest_particle_info :: ktmp(nk)) + end select + + select type(inserts) + class is (symba_particle_info) + allocate(symba_particle_info :: itmp(ni)) + class is (swiftest_particle_info) + allocate(swiftest_particle_info :: itmp(ni)) + end select + + ktmp(:) = unpack(ktmp(:), .not.lfill_list(:), ktmp(:)) + ktmp(:) = unpack(itmp(:), lfill_list(:), ktmp(:)) + + keeps(:) = ktmp(:) return end subroutine util_fill_arr_info diff --git a/src/util/util_sort.f90 b/src/util/util_sort.f90 index 4b96dd9d1..4364aa02a 100644 --- a/src/util/util_sort.f90 +++ b/src/util/util_sort.f90 @@ -466,15 +466,21 @@ module subroutine util_sort_rearrange_arr_info(arr, ind, n) implicit none ! Arguments class(swiftest_particle_info), dimension(:), allocatable, intent(inout) :: arr !! Destination array - integer(I4B), dimension(:), intent(in) :: ind !! Index to rearrange against + integer(I4B), dimension(:), intent(in) :: ind !! Index to rearrange against integer(I4B), intent(in) :: n !! Number of elements in arr and ind to rearrange ! Internals class(swiftest_particle_info), dimension(:), allocatable :: tmp !! Temporary copy of array used during rearrange operation + integer(I4B) :: i if (.not. allocated(arr) .or. n <= 0) return - allocate(tmp, source=arr) + select type(arr) + class is (symba_particle_info) + allocate(symba_particle_info :: tmp(n)) + class is (swiftest_particle_info) + allocate(swiftest_particle_info :: tmp(n)) + end select tmp(1:n) = arr(ind(1:n)) - call move_alloc(tmp, arr) + arr(1:n) = tmp(1:n) return end subroutine util_sort_rearrange_arr_info