Fixed memory dissues with the recursion level code

Makefile.Defines

@@ -71,7 +71,7 @@ GWARNINGS = -Wall -Warray-bounds  -Wimplicit-interface -Wextra -Warray-temporari
 #AR     = xiar
 
 FORTRAN = gfortran
-FFLAGS =  -ffree-line-length-none $(GDEBUG) #$(GMEM)
+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

@@ -31,7 +31,7 @@
 # Simulation start, stop, and output cadence times
 t_0	  = 0 # simulation start time
 deltaT	= 0.25 * swiftest.JD2S / TU2S   # simulation step size
-end_sim = deltaT #0.15
+end_sim = 0.15
 t_print = deltaT  #output interval to print results
 
 iout = int(np.ceil(t_print / deltaT))

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

@@ -1,6 +1,6 @@
 ! Swifter input file generated using init_cond.py
 T0            0 
-TSTOP         0.0006844626967830253
+TSTOP         0.15
 DT            0.0006844626967830253
 PL_IN         pl.swifter.in
 TP_IN         tp.swifter.in

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

@@ -1,6 +1,6 @@
 ! Swiftest input file generated using init_cond.py
 T0             0 
-TSTOP          0.0006844626967830253
+TSTOP          0.15
 DT             0.0006844626967830253
 CB_IN          cb.swiftest.in
 PL_IN          pl.swiftest.in

src/symba/symba_collision.f90

@@ -19,6 +19,8 @@ module subroutine symba_collision_check_plplenc(self, system, param, t, dt, irec
       real(DP),                   intent(in)    :: t      !! current time
       real(DP),                   intent(in)    :: dt     !! step size
       integer(I4B),               intent(in)    :: irec   !! Current recursion level
+
+      return
    end subroutine symba_collision_check_plplenc
 
 
@@ -43,6 +45,8 @@ module subroutine symba_collision_check_pltpenc(self, system, param, t, dt, irec
       real(DP), dimension(NDIM)                 :: xr, vr
       integer(I4B)                              :: k
 
+      if (self%nenc == 0) return
+
       select type(pl => system%pl)
       class is (symba_pl)
          select type(tp => system%tp)

src/symba/symba_step.f90

@@ -99,38 +99,44 @@ module subroutine symba_step_set_recur_levels_system(self, ireci)
       !!
       !! Resets pl, tp,and encounter structures at the start of a new step
       !! 
+      !! Adapted from David E. Kaufmann's Swifter routine: symba_step_recur.f90
+      !! Adapted from Hal Levison's Swift routine symba5_step_recur.f
       implicit none
       ! Arguments
       class(symba_nbody_system),  intent(inout) :: self  !! SyMBA nbody system object
       integer(I4B),               intent(in)    :: ireci !! Input recursion level 
       ! Internals
-      integer(I4B) :: i, irecp
+      integer(I4B) :: k, irecp
 
       associate(system => self, plplenc_list => self%plplenc_list, pltpenc_list => self%pltpenc_list)
          select type(pl => self%pl)
          class is (symba_pl)
             select type(tp => self%tp)
             class is (symba_tp)
-               associate (plind1 => plplenc_list%index1(1:plplenc_list%nenc), &
-                          plind2 => plplenc_list%index2(1:plplenc_list%nenc), &
-                          plind3 => pltpenc_list%index1(1:pltpenc_list%nenc), &
-                          tpind  => pltpenc_list%index2(1:pltpenc_list%nenc))
-
-                  irecp = ireci + 1
+               irecp = ireci + 1
 
-                  do i = 1, plplenc_list%nenc
-                     if (pl%levelg(plind1(i)) == irecp) pl%levelg(plind1(i)) = ireci
-                     if (pl%levelg(plind2(i)) == irecp) pl%levelg(plind2(i)) = ireci
+               if (plplenc_list%nenc > 0) then
+                  do k = 1, plplenc_list%nenc
+                     associate(i => plplenc_list%index1(k), j => plplenc_list%index2(k))
+                        if (pl%levelg(i) == irecp) pl%levelg(i) = ireci
+                        if (pl%levelg(j) == irecp) pl%levelg(j) = ireci
+                     end associate
                   end do
-                  do i = 1, pltpenc_list%nenc
-                     if (pl%levelg(plind3(i)) == irecp) pl%levelg(plind3(i)) = ireci
-                     if (tp%levelg(tpind(i))  == irecp) tp%levelg(tpind(i))  = ireci
+                  where(plplenc_list%level(1:plplenc_list%nenc) == irecp) plplenc_list%level(1:plplenc_list%nenc) = ireci
+               end if
+
+               if (pltpenc_list%nenc > 0) then
+                  do k = 1, pltpenc_list%nenc
+                     associate(i => pltpenc_list%index1(k), j => pltpenc_list%index2(k))
+                        if (pl%levelg(i) == irecp) pl%levelg(i) = ireci
+                        if (tp%levelg(j) == irecp) tp%levelg(j) = ireci
+                     end associate
                   end do
-               end associate
+                  where(pltpenc_list%level(1:pltpenc_list%nenc) == irecp) pltpenc_list%level(1:pltpenc_list%nenc) = ireci
+               end if
 
-               if (plplenc_list%nenc > 0) where(plplenc_list%level(1:plplenc_list%nenc) == irecp) plplenc_list%level(:) = ireci
-               if (pltpenc_list%nenc > 0) where(pltpenc_list%level(1:pltpenc_list%nenc) == irecp) pltpenc_list%level(:) = ireci
                system%irec = ireci
+
             end select
          end select
       end associate
@@ -223,6 +229,8 @@ module subroutine symba_step_reset_system(self)
       !!
       !! Resets pl, tp,and encounter structures at the start of a new step
       !! 
+      !! Adapted from David E. Kaufmann's Swifter routine: symba_step.f90
+      !! Adapted from Hal Levison's Swift routine symba5_step.f
       implicit none
       ! Arguments
       class(symba_nbody_system),  intent(inout) :: self !! SyMBA nbody system object