Improved encounter check memory requirements

examples/rmvs_swifter_comparison/mars_ejecta/config.swiftest.in

@@ -4,7 +4,7 @@
 !NPLMAX         -1                 ! not used
 !NTPMAX         -1                 ! not used
 T0             0.0e0
-TSTOP          365.25e5           ! simulation length in days
+TSTOP          365.25e3           ! simulation length in days
 DT             2e0            ! stepsize in days 
 ISTEP_OUT      73500                   ! output cadence
 ISTEP_DUMP     73500                   ! system dump cadence

examples/rmvs_swifter_comparison/mars_ejecta/param.swifter.in

@@ -4,7 +4,7 @@
 !NPLMAX         -1                 ! not used
 !NTPMAX         -1                 ! not used
 T0             0.0e0
-TSTOP          365.25e5           ! simulation length in days
+TSTOP          365.25e3           ! simulation length in days
 DT             2e0            ! stepsize in days 
 ISTEP_OUT      73500                   ! output cadence
 ISTEP_DUMP     73500             ! system dump cadence

examples/rmvs_swifter_comparison/mars_ejecta/profswifter.sh

@@ -0,0 +1,2 @@
+#!/bin/bash
+gprof ./swifter_rmvs | /home/daminton/git/gprof2dot/gprof2dot.py | dot -Tpng -o swifter_profile.png

examples/rmvs_swifter_comparison/mars_ejecta/tp.in

@@ -1,4 +1,4 @@
-         10 
+         200
      6000200
   0.321794727714005       -1.39371227734394      -3.555372224179648E-002
   1.263074812292130E-002  3.554726326213714E-003  5.297001107625803E-004

src/rmvs/rmvs_encounter_check.f90

@@ -21,26 +21,25 @@ module function rmvs_encounter_check_tp(self, cb, pl, dt, rts) result(lencounter
       real(DP)                                 :: r2, v2, vdotr
       real(DP), dimension(NDIM)                :: xr, vr
       real(DP), dimension(pl%nbody)            :: r2crit
+      logical                                  :: lflag
 
       associate(tp => self, ntp => self%nbody, npl => pl%nbody, rhill => pl%rhill, &
                 xht => self%xh, vht => self%vh, xbeg => self%xbeg, vbeg => self%vbeg)
-         if (.not.allocated(pl%encmask)) allocate(pl%encmask(ntp, npl))
-         pl%encmask(:,:) = .false.
          r2crit(:) = (rts * rhill(:))**2
-         do i = 1, ntp
-            if (tp%status(i) /= ACTIVE) cycle
-            tp%plencP(i) = 0
-            do j = 1, npl
+         tp%plencP(:) = 0
+         do j = 1, npl
+            do i = 1, ntp
+               if ((tp%status(i) /= ACTIVE).or.(tp%plencP(i) /=0)) cycle
                xr(:) = xht(:, i) - xbeg(:, j)
                vr(:) = vht(:, i) - vbeg(:, j)
                r2 = dot_product(xr(:), xr(:))
                v2 = dot_product(vr(:), vr(:))
                vdotr = dot_product(vr(:), xr(:))
-               pl%encmask(i,j) = rmvs_chk_ind(r2, v2, vdotr, dt, r2crit(j))
-               if (pl%encmask(i,j)) tp%plencP(i) = j
+               lflag = rmvs_chk_ind(r2, v2, vdotr, dt, r2crit(j))
+               if (lflag) tp%plencP(i) = j
             end do
+            pl%nenc(j) = count(tp%plencP(:) == j)
          end do
-         pl%nenc(:) = count(pl%encmask(:,:), dim=1)
          lencounter = any(pl%nenc(:) > 0)
       end associate
       return

src/rmvs/rmvs_step.f90

@@ -224,13 +224,17 @@ module subroutine rmvs_step_make_planetocentric(self, cb, tp, config)
       class(swiftest_configuration),  intent(in)     :: config !! Input collection of configuration parameters 
       ! Internals
       integer(I4B)                                   :: i, j
+      logical, dimension(:), allocatable             :: encmask
 
       associate(pl => self, npl => self%nbody, nenc => self%nenc, tpenc => self%tpenc, cbenc => self%cbenc, &
          plenc => self%plenc, GMpl => self%Gmass)
 
          do i = 1, npl
             if (nenc(i) == 0) cycle 
             ! There are inner encounters with this planet
+            if (allocated(encmask)) deallocate(encmask)
+            allocate(encmask(nenc(i)))
+            encmask(:) = tp%plencP(:) == i
 
             ! Save the index value of the planet corresponding to this encounter 
             tpenc(i)%ipleP = i
@@ -252,13 +256,13 @@ module subroutine rmvs_step_make_planetocentric(self, cb, tp, config)
             tpenc(i)%status(:) = ACTIVE
             tpenc(i)%lperi(:) = .false.
             tpenc(i)%plperP(:) = 0
-            tpenc(i)%name(:) = pack(tp%name(:), pl%encmask(:,i)) 
+            tpenc(i)%name(:) = pack(tp%name(:), encmask(:)) 
             !call tp%spill(tpenc(i), pl%encmask(:,i))
             ! Grab all the encountering test particles and convert them to a planetocentric frame
             do j = 1, NDIM 
-               tpenc(i)%xheliocen(j, :) = pack(tp%xh(j,:), pl%encmask(:,i)) 
+               tpenc(i)%xheliocen(j, :) = pack(tp%xh(j,:), encmask(:)) 
                tpenc(i)%xh(j, :) = tpenc(i)%xheliocen(j, :) - pl%xin(j, i, 0)
-               tpenc(i)%vh(j, :) = pack(tp%vh(j,:), pl%encmask(:,i)) - pl%vin(j, i, 0)
+               tpenc(i)%vh(j, :) = pack(tp%vh(j,:), encmask(:)) - pl%vin(j, i, 0)
             end do
 
             ! Make sure that the test particles get the planetocentric value of mu 
@@ -289,15 +293,19 @@ module subroutine rmvs_step_end_planetocentric(self, cb, tp)
       ! Internals
       integer(I4B) :: i, j
       integer(I4B), dimension(:), allocatable :: tpind
+      logical, dimension(:), allocatable :: encmask
 
       associate(pl => self, nenc => self%nenc, npl => self%nbody, ntp => tp%nbody, name => tp%name, &
-         tpenc => self%tpenc, plenc => self%plenc, cbenc => self%cbenc, encmask => self%encmask)
+         tpenc => self%tpenc, plenc => self%plenc, cbenc => self%cbenc)
 
          do i = 1, npl
             if (nenc(i) == 0) cycle
             allocate(tpind(nenc(i)))
             ! Index array of encountering test particles
-            tpind(:) = pack([(j,j=1,ntp)], encmask(1:ntp, i))
+            if (allocated(encmask)) deallocate(encmask)
+            allocate(encmask(nenc(i)))
+            encmask(:) = tp%plencP(:) == i
+            tpind(:) = pack([(j,j=1,ntp)], encmask(:))
 
             ! Copy the results of the integration back over and shift back to heliocentric reference
             tp%status(tpind(1:nenc(i))) = tpenc(i)%status(1:nenc(i)) 
@@ -319,7 +327,6 @@ module subroutine rmvs_step_end_planetocentric(self, cb, tp)
             deallocate(tpind)
          end do
       end associate
-      deallocate(self%encmask)
 
       return
    end subroutine rmvs_step_end_planetocentric