Merge branch 'OOPSymba' into OOPTides

Makefile.Defines

@@ -65,13 +65,13 @@ GPAR = -fopenmp -ftree-parallelize-loops=4
 GMEM = -fsanitize=undefined -fsanitize=address -fsanitize=leak 
 GWARNINGS = -Wall -Warray-bounds  -Wimplicit-interface -Wextra -Warray-temporaries
 
-#FFLAGS  = $(IDEBUG) $(HEAPARR)
+FFLAGS  = $(IDEBUG) $(HEAPARR)
 #FFLAGS  = -init=snan,arrays -no-wrap-margin -O3 $(STRICTREAL) $(SIMDVEC) $(PAR) 
-#FORTRAN = ifort
+FORTRAN = ifort
 #AR     = xiar
 
-FORTRAN = gfortran
-FFLAGS =  -ffree-line-length-none $(GDEBUG) $(GMEM)
+#FORTRAN = gfortran
+#FFLAGS =  -ffree-line-length-none $(GDEBUG) $(GMEM)
 AR     = ar
 
 # DO NOT include in CFLAGS the "-c" option to compile object only

examples/symba_swifter_comparison/1pl_1pl_encounter/init_cond.py

@@ -132,7 +132,7 @@
 plfile = FortranFile(swiftest_pl, 'w')
 plfile.write_record(npl)
 
-plfile.write_record(plid1)
+plfile.write_record(np.array([plid1, plid2]))
 plfile.write_record(np.vstack([p_pl1[0],p_pl2[0]]))
 plfile.write_record(np.vstack([p_pl1[1],p_pl2[1]]))
 plfile.write_record(np.vstack([p_pl1[2],p_pl2[2]]))

examples/symba_swifter_comparison/1pl_1pl_encounter/param.swiftest.in

@@ -5,7 +5,7 @@ DT             0.0006844626967830253
 CB_IN          cb.swiftest.in
 PL_IN          pl.swiftest.in
 TP_IN          tp.swiftest.in
-IN_TYPE        REAL8
+IN_TYPE        ASCII
 ISTEP_OUT      1
 ISTEP_DUMP     1
 BIN_OUT        bin.swiftest.dat

examples/symba_swifter_comparison/1pl_1pl_encounter/pl.swifter.in

@@ -1,5 +1,5 @@
 3 ! Planet input file generated using init_cond.py
-1 39.476926408897625196
+1 39.47692640889762629
 0.0 0.0 0.0
 0.0 0.0 0.0
 2 0.00012002693582795244940133 0.010044724833237892

examples/symba_swifter_comparison/1pl_1pl_encounter/swiftest_vs_swifter.ipynb

@@ -43,7 +43,7 @@
      "output_type": "stream",
      "text": [
       "Reading Swiftest file param.swiftest.in\n",
-      "Reading in time 6.845e-04\n",
+      "Reading in time 1.506e-01\n",
       "Creating Dataset\n"
      ]
     },

python/swiftest/tests/bin2xr/swiftest/bin2xr_swiftest.py

@@ -1,5 +1,4 @@
 import swiftest
-import swiftest.io as swio
-param_file = "param.swiftest.in"
 sim = swiftest.Simulation(param_file="param.swiftest.in")
-ds = swio.swiftest2xr(sim.param)
+sim.bin2xr()
+sim.ds

python/swiftest/tests/bin2xr/swiftest/param.swiftest.in

@@ -1,31 +1,31 @@
-! VERSION        Swiftest parameter file converted from Swift
-BIG_DISCARD      NO
-BIN_OUT          bin.dat
-CB_IN            cb.swiftest.in
-CHK_CLOSE        YES
-CHK_QMIN         0.00468
-CHK_QMIN_COORD   HELIO
-CHK_QMIN_RANGE   4.68e-03 100.0
-CHK_RMAX         100.0
-CHK_RMIN         0.00468
-DT               5.0
-DU2M             149597870700.0
-ENC_OUT          enc.dat
-ENERGY           NO
-EXTRA_FORCE      NO
-FRAGMENTATION    NO
-GR               NO
-ISTEP_DUMP       7305000
-ISTEP_OUT        7305000
-MU2KG            1.988409870698051e+30
-OUT_FORM         EL
-OUT_STAT         UNKNOWN
-PL_IN            pl.swiftest.in
-ROTATION         NO
-T0               0.0
-TIDES            NO
-TP_IN            tp.swiftest.in
-TSTOP            365250000000.0
-TU2S             86400
-YARKOVSKY        NO
-YORP             NO
+! Swiftest input file generated using init_cond.py
+T0             0 
+TSTOP          0.15
+DT             0.0006844626967830253
+CB_IN          cb.swiftest.in
+PL_IN          pl.swiftest.in
+TP_IN          tp.swiftest.in
+IN_TYPE        REAL8
+ISTEP_OUT      1
+ISTEP_DUMP     1
+BIN_OUT        bin.swiftest.dat
+OUT_TYPE       REAL8
+OUT_FORM       XV
+OUT_STAT       REPLACE
+CHK_CLOSE      yes
+CHK_RMIN       0.004650467260962157
+CHK_RMAX       1000.0
+CHK_EJECT      1000.0
+CHK_QMIN       0.004650467260962157
+CHK_QMIN_COORD HELIO
+CHK_QMIN_RANGE 0.004650467260962157 1000.0
+ENC_OUT        enc.swiftest.dat
+EXTRA_FORCE    no
+BIG_DISCARD    no
+ROTATION       no
+GR             no
+MU2KG          1.988409870698051e+30
+DU2M           149597870700.0
+TU2S           31557600.0
+RHILL_PRESENT  yes
+MTINY          1e-12

