Merge branch 'OOPSymba' into OOPTides

examples/rmvs_swifter_comparison/1pl_1tp_encounter/swiftest_vs_swifter.ipynb

@@ -81,8 +81,8 @@
     {
      "data": {
       "text/plain": [
-       "[<matplotlib.lines.Line2D at 0x2acdc3dd7450>,\n",
-       " <matplotlib.lines.Line2D at 0x2acdc3db7c10>]"
+       "[<matplotlib.lines.Line2D at 0x2b768e2c9d10>,\n",
+       " <matplotlib.lines.Line2D at 0x2b768e2b9c90>]"
       ]
      },
      "execution_count": 6,

examples/whm_gr_test/init_cond.py

@@ -13,7 +13,7 @@
 sim.param['DU2M'] = swiftest.AU2M
 sim.param['T0'] = 0.0
 sim.param['DT'] = 0.25 * swiftest.JD2S / swiftest.YR2S
-sim.param['TSTOP'] = 100.0
+sim.param['TSTOP'] = 1000.0
 sim.param['ISTEP_OUT']  = 1461
 sim.param['ISTEP_DUMP'] = 1461
 sim.param['CHK_QMIN_COORD'] = "HELIO"

examples/whm_gr_test/param.swifter.in

@@ -1,6 +1,6 @@
 ! VERSION        Swifter parameter file converted from Swiftest
 T0               0.0
-TSTOP            100.0
+TSTOP            1000.0
 DT               0.0006844626967830253
 ISTEP_OUT        1461
 ISTEP_DUMP       1461

examples/whm_gr_test/param.swiftest.in

@@ -1,6 +1,6 @@
 ! VERSION        Swiftest parameter input
 T0               0.0
-TSTOP            100.0
+TSTOP            1000.0
 DT               0.0006844626967830253
 ISTEP_OUT        1461
 ISTEP_DUMP       1461

examples/whm_gr_test/pl.swifter.in

@@ -2,35 +2,35 @@
 0 39.476926408897625196
 0.0 0.0 0.0
 0.0 0.0 0.0
-1 6.5537098095653139645e-06 0.0014751243077781048702
+1 6.5537098095653139645e-06 0.0014751234419554511911
 1.6306381826061645943e-05
-0.33206272695596028566 0.07436707001147663254 -0.02438290851908785084
--4.2340114788918336805 10.486553514018327622 1.2453138107251555947
-2 9.663313399581537916e-05 0.006759104275397271956
+0.13267502226188271353 0.2786606257975073886 0.010601098875389479426
+-11.331978934667442676 4.8184460126705647045 1.4332264599878684131
+2 9.663313399581537916e-05 0.00675908960945781479
 4.0453784346544178454e-05
--0.7188115337296047125 -0.0118554711069603201795 0.041316403191083782287
-0.07826338813583945357 -7.419533988988633545 -0.10634201014368884618
-3 0.000120026935827952453094 0.010044787321379672528
+-0.69398700025820403425 -0.19235393648106968723 0.03740673057980103272
+1.9245789988923785786 -7.1528261190002948057 -0.20922405362759749996
+3 0.000120026935827952453094 0.010044837538502923644
 4.25875607065040958e-05
-0.35677088372527121507 -0.95189300879814897627 4.4027442504036787155e-05
-5.7819217550992820422 2.18192814489641851 -0.00012230072278352209966
-4 1.2739802010675941456e-05 0.007246743835971885302
+0.49463573470256239073 -0.8874896493821613497 4.051630875713834232e-05
+5.386704768180099809 3.0357508899436080915 -0.00016218409216515533796
+4 1.2739802010675941456e-05 0.0072467236860282326973
 2.265740805092889601e-05
--1.5233712071242269115 0.6723825347339112968 0.051459143378398922164
--1.8728417739956807141 -4.239719661832373223 -0.042909557750301418264
-5 0.037692251088985676735 0.35527126534549128905
+-1.5655322071100350456 0.56626121192188216824 0.050269397991054412533
+-1.5477080637857006753 -4.370087697214287981 -0.05361768768801557225
+5 0.037692251088985676735 0.35527094075555771578
 0.00046732617030490929307
-4.049944927347420176 -2.9910878677758190314 -0.078187280837353656526
-1.6060801375519682711 2.349356876761497338 -0.045690062807172619064
-6 0.011285899820091272997 0.4376527512949726007
+4.0891378954287338487 -2.9329188614380639066 -0.07930573161132697946
+1.575024788882753283 2.3719591091996699917 -0.045089307261129988257
+6 0.011285899820091272997 0.43765464106459166412
 0.00038925687730393611812
-6.298929503477405767 -7.706413024510769816 -0.11669919842191249504
-1.4661378456572359413 1.2872251175075805794 -0.08070991686100478242
-7 0.0017236589478267730203 0.4695362423191493196
+6.3349788609660162564 -7.674600716671800882 -0.11868650931385750502
+1.4598618704191345578 1.2948691245181617393 -0.080593167691228835176
+7 0.0017236589478267730203 0.46956055286931676728
 0.00016953449859497231466
-14.856082147529010129 13.007589275314199284 -0.14417795763685259391
--0.9554310497290159123 1.0161753499437922057 0.016099529164307530124
-8 0.0020336100526728302319 0.7812870996943599397
+14.832516206189200858 13.032608531076540714 -0.14378102535616668622
+-0.9573374666934839659 1.014553546383260322 0.016118112341773867214
+8 0.0020336100526728302319 0.7813163071687303693
 0.000164587904124493665
-29.55744967800954015 -4.629377558152945049 -0.58590957207831262377
-0.17162147939801157335 1.1422848961108499101 -0.027445465472921385952
+29.561664938083289655 -4.6012285192418387325 -0.586585578731106283
+0.17051705220469790965 1.1424784769020628332 -0.027423757798549895085

examples/whm_gr_test/pl.swiftest.in

@@ -1,33 +1,33 @@
 8
 1 6.5537098095653139645e-06
 1.6306381826061645943e-05
-0.33206272695596028566 0.07436707001147663254 -0.02438290851908785084
--4.2340114788918336805 10.486553514018327622 1.2453138107251555947
+0.13267502226188271353 0.2786606257975073886 0.010601098875389479426
+-11.331978934667442676 4.8184460126705647045 1.4332264599878684131
 2 9.663313399581537916e-05
 4.0453784346544178454e-05
--0.7188115337296047125 -0.0118554711069603201795 0.041316403191083782287
-0.07826338813583945357 -7.419533988988633545 -0.10634201014368884618
+-0.69398700025820403425 -0.19235393648106968723 0.03740673057980103272
+1.9245789988923785786 -7.1528261190002948057 -0.20922405362759749996
 3 0.000120026935827952453094
 4.25875607065040958e-05
-0.35677088372527121507 -0.95189300879814897627 4.4027442504036787155e-05
-5.7819217550992820422 2.18192814489641851 -0.00012230072278352209966
+0.49463573470256239073 -0.8874896493821613497 4.051630875713834232e-05
+5.386704768180099809 3.0357508899436080915 -0.00016218409216515533796
 4 1.2739802010675941456e-05
 2.265740805092889601e-05
--1.5233712071242269115 0.6723825347339112968 0.051459143378398922164
--1.8728417739956807141 -4.239719661832373223 -0.042909557750301418264
+-1.5655322071100350456 0.56626121192188216824 0.050269397991054412533
+-1.5477080637857006753 -4.370087697214287981 -0.05361768768801557225
 5 0.037692251088985676735
 0.00046732617030490929307
-4.049944927347420176 -2.9910878677758190314 -0.078187280837353656526
-1.6060801375519682711 2.349356876761497338 -0.045690062807172619064
+4.0891378954287338487 -2.9329188614380639066 -0.07930573161132697946
+1.575024788882753283 2.3719591091996699917 -0.045089307261129988257
 6 0.011285899820091272997
 0.00038925687730393611812
-6.298929503477405767 -7.706413024510769816 -0.11669919842191249504
-1.4661378456572359413 1.2872251175075805794 -0.08070991686100478242
+6.3349788609660162564 -7.674600716671800882 -0.11868650931385750502
+1.4598618704191345578 1.2948691245181617393 -0.080593167691228835176
 7 0.0017236589478267730203
 0.00016953449859497231466
-14.856082147529010129 13.007589275314199284 -0.14417795763685259391
--0.9554310497290159123 1.0161753499437922057 0.016099529164307530124
+14.832516206189200858 13.032608531076540714 -0.14378102535616668622
+-0.9573374666934839659 1.014553546383260322 0.016118112341773867214
 8 0.0020336100526728302319
 0.000164587904124493665
-29.55744967800954015 -4.629377558152945049 -0.58590957207831262377
-0.17162147939801157335 1.1422848961108499101 -0.027445465472921385952
+29.561664938083289655 -4.6012285192418387325 -0.586585578731106283
+0.17051705220469790965 1.1424784769020628332 -0.027423757798549895085

src/drift/drift.f90

@@ -9,6 +9,81 @@
    integer(I2B), parameter :: NLAG2    = 40
 
 contains
+
+   module subroutine drift_body(self, system, param, dt, mask)
+      !! author: David A. Minton
+      !!
+      !! Loop bodies and call Danby drift routine on the heliocentric position and velocities.
+      !!
+      !! Adapted from Hal Levison's Swift routine drift_tp.f 
+      !! Adapted from David E. Kaufmann's Swifter routine whm_drift_tp.f90
+      implicit none
+      ! Arguments
+      class(swiftest_body),         intent(inout) :: self   !! Swiftest test particle data structure
+      class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object
+      class(swiftest_parameters),   intent(in)    :: param  !! Current run configuration parameters 
+      real(DP),                     intent(in)    :: dt     !! Stepsize
+      logical, dimension(:),        intent(in)    :: mask   !! Logical mask of size self%nbody that determines which bodies to drift.
+      ! Internals
+      integer(I4B)                              :: i   
+      integer(I4B), dimension(:), allocatable   :: iflag
+
+      associate(n => self%nbody)
+         allocate(iflag(n))
+         iflag(:) = 0
+         call drift_all(self%mu, self%xh, self%vh, self%nbody, param, dt, mask, iflag)
+         if (any(iflag(1:n) /= 0)) then
+            where(iflag(1:n) /= 0) self%status(1:n) = DISCARDED_DRIFTERR
+            do i = 1, n
+               if (iflag(i) /= 0) write(*, *) " Body ", self%id(i), " lost due to error in Danby drift"
+            end do
+         end if
+      end associate
+
+      return
+   end subroutine drift_body
+
+   module pure subroutine drift_all(mu, x, v, n, param, dt, mask, iflag)
+      !! author: David A. Minton
+      !!
+      !! Loop bodies and call Danby drift routine on all bodies for the given position and velocity vector.
+      !!
+      !! Adapted from Hal Levison's Swift routine drift_tp.f 
+      !! Adapted from David E. Kaufmann's Swifter routine whm_drift_tp.f9
+      implicit none
+      ! Arguments
+      real(DP), dimension(:),     intent(in)    :: mu    !! Vector of gravitational constants
+      real(DP), dimension(:,:),   intent(inout) :: x, v  !! Position and velocity vectors
+      integer(I4B),               intent(in)    :: n     !! number of bodies
+      class(swiftest_parameters), intent(in)    :: param  !! Current run configuration parameters
+      real(DP),                   intent(in)    :: dt    !! Stepsize
+      logical, dimension(:),      intent(in)    :: mask  !! Logical mask of size self%nbody that determines which bodies to drift.
+      integer(I4B), dimension(:), intent(out)   :: iflag !! Vector of error flags. 0 means no problem
+      ! Internals
+      integer(I4B)                              :: i   
+      real(DP)                                  :: energy, vmag2, rmag  !! Variables used in GR calculation
+      real(DP), dimension(:), allocatable       :: dtp
+
+      if (n == 0) return
+
+      allocate(dtp(n))
+      if (param%lgr) then
+         do concurrent(i = 1:n, mask(i))
+            rmag = norm2(x(:, i))
+            vmag2 = dot_product(v(:, i), v(:, i))
+            energy = 0.5_DP * vmag2 - mu(i) / rmag
+            dtp(i) = dt * (1.0_DP + 3 * param%inv_c2 * energy)
+         end do
+      else
+         where(mask(1:n)) dtp(1:n) = dt
+      end if 
+      do concurrent(i = 1:n, mask(i))
+         call drift_one(mu(i), x(1,i), x(2,i), x(3,i), v(1,i), v(2,i), v(3,i), dtp(i), iflag(i))
+      end do
+
+      return
+   end subroutine drift_all
+
    module pure elemental subroutine drift_one(mu, px, py, pz, vx, vy, vz, dt, iflag) 
       !! author: The Purdue Swiftest Team - David A. Minton, Carlisle A. Wishard, Jennifer L.L. Pouplin, and Jacob R. Elliott
       !!
@@ -406,4 +481,5 @@ pure subroutine drift_kepu_stumpff(x, c0, c1, c2, c3)
       return
    end subroutine drift_kepu_stumpff
 
+
 end submodule drift_implementation

src/eucl/eucl.f90

@@ -14,64 +14,24 @@ module subroutine eucl_dist_index_plpl(self)
       ! Arguments
       class(swiftest_pl),             intent(inout) :: self  !! Swiftest massive body objec
       ! Internals
-      integer(I4B) :: i, j, k, kp, p
+      integer(I8B) :: i, j, counter, npl
 
-      associate(npl => self%nbody, num_comparisons => self%num_comparisons, k_eucl => self%k_eucl)
-         num_comparisons = (npl * (npl - 1) / 2) ! number of entries in a strict lower triangle, nplm x npl, minus first column
-         if (allocated(self%k_eucl)) deallocate(self%k_eucl) ! Reset the index array if it's been set previously
-         if (allocated(self%irij3)) deallocate(self%irij3)  
-         allocate(self%k_eucl(2, num_comparisons))
-         allocate(self%irij3(num_comparisons))
-         !do concurrent(k = 1:num_comparisons) !shared(num_comparisons, k_eucl, npl) local(kp, i, j, p)
-         do k = 1, num_comparisons
-            kp = npl * (npl - 1) / 2 - k
-            p = floor((sqrt(1._DP + 8 * kp) - 1._DP) / 2._DP)
-            i = k - (npl - 1) * (npl - 2) / 2 + p * (p + 1) / 2 + 1
-            j = npl - 1 - p 
-            k_eucl(1, k) = i 
-            k_eucl(2, k) = j
+      npl = int(self%nbody, kind=I8B)
+      associate(nplpl => self%nplpl)
+         nplpl = (npl * (npl - 1) / 2) ! number of entries in a strict lower triangle, nplm x npl, minus first column
+         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 i = 1, npl
+            counter = (i - 1_I8B) * npl - i * (i - 1_I8B) / 2_I8B + 1_I8B
+            do j = i + 1_I8B, npl
+               self%k_plpl(1, counter) = i
+               self%k_plpl(2, counter) = j
+               counter = counter + 1_I8B
+            end do
          end do
       end associate
       return
 
    end subroutine eucl_dist_index_plpl
 
-   module subroutine eucl_dist_index_pltp(self, pl)
-      !! author: Jacob R. Elliott and David A. Minton
-      !!
-      !! Turns i,j indices into k index for use in the Euclidean distance matrix
-      !!
-      !! Reference:
-      !!
-      !!    Mélodie Angeletti, Jean-Marie Bonny, Jonas Koko. Parallel Euclidean distance matrix computation on big datasets *. 
-      !!       2019. hal-0204751      implicit none
-      class(swiftest_tp),             intent(inout) :: self  !! Swiftest test particle object
-      class(swiftest_pl),             intent(inout) :: pl    !! Swiftest massive body object
-   end subroutine eucl_dist_index_pltp
-
-   module subroutine eucl_irij3_plpl(self)
-      !! author: Jacob R. Elliott and David A. Minton
-      !!
-      !! Efficient parallel loop-blocking algrorithm for evaluating the Euclidean distance matrix for planet-planet
-      implicit none
-      ! Arguments
-      class(swiftest_pl),             intent(inout) :: self  !! Swiftest massive body object
-      ! Internals
-      integer(I4B) :: k, i, j
-      real(DP), dimension(NDIM) :: dx
-      real(DP) :: rji2
-
-      associate(k_eucl => self%k_eucl, xh => self%xh, irij3 => self%irij3, nk => self%num_comparisons)
-         do k = 1, nk
-            i = k_eucl(1, k)
-            j = k_eucl(2, k)
-            dx(:) = xh(:, j) - xh(:, i)
-            rji2  = dot_product(dx(:), dx(:))
-            irij3(k) = 1.0_DP / (rji2 * sqrt(rji2))
-         end do
-      end associate
-
-      return
-   end subroutine eucl_irij3_plpl
-
  end submodule s_eucl

src/gr/gr.f90

@@ -108,7 +108,7 @@ module pure subroutine gr_pv2vh_body(self, param)
       implicit none
       ! Arguments
       class(swiftest_body),       intent(inout) :: self  !! Swiftest particle object
-      class(swiftest_parameters), intent(in)    :: param !! Current run configuration parameters of on parameters 
+      class(swiftest_parameters), intent(in)    :: param !! Current run configuration parameters 
       ! Internals
       integer(I4B)                              :: i
       real(DP), dimension(:,:), allocatable     :: vh    !! Temporary holder of pseudovelocity for in-place conversion
@@ -207,7 +207,7 @@ module pure subroutine gr_vh2pv_body(self, param)
       implicit none
       ! Arguments
       class(swiftest_body),       intent(inout) :: self  !! Swiftest particle object
-      class(swiftest_parameters), intent(in)    :: param !! Current run configuration parameters of on parameters 
+      class(swiftest_parameters), intent(in)    :: param !! Current run configuration parameters 
       ! Internals
       integer(I4B)                                 :: i
       real(DP), dimension(:,:), allocatable        :: pv !! Temporary holder of pseudovelocity for in-place conversion