diff --git a/Makefile.Defines b/Makefile.Defines index 1346f4b09..70069bb71 100644 --- a/Makefile.Defines +++ b/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 diff --git a/examples/symba_swifter_comparison/1pl_1pl_encounter/init_cond.py b/examples/symba_swifter_comparison/1pl_1pl_encounter/init_cond.py index eeb2791d0..f63022624 100755 --- a/examples/symba_swifter_comparison/1pl_1pl_encounter/init_cond.py +++ b/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)) diff --git a/examples/symba_swifter_comparison/1pl_1pl_encounter/param.swifter.in b/examples/symba_swifter_comparison/1pl_1pl_encounter/param.swifter.in index 037d91c09..853815639 100644 --- a/examples/symba_swifter_comparison/1pl_1pl_encounter/param.swifter.in +++ b/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 diff --git a/examples/symba_swifter_comparison/1pl_1pl_encounter/param.swiftest.in b/examples/symba_swifter_comparison/1pl_1pl_encounter/param.swiftest.in index 3e8f808ce..a7f91ba33 100644 --- a/examples/symba_swifter_comparison/1pl_1pl_encounter/param.swiftest.in +++ b/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 diff --git a/src/symba/symba_collision.f90 b/src/symba/symba_collision.f90 index d307fa7a2..133190dc2 100644 --- a/src/symba/symba_collision.f90 +++ b/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) diff --git a/src/symba/symba_step.f90 b/src/symba/symba_step.f90 index c34f9fb0c..e8badd577 100644 --- a/src/symba/symba_step.f90 +++ b/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