src/symba/symba_encounter_check.f90

@@ -71,10 +71,11 @@ module function symba_encounter_check_pltpenc(self, system, dt, irec) result(lan
       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
+      integer(I4B)              :: k
       real(DP), dimension(NDIM) :: xr, vr
       logical                   :: lencounter, isplpl
       real(DP)                  :: rlim2, rji2
+      logical, dimension(:), allocatable :: lencmask
 
       lany_encounter = .false.
       if (self%nenc == 0) return
@@ -90,40 +91,43 @@ module function symba_encounter_check_pltpenc(self, system, dt, irec) result(lan
       class is (symba_pl)
          select type(tp => system%tp)
          class is (symba_tp)
-            do concurrent(i = 1:self%nenc, self%status(i) == ACTIVE .and. self%level(i) == irec - 1)
-               associate(index_i => self%index1(i), index_j => self%index2(i))
+            allocate(lencmask(self%nenc))
+            lencmask(:) = (self%status(1:self%nenc) == ACTIVE) .and. (self%level(1:self%nenc) == irec - 1)
+            if (.not.any(lencmask(:))) return
+            associate(ind1 => self%index1, ind2 => self%index2) 
+               do concurrent(k = 1:self%nenc, lencmask(k))
                   if (isplpl) then
-                     xr(:) = pl%xh(:,index_j) - pl%xh(:,index_i)
-                     vr(:) = pl%vb(:,index_j) - pl%vb(:,index_i)
-                     call symba_encounter_check_one(xr(1), xr(2), xr(3), vr(1), vr(2), vr(3), pl%rhill(index_i), pl%rhill(index_j), dt, irec, lencounter, self%lvdotr(i))
+                     xr(:) = pl%xh(:,ind2(k)) - pl%xh(:,ind1(k))
+                     vr(:) = pl%vb(:,ind2(k)) - pl%vb(:,ind1(k))
+                     call symba_encounter_check_one(xr(1), xr(2), xr(3), vr(1), vr(2), vr(3), pl%rhill(ind1(k)), pl%rhill(ind2(k)), dt, irec, lencounter, self%lvdotr(k))
                   else
-                     xr(:) = tp%xh(:,index_j) - pl%xh(:,index_i)
-                     vr(:) = tp%vb(:,index_j) - pl%vb(:,index_i)
-                     call symba_encounter_check_one(xr(1), xr(2), xr(3), vr(1), vr(2), vr(3), pl%rhill(index_i), 0.0_DP, dt, irec, lencounter, self%lvdotr(i))
+                     xr(:) = tp%xh(:,ind2(k)) - pl%xh(:,ind1(k))
+                     vr(:) = tp%vb(:,ind2(k)) - pl%vb(:,ind1(k))
+                     call symba_encounter_check_one(xr(1), xr(2), xr(3), vr(1), vr(2), vr(3), pl%rhill(ind1(k)), 0.0_DP, dt, irec, lencounter, self%lvdotr(k))
                   end if
                   if (lencounter) then
                      if (isplpl) then
-                        rlim2 = (pl%radius(index_i) + pl%radius(index_j))**2
+                        rlim2 = (pl%radius(ind1(k)) + pl%radius(ind2(k)))**2
                      else
-                        rlim2 = (pl%radius(index_i))**2
+                        rlim2 = (pl%radius(ind1(k)))**2
                      end if
                      rji2 = dot_product(xr(:), xr(:))! Check to see if these are physically overlapping bodies first, which we should ignore
                      if (rji2 > rlim2) then
                         lany_encounter = .true.
-                        pl%levelg(index_i) = irec
-                        pl%levelm(index_i) = MAX(irec, pl%levelm(index_i))
+                        pl%levelg(ind1(k)) = irec
+                        pl%levelm(ind1(k)) = MAX(irec, pl%levelm(ind1(k)))
                         if (isplpl) then
-                           pl%levelg(index_j) = irec
-                           pl%levelm(index_j) = MAX(irec, pl%levelm(index_j))
+                           pl%levelg(ind2(k)) = irec
+                           pl%levelm(ind2(k)) = MAX(irec, pl%levelm(ind2(k)))
                         else
-                           tp%levelg(index_j) = irec
-                           tp%levelm(index_j) = MAX(irec, tp%levelm(index_j))
+                           tp%levelg(ind2(k)) = irec
+                           tp%levelm(ind2(k)) = MAX(irec, tp%levelm(ind2(k)))
                         end if
-                        self%level(i) = irec
+                        self%level(k) = irec
                      end if
                   end if   
-               end associate
-            end do
+               end do
+            end associate
          end select
       end select