Skip to content
This repository was archived by the owner on Aug 28, 2024. It is now read-only.

Commit

Permalink
Browse files Browse the repository at this point in the history 
Moved associated index variables outside of loops (the new version of ifort doesn't play nicely with associate blocks inside of do concurrent, apparently
  • Loading branch information
@daminton
daminton committed Jul 31, 2021
1 parent 8f771a1 commit a935642
Show file tree
Hide file tree
Showing 2 changed files with 53 additions and 53 deletions.
25 changes: 15 additions & 10 deletions src/symba/symba_collision.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -41,7 +41,7 @@ module subroutine symba_collision_check_pltpenc(self, system, param, t, dt, irec
real(DP), intent(in) :: dt !! step size
integer(I4B), intent(in) :: irec !! Current recursion level
! Internals
logical, dimension(:), allocatable :: lcollision, mask
logical, dimension(:), allocatable :: lcollision, lmask
real(DP), dimension(NDIM) :: xr, vr
integer(I4B) :: k

Expand All @@ -51,16 +51,20 @@ module subroutine symba_collision_check_pltpenc(self, system, param, t, dt, irec
class is (symba_pl)
select type(tp => system%tp)
class is (symba_tp)
associate(pltpenc_list => self, npltpenc => self%nenc, plind => self%index1(1:self%nenc), tpind => self%index2(1:self%nenc))
allocate(lcollision(npltpenc), mask(npltpenc))
mask(:) = ((pltpenc_list%status(1:npltpenc) == ACTIVE) .and. (pl%levelg(plind) >= irec) .and. (tp%levelg(tpind) >= irec))
associate(pltpenc_list => self, npltpenc => self%nenc, plind => self%index1, tpind => self%index2 )
allocate(lmask(npltpenc))
lmask(:) = ((pltpenc_list%status(1:npltpenc) == ACTIVE) &
.and. (pl%levelg(plind(1:npltpenc)) >= irec) &
.and. (tp%levelg(tpind(1:npltpenc)) >= irec))
if (.not.any(lmask(:))) return

allocate(lcollision(npltpenc))
lcollision(:) = .false.
do concurrent(k = 1:npltpenc, mask(k))
associate(i => plind(k), j => tpind(k))
xr(:) = pl%xh(:, i) - tp%xh(:, j)
vr(:) = pl%vb(:, i) - tp%vb(:, j)
lcollision(k) = symba_collision_check_one(xr(1), xr(2), xr(3), vr(1), vr(2), vr(3), pl%Gmass(i), pl%radius(i), dt, pltpenc_list%lvdotr(k))
end associate

do concurrent(k = 1:npltpenc, lmask(k))
xr(:) = pl%xh(:, plind(k)) - tp%xh(:, tpind(k))
vr(:) = pl%vb(:, plind(k)) - tp%vb(:, tpind(k))
lcollision(k) = symba_collision_check_one(xr(1), xr(2), xr(3), vr(1), vr(2), vr(3), pl%Gmass(plind(k)), pl%radius(plind(k)), dt, pltpenc_list%lvdotr(k))
end do

if (any(lcollision(:))) then
Expand All @@ -69,6 +73,7 @@ module subroutine symba_collision_check_pltpenc(self, system, param, t, dt, irec
tp%status(tpind(1:npltpenc)) = DISCARDED_PLR
tp%ldiscard(tpind(1:npltpenc)) = .true.
end where

do k = 1, npltpenc
if (pltpenc_list%status(k) /= COLLISION) cycle
write(*,*) 'Test particle ',tp%id(tpind(k)), ' collided with massive body ',pl%id(plind(k)), ' at time ',t
Expand Down
81 changes: 38 additions & 43 deletions src/symba/symba_kick.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -123,38 +123,37 @@ module subroutine symba_kick_pltpenc(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(:) == ACTIVE
if (tp%nbody > 0) tp%lmask(:) = tp%status(:) == ACTIVE

if (pl%nbody > 0) pl%lmask(:) = pl%status(:) == ACTIVE
if (tp%nbody > 0) tp%lmask(:) = tp%status(:) == ACTIVE

irm1 = irec - 1
if (sgn < 0) then
irecl = irec - 1
else
irecl = irec
end if
do k = 1, self%nenc
associate(i => self%index1(k), j => self%index2(k))
irm1 = irec - 1
if (sgn < 0) then
irecl = irec - 1
else
irecl = irec
end if
do k = 1, self%nenc
if (isplpl) then
pl%ah(:,i) = 0.0_DP
pl%ah(:,j) = 0.0_DP
pl%ah(:,ind1(k)) = 0.0_DP
pl%ah(:,ind2(k)) = 0.0_DP
else
tp%ah(:,j) = 0.0_DP
tp%ah(:,ind2(k)) = 0.0_DP
end if
if (isplpl) then
lgoodlevel = (pl%levelg(i) >= irm1) .and. (pl%levelg(j) >= irm1)
lgoodlevel = (pl%levelg(ind1(k)) >= irm1) .and. (pl%levelg(ind2(k)) >= irm1)
else
lgoodlevel = (pl%levelg(i) >= irm1) .and. (tp%levelg(j) >= irm1)
lgoodlevel = (pl%levelg(ind1(k)) >= irm1) .and. (tp%levelg(ind2(k)) >= irm1)
end if
if ((self%status(k) == ACTIVE) .and. lgoodlevel) then
if (isplpl) then
ri = ((pl%rhill(i) + pl%rhill(j))**2) * (RHSCALE**2) * (RSHELL**(2*irecl))
ri = ((pl%rhill(ind1(k)) + pl%rhill(ind2(k)))**2) * (RHSCALE**2) * (RSHELL**(2*irecl))
rim1 = ri * (RSHELL**2)
dx(:) = pl%xh(:,j) - pl%xh(:,i)
dx(:) = pl%xh(:,ind2(k)) - pl%xh(:,ind1(k))
else
ri = ((pl%rhill(i))**2) * (RHSCALE**2) * (RSHELL**(2*irecl))
ri = ((pl%rhill(ind1(k)))**2) * (RHSCALE**2) * (RSHELL**(2*irecl))
rim1 = ri * (RSHELL**2)
dx(:) = tp%xh(:,j) - pl%xh(:,i)
dx(:) = tp%xh(:,ind2(k)) - pl%xh(:,ind1(k))
end if
r2 = dot_product(dx(:), dx(:))
if (r2 < rim1) then
Expand All @@ -168,34 +167,30 @@ module subroutine symba_kick_pltpenc(self, system, dt, irec, sgn)
ir3 = 1.0_DP / (r2 * sqrt(r2))
fac = ir3
end if
faci = fac * pl%Gmass(i)
faci = fac * pl%Gmass(ind1(k))
if (isplpl) then
facj = fac * pl%Gmass(j)
pl%ah(:,i) = pl%ah(:,i) + facj*dx(:)
pl%ah(:,j) = pl%ah(:,j) - faci*dx(:)
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(:,j) = tp%ah(:,j) - faci*dx(:)
tp%ah(:,ind2(k)) = tp%ah(:,ind2(k)) - faci * dx(:)
end if
end if
end associate
end do
if (isplpl) then
do k = 1, self%nenc
associate(i => self%index1(k), j => self%index2(k))
pl%vb(:,i) = pl%vb(:,i) + sgn * dt * pl%ah(:,i)
pl%vb(:,j) = pl%vb(:,j) + sgn * dt * pl%ah(:,j)
pl%ah(:,i) = 0.0_DP
pl%ah(:,j) = 0.0_DP
end associate
end do
else
where(tp%lmask(self%index2(1:self%nenc)))
tp%vb(1,self%index2(:)) = tp%vb(1,self%index2(:)) + sgn * dt * tp%ah(1,self%index2(:))
tp%vb(2,self%index2(:)) = tp%vb(2,self%index2(:)) + sgn * dt * tp%ah(2,self%index2(:))
tp%vb(3,self%index2(:)) = tp%vb(3,self%index2(:)) + sgn * dt * tp%ah(3,self%index2(:))
end where
tp%ah(:,self%index2(1:self%nenc)) = 0.0_DP
end if
if (isplpl) then
do k = 1, self%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(:,ind1(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
end do
end if
end associate
end select
end select

Expand Down

0 comments on commit a935642

Please sign in to comment.