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

Commit

Permalink
Merge branch 'debug' into IntelAdvisor
Browse files Browse the repository at this point in the history
  • Loading branch information
@daminton
daminton committed Oct 20, 2021
2 parents 6217e3a + 2cb3812 commit bdb9f0c
Show file tree
Hide file tree
Showing 6 changed files with 134 additions and 79 deletions.
22 changes: 16 additions & 6 deletions python/swiftest/swiftest/io.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -409,7 +409,8 @@ def make_swiftest_labels(param):
tlab.append('capm')
plab = tlab.copy()
plab.append('Gmass')
plab.append('radius')
if param['CHK_CLOSE'] == 'YES':
plab.append('radius')
if param['RHILL_PRESENT'] == 'YES':
plab.append('rhill')
clab = ['Gmass', 'radius', 'J_2', 'J_4']
Expand Down Expand Up @@ -817,7 +818,10 @@ def swiftest_xr2infile(ds, param, framenum=-1):
tp = frame.where(np.isnan(frame['Gmass']), drop=True).drop_vars(['Gmass', 'radius', 'J_2', 'J_4'])

GMSun = np.double(cb['Gmass'])
RSun = np.double(cb['radius'])
if param['CHK_CLOSE'] == 'YES':
RSun = np.double(cb['radius'])
else:
RSun = param['CHK_RMIN']
J2 = np.double(cb['J_2'])
J4 = np.double(cb['J_4'])
cbname = cb['name'].values[0]
Expand Down Expand Up @@ -851,7 +855,8 @@ def swiftest_xr2infile(ds, param, framenum=-1):
print(pli['name'].values, pli['Gmass'].values, pli['rhill'].values, file=plfile)
else:
print(pli['name'].values, pli['Gmass'].values, file=plfile)
print(pli['radius'].values, file=plfile)
if param['CHK_CLOSE'] == 'YES':
print(pli['radius'].values, file=plfile)
if param['IN_FORM'] == 'XV':
print(pli['xhx'].values, pli['xhy'].values, pli['xhz'].values, file=plfile)
print(pli['vhx'].values, pli['vhy'].values, pli['vhz'].values, file=plfile)
Expand Down Expand Up @@ -918,7 +923,8 @@ def swiftest_xr2infile(ds, param, framenum=-1):
else:
print(f"{param['IN_FORM']} is not a valid input format type.")
Gmass = pl['Gmass'].values
radius = pl['radius'].values
if param['CHK_CLOSE'] == 'YES':
radius = pl['radius'].values

plfile.write_record(npl)
plfile.write_record(plid)
Expand All @@ -931,7 +937,8 @@ def swiftest_xr2infile(ds, param, framenum=-1):
plfile.write_record(Gmass)
if param['RHILL_PRESENT'] == 'YES':
plfile.write_record(pl['rhill'].values)
plfile.write_record(radius)
if param['CHK_CLOSE'] == 'YES':
plfile.write_record(radius)
if param['ROTATION'] == 'YES':
plfile.write_record(pl['Ip1'].values)
plfile.write_record(pl['Ip2'].values)
Expand Down Expand Up @@ -995,7 +1002,10 @@ def swifter_xr2infile(ds, param, framenum=-1):
tp = frame.where(np.isnan(frame['Gmass']), drop=True).drop_vars(['Gmass', 'radius', 'J_2', 'J_4'])

GMSun = np.double(cb['Gmass'])
RSun = np.double(cb['radius'])
if param['CHK_CLOSE'] == 'YES':
RSun = np.double(cb['radius'])
else:
RSun = param['CHK_RMIN']
param['J2'] = np.double(cb['J_2'])
param['J4'] = np.double(cb['J_4'])

