From 55f745dff0ac43d1a398d439997750aabcd21510 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Wed, 18 Aug 2021 11:05:07 -0400 Subject: [PATCH 1/6] Changed whole array operations from (:) to (1:self%nbody) in order to only apply operations on valid array elements --- src/util/util_set.f90 | 21 +++++++++++---------- 1 file changed, 11 insertions(+), 10 deletions(-) diff --git a/src/util/util_set.f90 b/src/util/util_set.f90 index 86e021ab6..4e686eec8 100644 --- a/src/util/util_set.f90 +++ b/src/util/util_set.f90 @@ -77,7 +77,7 @@ module subroutine util_set_mu_pl(self, cb) class(swiftest_pl), intent(inout) :: self !! Swiftest massive body object class(swiftest_cb), intent(inout) :: cb !! Swiftest central body object - if (self%nbody > 0) self%mu(:) = cb%Gmass + self%Gmass(:) + if (self%nbody > 0) self%mu(1:self%nbody) = cb%Gmass + self%Gmass(1:self%nbody) return end subroutine util_set_mu_pl @@ -92,7 +92,8 @@ module subroutine util_set_mu_tp(self, cb) class(swiftest_tp), intent(inout) :: self !! Swiftest test particle object class(swiftest_cb), intent(inout) :: cb !! Swiftest central body object - if (self%nbody > 0) self%mu(:) = cb%Gmass + if (self%nbody == 0) return + self%mu(1:self%nbody) = cb%Gmass return end subroutine util_set_mu_tp @@ -107,10 +108,10 @@ module subroutine util_set_rhill(self,cb) class(swiftest_pl), intent(inout) :: self !! Swiftest massive body object class(swiftest_cb), intent(inout) :: cb !! Swiftest central body object - if (self%nbody > 0) then - call self%xv2el(cb) - self%rhill(:) = self%a(:) * (self%Gmass(:) / cb%Gmass / 3)**THIRD - end if + if (self%nbody == 0) return + + call self%xv2el(cb) + self%rhill(1:self%nbody) = self%a(1:self%nbody) * (self%Gmass(1:self%nbody) / cb%Gmass / 3)**THIRD return end subroutine util_set_rhill @@ -127,10 +128,10 @@ module subroutine util_set_rhill_approximate(self,cb) ! Internals real(DP), dimension(:), allocatable :: rh - if (self%nbody > 0) then - rh(:) = .mag. self%xh(:,:) - self%rhill(:) = rh(:) * (self%Gmass(:) / cb%Gmass / 3)**THIRD - end if + if (self%nbody == 0) return + + rh(1:self%nbody) = .mag. self%xh(:,1:self%nbody) + self%rhill(1:self%nbody) = rh(1:self%nbody) * (self%Gmass(1:self%nbody) / cb%Gmass / 3)**THIRD return end subroutine util_set_rhill_approximate From c6b7a00089a47c2440304eeed032a4f5e872606f Mon Sep 17 00:00:00 2001 From: David A Minton Date: Wed, 18 Aug 2021 14:50:33 -0400 Subject: [PATCH 2/6] Added explicit array range for a bunch of whole array operations. Also fixed a lingering bug in RMVS that was causing the planet accelerations to be miscalculated after an encounter was done if oblateness was turned on --- src/discard/discard.f90 | 6 +- src/helio/helio_kick.f90 | 8 +- src/kick/kick.f90 | 2 + src/rmvs/rmvs_encounter_check.f90 | 9 +- src/rmvs/rmvs_kick.f90 | 6 +- src/rmvs/rmvs_step.f90 | 244 ++++++++++++++-------------- src/symba/symba_collision.f90 | 40 ++--- src/symba/symba_drift.f90 | 30 ++-- src/symba/symba_encounter_check.f90 | 8 +- src/symba/symba_io.f90 | 52 +++--- src/symba/symba_kick.f90 | 22 +-- src/whm/whm_kick.f90 | 28 ++-- 12 files changed, 234 insertions(+), 221 deletions(-) diff --git a/src/discard/discard.f90 b/src/discard/discard.f90 index 9c8044c61..2998e8804 100644 --- a/src/discard/discard.f90 +++ b/src/discard/discard.f90 @@ -53,7 +53,7 @@ module subroutine discard_pl(self, system, param) class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameter if (self%nbody == 0) return - self%ldiscard(:) = .false. + self%ldiscard(1:self%nbody) = .false. return end subroutine discard_pl @@ -86,9 +86,9 @@ module subroutine discard_tp(self, system, param) if ((param%rmin >= 0.0_DP) .or. (param%rmax >= 0.0_DP) .or. (param%rmaxu >= 0.0_DP)) call discard_cb_tp(tp, system, param) if (param%qmin >= 0.0_DP) call discard_peri_tp(tp, system, param) if (param%lclose) call discard_pl_tp(tp, system, param) - if (any(tp%ldiscard)) then + if (any(tp%ldiscard(1:ntp))) then allocate(ldiscard, source=tp%ldiscard) - call tp%spill(system%tp_discards, ldiscard, ldestructive=.true.) + call tp%spill(system%tp_discards, ldiscard(1:ntp), ldestructive=.true.) end if end associate diff --git a/src/helio/helio_kick.f90 b/src/helio/helio_kick.f90 index eebd17f53..8e7571a0b 100644 --- a/src/helio/helio_kick.f90 +++ b/src/helio/helio_kick.f90 @@ -65,9 +65,9 @@ module subroutine helio_kick_getacch_tp(self, system, param, t, lbeg) associate(tp => self, cb => system%cb, pl => system%pl, npl => system%pl%nbody) system%lbeg = lbeg if (system%lbeg) then - call tp%accel_int(pl%Gmass(:), pl%xbeg(:,:), npl) + call tp%accel_int(pl%Gmass(1:npl), pl%xbeg(:,1:npl), npl) else - call tp%accel_int(pl%Gmass(:), pl%xend(:,:), npl) + call tp%accel_int(pl%Gmass(1:npl), pl%xend(:,1:npl), npl) end if if (param%loblatecb) call tp%accel_obl(system) if (param%lextra_force) call tp%accel_user(system, param, t, lbeg) @@ -99,7 +99,7 @@ module subroutine helio_kick_vb_pl(self, system, param, t, dt, lbeg) if (self%nbody == 0) return associate(pl => self, npl => self%nbody) - pl%ah(:,:) = 0.0_DP + pl%ah(:, 1:npl) = 0.0_DP call pl%accel(system, param, t, lbeg) if (lbeg) then call pl%set_beg_end(xbeg = pl%xh) @@ -136,7 +136,7 @@ module subroutine helio_kick_vb_tp(self, system, param, t, dt, lbeg) if (self%nbody == 0) return associate(tp => self, ntp => self%nbody) - tp%ah(:,:) = 0.0_DP + tp%ah(:, 1:ntp) = 0.0_DP call tp%accel(system, param, t, lbeg) do concurrent(i = 1:ntp, tp%lmask(i)) tp%vb(:, i) = tp%vb(:, i) + tp%ah(:, i) * dt diff --git a/src/kick/kick.f90 b/src/kick/kick.f90 index c41bf7649..af704c13c 100644 --- a/src/kick/kick.f90 +++ b/src/kick/kick.f90 @@ -64,6 +64,8 @@ module pure subroutine kick_getacch_int_tp(self, GMpl, xhp, npl) real(DP) :: rji2, irij3, fac, r2 real(DP), dimension(NDIM) :: dx + if ((self%nbody == 0) .or. (npl == 0)) return + associate(tp => self, ntp => self%nbody) do concurrent(i = 1:ntp, tp%lmask(i)) do j = 1, npl diff --git a/src/rmvs/rmvs_encounter_check.f90 b/src/rmvs/rmvs_encounter_check.f90 index e4c441472..4cb119e48 100644 --- a/src/rmvs/rmvs_encounter_check.f90 +++ b/src/rmvs/rmvs_encounter_check.f90 @@ -23,13 +23,14 @@ module function rmvs_encounter_check_tp(self, system, dt) result(lencounter) real(DP), dimension(system%pl%nbody) :: r2crit logical :: lflag + lencounter = .false. if (self%nbody == 0) return select type(pl => system%pl) class is (rmvs_pl) associate(tp => self, ntp => self%nbody, npl => pl%nbody, rts => system%rts) - r2crit(:) = (rts * pl%rhill(:))**2 - tp%plencP(:) = 0 + r2crit(1:npl) = (rts * pl%rhill(1:npl))**2 + tp%plencP(1:ntp) = 0 do j = 1, npl do i = 1, ntp if ((.not.tp%lmask(i)).or.(tp%plencP(i) /= 0)) cycle @@ -41,9 +42,9 @@ module function rmvs_encounter_check_tp(self, system, dt) result(lencounter) 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) + pl%nenc(j) = count(tp%plencP(1:ntp) == j) end do - lencounter = any(pl%nenc(:) > 0) + lencounter = any(pl%nenc(1:npl) > 0) end associate end select return diff --git a/src/rmvs/rmvs_kick.f90 b/src/rmvs/rmvs_kick.f90 index 018ada8f3..5d9514763 100644 --- a/src/rmvs/rmvs_kick.f90 +++ b/src/rmvs/rmvs_kick.f90 @@ -61,7 +61,9 @@ module subroutine rmvs_kick_getacch_tp(self, system, param, t, lbeg) end if ! Swap the planetocentric and heliocentric position vectors and central body masses - tp%xh(:,:) = tp%xheliocentric(:,:) + do concurrent(i = 1:ntp, tp%lmask(i)) + tp%xh(:, i) = tp%xheliocentric(:, i) + end do GMcb_original = cb%Gmass cb%Gmass = tp%cb_heliocentric%Gmass @@ -71,7 +73,7 @@ module subroutine rmvs_kick_getacch_tp(self, system, param, t, lbeg) if (param%lgr) call tp%accel_gr(param) ! Put everything back the way we found it - tp%xh(:,:) = xh_original(:,:) + call move_alloc(xh_original, tp%xh) cb%Gmass = GMcb_original end associate diff --git a/src/rmvs/rmvs_step.f90 b/src/rmvs/rmvs_step.f90 index c801ee16f..b8cfcd688 100644 --- a/src/rmvs/rmvs_step.f90 +++ b/src/rmvs/rmvs_step.f90 @@ -21,6 +21,8 @@ module subroutine rmvs_step_system(self, param, t, dt) real(DP), dimension(:,:), allocatable :: xbeg, xend, vbeg integer(I4B) :: i + if (self%tp%nbody == 0) call whm_step_system(self, param, t, dt) + select type(cb => self%cb) class is (rmvs_cb) select type(pl => self%pl) @@ -36,11 +38,11 @@ module subroutine rmvs_step_system(self, param, t, dt) lencounter = tp%encounter_check(system, dt) if (lencounter) then lfirstpl = pl%lfirst - pl%outer(0)%x(:,:) = xbeg(:,:) - pl%outer(0)%v(:,:) = vbeg(:,:) + pl%outer(0)%x(:, 1:npl) = xbeg(:, 1:npl) + pl%outer(0)%v(:, 1:npl) = vbeg(:, 1:npl) call pl%step(system, param, t, dt) - pl%outer(NTENC)%x(:,:) = pl%xh(:,:) - pl%outer(NTENC)%v(:,:) = pl%vh(:,:) + pl%outer(NTENC)%x(:, 1:npl) = pl%xh(:, 1:npl) + pl%outer(NTENC)%v(:, 1:npl) = pl%vh(:, 1:npl) call rmvs_interp_out(cb, pl, dt) call rmvs_step_out(cb, pl, tp, system, param, t, dt) tp%lmask(1:ntp) = .not. tp%lmask(1:ntp) @@ -89,10 +91,10 @@ subroutine rmvs_interp_out(cb, pl, dt) allocate(GMcb(npl)) allocate(dto(npl)) allocate(iflag(npl)) - dto(:) = dt / dntenc - GMcb(:) = cb%Gmass - xtmp(:,:) = pl%outer(0)%x(:, :) - vtmp(:,:) = pl%outer(0)%v(:, :) + dto(1:npl) = dt / dntenc + GMcb(1:npl) = cb%Gmass + xtmp(:,1:npl) = pl%outer(0)%x(:, 1:npl) + vtmp(:,1:npl) = pl%outer(0)%v(:, 1:npl) do outer_index = 1, NTENC - 1 call drift_one(GMcb(1:npl), xtmp(1,1:npl), xtmp(2,1:npl), xtmp(3,1:npl), & vtmp(1,1:npl), vtmp(2,1:npl), vtmp(3,1:npl), & @@ -110,11 +112,11 @@ subroutine rmvs_interp_out(cb, pl, dt) end do end if frac = 1.0_DP - outer_index / dntenc - pl%outer(outer_index)%x(:, :) = frac * xtmp(:,:) - pl%outer(outer_index)%v(:, :) = frac * vtmp(:,:) + pl%outer(outer_index)%x(:, 1:npl) = frac * xtmp(:, 1:npl) + pl%outer(outer_index)%v(:, 1:npl) = frac * vtmp(:, 1:npl) end do - xtmp(:,:) = pl%outer(NTENC)%x(:, :) - vtmp(:,:) = pl%outer(NTENC)%v(:, :) + xtmp(:, 1:npl) = pl%outer(NTENC)%x(:, 1:npl) + vtmp(:, 1:npl) = pl%outer(NTENC)%v(:, 1:npl) do outer_index = NTENC - 1, 1, -1 call drift_one(GMcb(1:npl), xtmp(1,1:npl), xtmp(2,1:npl), xtmp(3,1:npl), & vtmp(1,1:npl), vtmp(2,1:npl), vtmp(3,1:npl), & @@ -132,8 +134,8 @@ subroutine rmvs_interp_out(cb, pl, dt) end do end if frac = outer_index / dntenc - pl%outer(outer_index)%x(:, :) = pl%outer(outer_index)%x(:, :) + frac * xtmp(:,:) - pl%outer(outer_index)%v(:, :) = pl%outer(outer_index)%v(:, :) + frac * vtmp(:,:) + pl%outer(outer_index)%x(:, 1:npl) = pl%outer(outer_index)%x(:, 1:npl) + frac * xtmp(:, 1:npl) + pl%outer(outer_index)%v(:, 1:npl) = pl%outer(outer_index)%v(:, 1:npl) + frac * vtmp(:, 1:npl) end do end associate @@ -165,16 +167,16 @@ subroutine rmvs_step_out(cb, pl, tp, system, param, t, dt) associate(npl => pl%nbody, ntp => tp%nbody) dto = dt / NTENC - where(tp%plencP(:) == 0) - tp%lmask(:) = .false. + where(tp%plencP(1:ntp) == 0) + tp%lmask(1:ntp) = .false. elsewhere - tp%lperi(:) = .false. + tp%lperi(1:ntp) = .false. end where do outer_index = 1, NTENC outer_time = t + (outer_index - 1) * dto - call pl%set_beg_end(xbeg = pl%outer(outer_index - 1)%x(:, :), & - vbeg = pl%outer(outer_index - 1)%v(:, :), & - xend = pl%outer(outer_index )%x(:, :)) + call pl%set_beg_end(xbeg = pl%outer(outer_index - 1)%x(:, 1:npl), & + vbeg = pl%outer(outer_index - 1)%v(:, 1:npl), & + xend = pl%outer(outer_index )%x(:, 1:npl)) system%rts = RHPSCALE lencounter = tp%encounter_check(system, dto) if (lencounter) then @@ -192,8 +194,8 @@ subroutine rmvs_step_out(cb, pl, tp, system, param, t, dt) do j = 1, npl if (pl%nenc(j) == 0) cycle tp%lfirst = .true. - where((tp%plencP(:) == j) .and. (.not.tp%lmask(:))) - tp%lmask(:) = .true. + where((tp%plencP(1:ntp) == j) .and. (.not.tp%lmask(1:ntp))) + tp%lmask(1:ntp) = .true. end where end do end do @@ -222,7 +224,7 @@ subroutine rmvs_interp_in(cb, pl, system, param, dt, outer_index) ! Internals integer(I4B) :: i, inner_index real(DP) :: frac, dntphenc - real(DP), dimension(:,:), allocatable :: xtmp, vtmp, xh_original + real(DP), dimension(:,:), allocatable :: xtmp, vtmp, xh_original, ah_original real(DP), dimension(:), allocatable :: GMcb, dti integer(I4B), dimension(:), allocatable :: iflag @@ -230,32 +232,33 @@ subroutine rmvs_interp_in(cb, pl, system, param, dt, outer_index) dntphenc = real(NTPHENC, kind=DP) ! Set the endpoints of the inner region from the outer region values in the current outer step index - pl%inner(0)%x(:,:) = pl%outer(outer_index - 1)%x(:, :) - pl%inner(0)%v(:,:) = pl%outer(outer_index - 1)%v(:, :) - pl%inner(NTPHENC)%x(:,:) = pl%outer(outer_index)%x(:, :) - pl%inner(NTPHENC)%v(:,:) = pl%outer(outer_index)%v(:, :) + pl%inner(0)%x(:, 1:npl) = pl%outer(outer_index - 1)%x(:, 1:npl) + pl%inner(0)%v(:, 1:npl) = pl%outer(outer_index - 1)%v(:, 1:npl) + pl%inner(NTPHENC)%x(:, 1:npl) = pl%outer(outer_index)%x(:, 1:npl) + pl%inner(NTPHENC)%v(:, 1:npl) = pl%outer(outer_index)%v(:, 1:npl) allocate(xtmp,mold=pl%xh) allocate(vtmp,mold=pl%vh) allocate(GMcb(npl)) allocate(dti(npl)) allocate(iflag(npl)) - dti(:) = dt / dntphenc - GMcb(:) = cb%Gmass - xtmp(:, :) = pl%inner(0)%x(:, :) - vtmp(:, :) = pl%inner(0)%v(:, :) + dti(1:npl) = dt / dntphenc + GMcb(1:npl) = cb%Gmass + xtmp(:, 1:npl) = pl%inner(0)%x(:, 1:npl) + vtmp(:, 1:npl) = pl%inner(0)%v(:, 1:npl) if ((param%loblatecb) .or. (param%ltides)) then allocate(xh_original, source=pl%xh) - pl%xh(:, :) = xtmp(:, :) ! Temporarily replace heliocentric position with inner substep values to calculate the oblateness terms + allocate(ah_original, source=pl%ah) + pl%xh(:, 1:npl) = xtmp(:, 1:npl) ! Temporarily replace heliocentric position with inner substep values to calculate the oblateness terms end if if (param%loblatecb) then call pl%accel_obl(system) - pl%inner(0)%aobl(:, :) = pl%aobl(:, :) ! Save the oblateness acceleration on the planet for this substep + pl%inner(0)%aobl(:, 1:npl) = pl%aobl(:, 1:npl) ! Save the oblateness acceleration on the planet for this substep end if if (param%ltides) then call pl%accel_tides(system) - pl%inner(0)%atide(:, :) = pl%atide(:, :) ! Save the oblateness acceleration on the planet for this substep + pl%inner(0)%atide(:, 1:npl) = pl%atide(:, 1:npl) ! Save the oblateness acceleration on the planet for this substep end if do inner_index = 1, NTPHENC - 1 @@ -275,12 +278,12 @@ subroutine rmvs_interp_in(cb, pl, system, param, dt, outer_index) end do end if frac = 1.0_DP - inner_index / dntphenc - pl%inner(inner_index)%x(:, :) = frac * xtmp(:,:) - pl%inner(inner_index)%v(:, :) = frac * vtmp(:,:) + pl%inner(inner_index)%x(:, 1:npl) = frac * xtmp(:, 1:npl) + pl%inner(inner_index)%v(:, 1:npl) = frac * vtmp(:, 1:npl) end do - xtmp(:,:) = pl%inner(NTPHENC)%x(:, :) - vtmp(:,:) = pl%inner(NTPHENC)%v(:, :) + xtmp(:, 1:npl) = pl%inner(NTPHENC)%x(:, 1:npl) + vtmp(:, 1:npl) = pl%inner(NTPHENC)%v(:, 1:npl) do inner_index = NTPHENC - 1, 1, -1 call drift_one(GMcb(1:npl), xtmp(1,1:npl), xtmp(2,1:npl), xtmp(3,1:npl), & @@ -299,31 +302,32 @@ subroutine rmvs_interp_in(cb, pl, system, param, dt, outer_index) end do end if frac = inner_index / dntphenc - pl%inner(inner_index)%x(:, :) = pl%inner(inner_index)%x(:, :) + frac * xtmp(:, :) - pl%inner(inner_index)%v(:, :) = pl%inner(inner_index)%v(:, :) + frac * vtmp(:, :) + pl%inner(inner_index)%x(:, 1:npl) = pl%inner(inner_index)%x(:, 1:npl) + frac * xtmp(:, 1:npl) + pl%inner(inner_index)%v(:, 1:npl) = pl%inner(inner_index)%v(:, 1:npl) + frac * vtmp(:, 1:npl) if (param%loblatecb) then - pl%xh(:,:) = pl%inner(inner_index)%x(:, :) + pl%xh(:,1:npl) = pl%inner(inner_index)%x(:, 1:npl) call pl%accel_obl(system) - pl%inner(inner_index)%aobl(:, :) = pl%aobl(:, :) + pl%inner(inner_index)%aobl(:, 1:npl) = pl%aobl(:, 1:npl) end if if (param%ltides) then call pl%accel_tides(system) - pl%inner(inner_index)%atide(:, :) = pl%atide(:, :) + pl%inner(inner_index)%atide(:, 1:npl) = pl%atide(:, 1:npl) end if end do if (param%loblatecb) then ! Calculate the final value of oblateness accelerations at the final inner substep - pl%xh(:,:) = pl%inner(NTPHENC)%x(:, :) + pl%xh(:, 1:npl) = pl%inner(NTPHENC)%x(:, 1:npl) call pl%accel_obl(system) - pl%inner(NTPHENC)%aobl(:, :) = pl%aobl(:, :) + pl%inner(NTPHENC)%aobl(:, 1:npl) = pl%aobl(:, 1:npl) end if if (param%ltides) then call pl%accel_tides(system) - pl%inner(NTPHENC)%atide(:, :) = pl%atide(:, :) + pl%inner(NTPHENC)%atide(:, 1:npl) = pl%atide(:, 1:npl) end if - ! Put the planet positions back into place + ! Put the planet positions and accelerations back into place if (allocated(xh_original)) call move_alloc(xh_original, pl%xh) + if (allocated(ah_original)) call move_alloc(ah_original, pl%ah) end associate return @@ -371,7 +375,7 @@ subroutine rmvs_step_in(cb, pl, tp, param, outer_time, dto) ! now step the encountering test particles fully through the inner encounter lfirsttp = .true. do inner_index = 1, NTPHENC ! Integrate over the encounter region, using the "substitute" planetocentric systems at each level - plenci%xh(:,:) = plenci%inner(inner_index - 1)%x(:,:) + plenci%xh(:, 1:npl) = plenci%inner(inner_index - 1)%x(:, 1:npl) call plenci%set_beg_end(xbeg = plenci%inner(inner_index - 1)%x, & xend = plenci%inner(inner_index)%x) @@ -428,7 +432,7 @@ subroutine rmvs_make_planetocentric(param, cb, pl, tp) ! There are inner encounters with this planet if (allocated(encmask)) deallocate(encmask) allocate(encmask(ntp)) - encmask(:) = tp%plencP(:) == i + encmask(1:ntp) = tp%plencP(1:ntp) == i allocate(rmvs_tp :: pl%planetocentric(i)%tp) ! Create encountering test particle structure select type(cbenci => pl%planetocentric(i)%cb) @@ -438,51 +442,53 @@ subroutine rmvs_make_planetocentric(param, cb, pl, tp) select type(tpenci => pl%planetocentric(i)%tp) class is (rmvs_tp) tpenci%lplanetocentric = .true. - call tpenci%setup(pl%nenc(i), param) - tpenci%cb_heliocentric = cb - tpenci%ipleP = i - tpenci%lmask(:) = .true. - tpenci%status(:) = ACTIVE - ! Grab all the encountering test particles and convert them to a planetocentric frame - tpenci%id(:) = pack(tp%id(:), encmask(:)) - do j = 1, NDIM - tpenci%xheliocentric(j, :) = pack(tp%xh(j,:), encmask(:)) - tpenci%xh(j, :) = tpenci%xheliocentric(j, :) - pl%inner(0)%x(j, i) - tpenci%vh(j, :) = pack(tp%vh(j,:), encmask(:)) - pl%inner(0)%v(j, i) - end do - tpenci%lperi(:) = pack(tp%lperi(:), encmask(:)) - tpenci%plperP(:) = pack(tp%plperP(:), encmask(:)) - ! Make sure that the test particles get the planetocentric value of mu - allocate(cbenci%inner(0:NTPHENC)) - do inner_index = 0, NTPHENC - allocate(plenci%inner(inner_index)%x, mold=pl%inner(inner_index)%x) - allocate(plenci%inner(inner_index)%v, mold=pl%inner(inner_index)%x) - allocate(cbenci%inner(inner_index)%x(NDIM,1)) - allocate(cbenci%inner(inner_index)%v(NDIM,1)) - cbenci%inner(inner_index)%x(:,1) = pl%inner(inner_index)%x(:, i) - cbenci%inner(inner_index)%v(:,1) = pl%inner(inner_index)%v(:, i) - plenci%inner(inner_index)%x(:,1) = -cbenci%inner(inner_index)%x(:,1) - plenci%inner(inner_index)%v(:,1) = -cbenci%inner(inner_index)%v(:,1) - - if (param%loblatecb) then - allocate(plenci%inner(inner_index)%aobl, mold=pl%inner(inner_index)%aobl) - allocate(cbenci%inner(inner_index)%aobl(NDIM,1)) - cbenci%inner(inner_index)%aobl(:,1) = pl%inner(inner_index)%aobl(:, i) - end if - - if (param%ltides) then - allocate(plenci%inner(inner_index)%atide, mold=pl%inner(inner_index)%atide) - allocate(cbenci%inner(inner_index)%atide(NDIM,1)) - cbenci%inner(inner_index)%atide(:,1) = pl%inner(inner_index)%atide(:, i) - end if - - do j = 2, npl - ipc2hc = plenci%plind(j) - plenci%inner(inner_index)%x(:,j) = pl%inner(inner_index)%x(:, ipc2hc) - cbenci%inner(inner_index)%x(:,1) - plenci%inner(inner_index)%v(:,j) = pl%inner(inner_index)%v(:, ipc2hc) - cbenci%inner(inner_index)%v(:,1) + associate(nenci => pl%nenc(i)) + call tpenci%setup(nenci, param) + tpenci%cb_heliocentric = cb + tpenci%ipleP = i + tpenci%lmask(1:nenci) = .true. + tpenci%status(1:nenci) = ACTIVE + ! Grab all the encountering test particles and convert them to a planetocentric frame + tpenci%id(1:nenci) = pack(tp%id(1:ntp), encmask(1:ntp)) + do j = 1, NDIM + tpenci%xheliocentric(j, 1:nenci) = pack(tp%xh(j,1:ntp), encmask(:)) + tpenci%xh(j, 1:nenci) = tpenci%xheliocentric(j, 1:nenci) - pl%inner(0)%x(j, i) + tpenci%vh(j, 1:nenci) = pack(tp%vh(j, 1:ntp), encmask(1:ntp)) - pl%inner(0)%v(j, i) + end do + tpenci%lperi(1:nenci) = pack(tp%lperi(1:ntp), encmask(1:ntp)) + tpenci%plperP(1:nenci) = pack(tp%plperP(1:ntp), encmask(1:ntp)) + ! Make sure that the test particles get the planetocentric value of mu + allocate(cbenci%inner(0:NTPHENC)) + do inner_index = 0, NTPHENC + allocate(plenci%inner(inner_index)%x, mold=pl%inner(inner_index)%x) + allocate(plenci%inner(inner_index)%v, mold=pl%inner(inner_index)%x) + allocate(cbenci%inner(inner_index)%x(NDIM,1)) + allocate(cbenci%inner(inner_index)%v(NDIM,1)) + cbenci%inner(inner_index)%x(:,1) = pl%inner(inner_index)%x(:, i) + cbenci%inner(inner_index)%v(:,1) = pl%inner(inner_index)%v(:, i) + plenci%inner(inner_index)%x(:,1) = -cbenci%inner(inner_index)%x(:,1) + plenci%inner(inner_index)%v(:,1) = -cbenci%inner(inner_index)%v(:,1) + + if (param%loblatecb) then + allocate(plenci%inner(inner_index)%aobl, mold=pl%inner(inner_index)%aobl) + allocate(cbenci%inner(inner_index)%aobl(NDIM,1)) + cbenci%inner(inner_index)%aobl(:,1) = pl%inner(inner_index)%aobl(:, i) + end if + + if (param%ltides) then + allocate(plenci%inner(inner_index)%atide, mold=pl%inner(inner_index)%atide) + allocate(cbenci%inner(inner_index)%atide(NDIM,1)) + cbenci%inner(inner_index)%atide(:,1) = pl%inner(inner_index)%atide(:, i) + end if + + do j = 2, npl + ipc2hc = plenci%plind(j) + plenci%inner(inner_index)%x(:,j) = pl%inner(inner_index)%x(:, ipc2hc) - cbenci%inner(inner_index)%x(:,1) + plenci%inner(inner_index)%v(:,j) = pl%inner(inner_index)%v(:, ipc2hc) - cbenci%inner(inner_index)%v(:,1) + end do end do - end do - call tpenci%set_mu(cbenci) + call tpenci%set_mu(cbenci) + end associate end select end select end select @@ -599,31 +605,33 @@ subroutine rmvs_end_planetocentric(pl, tp) class is (rmvs_pl) select type(tpenci => pl%planetocentric(i)%tp) class is (rmvs_tp) - if (allocated(tpind)) deallocate(tpind) - allocate(tpind(pl%nenc(i))) - ! Index array of encountering test particles - if (allocated(encmask)) deallocate(encmask) - allocate(encmask(ntp)) - 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:pl%nenc(i))) = tpenci%status(1:pl%nenc(i)) - tp%lmask(tpind(1:pl%nenc(i))) = tpenci%lmask(1:pl%nenc(i)) - do j = 1, NDIM - tp%xh(j, tpind(1:pl%nenc(i))) = tpenci%xh(j,1:pl%nenc(i)) + pl%inner(NTPHENC)%x(j, i) - tp%vh(j, tpind(1:pl%nenc(i))) = tpenci%vh(j,1:pl%nenc(i)) + pl%inner(NTPHENC)%v(j, i) - end do - tp%lperi(tpind(1:pl%nenc(i))) = tpenci%lperi(1:pl%nenc(i)) - tp%plperP(tpind(1:pl%nenc(i))) = tpenci%plperP(1:pl%nenc(i)) - deallocate(pl%planetocentric(i)%tp) - deallocate(cbenci%inner) - do inner_index = 0, NTPHENC - deallocate(plenci%inner(inner_index)%x) - deallocate(plenci%inner(inner_index)%v) - if (allocated(plenci%inner(inner_index)%aobl)) deallocate(plenci%inner(inner_index)%aobl) - if (allocated(plenci%inner(inner_index)%atide)) deallocate(plenci%inner(inner_index)%atide) - end do + associate(nenci => pl%nenc(i)) + if (allocated(tpind)) deallocate(tpind) + allocate(tpind(nenci)) + ! Index array of encountering test particles + if (allocated(encmask)) deallocate(encmask) + allocate(encmask(ntp)) + encmask(1:ntp) = tp%plencP(1:ntp) == i + tpind(1:nenci) = pack([(j,j=1,ntp)], encmask(1:ntp)) + + ! Copy the results of the integration back over and shift back to heliocentric reference + tp%status(tpind(1:nenci)) = tpenci%status(1:nenci) + tp%lmask(tpind(1:nenci)) = tpenci%lmask(1:nenci) + do j = 1, NDIM + tp%xh(j, tpind(1:nenci)) = tpenci%xh(j,1:nenci) + pl%inner(NTPHENC)%x(j, i) + tp%vh(j, tpind(1:nenci)) = tpenci%vh(j,1:nenci) + pl%inner(NTPHENC)%v(j, i) + end do + tp%lperi(tpind(1:nenci)) = tpenci%lperi(1:nenci) + tp%plperP(tpind(1:nenci)) = tpenci%plperP(1:nenci) + deallocate(pl%planetocentric(i)%tp) + deallocate(cbenci%inner) + do inner_index = 0, NTPHENC + deallocate(plenci%inner(inner_index)%x) + deallocate(plenci%inner(inner_index)%v) + if (allocated(plenci%inner(inner_index)%aobl)) deallocate(plenci%inner(inner_index)%aobl) + if (allocated(plenci%inner(inner_index)%atide)) deallocate(plenci%inner(inner_index)%atide) + end do + end associate end select end select end select diff --git a/src/symba/symba_collision.f90 b/src/symba/symba_collision.f90 index 7c49009f2..b00a4a8e6 100644 --- a/src/symba/symba_collision.f90 +++ b/src/symba/symba_collision.f90 @@ -841,32 +841,32 @@ subroutine symba_collision_mergeaddsub(system, param, family, id_frag, Ip_frag, ibiggest = maxloc(pl%Gmass(family(:)), dim=1) ! Copy over identification, information, and physical properties of the new bodies from the fragment list - plnew%id(:) = id_frag(:) + plnew%id(1:nfrag) = id_frag(1:nfrag) system%maxid = system%maxid + nfrag - plnew%xb(:,:) = xb_frag(:, :) - plnew%vb(:,:) = vb_frag(:, :) + plnew%xb(:, 1:nfrag) = xb_frag(:, 1:nfrag) + plnew%vb(:, 1:nfrag) = vb_frag(:, 1:nfrag) do i = 1, nfrag plnew%xh(:,i) = xb_frag(:, i) - cb%xb(:) plnew%vh(:,i) = vb_frag(:, i) - cb%vb(:) end do - plnew%mass(:) = m_frag(:) - plnew%Gmass(:) = param%GU * m_frag(:) - plnew%radius(:) = rad_frag(:) - plnew%density(:) = m_frag(:) / rad_frag(:) + plnew%mass(1:nfrag) = m_frag(1:nfrag) + plnew%Gmass(1:nfrag) = param%GU * m_frag(1:nfrag) + plnew%radius(1:nfrag) = rad_frag(1:nfrag) + plnew%density(1:nfrag) = m_frag(1:nfrag) / rad_frag(1:nfrag) select case(status) case(DISRUPTION) - plnew%info(:)%origin_type = "Disruption" - plnew%status(:) = NEW_PARTICLE - plnew%info(:)%origin_time = param%t + plnew%info(1:nfrag)%origin_type = "Disruption" + plnew%status(1:nfrag) = NEW_PARTICLE + plnew%info(1:nfrag)%origin_time = param%t do i = 1, nfrag plnew%info(i)%origin_xh(:) = plnew%xh(:,i) plnew%info(i)%origin_vh(:) = plnew%vh(:,i) end do case(SUPERCATASTROPHIC) - plnew%info(:)%origin_type = "Supercatastrophic" - plnew%status(:) = NEW_PARTICLE - plnew%info(:)%origin_time = param%t + plnew%info(1:nfrag)%origin_type = "Supercatastrophic" + plnew%status(1:nfrag) = NEW_PARTICLE + plnew%info(1:nfrag)%origin_time = param%t do i = 1, nfrag plnew%info(i)%origin_xh(:) = plnew%xh(:,i) plnew%info(i)%origin_vh(:) = plnew%vh(:,i) @@ -887,8 +887,8 @@ subroutine symba_collision_mergeaddsub(system, param, family, id_frag, Ip_frag, end select if (param%lrotation) then - plnew%Ip(:,:) = Ip_frag(:,:) - plnew%rot(:,:) = rot_frag(:,:) + plnew%Ip(:, 1:nfrag) = Ip_frag(:, 1:nfrag) + plnew%rot(:, 1:nfrag) = rot_frag(:, 1:nfrag) end if if (param%ltides) then @@ -899,11 +899,11 @@ subroutine symba_collision_mergeaddsub(system, param, family, id_frag, Ip_frag, call plnew%set_mu(cb) !Copy over or set integration parameters for new bodies - plnew%lcollision(:) = .false. - plnew%ldiscard(:) = .false. - plnew%lmtiny(:) = plnew%Gmass(:) > param%GMTINY - plnew%levelg(:) = pl%levelg(ibiggest) - plnew%levelm(:) = pl%levelm(ibiggest) + plnew%lcollision(1:nfrag) = .false. + plnew%ldiscard(1:nfrag) = .false. + plnew%lmtiny(1:nfrag) = plnew%Gmass(1:nfrag) > param%GMTINY + plnew%levelg(1:nfrag) = pl%levelg(ibiggest) + plnew%levelm(1:nfrag) = pl%levelm(ibiggest) ! Append the new merged body to the list and record how many we made nstart = pl_adds%nbody + 1 diff --git a/src/symba/symba_drift.f90 b/src/symba/symba_drift.f90 index c4efee05f..e3bbc3660 100644 --- a/src/symba/symba_drift.f90 +++ b/src/symba/symba_drift.f90 @@ -14,13 +14,14 @@ module subroutine symba_drift_pl(self, system, param, dt) real(DP), intent(in) :: dt !! Stepsize if (self%nbody == 0) return - - select type(system) - class is (symba_nbody_system) - self%lmask(:) = self%status(:) /= INACTIVE .and. self%levelg(:) == system%irec - call helio_drift_body(self, system, param, dt) - self%lmask(:) = self%status(:) /= INACTIVE - end select + associate(pl => self, npl => self%nbody) + select type(system) + class is (symba_nbody_system) + pl%lmask(1:npl) = pl%status(1:npl) /= INACTIVE .and. pl%levelg(1:npl) == system%irec + call helio_drift_body(pl, system, param, dt) + pl%lmask(1:npl) = pl%status(1:npl) /= INACTIVE + end select + end associate return end subroutine symba_drift_pl @@ -38,13 +39,14 @@ module subroutine symba_drift_tp(self, system, param, dt) real(DP), intent(in) :: dt !! Stepsize if (self%nbody == 0) return - - select type(system) - class is (symba_nbody_system) - self%lmask(:) = self%status(:) /= INACTIVE .and. self%levelg(:) == system%irec - call helio_drift_body(self, system, param, dt) - self%lmask(:) = self%status(:) /= INACTIVE - end select + associate (tp => self, ntp => self%nbody) + select type(system) + class is (symba_nbody_system) + tp%lmask(1:ntp) = tp%status(1:ntp) /= INACTIVE .and. tp%levelg(1:ntp) == system%irec + call helio_drift_body(tp, system, param, dt) + tp%lmask(1:ntp) = tp%status(1:ntp) /= INACTIVE + end select + end associate return end subroutine symba_drift_tp diff --git a/src/symba/symba_encounter_check.f90 b/src/symba/symba_encounter_check.f90 index eb230b7e0..fe593bf7c 100644 --- a/src/symba/symba_encounter_check.f90 +++ b/src/symba/symba_encounter_check.f90 @@ -177,14 +177,14 @@ module function symba_encounter_check_tp(self, system, dt, irec) result(lany_enc call pltpenc_list%resize(nenc) 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=1, ncopies=ntp), lencounter(:,:)) - pltpenc_list%index2(1:nenc) = pack(spread([(i, i = 1, ntp)], dim=2, ncopies=npl), lencounter(:,:)) + pltpenc_list%lvdotr(1:nenc) = pack(loc_lvdotr(1:ntp, 1:npl), lencounter(1:ntp, 1:npl)) + pltpenc_list%index1(1:nenc) = pack(spread([(i, i = 1, npl)], dim=1, ncopies=ntp), lencounter(1:ntp, 1:npl)) + pltpenc_list%index2(1:nenc) = pack(spread([(i, i = 1, ntp)], dim=2, ncopies=npl), lencounter(1:ntp, 1:npl)) pltpenc_list%id1(1:nenc) = pl%id(pltpenc_list%index1(1:nenc)) pltpenc_list%id2(1:nenc) = tp%id(pltpenc_list%index2(1:nenc)) select type(pl) class is (symba_pl) - pl%lencounter(:) = .false. + pl%lencounter(1:npl) = .false. do k = 1, nenc pl%lencounter(pltpenc_list%index1(k)) = .true. end do diff --git a/src/symba/symba_io.f90 b/src/symba/symba_io.f90 index 775efa962..5ab525e23 100644 --- a/src/symba/symba_io.f90 +++ b/src/symba/symba_io.f90 @@ -273,34 +273,36 @@ module subroutine symba_io_read_particle(system, param) class is (symba_pl) select type(tp => system%tp) class is (symba_tp) - do - lmatch = .false. - read(LUN, err = 667, iomsg = errmsg) id - - if (idx == cb%id) then - read(LUN, err = 667, iomsg = errmsg) cb%info - lmatch = .true. - else - if (pl%nbody > 0) then - idx = findloc(pl%id(:), id, dim=1) - if (idx /= 0) then - read(LUN, err = 667, iomsg = errmsg) pl%info(idx) - lmatch = .true. + associate(npl => pl%nbody, ntp => tp%nbody) + do + lmatch = .false. + read(LUN, err = 667, iomsg = errmsg) id + + if (idx == cb%id) then + read(LUN, err = 667, iomsg = errmsg) cb%info + lmatch = .true. + else + if (npl > 0) then + idx = findloc(pl%id(1:npl), id, dim=1) + if (idx /= 0) then + read(LUN, err = 667, iomsg = errmsg) pl%info(idx) + lmatch = .true. + end if end if - end if - if (.not.lmatch .and. tp%nbody > 0) then - idx = findloc(tp%id(:), id, dim=1) - if (idx /= 0) then - read(LUN, err = 667, iomsg = errmsg) tp%info(idx) - lmatch = .true. + if (.not.lmatch .and. ntp > 0) then + idx = findloc(tp%id(1:ntp), id, dim=1) + if (idx /= 0) then + read(LUN, err = 667, iomsg = errmsg) tp%info(idx) + lmatch = .true. + end if end if end if - end if - if (.not.lmatch) then - write(*,*) 'Particle id ',id,' not found. Skipping' - read(LUN, err = 667, iomsg = errmsg) tmpinfo - end if - end do + if (.not.lmatch) then + write(*,*) 'Particle id ',id,' not found. Skipping' + read(LUN, err = 667, iomsg = errmsg) tmpinfo + end if + end do + end associate close(unit = LUN, err = 667, iomsg = errmsg) end select end select diff --git a/src/symba/symba_kick.f90 b/src/symba/symba_kick.f90 index ccf56e039..1eabf9b6d 100644 --- a/src/symba/symba_kick.f90 +++ b/src/symba/symba_kick.f90 @@ -126,9 +126,9 @@ module subroutine symba_kick_encounter(self, system, dt, irec, sgn) class is (symba_pl) select type(tp => system%tp) class is (symba_tp) - associate(ind1 => self%index1, ind2 => self%index2) - if (pl%nbody > 0) pl%lmask(:) = pl%status(:) /= INACTIVE - if (tp%nbody > 0) tp%lmask(:) = tp%status(:) /= INACTIVE + associate(ind1 => self%index1, ind2 => self%index2, npl => pl%nbody, ntp => tp%nbody, nenc => self%nenc) + if (npl > 0) pl%lmask(1:npl) = pl%status(1:npl) /= INACTIVE + if (ntp > 0) tp%lmask(1:ntp) = tp%status(1:ntp) /= INACTIVE irm1 = irec - 1 @@ -139,13 +139,13 @@ module subroutine symba_kick_encounter(self, system, dt, irec, sgn) end if if (isplpl) then - pl%ah(:,ind1(1:self%nenc)) = 0.0_DP - pl%ah(:,ind2(1:self%nenc)) = 0.0_DP + pl%ah(:,ind1(1:nenc)) = 0.0_DP + pl%ah(:,ind2(1:nenc)) = 0.0_DP else - tp%ah(:,ind2(1:self%nenc)) = 0.0_DP + tp%ah(:,ind2(1:nenc)) = 0.0_DP end if - do k = 1, self%nenc + do k = 1, nenc if (isplpl) then lgoodlevel = (pl%levelg(ind1(k)) >= irm1) .and. (pl%levelg(ind2(k)) >= irm1) else @@ -176,15 +176,15 @@ module subroutine symba_kick_encounter(self, system, dt, irec, sgn) faci = fac * pl%Gmass(ind1(k)) if (isplpl) then facj = fac * pl%Gmass(ind2(k)) - pl%ah(:,ind1(k)) = pl%ah(:,ind1(k)) + facj * dx(:) - pl%ah(:,ind2(k)) = pl%ah(:,ind2(k)) - faci * dx(:) + pl%ah(:, ind1(k)) = pl%ah(:, ind1(k)) + facj * dx(:) + pl%ah(:, ind2(k)) = pl%ah(:, ind2(k)) - faci * dx(:) else - tp%ah(:,ind2(k)) = tp%ah(:,ind2(k)) - faci * dx(:) + tp%ah(:, ind2(k)) = tp%ah(:, ind2(k)) - faci * dx(:) end if end if end do if (isplpl) then - do k = 1, self%nenc + do k = 1, nenc pl%vb(:,ind1(k)) = pl%vb(:,ind1(k)) + sgn * dt * pl%ah(:,ind1(k)) pl%vb(:,ind2(k)) = pl%vb(:,ind2(k)) + sgn * dt * pl%ah(:,ind2(k)) pl%ah(:,ind1(k)) = 0.0_DP diff --git a/src/whm/whm_kick.f90 b/src/whm/whm_kick.f90 index 2da00c332..f831ce15b 100644 --- a/src/whm/whm_kick.f90 +++ b/src/whm/whm_kick.f90 @@ -80,17 +80,17 @@ module subroutine whm_kick_getacch_tp(self, system, param, t, lbeg) system%lbeg = lbeg if (lbeg) then - ah0(:) = whm_kick_getacch_ah0(pl%Gmass(:), pl%xbeg(:,:), npl) + ah0(:) = whm_kick_getacch_ah0(pl%Gmass(1:npl), pl%xbeg(:, 1:npl), npl) do concurrent(i = 1:ntp, tp%lmask(i)) tp%ah(:, i) = tp%ah(:, i) + ah0(:) end do - call tp%accel_int(pl%Gmass(:), pl%xbeg(:,:), npl) + call tp%accel_int(pl%Gmass(1:npl), pl%xbeg(:, 1:npl), npl) else - ah0(:) = whm_kick_getacch_ah0(pl%Gmass(:), pl%xend(:,:), npl) + ah0(:) = whm_kick_getacch_ah0(pl%Gmass(1:npl), pl%xend(:, 1:npl), npl) do concurrent(i = 1:ntp, tp%lmask(i)) tp%ah(:, i) = tp%ah(:, i) + ah0(:) end do - call tp%accel_int(pl%Gmass(:), pl%xend(:,:), npl) + call tp%accel_int(pl%Gmass(1:npl), pl%xend(:, 1:npl), npl) end if if (param%loblatecb) call tp%accel_obl(system) @@ -213,13 +213,13 @@ module subroutine whm_kick_vh_pl(self, system, param, t, dt, lbeg) if (lbeg) then if (pl%lfirst) then call pl%h2j(cb) - pl%ah(:,:) = 0.0_DP + pl%ah(:, 1:npl) = 0.0_DP call pl%accel(system, param, t, lbeg) pl%lfirst = .false. end if call pl%set_beg_end(xbeg = pl%xh) else - pl%ah(:,:) = 0.0_DP + pl%ah(:, 1:npl) = 0.0_DP call pl%accel(system, param, t, lbeg) call pl%set_beg_end(xend = pl%xh) end if @@ -254,20 +254,16 @@ module subroutine whm_kick_vh_tp(self, system, param, t, dt, lbeg) associate(tp => self, ntp => self%nbody) if (tp%lfirst) then - where(tp%lmask(1:ntp)) - tp%ah(1,1:ntp) = 0.0_DP - tp%ah(2,1:ntp) = 0.0_DP - tp%ah(3,1:ntp) = 0.0_DP - end where + do concurrent(i = 1:ntp, tp%lmask(i)) + tp%ah(:, i) = 0.0_DP + end do call tp%accel(system, param, t, lbeg=.true.) tp%lfirst = .false. end if if (.not.lbeg) then - where(tp%lmask(1:ntp)) - tp%ah(1,1:ntp) = 0.0_DP - tp%ah(2,1:ntp) = 0.0_DP - tp%ah(3,1:ntp) = 0.0_DP - end where + do concurrent(i = 1:ntp, tp%lmask(i)) + tp%ah(:, i) = 0.0_DP + end do call tp%accel(system, param, t, lbeg) end if do concurrent(i = 1:ntp, tp%lmask(i)) From 70bfb25a8639475970bceac6787518850721c625 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Wed, 18 Aug 2021 16:48:59 -0400 Subject: [PATCH 3/6] Fixed some indexing issues with SyMBA recursion levels and tried to make the kick routine cleaner --- src/symba/symba_kick.f90 | 138 ++++++++++++++++++++++++--------------- src/symba/symba_step.f90 | 4 +- 2 files changed, 86 insertions(+), 56 deletions(-) diff --git a/src/symba/symba_kick.f90 b/src/symba/symba_kick.f90 index 1eabf9b6d..b531e447a 100644 --- a/src/symba/symba_kick.f90 +++ b/src/symba/symba_kick.f90 @@ -109,10 +109,12 @@ module subroutine symba_kick_encounter(self, system, dt, irec, sgn) integer(I4B), intent(in) :: irec !! Current recursion level integer(I4B), intent(in) :: sgn !! sign to be applied to acceleration ! Internals - integer(I4B) :: k, irm1, irecl + integer(I4B) :: i, irm1, irecl, ngood real(DP) :: r, rr, ri, ris, rim1, r2, ir3, fac, faci, facj real(DP), dimension(NDIM) :: dx - logical :: isplpl, lgoodlevel + logical :: isplpl + logical, dimension(self%nenc) :: lgoodlevel + integer(I4B), dimension(:), allocatable :: good_idx if (self%nenc == 0) return @@ -138,63 +140,91 @@ module subroutine symba_kick_encounter(self, system, dt, irec, sgn) irecl = irec end if - if (isplpl) then - pl%ah(:,ind1(1:nenc)) = 0.0_DP - pl%ah(:,ind2(1:nenc)) = 0.0_DP - else - tp%ah(:,ind2(1:nenc)) = 0.0_DP - end if - - do k = 1, nenc + do i = 1, nenc if (isplpl) then - lgoodlevel = (pl%levelg(ind1(k)) >= irm1) .and. (pl%levelg(ind2(k)) >= irm1) + lgoodlevel(i) = (pl%levelg(ind1(i)) >= irm1) .and. (pl%levelg(ind2(i)) >= irm1) else - lgoodlevel = (pl%levelg(ind1(k)) >= irm1) .and. (tp%levelg(ind2(k)) >= irm1) - end if - if ((self%status(k) /= INACTIVE) .and. lgoodlevel) then - if (isplpl) then - ri = ((pl%rhill(ind1(k)) + pl%rhill(ind2(k)))**2) * (RHSCALE**2) * (RSHELL**(2*irecl)) - rim1 = ri * (RSHELL**2) - dx(:) = pl%xh(:,ind2(k)) - pl%xh(:,ind1(k)) - else - ri = ((pl%rhill(ind1(k)))**2) * (RHSCALE**2) * (RSHELL**(2*irecl)) - rim1 = ri * (RSHELL**2) - dx(:) = tp%xh(:,ind2(k)) - pl%xh(:,ind1(k)) - end if - r2 = dot_product(dx(:), dx(:)) - if (r2 < rim1) then - fac = 0.0_DP - else if (r2 < ri) then - ris = sqrt(ri) - r = sqrt(r2) - rr = (ris - r) / (ris * (1.0_DP - RSHELL)) - fac = (r2**(-1.5_DP)) * (1.0_DP - 3 * (rr**2) + 2 * (rr**3)) - else - ir3 = 1.0_DP / (r2 * sqrt(r2)) - fac = ir3 - end if - faci = fac * pl%Gmass(ind1(k)) - if (isplpl) then - facj = fac * pl%Gmass(ind2(k)) - pl%ah(:, ind1(k)) = pl%ah(:, ind1(k)) + facj * dx(:) - pl%ah(:, ind2(k)) = pl%ah(:, ind2(k)) - faci * dx(:) - else - tp%ah(:, ind2(k)) = tp%ah(:, ind2(k)) - faci * dx(:) - end if + lgoodlevel(i) = (pl%levelg(ind1(i)) >= irm1) .and. (tp%levelg(ind2(i)) >= irm1) end if + lgoodlevel(i) = (self%status(i) == ACTIVE) .and. lgoodlevel(i) end do - if (isplpl) then - do k = 1, nenc - pl%vb(:,ind1(k)) = pl%vb(:,ind1(k)) + sgn * dt * pl%ah(:,ind1(k)) - pl%vb(:,ind2(k)) = pl%vb(:,ind2(k)) + sgn * dt * pl%ah(:,ind2(k)) - pl%ah(:,ind1(k)) = 0.0_DP - pl%ah(:,ind2(k)) = 0.0_DP - end do - else - do k = 1, self%nenc - tp%vb(:,ind2(k)) = tp%vb(:,ind2(k)) + sgn * dt * tp%ah(:,ind2(k)) - tp%ah(:,ind2(k)) = 0.0_DP + ngood = count(lgoodlevel(:)) + if (ngood > 0) then + allocate(good_idx(ngood)) + good_idx(:) = pack([(i, i = 1, nenc)], lgoodlevel(:)) + + if (isplpl) then + do i = 1, ngood + associate(i1 => ind1(good_idx(i)), i2 => ind2(good_idx(i))) + pl%ah(:,i1) = 0.0_DP + pl%ah(:,i2) = 0.0_DP + end associate + end do + else + do i = 1, ngood + associate(i2 => ind2(good_idx(i))) + tp%ah(:,i2) = 0.0_DP + end associate + end do + end if + + do i = 1, ngood + associate(k => good_idx(i), i1 => ind1(good_idx(i)), i2 => ind2(good_idx(i))) + if (isplpl) then + ri = ((pl%rhill(i1) + pl%rhill(i2))**2) * (RHSCALE**2) * (RSHELL**(2*irecl)) + rim1 = ri * (RSHELL**2) + dx(:) = pl%xh(:,i2) - pl%xh(:,i1) + else + ri = ((pl%rhill(i1))**2) * (RHSCALE**2) * (RSHELL**(2*irecl)) + rim1 = ri * (RSHELL**2) + dx(:) = tp%xh(:,i2) - pl%xh(:,i1) + end if + r2 = dot_product(dx(:), dx(:)) + if (r2 < rim1) then + fac = 0.0_DP + lgoodlevel(k) = .false. + cycle + end if + if (r2 < ri) then + ris = sqrt(ri) + r = sqrt(r2) + rr = (ris - r) / (ris * (1.0_DP - RSHELL)) + fac = (r2**(-1.5_DP)) * (1.0_DP - 3 * (rr**2) + 2 * (rr**3)) + else + ir3 = 1.0_DP / (r2 * sqrt(r2)) + fac = ir3 + end if + faci = fac * pl%Gmass(i1) + if (isplpl) then + facj = fac * pl%Gmass(i2) + pl%ah(:, i1) = pl%ah(:, i1) + facj * dx(:) + pl%ah(:, i2) = pl%ah(:, i2) - faci * dx(:) + else + tp%ah(:, i2) = tp%ah(:, i2) - faci * dx(:) + end if + end associate end do + ngood = count(lgoodlevel(:)) + if (ngood == 0) return + good_idx(1:ngood) = pack([(i, i = 1, nenc)], lgoodlevel(:)) + + if (isplpl) then + do i = 1, ngood + associate(i1 => ind1(good_idx(i)), i2 => ind2(good_idx(i))) + pl%vb(:,i1) = pl%vb(:,i1) + sgn * dt * pl%ah(:,i1) + pl%vb(:,i2) = pl%vb(:,i2) + sgn * dt * pl%ah(:,i2) + pl%ah(:,ind1(1:nenc)) = 0.0_DP + pl%ah(:,ind2(1:nenc)) = 0.0_DP + end associate + end do + else + do i = 1, ngood + associate(i2 => ind2(good_idx(i))) + tp%vb(:,i2) = tp%vb(:,i2) + sgn * dt * tp%ah(:,i2) + tp%ah(:,i2) = 0.0_DP + end associate + end do + end if end if end associate end select diff --git a/src/symba/symba_step.f90 b/src/symba/symba_step.f90 index 249c73a82..1988b4fbc 100644 --- a/src/symba/symba_step.f90 +++ b/src/symba/symba_step.f90 @@ -243,8 +243,8 @@ module subroutine symba_step_reset_system(self, param) end do pl%nplenc(:) = 0 pl%ntpenc(:) = 0 - pl%levelg(:) = 0 - pl%levelm(:) = 0 + pl%levelg(:) = -1 + pl%levelm(:) = -1 pl%lencounter = .false. pl%lcollision = .false. system%plplenc_list%nenc = 0 From 5d8588061616d8e3def9b585e7f3616399166fd7 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Wed, 18 Aug 2021 17:59:32 -0400 Subject: [PATCH 4/6] Fixed recursion level problem in SyMBA --- src/symba/symba_encounter_check.f90 | 26 +++++++++++++++++++++----- src/symba/symba_step.f90 | 2 +- 2 files changed, 22 insertions(+), 6 deletions(-) diff --git a/src/symba/symba_encounter_check.f90 b/src/symba/symba_encounter_check.f90 index fe593bf7c..372a40996 100644 --- a/src/symba/symba_encounter_check.f90 +++ b/src/symba/symba_encounter_check.f90 @@ -46,10 +46,18 @@ module function symba_encounter_check_pl(self, system, dt, irec) result(lany_enc plplenc_list%id1(1:nenc) = pl%id(plplenc_list%index1(1:nenc)) plplenc_list%id2(1:nenc) = pl%id(plplenc_list%index2(1:nenc)) do k = 1, nenc - plplenc_list%status(k) = ACTIVE - plplenc_list%level(k) = irec - pl%lencounter(plplenc_list%index1(k)) = .true. - pl%lencounter(plplenc_list%index2(k)) = .true. + associate(i => plplenc_list%index1(k), j => plplenc_list%index2(k)) + plplenc_list%status(k) = ACTIVE + plplenc_list%level(k) = irec + pl%lencounter(i) = .true. + pl%lencounter(j) = .true. + pl%levelg(i) = irec + pl%levelm(i) = irec + pl%levelg(j) = irec + pl%levelm(j) = irec + pl%nplenc(i) = pl%nplenc(i) + 1 + pl%nplenc(j) = pl%nplenc(j) + 1 + end associate end do end associate end if @@ -186,7 +194,15 @@ module function symba_encounter_check_tp(self, system, dt, irec) result(lany_enc class is (symba_pl) pl%lencounter(1:npl) = .false. do k = 1, nenc - pl%lencounter(pltpenc_list%index1(k)) = .true. + associate(plind => pltpenc_list%index1(k), tpind => pltpenc_list%index2(k)) + pl%lencounter(plind) = .true. + pl%levelg(plind) = irec + pl%levelm(plind) = irec + tp%levelg(tpind) = irec + tp%levelm(tpind) = irec + pl%ntpenc(plind) = pl%ntpenc(plind) + 1 + tp%nplenc(tpind) = tp%nplenc(tpind) + 1 + end associate end do end select end associate diff --git a/src/symba/symba_step.f90 b/src/symba/symba_step.f90 index 1988b4fbc..2b68c1f1d 100644 --- a/src/symba/symba_step.f90 +++ b/src/symba/symba_step.f90 @@ -255,7 +255,7 @@ module subroutine symba_step_reset_system(self, param) tp%nplenc(:) = 0 tp%levelg(:) = 0 tp%levelm(:) = 0 - system%pltpenc_list%nenc = 0 + system%pltpenc_list%nenc = 0 end if call system%pl_adds%setup(0, param) From b63c173bcf62014ee7439aa632bef82aa0a94211 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Wed, 18 Aug 2021 18:15:27 -0400 Subject: [PATCH 5/6] Fixed more encounter level depth index values --- src/symba/symba_step.f90 | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/symba/symba_step.f90 b/src/symba/symba_step.f90 index 2b68c1f1d..cde5aff29 100644 --- a/src/symba/symba_step.f90 +++ b/src/symba/symba_step.f90 @@ -33,6 +33,7 @@ module subroutine symba_step_system(self, param, t, dt) pl%lfirst = .true. tp%lfirst = .true. else + self%irec = -1 call helio_step_system(self, param, t, dt) end if end select @@ -253,8 +254,8 @@ module subroutine symba_step_reset_system(self, param) if (tp%nbody > 0) then tp%nplenc(:) = 0 - tp%levelg(:) = 0 - tp%levelm(:) = 0 + tp%levelg(:) = -1 + tp%levelm(:) = -1 system%pltpenc_list%nenc = 0 end if From d709d55591586909fba5a55d16d6459b553d2b5d Mon Sep 17 00:00:00 2001 From: David A Minton Date: Wed, 18 Aug 2021 20:29:49 -0400 Subject: [PATCH 6/6] Maintenence on SyMBA to try to track down an auto-paralellization bug --- Makefile.Defines | 12 +++---- src/modules/symba_classes.f90 | 4 +-- src/symba/symba_step.f90 | 60 +++++++++++++++++++---------------- 3 files changed, 41 insertions(+), 35 deletions(-) diff --git a/Makefile.Defines b/Makefile.Defines index 4a283b655..7bb6ae32e 100644 --- a/Makefile.Defines +++ b/Makefile.Defines @@ -54,24 +54,24 @@ VTUNE_FLAGS = -g -O2 -vec -simd -shared-intel -qopenmp -debug inline-debug-info IDEBUG = -O0 -nogen-interfaces -no-pie -no-ftz -fpe-all=0 -g -traceback -mp1 -fp-model strict -fpe0 -debug all -align all -pad -ip -prec-div -prec-sqrt -assume protect-parens -CB -no-wrap-margin -init=snan,arrays STRICTREAL = -fp-model strict -fp-model no-except -prec-div -prec-sqrt -assume protect-parens SIMDVEC = -simd -xhost -align all -assume contiguous_assumed_shape -vecabi=cmdtarget -prec-div -prec-sqrt -assume protect-parens -PAR = -parallel -qopenmp +PAR = -qopenmp #-parallel Something goes wrong in SyMBA at the moment with auto-paralellization enabled HEAPARR = -heap-arrays 1048576 -OPTIMIZE = -qopt-report=5 +OPTREPORT = -qopt-report=5 #gfortran flags -GDEBUG = -g -O0 -fbacktrace -fbounds-check +GDEBUG = -g -Og -fbacktrace -fbounds-check GPRODUCTION = -O3 GPAR = -fopenmp -ftree-parallelize-loops=4 GMEM = -fsanitize=undefined -fsanitize=address -fsanitize=leak GWARNINGS = -Wall -Warray-bounds -Wimplicit-interface -Wextra -Warray-temporaries -#FFLAGS = $(IDEBUG) $(HEAPARR) -FFLAGS = -init=snan,arrays -no-wrap-margin -O3 $(STRICTREAL) $(SIMDVEC) $(PAR) $(HEAPARR) +#FFLAGS = $(IDEBUG) $(HEAPARR) $(SIMDVEC) $(PAR) +FFLAGS = -init=snan,arrays -no-wrap-margin -O3 $(STRICTREAL) $(PAR) $(SIMDVEC) $(HEAPARR) FORTRAN = ifort #AR = xiar #FORTRAN = gfortran -#FFLAGS = -ffree-line-length-none $(GDEBUG) #$(GMEM) +#FFLAGS = -ffree-line-length-none $(GPRODUCTION) $(GPAR) #$(GDEBUG) $(GMEM) AR = ar # DO NOT include in CFLAGS the "-c" option to compile object only diff --git a/src/modules/symba_classes.f90 b/src/modules/symba_classes.f90 index 05698d44c..056e6f390 100644 --- a/src/modules/symba_classes.f90 +++ b/src/modules/symba_classes.f90 @@ -497,8 +497,8 @@ module recursive subroutine symba_step_recur_system(self, param, t, ireci) implicit none class(symba_nbody_system), intent(inout) :: self !! SyMBA nbody system object class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters - real(DP), value :: t - integer(I4B), value :: ireci !! input recursion level + real(DP), intent(in) :: t + integer(I4B), intent(in) :: ireci !! input recursion level end subroutine symba_step_recur_system module subroutine symba_step_reset_system(self, param) diff --git a/src/symba/symba_step.f90 b/src/symba/symba_step.f90 index cde5aff29..ed5985c0b 100644 --- a/src/symba/symba_step.f90 +++ b/src/symba/symba_step.f90 @@ -146,8 +146,8 @@ module recursive subroutine symba_step_recur_system(self, param, t, ireci) ! Arguments class(symba_nbody_system), intent(inout) :: self !! SyMBA nbody system object class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters - real(DP), value :: t - integer(I4B), value :: ireci !! input recursion level + real(DP), intent(in) :: t + integer(I4B), intent(in) :: ireci !! input recursion level ! Internals integer(I4B) :: i, j, irecp, nloops real(DP) :: dtl, dth @@ -235,32 +235,38 @@ module subroutine symba_step_reset_system(self, param) class is (symba_pl) select type(tp => system%tp) class is (symba_tp) - if (pl%nbody > 0) then - pl%lcollision(:) = .false. - pl%kin(:)%parent = [(i, i=1, pl%nbody)] - pl%kin(:)%nchild = 0 - do i = 1, pl%nbody - if (allocated(pl%kin(i)%child)) deallocate(pl%kin(i)%child) - end do - pl%nplenc(:) = 0 - pl%ntpenc(:) = 0 - pl%levelg(:) = -1 - pl%levelm(:) = -1 - pl%lencounter = .false. - pl%lcollision = .false. - system%plplenc_list%nenc = 0 - system%plplcollision_list%nenc = 0 - end if - - if (tp%nbody > 0) then - tp%nplenc(:) = 0 - tp%levelg(:) = -1 - tp%levelm(:) = -1 - system%pltpenc_list%nenc = 0 - end if + associate(npl => pl%nbody, ntp => tp%nbody) + if (npl > 0) then + pl%lcollision(1:npl) = .false. + pl%kin(1:npl)%parent = [(i, i=1, npl)] + pl%kin(1:npl)%nchild = 0 + do i = 1, npl + if (allocated(pl%kin(i)%child)) deallocate(pl%kin(i)%child) + end do + pl%nplenc(1:npl) = 0 + pl%ntpenc(1:npl) = 0 + pl%levelg(1:npl) = -1 + pl%levelm(1:npl) = -1 + pl%lencounter(1:npl) = .false. + pl%lcollision(1:npl) = .false. + pl%ldiscard(1:npl) = .false. + pl%lmask(1:npl) = .true. + system%plplenc_list%nenc = 0 + system%plplcollision_list%nenc = 0 + end if + + if (ntp > 0) then + tp%nplenc(1:ntp) = 0 + tp%levelg(1:ntp) = -1 + tp%levelm(1:ntp) = -1 + tp%lmask(1:ntp) = .true. + pl%ldiscard(1:npl) = .false. + system%pltpenc_list%nenc = 0 + end if - call system%pl_adds%setup(0, param) - call system%pl_discards%setup(0, param) + call system%pl_adds%setup(0, param) + call system%pl_discards%setup(0, param) + end associate end select end select end associate