From ea58038b6f7f08c27f9340e3b1d528ed31ac4ded Mon Sep 17 00:00:00 2001 From: David A Minton Date: Wed, 28 Jul 2021 16:45:59 -0400 Subject: [PATCH] Fixed encounter check code for particles so that it overwrites the encounter list rather than appending --- src/symba/symba_encounter_check.f90 | 36 ++++++++++++++--------------- src/symba/symba_step.f90 | 6 ++++- src/symba/symba_util.f90 | 2 ++ 3 files changed, 24 insertions(+), 20 deletions(-) diff --git a/src/symba/symba_encounter_check.f90 b/src/symba/symba_encounter_check.f90 index 796df5d4a..1a202ba60 100644 --- 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 index 4e7082c4b..09102d82c 100644 --- 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 index 71db85cc3..a19886853 100644 --- 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