Expand Down
6 changes: 3 additions & 3 deletions src/io/io.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1525,7 +1525,7 @@ module function io_read_frame_body(self, iu, param) result(ierr)
read(iu, err = 667, iomsg = errmsg) pl%Gmass(:)
pl%mass(:) = pl%Gmass(:) / param%GU
if (param%lrhill_present) read(iu, err = 667, iomsg = errmsg) pl%rhill(:)
read(iu, err = 667, iomsg = errmsg) pl%radius(:)
if (param%lclose) read(iu, err = 667, iomsg = errmsg) pl%radius(:)
if (param%lrotation) then
read(iu, err = 667, iomsg = errmsg) pl%Ip(1, :)
read(iu, err = 667, iomsg = errmsg) pl%Ip(2, :)
Expand Down Expand Up @@ -1553,7 +1553,7 @@ module function io_read_frame_body(self, iu, param) result(ierr)
end if
self%Gmass(i) = real(val, kind=DP)
self%mass(i) = real(val / param%GU, kind=DP)
read(iu, *, err = 667, iomsg = errmsg) self%radius(i)
if (param%lclose) read(iu, *, err = 667, iomsg = errmsg) self%radius(i)
class is (swiftest_tp)
read(iu, *, err = 667, iomsg = errmsg) name(i)
end select
Expand Down Expand Up @@ -2027,7 +2027,7 @@ module subroutine io_write_frame_body(self, iu, param)
class is (swiftest_pl) ! Additional output if the passed polymorphic object is a massive body
write(iu, err = 667, iomsg = errmsg) pl%Gmass(1:n)
if (param%lrhill_present) write(iu, err = 667, iomsg = errmsg) pl%rhill(1:n)
write(iu, err = 667, iomsg = errmsg) pl%radius(1:n)
if (param%lclose) write(iu, err = 667, iomsg = errmsg) pl%radius(1:n)
if (param%lrotation) then
write(iu, err = 667, iomsg = errmsg) pl%Ip(1, 1:n)
write(iu, err = 667, iomsg = errmsg) pl%Ip(2, 1:n)
Expand Down
142 changes: 94 additions & 48 deletions src/kick/kick.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -33,9 +33,17 @@ module subroutine kick_getacch_int_pl(self, param)
end if

if (param%lflatten_interactions) then
call kick_getacch_int_all_flat_pl(self%nbody, self%nplpl, self%k_plpl, self%xh, self%Gmass, self%radius, self%ah)
if (param%lclose) then
call kick_getacch_int_all_flat_pl(self%nbody, self%nplpl, self%k_plpl, self%xh, self%Gmass, self%radius, self%ah)
else
call kick_getacch_int_all_flat_pl(self%nbody, self%nplpl, self%k_plpl, self%xh, self%Gmass, acc=self%ah)
end if
else
call kick_getacch_int_all_triangular_pl(self%nbody, self%nbody, self%xh, self%Gmass, self%radius, self%ah)
if (param%lclose) then
call kick_getacch_int_all_triangular_pl(self%nbody, self%nbody, self%xh, self%Gmass, self%radius, self%ah)
else
call kick_getacch_int_all_triangular_pl(self%nbody, self%nbody, self%xh, self%Gmass, acc=self%ah)
end if
end if

