diff --git a/Makefile.Defines b/Makefile.Defines index 6effd7332..9e06d56ba 100644 --- a/Makefile.Defines +++ b/Makefile.Defines @@ -53,11 +53,11 @@ VTUNE_FLAGS = -g -O2 -qopt-report=5 -simd -shared-intel -qopenmp -debug inline-d IDEBUG = -O0 -init=snan,arrays -nogen-interfaces -no-pie -no-ftz -fpe-all=0 -g -traceback -mp1 -fp-model strict -fpe0 -debug all -align all -pad -ip -prec-div -prec-sqrt -assume protect-parens -CB -no-wrap-margin STRICTREAL = -fp-model strict -prec-div -prec-sqrt -assume protect-parens -SIMDVEC = -simd -xhost -align all -assume contiguous_assumed_shape -vecabi=cmdtarget -fp-model no-except +SIMDVEC = -simd -xhost -align all -assume contiguous_assumed_shape -vecabi=cmdtarget -fp-model no-except -fma PAR = -qopenmp -parallel HEAPARR = -heap-arrays 4194304 OPTREPORT = -qopt-report=5 -IPRODUCTION = -no-wrap-margin -O3 $(PAR) $(SIMDVEC) #$(HEAPARR) +IPRODUCTION = -no-wrap-margin -O3 -ipo -qopt-prefetch=0 -sox $(PAR) $(SIMDVEC) #$(HEAPARR) #gfortran flags GDEBUG = -g -Og -fbacktrace -fbounds-check -ffree-line-length-none @@ -68,10 +68,11 @@ GPRODUCTION = -O2 -ffree-line-length-none $(GPAR) #FFLAGS = $(IDEBUG) $(SIMDVEC) $(PAR) -FFLAGS = $(IPRODUCTION) $(STRICTREAL) $(OPTREPORT) -FFASTFLAGS = $(IPRODUCTION) -fp-model fast $(OPTREPORT) +#FFASTFLAGS = $(IDEBUG) $(SIMDVEC) $(PAR) +FFLAGS = $(IPRODUCTION) $(STRICTREAL) +FFASTFLAGS = $(IPRODUCTION) -fp-model fast FORTRAN = ifort -#AR = xiar +AR = xiar #FORTRAN = gfortran #FFLAGS = $(GDEBUG) $(GMEM) $(GPAR) diff --git a/examples/symba_mars_disk/Untitled.ipynb b/examples/symba_mars_disk/Untitled.ipynb deleted file mode 100644 index 7fec51502..000000000 --- a/examples/symba_mars_disk/Untitled.ipynb +++ /dev/null @@ -1,6 +0,0 @@ -{ - "cells": [], - "metadata": {}, - "nbformat": 4, - "nbformat_minor": 4 -} diff --git a/examples/symba_mars_disk/param.in b/examples/symba_mars_disk/param.in index 1769c6c74..b75341794 100644 --- a/examples/symba_mars_disk/param.in +++ b/examples/symba_mars_disk/param.in @@ -1,17 +1,17 @@ !Parameter file for the SyMBA-RINGMOONS test T0 0.0 -TSTOP 60000.0 +TSTOP 12000.0 DT 600.0 CB_IN cb.in PL_IN mars.in TP_IN tp.in IN_TYPE ASCII -ISTEP_OUT 100 -ISTEP_DUMP 100 -BIN_OUT bin.nc +ISTEP_OUT 1 +ISTEP_DUMP 1 +BIN_OUT bin.nc PARTICLE_OUT particle.dat -OUT_TYPE REAL8 -OUT_FORM XV +OUT_TYPE NETCDF_DOUBLE +OUT_FORM XVEL OUT_STAT REPLACE CHK_CLOSE yes CHK_RMIN 3389500.0 @@ -31,3 +31,4 @@ MU2KG 1.0 DU2M 1.0 TU2S 1.0 SEED 2 3080983 2220830 +FLATTEN_INTERACTIONS yes diff --git a/examples/symba_mars_disk/param.real8.in b/examples/symba_mars_disk/param.real8.in new file mode 100644 index 000000000..6ebc9197a --- /dev/null +++ b/examples/symba_mars_disk/param.real8.in @@ -0,0 +1,33 @@ +!Parameter file for the SyMBA-RINGMOONS test +T0 0.0 +TSTOP 6000.0 +DT 600.0 +CB_IN cb.in +PL_IN mars.in +TP_IN tp.in +IN_TYPE ASCII +ISTEP_OUT 1 +ISTEP_DUMP 1 +BIN_OUT bin.dat +PARTICLE_OUT particle.dat +OUT_TYPE REAL8 +OUT_FORM EL +OUT_STAT REPLACE +CHK_CLOSE yes +CHK_RMIN 3389500.0 +CHK_RMAX 3389500000.0 +CHK_EJECT 3389500000.0 +CHK_QMIN 3389500.0 +CHK_QMIN_COORD HELIO +CHK_QMIN_RANGE 3389500.0 338950000000.0 +EXTRA_FORCE no +RHILL_PRESENT yes +GMTINY 10000.0 +MIN_GMFRAG 1000.0 +ENERGY yes +FRAGMENTATION yes +ROTATION yes +MU2KG 1.0 +DU2M 1.0 +TU2S 1.0 +SEED 2 3080983 2220830 diff --git a/examples/symba_mars_disk/testnetcdf.ipynb b/examples/symba_mars_disk/testnetcdf.ipynb index be4bf6a11..a0dd92fe6 100644 --- a/examples/symba_mars_disk/testnetcdf.ipynb +++ b/examples/symba_mars_disk/testnetcdf.ipynb @@ -2,31 +2,42 @@ "cells": [ { "cell_type": "code", - "execution_count": 1, + "execution_count": 6, "metadata": {}, "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "The autoreload extension is already loaded. To reload it, use:\n", + " %reload_ext autoreload\n" + ] + }, { "data": { "text/plain": [ - "'/home/daminton/git/swiftest/examples/symba_mars_disk'" + "'/home/daminton/gittmp/swiftest/examples/symba_mars_disk'" ] }, - "execution_count": 1, + "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ + "%load_ext autoreload\n", + "%autoreload 2\n", "import swiftest\n", "import os\n", "import xarray as xr\n", "import numpy as np\n", + "from netCDF4 import Dataset\n", "os.getcwd()" ] }, { "cell_type": "code", - "execution_count": 2, + "execution_count": 7, "metadata": {}, "outputs": [ { @@ -35,20 +46,20 @@ "text": [ "Reading Swiftest file param.in\n", "\n", - "Creating Dataset\n", - "Successfully converted 6 output frames.\n", + "Creating Dataset from NetCDF file\n", + "Successfully converted 11 output frames.\n", "Swiftest simulation data stored as xarray DataSet .ds\n" ] } ], "source": [ "sim = swiftest.Simulation(param_file=\"param.in\")\n", - "sim.bin2xr()\n" + "sim.bin2xr()" ] }, { "cell_type": "code", - "execution_count": 3, + "execution_count": 8, "metadata": {}, "outputs": [ { @@ -405,43 +416,35 @@ " stroke: currentColor;\n", " fill: currentColor;\n", "}\n", - "
<xarray.DataArray 'status' (id: 1521)>\n",
-       "array(['ACTIVE', 'ACTIVE', 'ACTIVE', ..., 'ACTIVE', 'ACTIVE', 'ACTIVE'],\n",
-       "      dtype='<U17')\n",
+       "
<xarray.DataArray 'id' (id: 1525)>\n",
+       "array([   0,    1,    2, ..., 1522, 1523, 1524], dtype=int32)\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 ... 1518 1519 1520 1521 1522 1523 1524\n",
+       "Attributes:\n",
+       "    _FillValue:  -2147483647
" ], "text/plain": [ - "\n", - "array(['ACTIVE', 'ACTIVE', 'ACTIVE', ..., 'ACTIVE', 'ACTIVE', 'ACTIVE'],\n", - " dtype='\n", + "array([ 0, 1, 2, ..., 1522, 1523, 1524], dtype=int32)\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 ... 1518 1519 1520 1521 1522 1523 1524\n", + "Attributes:\n", + " _FillValue: -2147483647" ] }, - "execution_count": 3, + "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ - "sim.ds['status']" + "sim.ds.id" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, - "outputs": [], - "source": [ - "ds = xr.open_dataset('bin.nc')" - ] - }, - { - "cell_type": "code", - "execution_count": 5, - "metadata": {}, "outputs": [ { "data": { @@ -797,5372 +800,99 @@ " stroke: currentColor;\n", " fill: currentColor;\n", "}\n", - "
<xarray.DataArray 'status' ()>\n",
-       "array('Supercatastrophic', dtype='<U17')\n",
+       "
<xarray.DataArray 'id' (id: 1525)>\n",
+       "array([0.000e+00, 1.000e+00, 2.000e+00, ..., 1.522e+03, 1.523e+03, 1.524e+03])\n",
        "Coordinates:\n",
-       "    id       int32 726
"
+       "  * id       (id) float64 0.0 1.0 2.0 3.0 ... 1.522e+03 1.523e+03 1.524e+03
" ], "text/plain": [ - "\n", - "array('Supercatastrophic', dtype='\n", + "array([0.000e+00, 1.000e+00, 2.000e+00, ..., 1.522e+03, 1.523e+03, 1.524e+03])\n", "Coordinates:\n", - " id int32 726" + " * id (id) float64 0.0 1.0 2.0 3.0 ... 1.522e+03 1.523e+03 1.524e+03" ] }, - "execution_count": 5, + "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ - "sim.ds['status'].sel(id=726)" + "sim.ds.id" ] }, { "cell_type": "code", - "execution_count": 6, + "execution_count": 5, "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: 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" + "ename": "NameError", + "evalue": "name 'osd' is not defined", + "output_type": "error", + "traceback": [ + "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", + "\u001b[0;31mNameError\u001b[0m Traceback (most recent call last)", + "\u001b[0;32m\u001b[0m in \u001b[0;36m\u001b[0;34m\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mosd\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'id'\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m", + "\u001b[0;31mNameError\u001b[0m: name 'osd' is not defined" + ] } ], "source": [ - "sim.ds['discard_time']" + "osd['id']" ] }, { "cell_type": "code", - "execution_count": 7, + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, "metadata": {}, "outputs": [], "source": [ - "dslost = sim.ds.where(sim.ds['status'] != \"ACTIVE\", drop=True)" + "ds = xr.open_dataset(sim.param['BIN_OUT'], mask_and_scale=False)\n", + "ds = swiftest.io.clean_string_values(sim.param, ds)" ] }, { "cell_type": "code", - "execution_count": 8, + "execution_count": null, "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" - } - ], + "outputs": [], "source": [ - "dslost['discard_time']" + "ods = xr.open_dataset('bin.nc', mask_and_scale=False)" ] }, { "cell_type": "code", - "execution_count": 9, + "execution_count": null, "metadata": {}, "outputs": [], "source": [ - "lastloss = dslost.where(dslost['discard_time'] == 98400.0, drop=True)" + "ods['a']" ] }, { "cell_type": "code", - "execution_count": 10, + "execution_count": null, "metadata": {}, "outputs": [], "source": [ - "dsactive = sim.ds.where(sim.ds['status'] == \"ACTIVE\", drop=True)" + "ods.id" ] }, { "cell_type": "code", - "execution_count": 11, + "execution_count": null, "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" - } - ], + "outputs": [], "source": [ - "lastframe.sel(id=2)" + "sim.ds.sel(id=1501)['a']" ] }, { @@ -6170,7 +900,9 @@ "execution_count": null, "metadata": {}, "outputs": [], - "source": [] + "source": [ + "cds.sel(id=1501)['a']" + ] } ], "metadata": { diff --git a/src/helio/helio_kick.f90 b/src/helio/helio_kick.f90 index b06abd581..928e7c197 100644 --- a/src/helio/helio_kick.f90 +++ b/src/helio/helio_kick.f90 @@ -20,7 +20,7 @@ module subroutine helio_kick_getacch_pl(self, system, param, t, lbeg) if (self%nbody == 0) return associate(cb => system%cb, pl => self, npl => self%nbody) - call pl%accel_int() + call pl%accel_int(param) if (param%loblatecb) then call pl%accel_obl(system) if (lbeg) then @@ -65,9 +65,9 @@ module subroutine helio_kick_getacch_tp(self, system, param, t, lbeg) associate(tp => self, cb => system%cb, pl => system%pl, npl => system%pl%nbody) system%lbeg = lbeg if (system%lbeg) then - call tp%accel_int(pl%Gmass(1:npl), pl%xbeg(:,1:npl), npl) + call tp%accel_int(param, pl%Gmass(1:npl), pl%xbeg(:,1:npl), npl) else - call tp%accel_int(pl%Gmass(1:npl), pl%xend(:,1:npl), npl) + call tp%accel_int(param, pl%Gmass(1:npl), pl%xend(:,1:npl), npl) end if if (param%loblatecb) call tp%accel_obl(system) if (param%lextra_force) call tp%accel_user(system, param, t, lbeg) diff --git a/src/io/io.f90 b/src/io/io.f90 index f04c07e4f..51d10561d 100644 --- a/src/io/io.f90 +++ b/src/io/io.f90 @@ -601,6 +601,9 @@ module subroutine io_param_reader(self, unit, iotype, v_list, iostat, iomsg) case ("TIDES") call io_toupper(param_value) if (param_value == "YES" .or. param_value == 'T') param%ltides = .true. + case ("FLATTEN_INTERACTIONS") + call io_toupper(param_value) + if (param_value == "NO" .or. param_value == 'F') param%lflatten_interactions = .false. case ("FIRSTKICK") call io_toupper(param_value) if (param_value == "NO" .or. param_value == 'F') param%lfirstkick = .false. @@ -755,9 +758,10 @@ module subroutine io_param_reader(self, unit, iotype, v_list, iostat, iomsg) if ((param%qmin_alo > 0.0_DP) .or. (param%qmin_ahi > 0.0_DP)) write(*,*) "CHK_QMIN_RANGE = ",param%qmin_alo, param%qmin_ahi write(*,*) "EXTRA_FORCE = ",param%lextra_force write(*,*) "RHILL_PRESENT = ",param%lrhill_present - write(*,*) "ROTATION = ", param%lrotation - write(*,*) "TIDES = ", param%ltides + write(*,*) "ROTATION = ", param%lrotation + write(*,*) "TIDES = ", param%ltides write(*,*) "ENERGY = ",param%lenergy + write(*,*) "FLATTEN_INTERACTIONS = ",param%lflatten_interactions if (param%lenergy) write(*,*) "ENERGY_OUT = ",trim(adjustl(param%energy_out)) write(*,*) "MU2KG = ",param%MU2KG write(*,*) "TU2S = ",param%TU2S @@ -899,6 +903,7 @@ module subroutine io_param_writer(self, unit, iotype, v_list, iostat, iomsg) write(param_name, *) "GR"; write(param_value, Lfmt) param%lgr; write(unit, *, err = 667, iomsg = iomsg) adjustl(param_name) // trim(adjustl(param_value)) write(param_name, *) "ROTATION"; write(param_value, Lfmt) param%lrotation; write(unit, *, err = 667, iomsg = iomsg) adjustl(param_name) // trim(adjustl(param_value)) write(param_name, *) "TIDES"; write(param_value, Lfmt) param%ltides; write(unit, *, err = 667, iomsg = iomsg) adjustl(param_name) // trim(adjustl(param_value)) + write(param_name, *) "FLATTEN_INTERACTIONS"; write(param_value, Lfmt) param%lflatten_interactions; write(unit, *, err = 667, iomsg = iomsg) adjustl(param_name) // trim(adjustl(param_value)) if (param%lenergy) then write(param_name, *) "FIRSTENERGY"; write(param_value, Lfmt) param%lfirstenergy; write(unit, *, err = 667, iomsg = iomsg) adjustl(param_name) // trim(adjustl(param_value)) diff --git a/src/kick/kick.f90 b/src/kick/kick.f90 index 45da5cb31..766e829d6 100644 --- a/src/kick/kick.f90 +++ b/src/kick/kick.f90 @@ -2,7 +2,7 @@ use swiftest contains - module subroutine kick_getacch_int_pl(self) + module subroutine kick_getacch_int_pl(self, param) !! author: David A. Minton !! !! Compute direct cross (third) term heliocentric accelerations of massive bodies @@ -11,15 +11,20 @@ module subroutine kick_getacch_int_pl(self) !! Adapted from David E. Kaufmann's Swifter routine whm_kick_getacch_ah3.f90 and helio_kick_getacch_int.f90 implicit none ! Arguments - class(swiftest_pl), intent(inout) :: self !! Swiftest massive body object + class(swiftest_pl), intent(inout) :: self !! Swiftest massive body object + class(swiftest_parameters), intent(in) :: param !! Current swiftest run configuration parameters - call kick_getacch_int_all_pl(self%nbody, self%nplpl, self%k_plpl, self%xh, self%Gmass, self%radius, self%ah) + if (param%lflatten_interactions) then + call kick_getacch_int_all_flat_pl(self%nbody, self%nplpl, self%k_plpl, self%xh, self%Gmass, self%radius, self%ah) + else + call kick_getacch_int_all_triangular_pl(self%nbody, self%nbody, self%xh, self%Gmass, self%radius, self%ah) + end if return end subroutine kick_getacch_int_pl - module subroutine kick_getacch_int_tp(self, GMpl, xhp, npl) + module subroutine kick_getacch_int_tp(self, param, GMpl, xhp, npl) !! author: David A. Minton !! !! Compute direct cross (third) term heliocentric accelerations of test particles by massive bodies @@ -28,10 +33,11 @@ module subroutine kick_getacch_int_tp(self, GMpl, xhp, npl) !! Adapted from David E. Kaufmann's Swifter routine whm_kick_getacch_ah3.f90 and helio_kick_getacch_int_tp.f90 implicit none ! Arguments - class(swiftest_tp), intent(inout) :: self !! Swiftest test particle object - real(DP), dimension(:), intent(in) :: GMpl !! Massive body masses - real(DP), dimension(:,:), intent(in) :: xhp !! Massive body position vectors - integer(I4B), intent(in) :: npl !! Number of active massive bodies + class(swiftest_tp), intent(inout) :: self !! Swiftest test particle object + class(swiftest_parameters), intent(in) :: param !! Current swiftest run configuration parameters + real(DP), dimension(:), intent(in) :: GMpl !! Massive body masses + real(DP), dimension(:,:), intent(in) :: xhp !! Massive body position vectors + integer(I4B), intent(in) :: npl !! Number of active massive bodies if ((self%nbody == 0) .or. (npl == 0)) return @@ -41,10 +47,11 @@ module subroutine kick_getacch_int_tp(self, GMpl, xhp, npl) end subroutine kick_getacch_int_tp - module subroutine kick_getacch_int_all_pl(npl, nplpl, k_plpl, x, Gmass, radius, acc) + module subroutine kick_getacch_int_all_flat_pl(npl, nplpl, k_plpl, x, Gmass, radius, acc) !! author: David A. Minton !! - !! Compute direct cross (third) term heliocentric accelerations for massive bodies, with parallelization + !! Compute direct cross (third) term heliocentric accelerations for massive bodies, with parallelization. + !! This is the flattened (single loop) version that uses the k_plpl interaction pair index array !! !! Adapted from Hal Levison's Swift routine getacch_ah3.f !! Adapted from David E. Kaufmann's Swifter routine whm_kick_getacch_ah3.f90 and helio_kick_getacch_int.f9 @@ -88,7 +95,56 @@ module subroutine kick_getacch_int_all_pl(npl, nplpl, k_plpl, x, Gmass, radius, end do return - end subroutine kick_getacch_int_all_pl + end subroutine kick_getacch_int_all_flat_pl + + + module subroutine kick_getacch_int_all_triangular_pl(npl, nplm, x, Gmass, radius, acc) + !! author: David A. Minton + !! + !! Compute direct cross (third) term heliocentric accelerations for massive bodies, with parallelization. + !! This is the upper triangular matrix (double loop) version. + !! + !! Adapted from Hal Levison's Swift routine getacch_ah3.f + !! Adapted from David E. Kaufmann's Swifter routine whm_kick_getacch_ah3.f90 and helio_kick_getacch_int.f9 + implicit none + integer(I4B), intent(in) :: npl !! Total number of massive bodies + integer(I4B), intent(in) :: nplm !! Number of fully interacting massive bodies + real(DP), dimension(:,:), intent(in) :: x !! Position vector array + real(DP), dimension(:), intent(in) :: Gmass !! Array of massive body G*mass + real(DP), dimension(:), intent(in) :: radius !! Array of massive body radii + real(DP), dimension(:,:), intent(inout) :: acc !! Acceleration vector array + ! Internals + real(DP), dimension(NDIM,npl) :: ahi, ahj + integer(I4B) :: i, j + real(DP) :: rji2, rlim2 + real(DP) :: xr, yr, zr + + ahi(:,:) = 0.0_DP + ahj(:,:) = 0.0_DP + + !$omp parallel do default(private) schedule(static)& + !$omp shared(npl, nplm, x, Gmass, radius) & + !$omp lastprivate(rji2, rlim2, xr, yr, zr) & + !$omp reduction(+:ahi) & + !$omp reduction(-:ahj) + do i = 1, nplm + do concurrent(j = i+1:npl) + xr = x(1, j) - x(1, i) + yr = x(2, j) - x(2, i) + zr = x(3, j) - x(3, i) + rji2 = xr**2 + yr**2 + zr**2 + rlim2 = (radius(i) + radius(j))**2 + if (rji2 > rlim2) call kick_getacch_int_one_pl(rji2, xr, yr, zr, Gmass(i), Gmass(j), ahi(1,i), ahi(2,i), ahi(3,i), ahj(1,j), ahj(2,j), ahj(3,j)) + end do + end do + !$omp end parallel do + + do concurrent(i = 1:npl) + acc(:,i) = acc(:,i) + ahi(:,i) + ahj(:,i) + end do + + return + end subroutine kick_getacch_int_all_triangular_pl module subroutine kick_getacch_int_all_tp(ntp, npl, xtp, xpl, GMpl, lmask, acc) @@ -112,7 +168,7 @@ module subroutine kick_getacch_int_all_tp(ntp, npl, xtp, xpl, GMpl, lmask, acc) integer(I4B) :: i, j !$omp parallel do default(private) schedule(static)& - !$omp shared(npl, ntp, lmask, xtp, xpl) & + !$omp shared(npl, ntp, lmask, xtp, xpl, GMpl) & !$omp reduction(-:acc) do i = 1, ntp if (lmask(i)) then diff --git a/src/modules/rmvs_classes.f90 b/src/modules/rmvs_classes.f90 index b29cd02c4..680612de5 100644 --- a/src/modules/rmvs_classes.f90 +++ b/src/modules/rmvs_classes.f90 @@ -122,12 +122,13 @@ module subroutine rmvs_discard_tp(self, system, param) class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters end subroutine rmvs_discard_tp - module function rmvs_encounter_check_tp(self, system, dt) result(lencounter) + module function rmvs_encounter_check_tp(self, param, system, dt) result(lencounter) implicit none - class(rmvs_tp), intent(inout) :: self !! RMVS test particle object - class(rmvs_nbody_system), intent(inout) :: system !! RMVS nbody system object - real(DP), intent(in) :: dt !! step size - logical :: lencounter !! Returns true if there is at least one close encounter + class(rmvs_tp), intent(inout) :: self !! RMVS test particle object + class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters + class(rmvs_nbody_system), intent(inout) :: system !! RMVS nbody system object + real(DP), intent(in) :: dt !! step size + logical :: lencounter !! Returns true if there is at least one close encounter end function rmvs_encounter_check_tp module subroutine rmvs_io_write_encounter(t, id1, id2, Gmass1, Gmass2, radius1, radius2, xh1, xh2, vh1, vh2, enc_out) diff --git a/src/modules/swiftest_classes.f90 b/src/modules/swiftest_classes.f90 index 598cb4811..a1f60c0dd 100644 --- a/src/modules/swiftest_classes.f90 +++ b/src/modules/swiftest_classes.f90 @@ -133,6 +133,7 @@ module swiftest_classes logical :: loblatecb = .false. !! Calculate acceleration from oblate central body (automatically turns true if nonzero J2 is input) logical :: lrotation = .false. !! Include rotation states of big bodies logical :: ltides = .false. !! Include tidal dissipation + logical :: lflatten_interactions = .true. !! Use the flattened upper triangular matrix for pl-pl interactions (turning this on improves the speed but uses more memory) ! Initial values to pass to the energy report subroutine (usually only used in the case of a restart, otherwise these will be updated with initial conditions values) real(DP) :: Eorbit_orig = 0.0_DP !! Initial orbital energy @@ -432,7 +433,6 @@ module swiftest_classes integer(I4B) :: nenc !! Total number of encounters logical, dimension(:), allocatable :: lvdotr !! relative vdotr flag integer(I4B), dimension(:), allocatable :: status !! status of the interaction - integer(I8B), dimension(:), allocatable :: kidx !! index value of the encounter from the master k_plpl encounter list integer(I4B), dimension(:), allocatable :: index1 !! position of the first body in the encounter integer(I4B), dimension(:), allocatable :: index2 !! position of the second body in the encounter integer(I4B), dimension(:), allocatable :: id1 !! id of the first body in the encounter @@ -840,20 +840,22 @@ module subroutine io_write_hdr_system(self, iu, param) class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters end subroutine io_write_hdr_system - module subroutine kick_getacch_int_pl(self) + module subroutine kick_getacch_int_pl(self, param) implicit none - class(swiftest_pl), intent(inout) :: self !! Swiftest massive body object + class(swiftest_pl), intent(inout) :: self !! Swiftest massive body object + class(swiftest_parameters), intent(in) :: param !! Current swiftest run configuration parameters end subroutine kick_getacch_int_pl - module subroutine kick_getacch_int_tp(self, GMpl, xhp, npl) + module subroutine kick_getacch_int_tp(self, param, GMpl, xhp, npl) implicit none - class(swiftest_tp), intent(inout) :: self !! Swiftest test particle - real(DP), dimension(:), intent(in) :: GMpl !! Massive body masses - real(DP), dimension(:,:), intent(in) :: xhp !! Massive body position vectors - integer(I4B), intent(in) :: npl !! Number of active massive bodies + class(swiftest_tp), intent(inout) :: self !! Swiftest test particle object + class(swiftest_parameters), intent(in) :: param !! Current swiftest run configuration parameters + real(DP), dimension(:), intent(in) :: GMpl !! Massive body masses + real(DP), dimension(:,:), intent(in) :: xhp !! Massive body position vectors + integer(I4B), intent(in) :: npl !! Number of active massive bodies end subroutine kick_getacch_int_tp - module subroutine kick_getacch_int_all_pl(npl, nplpl, k_plpl, x, Gmass, radius, acc) + module subroutine kick_getacch_int_all_flat_pl(npl, nplpl, k_plpl, x, Gmass, radius, acc) implicit none integer(I4B), intent(in) :: npl !! Number of massive bodies integer(I8B), intent(in) :: nplpl !! Number of massive body interactions to compute @@ -862,7 +864,17 @@ module subroutine kick_getacch_int_all_pl(npl, nplpl, k_plpl, x, Gmass, radius, real(DP), dimension(:), intent(in) :: Gmass !! Array of massive body G*mass real(DP), dimension(:), intent(in) :: radius !! Array of massive body radii real(DP), dimension(:,:), intent(inout) :: acc !! Acceleration vector array - end subroutine kick_getacch_int_all_pl + end subroutine kick_getacch_int_all_flat_pl + + module subroutine kick_getacch_int_all_triangular_pl(npl, nplm, x, Gmass, radius, acc) + implicit none + integer(I4B), intent(in) :: npl !! Total number of massive bodies + integer(I4B), intent(in) :: nplm !! Number of fully interacting massive bodies + real(DP), dimension(:,:), intent(in) :: x !! Position vector array + real(DP), dimension(:), intent(in) :: Gmass !! Array of massive body G*mass + real(DP), dimension(:), intent(in) :: radius !! Array of massive body radii + real(DP), dimension(:,:), intent(inout) :: acc !! Acceleration vector array + end subroutine kick_getacch_int_all_triangular_pl module subroutine kick_getacch_int_all_tp(ntp, npl, xtp, xpl, GMpl, lmask, acc) implicit none diff --git a/src/modules/symba_classes.f90 b/src/modules/symba_classes.f90 index 40267e298..0faefc3d1 100644 --- a/src/modules/symba_classes.f90 +++ b/src/modules/symba_classes.f90 @@ -256,40 +256,34 @@ module subroutine symba_drift_tp(self, system, param, dt) real(DP), intent(in) :: dt !! Stepsize end subroutine symba_drift_tp - module pure subroutine symba_encounter_check_one(xr, yr, zr, vxr, vyr, vzr, rhill1, rhill2, dt, irec, lencounter, lvdotr) - !$omp declare simd(symba_encounter_check_one) - implicit none - real(DP), intent(in) :: xr, yr, zr, vxr, vyr, vzr - real(DP), intent(in) :: rhill1, rhill2, dt - integer(I4B), intent(in) :: irec - logical, intent(out) :: lencounter, lvdotr - end subroutine symba_encounter_check_one - - module function symba_encounter_check_pl(self, system, dt, irec) result(lany_encounter) - implicit none - class(symba_pl), intent(inout) :: self !! SyMBA test particle object - class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object - real(DP), intent(in) :: dt !! step size - integer(I4B), intent(in) :: irec !! Current recursion level - logical :: lany_encounter !! Returns true if there is at least one close encounter + module function symba_encounter_check_pl(self, param, system, dt, irec) result(lany_encounter) + implicit none + class(symba_pl), intent(inout) :: self !! SyMBA test particle object + class(swiftest_parameters), intent(in) :: param !! Current swiftest run configuration parameters + class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object + real(DP), intent(in) :: dt !! step size + integer(I4B), intent(in) :: irec !! Current recursion level + logical :: lany_encounter !! Returns true if there is at least one close encounter end function symba_encounter_check_pl - module function symba_encounter_check(self, system, dt, irec) result(lany_encounter) + module function symba_encounter_check(self, param, system, dt, irec) result(lany_encounter) implicit none - class(symba_encounter), intent(inout) :: self !! SyMBA pl-pl encounter list object - class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object - real(DP), intent(in) :: dt !! step size - integer(I4B), intent(in) :: irec !! Current recursion level - logical :: lany_encounter !! Returns true if there is at least one close encounter + class(symba_encounter), intent(inout) :: self !! SyMBA pl-pl encounter list object + class(swiftest_parameters), intent(in) :: param !! Current swiftest run configuration parameters + class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object + real(DP), intent(in) :: dt !! step size + integer(I4B), intent(in) :: irec !! Current recursion level + logical :: lany_encounter !! Returns true if there is at least one close encounter end function symba_encounter_check - module function symba_encounter_check_tp(self, system, dt, irec) result(lany_encounter) + module function symba_encounter_check_tp(self, param, system, dt, irec) result(lany_encounter) implicit none - class(symba_tp), intent(inout) :: self !! SyMBA test particle object - class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object - real(DP), intent(in) :: dt !! step size - integer(I4B), intent(in) :: irec !! Current recursion level - logical :: lany_encounter !! Returns true if there is at least one close encounter + class(symba_tp), intent(inout) :: self !! SyMBA test particle object + class(swiftest_parameters), intent(in) :: param !! Current swiftest run configuration parameters + class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object + real(DP), intent(in) :: dt !! step size + integer(I4B), intent(in) :: irec !! Current recursion level + logical :: lany_encounter !! Returns true if there is at least one close encounter end function symba_encounter_check_tp module function symba_collision_casedisruption(system, param, colliders, frag) result(status) @@ -358,9 +352,10 @@ module subroutine symba_io_write_discard(self, param) class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters end subroutine symba_io_write_discard - module subroutine symba_kick_getacch_int_pl(self) + module subroutine symba_kick_getacch_int_pl(self, param) implicit none - class(symba_pl), intent(inout) :: self + class(symba_pl), intent(inout) :: self !! SyMBA massive body object + class(swiftest_parameters), intent(in) :: param !! Current swiftest run configuration parameters end subroutine symba_kick_getacch_int_pl module subroutine symba_kick_getacch_pl(self, system, param, t, lbeg) diff --git a/src/rmvs/rmvs_encounter_check.f90 b/src/rmvs/rmvs_encounter_check.f90 index 4c59f0a15..dd492032b 100644 --- a/src/rmvs/rmvs_encounter_check.f90 +++ b/src/rmvs/rmvs_encounter_check.f90 @@ -2,7 +2,7 @@ use swiftest contains - module function rmvs_encounter_check_tp(self, system, dt) result(lencounter) + module function rmvs_encounter_check_tp(self, param, system, dt) result(lencounter) !! author: David A. Minton !! !! Determine whether a test particle and planet are having or will have an encounter within the next time step @@ -11,9 +11,10 @@ module function rmvs_encounter_check_tp(self, system, dt) result(lencounter) !! Adapted from Hal Levison's Swift routine rmvs3_chk.f implicit none ! Arguments - class(rmvs_tp), intent(inout) :: self !! RMVS test particle object - class(rmvs_nbody_system), intent(inout) :: system !! RMVS nbody system object - real(DP), intent(in) :: dt !! step size + class(rmvs_tp), intent(inout) :: self !! RMVS test particle object + class(swiftest_parameters), intent(in) :: param !! Current swiftest run configuration parameters + class(rmvs_nbody_system), intent(inout) :: system !! RMVS nbody system object + real(DP), intent(in) :: dt !! step size ! Result logical :: lencounter !! Returns true if there is at least one close encounter ! Internals diff --git a/src/rmvs/rmvs_step.f90 b/src/rmvs/rmvs_step.f90 index b8cfcd688..00b66a9f8 100644 --- a/src/rmvs/rmvs_step.f90 +++ b/src/rmvs/rmvs_step.f90 @@ -35,7 +35,7 @@ module subroutine rmvs_step_system(self, param, t, dt) call pl%set_beg_end(xbeg = xbeg, vbeg = vbeg) ! ****** Check for close encounters ***** ! system%rts = RHSCALE - lencounter = tp%encounter_check(system, dt) + lencounter = tp%encounter_check(param, system, dt) if (lencounter) then lfirstpl = pl%lfirst pl%outer(0)%x(:, 1:npl) = xbeg(:, 1:npl) @@ -178,7 +178,7 @@ subroutine rmvs_step_out(cb, pl, tp, system, param, t, dt) vbeg = pl%outer(outer_index - 1)%v(:, 1:npl), & xend = pl%outer(outer_index )%x(:, 1:npl)) system%rts = RHPSCALE - lencounter = tp%encounter_check(system, dto) + lencounter = tp%encounter_check(param, system, dto) if (lencounter) then ! Interpolate planets in inner encounter region call rmvs_interp_in(cb, pl, system, param, dto, outer_index) diff --git a/src/setup/setup.f90 b/src/setup/setup.f90 index 41c0d9fba..18d8f2189 100644 --- a/src/setup/setup.f90 +++ b/src/setup/setup.f90 @@ -85,7 +85,6 @@ module subroutine setup_encounter(self, n) if (allocated(self%lvdotr)) deallocate(self%lvdotr) if (allocated(self%status)) deallocate(self%status) - if (allocated(self%kidx)) deallocate(self%kidx) if (allocated(self%index1)) deallocate(self%index1) if (allocated(self%index2)) deallocate(self%index2) if (allocated(self%id1)) deallocate(self%id1) @@ -101,7 +100,6 @@ module subroutine setup_encounter(self, n) allocate(self%lvdotr(n)) allocate(self%status(n)) - allocate(self%kidx(n)) allocate(self%index1(n)) allocate(self%index2(n)) allocate(self%id1(n)) @@ -114,7 +112,6 @@ module subroutine setup_encounter(self, n) self%lvdotr(:) = .false. self%status(:) = INACTIVE - self%kidx(:) = 0_I8B self%index1(:) = 0 self%index2(:) = 0 self%id1(:) = 0 diff --git a/src/symba/symba_collision.f90 b/src/symba/symba_collision.f90 index a3990cabd..864eaa723 100644 --- a/src/symba/symba_collision.f90 +++ b/src/symba/symba_collision.f90 @@ -291,7 +291,6 @@ module function symba_collision_check_encounter(self, system, param, t, dt, irec class is (symba_pl) select type(tp => system%tp) class is (symba_tp) - nenc = self%nenc nenc = self%nenc allocate(lmask(nenc)) lmask(:) = ((self%status(1:nenc) == ACTIVE) .and. (pl%levelg(self%index1(1:nenc)) >= irec)) diff --git a/src/symba/symba_encounter_check.f90 b/src/symba/symba_encounter_check.f90 index e20e51ea9..3ba6b4adf 100644 --- a/src/symba/symba_encounter_check.f90 +++ b/src/symba/symba_encounter_check.f90 @@ -2,18 +2,21 @@ use swiftest contains - subroutine symba_encounter_check_all(nplplm, k_plpl, x, v, rhill, dt, irec, lencounter, loc_lvdotr) + subroutine symba_encounter_check_all_flat(nplplm, k_plpl, x, v, rhill, dt, irec, lencounter, loc_lvdotr) !! author: David A. Minton !! - !! Check for encounters between massive bodies. Split off from the main subroutine for performance + !! Check for encounters between massive bodies. Split off from the main subroutine for performance. + !! This is the flat (single loop) version. implicit none - integer(I8B), intent(in) :: nplplm - integer(I4B), dimension(:,:), intent(in) :: k_plpl - real(DP), dimension(:,:), intent(in) :: x, v - real(DP), dimension(:), intent(in) :: rhill - real(DP), intent(in) :: dt - integer(I4B), intent(in) :: irec - logical, dimension(:), intent(out) :: lencounter, loc_lvdotr + integer(I8B), intent(in) :: nplplm !! Total number of plm-pl interactions + integer(I4B), dimension(:,:), intent(in) :: k_plpl !! List of all pl-pl interactions + real(DP), dimension(:,:), intent(in) :: x !! Position vectors of massive bodies + real(DP), dimension(:,:), intent(in) :: v !! Velocity vectors of massive bodies + real(DP), dimension(:), intent(in) :: rhill !! Hill's radii of massive bodies + real(DP), intent(in) :: dt !! Step size + integer(I4B), intent(in) :: irec !! Current recursion depth + logical, dimension(:), intent(out) :: lencounter !! Logical array indicating which pair is in an encounter state + logical, dimension(:), intent(out) :: loc_lvdotr !! Logical array indicating the sign of v .dot. x for each encounter ! Internals integer(I8B) :: k integer(I4B) :: i, j @@ -38,43 +41,123 @@ subroutine symba_encounter_check_all(nplplm, k_plpl, x, v, rhill, dt, irec, len !$omp end parallel do simd return - end subroutine symba_encounter_check_all + end subroutine symba_encounter_check_all_flat + + + subroutine symba_encounter_check_all_triangular(npl, nplm, x, v, rhill, dt, irec, lvdotr, index1, index2, nenc) + !! author: David A. Minton + !! + !! Check for encounters between massive bodies. Split off from the main subroutine for performance + !! This is the upper triangular (double loop) version. + implicit none + integer(I4B), intent(in) :: npl !! Total number of massive bodies + integer(I4B), intent(in) :: nplm !! Number of fully interacting massive bodies + real(DP), dimension(:,:), intent(in) :: x !! Position vectors of massive bodies + real(DP), dimension(:,:), intent(in) :: v !! Velocity vectors of massive bodies + real(DP), dimension(:), intent(in) :: rhill !! Hill's radii of massive bodies + real(DP), intent(in) :: dt !! Step size + integer(I4B), intent(in) :: irec !! Current recursion depth + logical, dimension(:), allocatable, intent(out) :: lvdotr !! Logical flag indicating the sign of v .dot. x + integer(I4B), dimension(:), allocatable, intent(out) :: index1 !! List of indices for body 1 in each interaction + integer(I4B), dimension(:), allocatable, intent(out) :: index2 !! List of indices for body 2 in each interaction + integer(I4B), intent(out) :: nenc !! Total number of interactions + ! Internals + integer(I4B) :: i, j, nenci, j0, j1 + real(DP) :: xr, yr, zr, vxr, vyr, vzr, rhill1, rhill2 + logical, dimension(npl) :: lencounteri, lvdotri + integer(I4B), dimension(npl) :: ind_arr + type lenctype + logical, dimension(:), allocatable :: lvdotr + integer(I4B), dimension(:), allocatable :: index2 + integer(I4B) :: nenc + end type + type(lenctype), dimension(nplm) :: lenc + + ind_arr(:) = [(i, i = 1, npl)] + !$omp parallel do default(private) schedule(static)& + !$omp shared(npl, nplm, x, v, rhill, dt, irec, lenc, ind_arr) & + !$omp lastprivate(xr, yr, zr, vxr, vyr, vzr, rhill1, rhill2, lencounteri, lvdotri) + do i = 1, nplm + do concurrent(j = i+1:npl) + xr = x(1, j) - x(1, i) + yr = x(2, j) - x(2, i) + zr = x(3, j) - x(3, i) + vxr = v(1, j) - v(1, i) + vyr = v(2, j) - v(2, i) + vzr = v(3, j) - v(3, i) + rhill1 = rhill(i) + rhill2 = rhill(j) + call symba_encounter_check_one(xr, yr, zr, vxr, vyr, vzr, rhill1, rhill2, dt, irec, lencounteri(j), lvdotri(j)) + end do + nenci = count(lencounteri(i+1:npl)) + if (nenci > 0) then + allocate(lenc(i)%lvdotr(nenci), lenc(i)%index2(nenci)) + lenc(i)%nenc = nenci + lenc(i)%lvdotr(:) = pack(lvdotri(i+1:npl), lencounteri(i+1:npl)) + lenc(i)%index2(:) = pack(ind_arr(i+1:npl), lencounteri(i+1:npl)) + end if + end do + !$omp end parallel do + + associate(nenc_arr => lenc(:)%nenc) + nenc = sum(nenc_arr(1:nplm)) + end associate + if (nenc > 0) then + allocate(lvdotr(nenc)) + allocate(index1(nenc)) + allocate(index2(nenc)) + j0 = 1 + do i = 1, nplm + if (lenc(i)%nenc > 0) then + j1 = j0 + lenc(i)%nenc - 1 + lvdotr(j0:j1) = lenc(i)%lvdotr(:) + index1(j0:j1) = i + index2(j0:j1) = lenc(i)%index2(:) + j0 = j1 + 1 + end if + end do + end if + + return + end subroutine symba_encounter_check_all_triangular - module function symba_encounter_check_pl(self, system, dt, irec) result(lany_encounter) + module function symba_encounter_check_pl(self, param, system, dt, irec) result(lany_encounter) !! author: David A. Minton !! !! Check for an encounter between massive bodies. !! implicit none ! Arguments - class(symba_pl), intent(inout) :: self !! SyMBA test particle object - class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object - real(DP), intent(in) :: dt !! step size - integer(I4B), intent(in) :: irec !! Current recursion level + class(symba_pl), intent(inout) :: self !! SyMBA test particle object + class(swiftest_parameters), intent(in) :: param !! Current swiftest run configuration parameters + class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object + real(DP), intent(in) :: dt !! step size + integer(I4B), intent(in) :: irec !! Current recursion level ! Result logical :: lany_encounter !! Returns true if there is at least one close encounter ! Internals integer(I8B) :: k, nplplm, kenc - integer(I4B) :: i, j, nenc, npl + integer(I4B) :: i, j, nenc, npl, nplm logical, dimension(:), allocatable :: lencounter, loc_lvdotr, lvdotr integer(I4B), dimension(:), allocatable :: index1, index2 + integer(I4B), dimension(:,:), allocatable :: k_plpl_enc if (self%nbody == 0) return - associate(pl => self) - nplplm = pl%nplplm + associate(pl => self, plplenc_list => system%plplenc_list) npl = pl%nbody - allocate(lencounter(nplplm)) - allocate(loc_lvdotr(nplplm)) + if (param%lflatten_interactions) then + nplplm = pl%nplplm + allocate(lencounter(nplplm)) + allocate(loc_lvdotr(nplplm)) - call symba_encounter_check_all(nplplm, pl%k_plpl, pl%xh, pl%vh, pl%rhill, dt, irec, lencounter, loc_lvdotr) + call symba_encounter_check_all_flat(nplplm, pl%k_plpl, pl%xh, pl%vh, pl%rhill, dt, irec, lencounter, loc_lvdotr) - nenc = count(lencounter(:)) + nenc = count(lencounter(:)) - lany_encounter = nenc > 0 - if (lany_encounter) then - associate(plplenc_list => system%plplenc_list) + lany_encounter = nenc > 0 + if (lany_encounter) then call plplenc_list%resize(nenc) allocate(lvdotr(nenc)) allocate(index1(nenc)) @@ -85,33 +168,45 @@ module function symba_encounter_check_pl(self, system, dt, irec) result(lany_enc deallocate(lencounter, loc_lvdotr) call move_alloc(lvdotr, plplenc_list%lvdotr) call move_alloc(index1, plplenc_list%index1) - call move_alloc(index2, plplenc_list%index2) - do k = 1, nenc - i = plplenc_list%index1(k) - j = plplenc_list%index2(k) - call util_index_eucl_ij_to_k(npl, i, j, kenc) - plplenc_list%kidx(k) = kenc - plplenc_list%id1(k) = pl%id(plplenc_list%index1(k)) - plplenc_list%id2(k) = pl%id(plplenc_list%index2(k)) - plplenc_list%status(k) = ACTIVE - plplenc_list%level(k) = irec - pl%lencounter(i) = .true. - pl%lencounter(j) = .true. - pl%levelg(i) = irec - pl%levelm(i) = irec - pl%levelg(j) = irec - pl%levelm(j) = irec - pl%nplenc(i) = pl%nplenc(i) + 1 - pl%nplenc(j) = pl%nplenc(j) + 1 - end do - end associate + call move_alloc(index2, plplenc_list%index2) + end if + else + nplm = pl%nplm + call symba_encounter_check_all_triangular(npl, nplm, pl%xh, pl%vh, pl%rhill, dt, irec, lvdotr, index1, index2, nenc) + lany_encounter = nenc > 0 + if (lany_encounter) then + call plplenc_list%resize(nenc) + call move_alloc(lvdotr, plplenc_list%lvdotr) + call move_alloc(index1, plplenc_list%index1) + call move_alloc(index2, plplenc_list%index2) + end if + end if + + if (lany_encounter) then + do k = 1, nenc + i = plplenc_list%index1(k) + j = plplenc_list%index2(k) + plplenc_list%id1(k) = pl%id(i) + plplenc_list%id2(k) = pl%id(j) + plplenc_list%status(k) = ACTIVE + plplenc_list%level(k) = irec + pl%lencounter(i) = .true. + pl%lencounter(j) = .true. + pl%levelg(i) = irec + pl%levelm(i) = irec + pl%levelg(j) = irec + pl%levelm(j) = irec + pl%nplenc(i) = pl%nplenc(i) + 1 + pl%nplenc(j) = pl%nplenc(j) + 1 + end do end if end associate + return end function symba_encounter_check_pl - module function symba_encounter_check(self, system, dt, irec) result(lany_encounter) + module function symba_encounter_check(self, param, system, dt, irec) result(lany_encounter) !! author: David A. Minton !! !! Check for an encounter between test particles and massive bodies in the pltpenc list. @@ -120,11 +215,12 @@ module function symba_encounter_check(self, system, dt, irec) result(lany_encoun !! Adapted from portions of David E. Kaufmann's Swifter routine: symba_step_recur.f90 implicit none ! Arguments - class(symba_encounter), intent(inout) :: self !! SyMBA pl-pl encounter list object - class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object - real(DP), intent(in) :: dt !! step size - integer(I4B), intent(in) :: irec !! Current recursion level - logical :: lany_encounter !! Returns true if there is at least one close encounter + class(symba_encounter), intent(inout) :: self !! SyMBA pl-pl encounter list object + class(swiftest_parameters), intent(in) :: param !! Current swiftest run configuration parameters + class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object + real(DP), intent(in) :: dt !! step size + integer(I4B), intent(in) :: irec !! Current recursion level + logical :: lany_encounter !! Returns true if there is at least one close encounter ! Internals integer(I4B) :: i, j, k, lidx, nenc_enc real(DP), dimension(NDIM) :: xr, vr @@ -197,17 +293,18 @@ module function symba_encounter_check(self, system, dt, irec) result(lany_encoun end function symba_encounter_check - module function symba_encounter_check_tp(self, system, dt, irec) result(lany_encounter) + module function symba_encounter_check_tp(self, param, system, dt, irec) result(lany_encounter) !! author: David A. Minton !! !! Check for an encounter between test particles and massive bodies. !! implicit none ! Arguments - class(symba_tp), intent(inout) :: self !! SyMBA test particle object - class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object - real(DP), intent(in) :: dt !! step size - integer(I4B), intent(in) :: irec !! Current recursion level + class(symba_tp), intent(inout) :: self !! SyMBA test particle object + class(swiftest_parameters), intent(in) :: param !! Current swiftest run configuration parameters + class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object + real(DP), intent(in) :: dt !! step size + integer(I4B), intent(in) :: irec !! Current recursion level ! Result logical :: lany_encounter !! Returns true if there is at least one close encounter ! Internals @@ -269,7 +366,7 @@ module function symba_encounter_check_tp(self, system, dt, irec) result(lany_enc end function symba_encounter_check_tp - module pure subroutine symba_encounter_check_one(xr, yr, zr, vxr, vyr, vzr, rhill1, rhill2, dt, irec, lencounter, lvdotr) + pure subroutine symba_encounter_check_one(xr, yr, zr, vxr, vyr, vzr, rhill1, rhill2, dt, irec, lencounter, lvdotr) !$omp declare simd(symba_encounter_check_one) !! author: David A. Minton !! diff --git a/src/symba/symba_kick.f90 b/src/symba/symba_kick.f90 index 325ef5a45..fd0354612 100644 --- a/src/symba/symba_kick.f90 +++ b/src/symba/symba_kick.f90 @@ -2,7 +2,7 @@ use swiftest contains - module subroutine symba_kick_getacch_int_pl(self) + module subroutine symba_kick_getacch_int_pl(self, param) !! author: David A. Minton !! !! Compute direct cross (third) term heliocentric accelerations of massive bodies, with no mutual interactions between bodies below GMTINY @@ -11,9 +11,14 @@ module subroutine symba_kick_getacch_int_pl(self) !! Adapted from David E. Kaufmann's Swifter routine helio_kick_getacch_int.f90 implicit none ! Arguments - class(symba_pl), intent(inout) :: self + class(symba_pl), intent(inout) :: self !! SyMBA massive body object + class(swiftest_parameters), intent(in) :: param !! Current swiftest run configuration parameter - call kick_getacch_int_all_pl(self%nbody, self%nplplm, self%k_plpl, self%xh, self%Gmass, self%radius, self%ah) + if (param%lflatten_interactions) then + call kick_getacch_int_all_flat_pl(self%nbody, self%nplplm, self%k_plpl, self%xh, self%Gmass, self%radius, self%ah) + else + call kick_getacch_int_all_triangular_pl(self%nbody, self%nplm, self%xh, self%Gmass, self%radius, self%ah) + end if return end subroutine symba_kick_getacch_int_pl @@ -43,16 +48,17 @@ module subroutine symba_kick_getacch_pl(self, system, param, t, lbeg) if (self%nbody == 0) return select type(system) class is (symba_nbody_system) - associate(pl => self, npl => self%nbody, plplenc_list => system%plplenc_list, radius => self%radius) + associate(pl => self, npl => self%nbody, nplm => self%nplm, plplenc_list => system%plplenc_list, radius => self%radius) ! Apply kicks to all bodies (including those in the encounter list) call helio_kick_getacch_pl(pl, system, param, t, lbeg) if (plplenc_list%nenc > 0) then ! Remove kicks from bodies involved currently in the encounter list, as these are dealt with separately. + ah_enc(:,:) = 0.0_DP nplplenc = int(plplenc_list%nenc, kind=I8B) allocate(k_plpl_enc(2,nplplenc)) - k_plpl_enc(:,1:nplplenc) = pl%k_plpl(:,plplenc_list%kidx(1:nplplenc)) - ah_enc(:,:) = 0.0_DP - call kick_getacch_int_all_pl(npl, nplplenc, k_plpl_enc, pl%xh, pl%Gmass, pl%radius, ah_enc) + k_plpl_enc(1,1:nplplenc) = plplenc_list%index1(1:nplplenc) + k_plpl_enc(2,1:nplplenc) = plplenc_list%index2(1:nplplenc) + call kick_getacch_int_all_flat_pl(npl, nplplenc, k_plpl_enc, pl%xh, pl%Gmass, pl%radius, ah_enc) pl%ah(:,1:npl) = pl%ah(:,1:npl) - ah_enc(:,1:npl) end if @@ -62,7 +68,6 @@ module subroutine symba_kick_getacch_pl(self, system, param, t, lbeg) return end subroutine symba_kick_getacch_pl - module subroutine symba_kick_getacch_tp(self, system, param, t, lbeg) !! author: David A. Minton !! diff --git a/src/symba/symba_step.f90 b/src/symba/symba_step.f90 index 9f4508550..dc07bf561 100644 --- a/src/symba/symba_step.f90 +++ b/src/symba/symba_step.f90 @@ -26,7 +26,7 @@ module subroutine symba_step_system(self, param, t, dt) select type(param) class is (symba_parameters) call self%reset(param) - lencounter = pl%encounter_check(self, dt, 0) .or. tp%encounter_check(self, dt, 0) + lencounter = pl%encounter_check(param, self, dt, 0) .or. tp%encounter_check(param, self, dt, 0) if (lencounter) then call self%interp(param, t, dt) param%lfirstkick = .true. @@ -173,7 +173,7 @@ module recursive subroutine symba_step_recur_system(self, param, t, ireci) nloops = NTENC end if do j = 1, nloops - lencounter = plplenc_list%encounter_check(system, dtl, irecp) .or. pltpenc_list%encounter_check(system, dtl, irecp) + lencounter = plplenc_list%encounter_check(param, system, dtl, irecp) .or. pltpenc_list%encounter_check(param, system, dtl, irecp) call plplenc_list%kick(system, dth, irecp, 1) call pltpenc_list%kick(system, dth, irecp, 1) diff --git a/src/symba/symba_util.f90 b/src/symba/symba_util.f90 index a59d1f0b8..0d062894f 100644 --- a/src/symba/symba_util.f90 +++ b/src/symba/symba_util.f90 @@ -292,18 +292,19 @@ module subroutine symba_util_index_eucl_plpl(self, param) npl = int(self%nbody, kind=I8B) nplm = count(.not. pl%lmtiny(1:npl)) pl%nplm = int(nplm, kind=I4B) - pl%nplpl = (npl * (npl - 1) / 2) ! number of entries in a strict lower triangle, npl x npl, minus first column pl%nplplm = nplm * npl - nplm * (nplm + 1) / 2 ! number of entries in a strict lower triangle, npl x npl, minus first column including only mutually interacting bodies - if (allocated(self%k_plpl)) deallocate(self%k_plpl) ! Reset the index array if it's been set previously - allocate(self%k_plpl(2, pl%nplpl)) - do concurrent (i = 1:npl) - do concurrent (j = i+1:npl) - call util_index_eucl_ij_to_k(npl, i, j, k) - self%k_plpl(1, k) = i - self%k_plpl(2, k) = j + if (param%lflatten_interactions) then + if (allocated(self%k_plpl)) deallocate(self%k_plpl) ! Reset the index array if it's been set previously + allocate(self%k_plpl(2, pl%nplpl)) + do concurrent (i = 1:npl) + do concurrent (j = i+1:npl) + call util_index_eucl_ij_to_k(npl, i, j, k) + self%k_plpl(1, k) = i + self%k_plpl(2, k) = j + end do end do - end do + end if end associate return @@ -486,7 +487,7 @@ module subroutine symba_util_rearray_pl(self, system, param) allocate(levelg_orig_pl, source=pl%levelg) allocate(levelm_orig_pl, source=pl%levelm) allocate(nplenc_orig_pl, source=pl%nplenc) - lencounter = pl%encounter_check(system, param%dt, 0) + lencounter = pl%encounter_check(param, system, param%dt, 0) if (system%tp%nbody > 0) then select type(tp => system%tp) class is (symba_tp) @@ -494,7 +495,7 @@ module subroutine symba_util_rearray_pl(self, system, param) allocate(levelg_orig_tp, source=tp%levelg) allocate(levelm_orig_tp, source=tp%levelm) allocate(nplenc_orig_tp, source=tp%nplenc) - lencounter = tp%encounter_check(system, param%dt, 0) + lencounter = tp%encounter_check(param, system, param%dt, 0) call move_alloc(levelg_orig_tp, tp%levelg) call move_alloc(levelm_orig_tp, tp%levelm) call move_alloc(nplenc_orig_tp, tp%nplenc) diff --git a/src/util/util_copy.f90 b/src/util/util_copy.f90 index 60d93e45e..f271c1c2e 100644 --- a/src/util/util_copy.f90 +++ b/src/util/util_copy.f90 @@ -15,7 +15,6 @@ module subroutine util_copy_encounter(self, source) self%nenc = n self%lvdotr(1:n) = source%lvdotr(1:n) self%status(1:n) = source%status(1:n) - self%kidx(1:n) = source%kidx(1:n) self%index1(1:n) = source%index1(1:n) self%index2(1:n) = source%index2(1:n) self%id1(1:n) = source%id1(1:n) diff --git a/src/util/util_get_energy_momentum.f90 b/src/util/util_get_energy_momentum.f90 index 677de5646..3474e2921 100644 --- a/src/util/util_get_energy_momentum.f90 +++ b/src/util/util_get_energy_momentum.f90 @@ -78,8 +78,12 @@ module subroutine util_get_energy_momentum_system(self, param) kespincb = 0.0_DP kespinpl(:) = 0.0_DP end if - - call util_get_energy_potential(npl, nplpl, pl%k_plpl, pl%lmask, cb%Gmass, pl%Gmass, pl%mass, pl%xb, system%pe) + + if (param%lflatten_interactions) then + call util_get_energy_potential_flat(npl, nplpl, pl%k_plpl, pl%lmask, cb%Gmass, pl%Gmass, pl%mass, pl%xb, system%pe) + else + call util_get_energy_potential_triangular(npl, pl%lmask, cb%Gmass, pl%Gmass, pl%mass, pl%xb, system%pe) + end if ! Potential energy from the oblateness term if (param%loblatecb) then @@ -108,7 +112,7 @@ module subroutine util_get_energy_momentum_system(self, param) end subroutine util_get_energy_momentum_system - subroutine util_get_energy_potential(npl, nplpl, k_plpl, lmask, GMcb, Gmass, mass, xb, pe) + subroutine util_get_energy_potential_flat(npl, nplpl, k_plpl, lmask, GMcb, Gmass, mass, xb, pe) !! author: David A. Minton !! !! Compute total system potential energy @@ -156,6 +160,41 @@ subroutine util_get_energy_potential(npl, nplpl, k_plpl, lmask, GMcb, Gmass, mas pe = sum(pepl(:), lstatpl(:)) + sum(pecb(1:npl), lmask(1:npl)) return - end subroutine util_get_energy_potential + end subroutine util_get_energy_potential_flat + + + subroutine util_get_energy_potential_triangular(npl, lmask, GMcb, Gmass, mass, xb, pe) + !! author: David A. Minton + !! + !! Compute total system potential energy + implicit none + ! Arguments + integer(I4B), intent(in) :: npl + logical, dimension(:), intent(in) :: lmask + real(DP), intent(in) :: GMcb + real(DP), dimension(:), intent(in) :: Gmass + real(DP), dimension(:), intent(in) :: mass + real(DP), dimension(:,:), intent(in) :: xb + real(DP), intent(out) :: pe + ! Internals + integer(I4B) :: i, j + real(DP), dimension(npl) :: pecb + + ! Do the central body potential energy component first + where(.not. lmask(1:npl)) + pecb(1:npl) = 0.0_DP + end where + + do concurrent(i = 1:npl, lmask(i)) + pecb(i) = -GMcb * mass(i) / norm2(xb(:,i)) + end do + + pe = sum(pecb(1:npl), lmask(1:npl)) + do concurrent(i = 1:npl, j = 1:npl, (j > i) .and. lmask(i) .and. lmask(j)) + pe = pe - (Gmass(i) * mass(j)) / norm2(xb(:, i) - xb(:, j)) + end do + + return + end subroutine util_get_energy_potential_triangular end submodule s_util_get_energy_momentum diff --git a/src/util/util_index.f90 b/src/util/util_index.f90 index 1ee60e400..d0dd83896 100644 --- a/src/util/util_index.f90 +++ b/src/util/util_index.f90 @@ -82,13 +82,15 @@ module subroutine util_index_eucl_plpl(self, param) npl = int(self%nbody, kind=I8B) associate(nplpl => self%nplpl) nplpl = (npl * (npl - 1) / 2) ! number of entries in a strict lower triangle, npl x npl - if (allocated(self%k_plpl)) deallocate(self%k_plpl) ! Reset the index array if it's been set previously - allocate(self%k_plpl(2, nplpl)) - do concurrent (i=1:npl, j=1:npl, j>i) - call util_index_eucl_ij_to_k(npl, i, j, k) - self%k_plpl(1, k) = i - self%k_plpl(2, k) = j - end do + if (param%lflatten_interactions) then + if (allocated(self%k_plpl)) deallocate(self%k_plpl) ! Reset the index array if it's been set previously + allocate(self%k_plpl(2, nplpl)) + do concurrent (i=1:npl, j=1:npl, j>i) + call util_index_eucl_ij_to_k(npl, i, j, k) + self%k_plpl(1, k) = i + self%k_plpl(2, k) = j + end do + end if end associate return diff --git a/src/util/util_spill.f90 b/src/util/util_spill.f90 index 16039c4da..e9d243cb6 100644 --- a/src/util/util_spill.f90 +++ b/src/util/util_spill.f90 @@ -373,7 +373,6 @@ module subroutine util_spill_encounter(self, discards, lspill_list, ldestructive associate(keeps => self) call util_spill(keeps%lvdotr, discards%lvdotr, lspill_list, ldestructive) call util_spill(keeps%status, discards%status, lspill_list, ldestructive) - call util_spill(keeps%kidx, discards%kidx, lspill_list, ldestructive) call util_spill(keeps%index1, discards%index1, lspill_list, ldestructive) call util_spill(keeps%index2, discards%index2, lspill_list, ldestructive) call util_spill(keeps%id1, discards%id1, lspill_list, ldestructive) diff --git a/src/whm/whm_kick.f90 b/src/whm/whm_kick.f90 index f831ce15b..c36c1ce23 100644 --- a/src/whm/whm_kick.f90 +++ b/src/whm/whm_kick.f90 @@ -32,7 +32,7 @@ module subroutine whm_kick_getacch_pl(self, system, param, t, lbeg) call whm_kick_getacch_ah1(cb, pl) call whm_kick_getacch_ah2(cb, pl) - call pl%accel_int() + call pl%accel_int(param) if (param%loblatecb) then call pl%accel_obl(system) @@ -84,13 +84,13 @@ module subroutine whm_kick_getacch_tp(self, system, param, t, lbeg) do concurrent(i = 1:ntp, tp%lmask(i)) tp%ah(:, i) = tp%ah(:, i) + ah0(:) end do - call tp%accel_int(pl%Gmass(1:npl), pl%xbeg(:, 1:npl), npl) + call tp%accel_int(param, pl%Gmass(1:npl), pl%xbeg(:, 1:npl), npl) else ah0(:) = whm_kick_getacch_ah0(pl%Gmass(1:npl), pl%xend(:, 1:npl), npl) do concurrent(i = 1:ntp, tp%lmask(i)) tp%ah(:, i) = tp%ah(:, i) + ah0(:) end do - call tp%accel_int(pl%Gmass(1:npl), pl%xend(:, 1:npl), npl) + call tp%accel_int(param, pl%Gmass(1:npl), pl%xend(:, 1:npl), npl) end if if (param%loblatecb) call tp%accel_obl(system)