From 25385d261d38befac20069309c53e90d5c6b95bc Mon Sep 17 00:00:00 2001 From: David A Minton Date: Thu, 16 Sep 2021 11:11:30 -0400 Subject: [PATCH 01/13] Setting the stage to have the option of using either the flat or triangular versions of the two big interaction loops. --- src/kick/kick.f90 | 57 ++++++++++++++++++-- src/modules/swiftest_classes.f90 | 13 ++++- src/symba/symba_encounter_check.f90 | 82 +++++++++++++++++++++++++++-- src/symba/symba_kick.f90 | 4 +- 4 files changed, 144 insertions(+), 12 deletions(-) diff --git a/src/kick/kick.f90 b/src/kick/kick.f90 index 45da5cb31..05603143b 100644 --- a/src/kick/kick.f90 +++ b/src/kick/kick.f90 @@ -13,7 +13,7 @@ module subroutine kick_getacch_int_pl(self) ! Arguments class(swiftest_pl), intent(inout) :: self !! Swiftest massive body object - call kick_getacch_int_all_pl(self%nbody, self%nplpl, self%k_plpl, self%xh, self%Gmass, self%radius, self%ah) + call kick_getacch_int_all_flat_pl(self%nbody, self%nplpl, self%k_plpl, self%xh, self%Gmass, self%radius, self%ah) return end subroutine kick_getacch_int_pl @@ -41,10 +41,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 +89,55 @@ 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, 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 !! Number of 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, x, Gmass, radius) & + !$omp lastprivate(rji2, rlim2, xr, yr, zr) & + !$omp reduction(+:ahi) & + !$omp reduction(-:ahj) + do i = 1, npl + 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) diff --git a/src/modules/swiftest_classes.f90 b/src/modules/swiftest_classes.f90 index 598cb4811..61f0de8ef 100644 --- a/src/modules/swiftest_classes.f90 +++ b/src/modules/swiftest_classes.f90 @@ -853,7 +853,7 @@ module subroutine kick_getacch_int_tp(self, GMpl, xhp, npl) 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 +862,16 @@ 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, x, Gmass, radius, acc) + implicit none + integer(I4B), intent(in) :: npl !! Number of 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/symba/symba_encounter_check.f90 b/src/symba/symba_encounter_check.f90 index e20e51ea9..fa5913507 100644 --- a/src/symba/symba_encounter_check.f90 +++ b/src/symba/symba_encounter_check.f90 @@ -2,10 +2,11 @@ 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 @@ -38,7 +39,80 @@ 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, loc_lvdotr, k_plplenc, 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, nplm + 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(:), allocatable, intent(out) :: loc_lvdotr + integer(I4B), dimension(:,:), allocatable, intent(out) :: k_plplenc + integer(I4B), intent(out) :: nenc + ! Internals + integer(I4B) :: i, j, nenci, j0, j1 + real(DP) :: xr, yr, zr, vxr, vyr, vzr, rhill1, rhill2 + logical, dimension(npl) :: lencounteri, loc_lvdotri + integer(I4B), dimension(npl) :: ind_arr + type lenctype + logical, dimension(:), allocatable :: loc_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 simd default(private) schedule(static)& + !$omp shared(npl, nplm, x, v, rhill, dt, irec, lenc) & + !$omp lastprivate(xr, yr, zr, vxr, vyr, vzr, rhill1, rhill2, lencounteri, loc_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), loc_lvdotri(j)) + end do + lenc(i)%nenc = count(lencounteri(i+1:npl)) + if (lenc(i)%nenc > 0) then + allocate(lenc(i)%loc_lvdotr(nenci), lenc(i)%index2(nenci)) + lenc(i)%loc_lvdotr(:) = pack(loc_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 simd + + associate(nenc_arr => lenc(:)%nenc) + nenc = sum(nenc_arr(1:nplm)) + end associate + if (nenc > 0) then + allocate(loc_lvdotr(nenc)) + allocate(k_plplenc(2,nenc)) + j0 = 1 + do i = 1, nplm + if (lenc(i)%nenc > 0) then + j1 = j0 + lenc(i)%nenc - 1 + loc_lvdotr(j0:j1) = lenc(i)%loc_lvdotr(:) + k_plplenc(1,j0:j1) = i + k_plplenc(2,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) @@ -68,7 +142,7 @@ module function symba_encounter_check_pl(self, system, dt, irec) result(lany_enc 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(:)) diff --git a/src/symba/symba_kick.f90 b/src/symba/symba_kick.f90 index 325ef5a45..97b52e96e 100644 --- a/src/symba/symba_kick.f90 +++ b/src/symba/symba_kick.f90 @@ -13,7 +13,7 @@ module subroutine symba_kick_getacch_int_pl(self) ! Arguments class(symba_pl), intent(inout) :: self - call kick_getacch_int_all_pl(self%nbody, self%nplplm, self%k_plpl, self%xh, self%Gmass, self%radius, self%ah) + call kick_getacch_int_all_flat_pl(self%nbody, self%nplplm, self%k_plpl, self%xh, self%Gmass, self%radius, self%ah) return end subroutine symba_kick_getacch_int_pl @@ -52,7 +52,7 @@ module subroutine symba_kick_getacch_pl(self, system, param, t, lbeg) 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) + 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 From b2401e0104ad92cf0827a3207df20e54648a52e8 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Thu, 16 Sep 2021 13:45:03 -0400 Subject: [PATCH 02/13] Added user parameter option to turn on or off the flattened interaction loops --- src/io/io.f90 | 9 +++++++-- src/modules/swiftest_classes.f90 | 1 + 2 files changed, 8 insertions(+), 2 deletions(-) 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/modules/swiftest_classes.f90 b/src/modules/swiftest_classes.f90 index 61f0de8ef..287654614 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 From 7b784f6aeb71d6714fe0e90f19c4063654d176f9 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Thu, 16 Sep 2021 13:53:50 -0400 Subject: [PATCH 03/13] Added param to the argument list of the interaction acceleration methods so that the choice of whether to use the flattened matrix or not is availble --- examples/symba_mars_disk/Untitled.ipynb | 6 - examples/symba_mars_disk/param.in | 13 +- examples/symba_mars_disk/param.real8.in | 33 + examples/symba_mars_disk/testnetcdf.ipynb | 5414 +-------------------- src/helio/helio_kick.f90 | 6 +- src/kick/kick.f90 | 16 +- src/modules/swiftest_classes.f90 | 16 +- src/modules/symba_classes.f90 | 5 +- src/symba/symba_kick.f90 | 5 +- src/whm/whm_kick.f90 | 6 +- 10 files changed, 143 insertions(+), 5377 deletions(-) delete mode 100644 examples/symba_mars_disk/Untitled.ipynb create mode 100644 examples/symba_mars_disk/param.real8.in 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/kick/kick.f90 b/src/kick/kick.f90 index 05603143b..de99a9bb0 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,7 +11,8 @@ 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_flat_pl(self%nbody, self%nplpl, self%k_plpl, self%xh, self%Gmass, self%radius, self%ah) @@ -19,7 +20,7 @@ module subroutine kick_getacch_int_pl(self) 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 +29,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 diff --git a/src/modules/swiftest_classes.f90 b/src/modules/swiftest_classes.f90 index 287654614..a88fb7d7a 100644 --- a/src/modules/swiftest_classes.f90 +++ b/src/modules/swiftest_classes.f90 @@ -841,17 +841,19 @@ 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_flat_pl(npl, nplpl, k_plpl, x, Gmass, radius, acc) diff --git a/src/modules/symba_classes.f90 b/src/modules/symba_classes.f90 index 40267e298..f805aa2c4 100644 --- a/src/modules/symba_classes.f90 +++ b/src/modules/symba_classes.f90 @@ -358,9 +358,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/symba/symba_kick.f90 b/src/symba/symba_kick.f90 index 97b52e96e..2011ee0c2 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,7 +11,8 @@ 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_flat_pl(self%nbody, self%nplplm, self%k_plpl, self%xh, self%Gmass, self%radius, self%ah) 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) From ea1f5f87369a325278c56959b95fb7d05b6a1165 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Thu, 16 Sep 2021 14:08:13 -0400 Subject: [PATCH 04/13] Added param to argument list of encounter methods so that the choice of whether to use the flattened matrix is available --- src/modules/rmvs_classes.f90 | 11 ++++--- src/modules/symba_classes.f90 | 50 +++++++++++++---------------- src/rmvs/rmvs_encounter_check.f90 | 9 +++--- src/rmvs/rmvs_step.f90 | 4 +-- src/symba/symba_encounter_check.f90 | 37 +++++++++++---------- src/symba/symba_step.f90 | 4 +-- src/symba/symba_util.f90 | 4 +-- 7 files changed, 59 insertions(+), 60 deletions(-) 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/symba_classes.f90 b/src/modules/symba_classes.f90 index f805aa2c4..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) 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/symba/symba_encounter_check.f90 b/src/symba/symba_encounter_check.f90 index fa5913507..82c4c1890 100644 --- a/src/symba/symba_encounter_check.f90 +++ b/src/symba/symba_encounter_check.f90 @@ -115,17 +115,18 @@ subroutine symba_encounter_check_all_triangular(npl, nplm, x, v, rhill, dt, ire 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 @@ -185,7 +186,7 @@ module function symba_encounter_check_pl(self, system, dt, irec) result(lany_enc 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. @@ -194,11 +195,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 @@ -271,17 +273,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 @@ -343,7 +346,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_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..752028cf7 100644 --- a/src/symba/symba_util.f90 +++ b/src/symba/symba_util.f90 @@ -486,7 +486,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 +494,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) From cc975089f89a19c3bb504fb0f67a57b166bf4a30 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Thu, 16 Sep 2021 15:01:05 -0400 Subject: [PATCH 05/13] Added ability to switch between flat and triangular versions of double loops --- src/kick/kick.f90 | 15 ++- src/modules/swiftest_classes.f90 | 5 +- src/symba/symba_encounter_check.f90 | 136 +++++++++++++++----------- src/symba/symba_kick.f90 | 11 ++- src/symba/symba_util.f90 | 19 ++-- src/util/util_get_energy_momentum.f90 | 47 ++++++++- src/util/util_index.f90 | 16 +-- 7 files changed, 159 insertions(+), 90 deletions(-) diff --git a/src/kick/kick.f90 b/src/kick/kick.f90 index de99a9bb0..ea21c42e8 100644 --- a/src/kick/kick.f90 +++ b/src/kick/kick.f90 @@ -14,7 +14,11 @@ module subroutine kick_getacch_int_pl(self, param) 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_flat_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 @@ -94,7 +98,7 @@ module subroutine kick_getacch_int_all_flat_pl(npl, nplpl, k_plpl, x, Gmass, rad end subroutine kick_getacch_int_all_flat_pl - module subroutine kick_getacch_int_all_triangular_pl(npl, x, Gmass, radius, acc) + 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. @@ -103,7 +107,8 @@ module subroutine kick_getacch_int_all_triangular_pl(npl, x, Gmass, radius, acc) !! 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 !! Number of massive bodies + 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 @@ -118,11 +123,11 @@ module subroutine kick_getacch_int_all_triangular_pl(npl, x, Gmass, radius, acc) ahj(:,:) = 0.0_DP !$omp parallel do default(private) schedule(static)& - !$omp shared(npl, x, Gmass, radius) & + !$omp shared(npl, nplm, x, Gmass, radius) & !$omp lastprivate(rji2, rlim2, xr, yr, zr) & !$omp reduction(+:ahi) & !$omp reduction(-:ahj) - do i = 1, npl + do i = 1, nplm do concurrent(j = i+1:npl) xr = x(1, j) - x(1, i) yr = x(2, j) - x(2, i) diff --git a/src/modules/swiftest_classes.f90 b/src/modules/swiftest_classes.f90 index a88fb7d7a..f4f9b2698 100644 --- a/src/modules/swiftest_classes.f90 +++ b/src/modules/swiftest_classes.f90 @@ -867,9 +867,10 @@ module subroutine kick_getacch_int_all_flat_pl(npl, nplpl, k_plpl, x, Gmass, rad real(DP), dimension(:,:), intent(inout) :: acc !! Acceleration vector array end subroutine kick_getacch_int_all_flat_pl - module subroutine kick_getacch_int_all_triangular_pl(npl, x, Gmass, radius, acc) + module subroutine kick_getacch_int_all_triangular_pl(npl, nplm, x, Gmass, radius, acc) implicit none - integer(I4B), intent(in) :: npl !! Number of massive bodies + 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 diff --git a/src/symba/symba_encounter_check.f90 b/src/symba/symba_encounter_check.f90 index 82c4c1890..1df9ca03c 100644 --- a/src/symba/symba_encounter_check.f90 +++ b/src/symba/symba_encounter_check.f90 @@ -8,13 +8,15 @@ subroutine symba_encounter_check_all_flat(nplplm, k_plpl, x, v, rhill, dt, irec !! 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 @@ -42,36 +44,40 @@ subroutine symba_encounter_check_all_flat(nplplm, k_plpl, x, v, rhill, dt, irec end subroutine symba_encounter_check_all_flat - subroutine symba_encounter_check_all_triangular(npl, nplm, x, v, rhill, dt, irec, loc_lvdotr, k_plplenc, nenc) + 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, nplm - 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(:), allocatable, intent(out) :: loc_lvdotr - integer(I4B), dimension(:,:), allocatable, intent(out) :: k_plplenc - integer(I4B), intent(out) :: nenc + 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, loc_lvdotri + logical, dimension(npl) :: lencounteri, lvdotri integer(I4B), dimension(npl) :: ind_arr type lenctype - logical, dimension(:), allocatable :: loc_lvdotr + 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 simd default(private) schedule(static)& + !$omp parallel do default(private) schedule(static)& !$omp shared(npl, nplm, x, v, rhill, dt, irec, lenc) & - !$omp lastprivate(xr, yr, zr, vxr, vyr, vzr, rhill1, rhill2, lencounteri, loc_lvdotri) + !$omp firstprivate(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) @@ -82,30 +88,31 @@ subroutine symba_encounter_check_all_triangular(npl, nplm, x, v, rhill, dt, ire 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), loc_lvdotri(j)) + call symba_encounter_check_one(xr, yr, zr, vxr, vyr, vzr, rhill1, rhill2, dt, irec, lencounteri(j), lvdotri(j)) end do lenc(i)%nenc = count(lencounteri(i+1:npl)) if (lenc(i)%nenc > 0) then - allocate(lenc(i)%loc_lvdotr(nenci), lenc(i)%index2(nenci)) - lenc(i)%loc_lvdotr(:) = pack(loc_lvdotri(i+1:npl), lencounteri(i+1:npl)) + allocate(lenc(i)%lvdotr(nenci), lenc(i)%index2(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 simd + !$omp end parallel do associate(nenc_arr => lenc(:)%nenc) nenc = sum(nenc_arr(1:nplm)) end associate if (nenc > 0) then - allocate(loc_lvdotr(nenc)) - allocate(k_plplenc(2,nenc)) + 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 - loc_lvdotr(j0:j1) = lenc(i)%loc_lvdotr(:) - k_plplenc(1,j0:j1) = i - k_plplenc(2,j0:j1) = lenc(i)%index2(:) + lvdotr(j0:j1) = lenc(i)%lvdotr(:) + index1(j0:j1) = i + index2(j0:j1) = lenc(i)%index2(:) j0 = j1 + 1 end if end do @@ -131,25 +138,26 @@ module function symba_encounter_check_pl(self, param, system, dt, irec) result(l 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_flat(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)) @@ -160,28 +168,38 @@ module function symba_encounter_check_pl(self, param, system, dt, irec) result(l 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) + call move_alloc(lvdotr, plplenc_list%lvdotr) + call move_alloc(index1, plplenc_list%index1) + call move_alloc(index2, plplenc_list%index2) + end if + + if (lany_encounter) then + 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 if end associate + return end function symba_encounter_check_pl diff --git a/src/symba/symba_kick.f90 b/src/symba/symba_kick.f90 index 2011ee0c2..414c89805 100644 --- a/src/symba/symba_kick.f90 +++ b/src/symba/symba_kick.f90 @@ -14,7 +14,11 @@ module subroutine symba_kick_getacch_int_pl(self, param) 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_flat_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 @@ -44,15 +48,15 @@ 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_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 @@ -63,7 +67,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_util.f90 b/src/symba/symba_util.f90 index 752028cf7..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 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 From 1e9794aa67349497e38483e18fe4c8ceb44cd901 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Thu, 16 Sep 2021 22:56:59 -0400 Subject: [PATCH 06/13] Added missing variable to openmp list --- src/kick/kick.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/kick/kick.f90 b/src/kick/kick.f90 index ea21c42e8..766e829d6 100644 --- a/src/kick/kick.f90 +++ b/src/kick/kick.f90 @@ -168,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 From d701f5d1afa626553c08881fe19530c5c7f0911d Mon Sep 17 00:00:00 2001 From: David A Minton Date: Thu, 16 Sep 2021 23:12:35 -0400 Subject: [PATCH 07/13] Set flatten to false by default --- src/io/io.f90 | 2 +- src/modules/swiftest_classes.f90 | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/io/io.f90 b/src/io/io.f90 index 51d10561d..b951130d1 100644 --- a/src/io/io.f90 +++ b/src/io/io.f90 @@ -603,7 +603,7 @@ module subroutine io_param_reader(self, unit, iotype, v_list, iostat, iomsg) 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. + if (param_value == "YES" .or. param_value == 'T') param%lflatten_interactions = .false. case ("FIRSTKICK") call io_toupper(param_value) if (param_value == "NO" .or. param_value == 'F') param%lfirstkick = .false. diff --git a/src/modules/swiftest_classes.f90 b/src/modules/swiftest_classes.f90 index f4f9b2698..92d671e9b 100644 --- a/src/modules/swiftest_classes.f90 +++ b/src/modules/swiftest_classes.f90 @@ -133,7 +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) + logical :: lflatten_interactions = .false. !! 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 From d22e6074e3e6fc4604f3061ef510278c357f8355 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Thu, 16 Sep 2021 23:18:53 -0400 Subject: [PATCH 08/13] Fixed issue where encounters were not being recorded correctly in triangular mode --- src/modules/swiftest_classes.f90 | 1 - src/setup/setup.f90 | 3 --- src/symba/symba_encounter_check.f90 | 10 ++++++---- src/symba/symba_kick.f90 | 3 ++- src/util/util_copy.f90 | 1 - src/util/util_spill.f90 | 1 - 6 files changed, 8 insertions(+), 11 deletions(-) diff --git a/src/modules/swiftest_classes.f90 b/src/modules/swiftest_classes.f90 index 92d671e9b..0842ab154 100644 --- a/src/modules/swiftest_classes.f90 +++ b/src/modules/swiftest_classes.f90 @@ -433,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 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_encounter_check.f90 b/src/symba/symba_encounter_check.f90 index 1df9ca03c..cd6b81577 100644 --- a/src/symba/symba_encounter_check.f90 +++ b/src/symba/symba_encounter_check.f90 @@ -173,9 +173,12 @@ module function symba_encounter_check_pl(self, param, system, dt, irec) result(l else nplm = pl%nplm call symba_encounter_check_all_triangular(npl, nplm, pl%xh, pl%vh, pl%rhill, dt, irec, lvdotr, index1, index2, nenc) - call move_alloc(lvdotr, plplenc_list%lvdotr) - call move_alloc(index1, plplenc_list%index1) - call move_alloc(index2, plplenc_list%index2) + lany_encounter = nenc > 0 + if (lany_encounter) then + 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 @@ -183,7 +186,6 @@ module function symba_encounter_check_pl(self, param, system, dt, irec) result(l 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 diff --git a/src/symba/symba_kick.f90 b/src/symba/symba_kick.f90 index 414c89805..fd0354612 100644 --- a/src/symba/symba_kick.f90 +++ b/src/symba/symba_kick.f90 @@ -56,7 +56,8 @@ module subroutine symba_kick_getacch_pl(self, system, param, t, lbeg) 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)) + 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 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_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) From 22681730a32014d8d5f6b86665eb54bc02baa0aa Mon Sep 17 00:00:00 2001 From: David A Minton Date: Thu, 16 Sep 2021 23:25:39 -0400 Subject: [PATCH 09/13] Switched back to default on for flatten_interactions --- src/io/io.f90 | 2 +- src/modules/swiftest_classes.f90 | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/io/io.f90 b/src/io/io.f90 index b951130d1..51d10561d 100644 --- a/src/io/io.f90 +++ b/src/io/io.f90 @@ -603,7 +603,7 @@ module subroutine io_param_reader(self, unit, iotype, v_list, iostat, iomsg) if (param_value == "YES" .or. param_value == 'T') param%ltides = .true. case ("FLATTEN_INTERACTIONS") call io_toupper(param_value) - if (param_value == "YES" .or. param_value == 'T') param%lflatten_interactions = .false. + 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. diff --git a/src/modules/swiftest_classes.f90 b/src/modules/swiftest_classes.f90 index 0842ab154..a1f60c0dd 100644 --- a/src/modules/swiftest_classes.f90 +++ b/src/modules/swiftest_classes.f90 @@ -133,7 +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 = .false. !! Use the flattened upper triangular matrix for pl-pl interactions (turning this on improves the speed but uses more memory) + 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 From 85b296b3b30e867cb8d0d1cae1af212c8a66dc09 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Fri, 17 Sep 2021 00:17:37 -0400 Subject: [PATCH 10/13] simplified encounter saves --- src/symba/symba_encounter_check.f90 | 10 ++++------ 1 file changed, 4 insertions(+), 6 deletions(-) diff --git a/src/symba/symba_encounter_check.f90 b/src/symba/symba_encounter_check.f90 index cd6b81577..1b107d6a2 100644 --- a/src/symba/symba_encounter_check.f90 +++ b/src/symba/symba_encounter_check.f90 @@ -175,8 +175,8 @@ module function symba_encounter_check_pl(self, param, system, dt, irec) result(l 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 move_alloc(lvdotr, plplenc_list%lvdotr) - call move_alloc(index1, plplenc_list%index1) + 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 @@ -185,10 +185,8 @@ module function symba_encounter_check_pl(self, param, system, dt, irec) result(l 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%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%id1(k) = pl%id(i) + plplenc_list%id2(k) = pl%id(j) plplenc_list%level(k) = irec pl%lencounter(i) = .true. pl%lencounter(j) = .true. From 8ae52155f67390f8aa4f2287b917de1080d04fc5 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Fri, 17 Sep 2021 11:19:45 -0400 Subject: [PATCH 11/13] Fixed bugs that prevented the encounter list from being properly built in triangular mode --- src/symba/symba_encounter_check.f90 | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/src/symba/symba_encounter_check.f90 b/src/symba/symba_encounter_check.f90 index 1b107d6a2..743ef677e 100644 --- a/src/symba/symba_encounter_check.f90 +++ b/src/symba/symba_encounter_check.f90 @@ -75,8 +75,7 @@ subroutine symba_encounter_check_all_triangular(npl, nplm, x, v, rhill, dt, ire ind_arr(:) = [(i, i = 1, npl)] !$omp parallel do default(private) schedule(static)& - !$omp shared(npl, nplm, x, v, rhill, dt, irec, lenc) & - !$omp firstprivate(ind_arr) & + !$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) @@ -90,9 +89,10 @@ subroutine symba_encounter_check_all_triangular(npl, nplm, x, v, rhill, dt, ire rhill2 = rhill(j) call symba_encounter_check_one(xr, yr, zr, vxr, vyr, vzr, rhill1, rhill2, dt, irec, lencounteri(j), lvdotri(j)) end do - lenc(i)%nenc = count(lencounteri(i+1:npl)) - if (lenc(i)%nenc > 0) then + 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 @@ -175,6 +175,7 @@ module function symba_encounter_check_pl(self, param, system, dt, irec) result(l 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) From 9bec1497ac7b6ea0090e50939194d524ed688097 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Fri, 17 Sep 2021 18:21:24 -0400 Subject: [PATCH 12/13] Fixed bug resulting from accidentally removing the status flag setting for encounters --- src/symba/symba_collision.f90 | 1 - src/symba/symba_encounter_check.f90 | 1 + 2 files changed, 1 insertion(+), 1 deletion(-) 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 743ef677e..3ba6b4adf 100644 --- a/src/symba/symba_encounter_check.f90 +++ b/src/symba/symba_encounter_check.f90 @@ -188,6 +188,7 @@ module function symba_encounter_check_pl(self, param, system, dt, irec) result(l 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. From 41265b7ad5ecff1f74ca7741325fe3b638666ac3 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Mon, 20 Sep 2021 07:47:35 -0400 Subject: [PATCH 13/13] Updated flags with more optimization parameters --- Makefile.Defines | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) 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)