if (param%ladaptive_interactions .and. self%nplpl > 0) then
Expand Down Expand Up @@ -78,13 +86,13 @@ module subroutine kick_getacch_int_all_flat_pl(npl, nplpl, k_plpl, x, Gmass, rad
!! Adapted from Hal Levison's Swift routine getacch_ah3.f
!! Adapted from David E. Kaufmann's Swifter routine whm_kick_getacch_ah3.f90 and helio_kick_getacch_int.f9
implicit none
integer(I4B), intent(in) :: npl !! Number of massive bodies
integer(I8B), intent(in) :: nplpl !! Number of massive body interactions to compute
integer(I4B), dimension(:,:), intent(in) :: k_plpl !! Array of interaction pair indices (flattened upper triangular matrix)
real(DP), dimension(:,:), intent(in) :: x !! Position vector array
real(DP), dimension(:), intent(in) :: Gmass !! Array of massive body G*mass
real(DP), dimension(:), intent(in) :: radius !! Array of massive body radii
real(DP), dimension(:,:), intent(inout) :: acc !! Acceleration vector array
integer(I4B), intent(in) :: npl !! Number of massive bodies
integer(I8B), intent(in) :: nplpl !! Number of massive body interactions to compute
integer(I4B), dimension(:,:), intent(in) :: k_plpl !! Array of interaction pair indices (flattened upper triangular matrix)
real(DP), dimension(:,:), intent(in) :: x !! Position vector array
real(DP), dimension(:), intent(in) :: Gmass !! Array of massive body G*mass
real(DP), dimension(:), intent(in), optional :: radius !! Array of massive body radii
real(DP), dimension(:,:), intent(inout) :: acc !! Acceleration vector array
! Internals
integer(I8B) :: k
real(DP), dimension(NDIM,npl) :: ahi, ahj
Expand All @@ -95,23 +103,42 @@ module subroutine kick_getacch_int_all_flat_pl(npl, nplpl, k_plpl, x, Gmass, rad
ahi(:,:) = 0.0_DP
ahj(:,:) = 0.0_DP

!$omp parallel do default(private) schedule(static)&
!$omp shared(nplpl, k_plpl, x, Gmass, radius) &
!$omp lastprivate(rji2, rlim2, xr, yr, zr) &
!$omp reduction(+:ahi) &
!$omp reduction(-:ahj)
do k = 1_I8B, nplpl
i = k_plpl(1, k)
j = k_plpl(2, k)
xr = x(1, j) - x(1, i)
yr = x(2, j) - x(2, i)
zr = x(3, j) - x(3, i)
rji2 = xr**2 + yr**2 + zr**2
rlim2 = (radius(i) + radius(j))**2
if (rji2 > rlim2) call kick_getacch_int_one_pl(rji2, xr, yr, zr, Gmass(i), Gmass(j), &
ahi(1,i), ahi(2,i), ahi(3,i), ahj(1,j), ahj(2,j), ahj(3,j))
end do
!$omp end parallel do
if (present(radius)) then
!$omp parallel do default(private) schedule(static)&
!$omp shared(nplpl, k_plpl, x, Gmass, radius) &
!$omp lastprivate(rji2, rlim2, xr, yr, zr) &
!$omp reduction(+:ahi) &
!$omp reduction(-:ahj)
do k = 1_I8B, nplpl
i = k_plpl(1, k)
j = k_plpl(2, k)
xr = x(1, j) - x(1, i)
yr = x(2, j) - x(2, i)
zr = x(3, j) - x(3, i)
rji2 = xr**2 + yr**2 + zr**2
rlim2 = (radius(i) + radius(j))**2
if (rji2 > rlim2) call kick_getacch_int_one_pl(rji2, xr, yr, zr, Gmass(i), Gmass(j), &
ahi(1,i), ahi(2,i), ahi(3,i), ahj(1,j), ahj(2,j), ahj(3,j))
end do
!$omp end parallel do
else
!$omp parallel do default(private) schedule(static)&
!$omp shared(nplpl, k_plpl, x, Gmass, radius) &
!$omp lastprivate(rji2, xr, yr, zr) &
!$omp reduction(+:ahi) &
!$omp reduction(-:ahj)
do k = 1_I8B, nplpl
i = k_plpl(1, k)
j = k_plpl(2, k)
xr = x(1, j) - x(1, i)
yr = x(2, j) - x(2, i)
zr = x(3, j) - x(3, i)
rji2 = xr**2 + yr**2 + zr**2
call kick_getacch_int_one_pl(rji2, xr, yr, zr, Gmass(i), Gmass(j), &
ahi(1,i), ahi(2,i), ahi(3,i), ahj(1,j), ahj(2,j), ahj(3,j))
end do
!$omp end parallel do
end if

do concurrent(i = 1:npl)
acc(:,i) = acc(:,i) + ahi(:,i) + ahj(:,i)
Expand All @@ -130,12 +157,12 @@ module subroutine kick_getacch_int_all_triangular_pl(npl, nplm, x, Gmass, radius
!! Adapted from Hal Levison's Swift routine getacch_ah3.f
!! Adapted from David E. Kaufmann's Swifter routine whm_kick_getacch_ah3.f90 and helio_kick_getacch_int.f9
implicit none
integer(I4B), intent(in) :: npl !! Total number of massive bodies
integer(I4B), intent(in) :: nplm !! Number of fully interacting massive bodies
real(DP), dimension(:,:), intent(in) :: x !! Position vector array
real(DP), dimension(:), intent(in) :: Gmass !! Array of massive body G*mass
real(DP), dimension(:), intent(in) :: radius !! Array of massive body radii
real(DP), dimension(:,:), intent(inout) :: acc !! Acceleration vector array
integer(I4B), intent(in) :: npl !! Total number of massive bodies
integer(I4B), intent(in) :: nplm !! Number of fully interacting massive bodies
real(DP), dimension(:,:), intent(in) :: x !! Position vector array
real(DP), dimension(:), intent(in) :: Gmass !! Array of massive body G*mass
real(DP), dimension(:), intent(in), optional :: radius !! Array of massive body radii
real(DP), dimension(:,:), intent(inout) :: acc !! Acceleration vector array
! Internals
real(DP), dimension(NDIM,npl) :: ahi, ahj
integer(I4B) :: i, j
Expand All @@ -145,23 +172,42 @@ module subroutine kick_getacch_int_all_triangular_pl(npl, nplm, x, Gmass, radius
ahi(:,:) = 0.0_DP
ahj(:,:) = 0.0_DP

!$omp parallel do default(private) schedule(static)&
!$omp shared(npl, nplm, x, Gmass, radius) &
!$omp lastprivate(rji2, rlim2, xr, yr, zr) &
!$omp reduction(+:ahi) &
!$omp reduction(-:ahj)
do i = 1, nplm
do concurrent(j = i+1:npl)
xr = x(1, j) - x(1, i)
yr = x(2, j) - x(2, i)
zr = x(3, j) - x(3, i)
rji2 = xr**2 + yr**2 + zr**2
rlim2 = (radius(i) + radius(j))**2
if (rji2 > rlim2) call kick_getacch_int_one_pl(rji2, xr, yr, zr, Gmass(i), Gmass(j), &
ahi(1,i), ahi(2,i), ahi(3,i), ahj(1,j), ahj(2,j), ahj(3,j))
if (present(radius)) then
!$omp parallel do default(private) schedule(static)&
!$omp shared(npl, nplm, x, Gmass, radius) &
!$omp lastprivate(rji2, rlim2, xr, yr, zr) &
!$omp reduction(+:ahi) &
!$omp reduction(-:ahj)
do i = 1, nplm
do concurrent(j = i+1:npl)
xr = x(1, j) - x(1, i)
yr = x(2, j) - x(2, i)
zr = x(3, j) - x(3, i)
rji2 = xr**2 + yr**2 + zr**2
rlim2 = (radius(i) + radius(j))**2
if (rji2 > rlim2) call kick_getacch_int_one_pl(rji2, xr, yr, zr, Gmass(i), Gmass(j), &
ahi(1,i), ahi(2,i), ahi(3,i), ahj(1,j), ahj(2,j), ahj(3,j))
end do
end do
end do
!$omp end parallel do
!$omp end parallel do
else
!$omp parallel do default(private) schedule(static)&
!$omp shared(npl, nplm, x, Gmass, radius) &
!$omp lastprivate(rji2, xr, yr, zr) &
!$omp reduction(+:ahi) &
!$omp reduction(-:ahj)
do i = 1, nplm
do concurrent(j = i+1:npl)
xr = x(1, j) - x(1, i)
yr = x(2, j) - x(2, i)
zr = x(3, j) - x(3, i)
rji2 = xr**2 + yr**2 + zr**2
call kick_getacch_int_one_pl(rji2, xr, yr, zr, Gmass(i), Gmass(j), &
ahi(1,i), ahi(2,i), ahi(3,i), ahj(1,j), ahj(2,j), ahj(3,j))
end do
end do
!$omp end parallel do
end if

do concurrent(i = 1:npl)
acc(:,i) = acc(:,i) + ahi(:,i) + ahj(:,i)
Expand Down
26 changes: 13 additions & 13 deletions src/modules/swiftest_classes.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -906,23 +906,23 @@ end subroutine kick_getacch_int_tp

module subroutine kick_getacch_int_all_flat_pl(npl, nplpl, k_plpl, x, Gmass, radius, acc)
implicit none
integer(I4B), intent(in) :: npl !! Number of massive bodies
integer(I8B), intent(in) :: nplpl !! Number of massive body interactions to compute
integer(I4B), dimension(:,:), intent(in) :: k_plpl !! Array of interaction pair indices (flattened upper triangular matrix)
real(DP), dimension(:,:), intent(in) :: x !! Position vector array
real(DP), dimension(:), intent(in) :: Gmass !! Array of massive body G*mass
real(DP), dimension(:), intent(in) :: radius !! Array of massive body radii
real(DP), dimension(:,:), intent(inout) :: acc !! Acceleration vector array
integer(I4B), intent(in) :: npl !! Number of massive bodies
integer(I8B), intent(in) :: nplpl !! Number of massive body interactions to compute
integer(I4B), dimension(:,:), intent(in) :: k_plpl !! Array of interaction pair indices (flattened upper triangular matrix)
real(DP), dimension(:,:), intent(in) :: x !! Position vector array
real(DP), dimension(:), intent(in) :: Gmass !! Array of massive body G*mass
real(DP), dimension(:), intent(in), optional :: radius !! Array of massive body radii
real(DP), dimension(:,:), intent(inout) :: acc !! Acceleration vector array
end subroutine kick_getacch_int_all_flat_pl

module subroutine kick_getacch_int_all_triangular_pl(npl, nplm, x, Gmass, radius, acc)
implicit none
integer(I4B), intent(in) :: npl !! Total number of massive bodies
integer(I4B), intent(in) :: nplm !! Number of fully interacting massive bodies
real(DP), dimension(:,:), intent(in) :: x !! Position vector array
real(DP), dimension(:), intent(in) :: Gmass !! Array of massive body G*mass
real(DP), dimension(:), intent(in) :: radius !! Array of massive body radii
real(DP), dimension(:,:), intent(inout) :: acc !! Acceleration vector array
integer(I4B), intent(in) :: npl !! Total number of massive bodies
integer(I4B), intent(in) :: nplm !! Number of fully interacting massive bodies
real(DP), dimension(:,:), intent(in) :: x !! Position vector array
real(DP), dimension(:), intent(in) :: Gmass !! Array of massive body G*mass
real(DP), dimension(:), intent(in), optional :: radius !! Array of massive body radii
real(DP), dimension(:,:), intent(inout) :: acc !! Acceleration vector array
end subroutine kick_getacch_int_all_triangular_pl

module subroutine kick_getacch_int_all_tp(ntp, npl, xtp, xpl, GMpl, lmask, acc)
Expand Down
13 changes: 6 additions & 7 deletions src/netcdf/netcdf.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -567,6 +567,9 @@ module function netcdf_read_frame_system(self, iu, param) result(ierr)
call check( nf90_get_var(iu%ncid, iu%radius_varid, rtemp, start=[1, tslot]) )
cb%radius = rtemp(1)
pl%radius(:) = pack(rtemp, plmask)
else
cb%radius = param%rmin
pl%radius(:) = 0.0_DP
end if

if (param%lrotation) then
Expand Down Expand Up @@ -896,12 +899,8 @@ module subroutine netcdf_write_frame_base(self, iu, param)
select type(self)
class is (swiftest_pl) ! Additional output if the passed polymorphic object is a massive body
call check( nf90_put_var(iu%ncid, iu%Gmass_varid, self%Gmass(j), start=[idslot, tslot]) )
if (param%lrhill_present) then
call check( nf90_put_var(iu%ncid, iu%rhill_varid, self%rhill(j), start=[idslot, tslot]) )
end if
if (param%lclose) then
call check( nf90_put_var(iu%ncid, iu%radius_varid, self%radius(j), start=[idslot, tslot]) )
end if
if (param%lrhill_present) call check( nf90_put_var(iu%ncid, iu%rhill_varid, self%rhill(j), start=[idslot, tslot]) )
if (param%lclose) call check( nf90_put_var(iu%ncid, iu%radius_varid, self%radius(j), start=[idslot, tslot]) )
if (param%lrotation) then
call check( nf90_put_var(iu%ncid, iu%Ip1_varid, self%Ip(1, j), start=[idslot, tslot]) )
call check( nf90_put_var(iu%ncid, iu%Ip2_varid, self%Ip(2, j), start=[idslot, tslot]) )
Expand All @@ -923,7 +922,7 @@ module subroutine netcdf_write_frame_base(self, iu, param)
call check( nf90_put_var(iu%ncid, iu%id_varid, self%id, start=[idslot]) )

call check( nf90_put_var(iu%ncid, iu%Gmass_varid, self%Gmass, start=[idslot, tslot]) )
call check( nf90_put_var(iu%ncid, iu%radius_varid, self%radius, start=[idslot, tslot]) )
if (param%lclose) call check( nf90_put_var(iu%ncid, iu%radius_varid, self%radius, start=[idslot, tslot]) )
if (param%lrotation) then
call check( nf90_put_var(iu%ncid, iu%Ip1_varid, self%Ip(1), start=[idslot, tslot]) )
call check( nf90_put_var(iu%ncid, iu%Ip2_varid, self%Ip(2), start=[idslot, tslot]) )
Expand Down
4 changes: 2 additions & 2 deletions src/util/util_dealloc.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -22,8 +22,8 @@ module subroutine util_dealloc_body(self)
if (allocated(self%vb)) deallocate(self%vb)
if (allocated(self%ah)) deallocate(self%ah)
if (allocated(self%aobl)) deallocate(self%aobl)
if (allocated(self%agr)) deallocate(self%lmask)
if (allocated(self%atide)) deallocate(self%lmask)
if (allocated(self%agr)) deallocate(self%agr)
if (allocated(self%atide)) deallocate(self%atide)
if (allocated(self%ir3h)) deallocate(self%ir3h)
if (allocated(self%a)) deallocate(self%a)
if (allocated(self%e)) deallocate(self%e)
Expand Down

0 comments on commit bdb9f0c

Please sign in to comment.