Fixed encounter check code for particles so that it overwrites the en…

…counter list rather than appending
MintonGroup · Jul 28, 2021 · ea58038 · ea58038
1 parent 7694a02
commit ea58038
Show file tree

Hide file tree

Showing 3 changed files with 24 additions and 20 deletions.
diff --git a/src/symba/symba_encounter_check.f90 b/src/symba/symba_encounter_check.f90
@@ -15,7 +15,6 @@ module function symba_encounter_check_pl(self, system, dt, irec) result(lany_enc
       ! Result
       logical                                   :: lany_encounter !! Returns true if there is at least one close encounter      
       ! Internals
-      integer(I4B)                              :: nenc_old
       integer(I8B)                              :: k
       real(DP),     dimension(NDIM)             :: xr, vr
       logical,      dimension(:),   allocatable :: lencounter, loc_lvdotr
@@ -35,15 +34,14 @@ module function symba_encounter_check_pl(self, system, dt, irec) result(lany_enc
          lany_encounter = any(lencounter(:))
          if (lany_encounter) then 
             associate(plplenc_list => system%plplenc_list, nenc => system%plplenc_list%nenc)
-               nenc_old = nenc
-               call plplenc_list%resize(nenc_old + count(lencounter(:)))
-               plplenc_list%status(nenc_old+1:nenc) = ACTIVE
-               plplenc_list%level(nenc_old+1:nenc) = irec
-               plplenc_list%lvdotr(nenc_old+1:nenc) = pack(loc_lvdotr(:), lencounter(:))
-               plplenc_list%index1(nenc_old+1:nenc) = pack(pl%k_plpl(1,:), lencounter(:))
-               plplenc_list%index2(nenc_old+1:nenc) = pack(pl%k_plpl(2,:), lencounter(:))
-               pl%lencounter(plplenc_list%index1(nenc_old+1:nenc)) = .true.
-               pl%lencounter(plplenc_list%index2(nenc_old+1:nenc)) = .true.
+               call plplenc_list%resize(count(lencounter(:)))
+               plplenc_list%status(1:nenc) = ACTIVE
+               plplenc_list%level(1:nenc) = irec
+               plplenc_list%lvdotr(1:nenc) = pack(loc_lvdotr(:), lencounter(:))
+               plplenc_list%index1(1:nenc) = pack(pl%k_plpl(1,:), lencounter(:))
+               plplenc_list%index2(1:nenc) = pack(pl%k_plpl(2,:), lencounter(:))
+               pl%lencounter(plplenc_list%index1(1:nenc)) = .true.
+               pl%lencounter(plplenc_list%index2(1:nenc)) = .true.
             end associate
          end if
       end associate
@@ -135,7 +133,7 @@ module function symba_encounter_check_tp(self, system, dt, irec) result(lany_enc
       logical                                   :: lany_encounter !! Returns true if there is at least one close encounter      
       ! Internals
       real(DP)                                  :: r2crit, vdotr, r2, v2, tmin, r2min, term2
-      integer(I4B)                              :: i, j, nenc_old
+      integer(I4B)                              :: i, j
       real(DP),     dimension(NDIM)             :: xr, vr
       logical,      dimension(:,:), allocatable :: lencounter, loc_lvdotr
 
@@ -154,16 +152,16 @@ module function symba_encounter_check_tp(self, system, dt, irec) result(lany_enc
          lany_encounter = any(lencounter(:,:))
          if (lany_encounter) then 
             associate(pltpenc_list => system%pltpenc_list, nenc => system%pltpenc_list%nenc)
-               nenc_old = nenc
-               call pltpenc_list%resize(nenc_old + count(lencounter(:,:)))
-               pltpenc_list%status(nenc_old+1:nenc) = ACTIVE
-               pltpenc_list%level(nenc_old+1:nenc) = irec
-               pltpenc_list%lvdotr(nenc_old+1:nenc) = pack(loc_lvdotr(:,:), lencounter(:,:))
-               pltpenc_list%index1(nenc_old+1:nenc) = pack(spread([(i, i = 1, npl)], dim=2, ncopies=ntp), lencounter(:,:)) 
-               pltpenc_list%index2(nenc_old+1:nenc) = pack(spread([(j, j = 1, ntp)], dim=1, ncopies=npl), lencounter(:,:))
+               call pltpenc_list%resize(count(lencounter(:,:)))
+               pltpenc_list%status(1:nenc) = ACTIVE
+               pltpenc_list%level(1:nenc) = irec
+               pltpenc_list%lvdotr(1:nenc) = pack(loc_lvdotr(:,:), lencounter(:,:))
+               pltpenc_list%index1(1:nenc) = pack(spread([(i, i = 1, npl)], dim=2, ncopies=ntp), lencounter(:,:)) 
+               pltpenc_list%index2(1:nenc) = pack(spread([(j, j = 1, ntp)], dim=1, ncopies=npl), lencounter(:,:))
                select type(pl)
                class is (symba_pl)
-                  pl%lencounter(pltpenc_list%index1(nenc_old+1:nenc)) = .true.
+                  pl%lencounter(:) = .false.
+                  pl%lencounter(pltpenc_list%index1(1:nenc)) = .true.
                end select
             end associate
          end if

diff --git a/src/symba/symba_step.f90 b/src/symba/symba_step.f90
@@ -137,7 +137,11 @@ module recursive subroutine symba_step_recur_system(self, param, t, ireci)
 
                   call pl%drift(system, param, dtl, mask=(pl%status(:) == ACTIVE .and. pl%levelg(:) == ireci))
                   call tp%drift(system, param, dtl, mask=(tp%status(:) == ACTIVE .and. tp%levelg(:) == ireci))
-
+                  if (ireci /=0) then
+                     if (pl%levelg(2) == ireci) then
+                        write(14,*) ireci,j,pl%xh(1,2)
+                     end if
+                  end if
                   if (lencounter) call system%recursive_step(param, t+dth,irecp)
 
                   call plplenc_list%kick(system, dth, irecp, 1)

diff --git a/src/symba/symba_util.f90 b/src/symba/symba_util.f90
@@ -60,6 +60,8 @@ module subroutine symba_util_resize_pltpenc(self, nrequested)
          call self%setup(2 * nrequested)
          call self%copy(enc_temp)
          deallocate(enc_temp)
+      else
+         self%status(nrequested+1:nold) = INACTIVE
       end if
       self%nenc = nrequested
       return