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

Commit

Permalink
Merge branch 'debug'
Browse files Browse the repository at this point in the history
  • Loading branch information
@daminton
daminton committed Aug 23, 2021
2 parents ca3d253 + e43ae9b commit 58c2051
Show file tree
Hide file tree
Showing 6 changed files with 139 additions and 133 deletions.
17 changes: 9 additions & 8 deletions src/io/io.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -189,14 +189,14 @@ module subroutine io_dump_system(self, param)
dump_param%incbfile = trim(adjustl(DUMP_CB_FILE(idx)))
dump_param%inplfile = trim(adjustl(DUMP_PL_FILE(idx)))
dump_param%intpfile = trim(adjustl(DUMP_TP_FILE(idx)))
dump_param%out_form = XV
dump_param%in_form = XV
dump_param%out_stat = 'APPEND'
dump_param%in_type = REAL8_TYPE
dump_param%T0 = param%t

call dump_param%dump(param_file_name)

dump_param%out_form = XV
call self%cb%dump(dump_param)
call self%pl%dump(dump_param)
call self%tp%dump(dump_param)
Expand Down Expand Up @@ -454,7 +454,7 @@ module subroutine io_param_reader(self, unit, iotype, v_list, iostat, iomsg)
call io_toupper(param_value)
param%out_stat = param_value
case ("ISTEP_DUMP")
read(param_value, *) param%istep_dump
read(param_value, *, err = 667, iomsg = iomsg) param%istep_dump
case ("CHK_CLOSE")
call io_toupper(param_value)
if (param_value == "YES" .or. param_value == 'T') param%lclose = .true.
Expand Down Expand Up @@ -551,6 +551,8 @@ module subroutine io_param_reader(self, unit, iotype, v_list, iostat, iomsg)
read(param_value, *, err = 667, iomsg = iomsg) param%Ecollisions
case("EUNTRACKED")
read(param_value, *, err = 667, iomsg = iomsg) param%Euntracked
case ("MAXID")
read(param_value, *, err = 667, iomsg = iomsg) param%maxid
case ("NPLMAX", "NTPMAX", "GMTINY", "MIN_GMFRAG", "PARTICLE_OUT", "FRAGMENTATION", "SEED", "YARKOVSKY", "YORP") ! Ignore SyMBA-specific, not-yet-implemented, or obsolete input parameters
case default
write(iomsg,*) "Unknown parameter -> ",param_name
Expand Down Expand Up @@ -760,6 +762,7 @@ module subroutine io_param_writer(self, unit, iotype, v_list, iostat, iomsg)
write(param_name, Afmt) "TP_IN"; write(param_value, Afmt) trim(adjustl(param%intpfile)); write(unit, Afmt, err = 667, iomsg = iomsg) adjustl(param_name), adjustl(param_value)
write(param_name, Afmt) "IN_TYPE"; write(param_value, Afmt) trim(adjustl(param%in_type)); write(unit, Afmt, err = 667, iomsg = iomsg) adjustl(param_name), adjustl(param_value)
write(param_name, Afmt) "IN_FORM"; write(param_value, Afmt) trim(adjustl(param%in_form)); write(unit, Afmt, err = 667, iomsg = iomsg) adjustl(param_name), adjustl(param_value)
if (param%istep_dump > 0) write(param_name, Afmt) "ISTEP_DUMP"; write(param_value, Ifmt) param%istep_dump; write(unit, Afmt, err = 667, iomsg = iomsg) adjustl(param_name), adjustl(param_value)
if (param%istep_out > 0) then
write(param_name, Afmt) "ISTEP_OUT"; write(param_value, Ifmt) param%istep_out; write(unit, Afmt, err = 667, iomsg = iomsg) adjustl(param_name), adjustl(param_value)
write(param_name, Afmt) "BIN_OUT"; write(param_value, Afmt) trim(adjustl(param%outfile)); write(unit, Afmt, err = 667, iomsg = iomsg) adjustl(param_name), adjustl(param_value)
Expand All @@ -768,9 +771,6 @@ module subroutine io_param_writer(self, unit, iotype, v_list, iostat, iomsg)
write(param_name, Afmt) "OUT_STAT"; write(param_value, Afmt) "APPEND"; write(unit, Afmt, err = 667, iomsg = iomsg) adjustl(param_name), adjustl(param_value)
end if
write(param_name, Afmt) "ENC_OUT"; write(param_value, Afmt) trim(adjustl(param%enc_out)); write(unit, Afmt, err = 667, iomsg = iomsg) adjustl(param_name), adjustl(param_value)
if (param%istep_dump > 0) then
write(param_name, Afmt) "ISTEP_DUMP"; write(param_value, Ifmt) param%istep_dump; write(unit, Afmt, err = 667, iomsg = iomsg) adjustl(param_name), adjustl(param_value)
end if
write(param_name, Afmt) "CHK_RMIN"; write(param_value, Rfmt) param%rmin; write(unit, Afmt, err = 667, iomsg = iomsg) adjustl(param_name), adjustl(param_value)
write(param_name, Afmt) "CHK_RMAX"; write(param_value, Rfmt) param%rmax; write(unit, Afmt, err = 667, iomsg = iomsg) adjustl(param_name), adjustl(param_value)
write(param_name, Afmt) "CHK_EJECT"; write(param_value, Rfmt) param%rmaxu; write(unit, Afmt, err = 667, iomsg = iomsg) adjustl(param_name), adjustl(param_value)
Expand Down Expand Up @@ -810,6 +810,7 @@ module subroutine io_param_writer(self, unit, iotype, v_list, iostat, iomsg)
write(param_name, Afmt) "EUNTRACKED"; write(param_value, Rfmt) param%Euntracked; write(unit, Afmt, err = 667, iomsg = iomsg) adjustl(param_name), adjustl(param_value)
end if
write(param_name, Afmt) "FIRSTKICK"; write(param_value, Lfmt) param%lfirstkick; write(unit, Afmt, err = 667, iomsg = iomsg) adjustl(param_name), adjustl(param_value)
write(param_name, Afmt) "MAXID"; write(param_value, Ifmt) param%maxid; write(unit, Afmt, err = 667, iomsg = iomsg) adjustl(param_name), adjustl(param_value)

iostat = 0
iomsg = "UDIO not implemented"
Expand Down Expand Up @@ -1053,7 +1054,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(:)
if (param%lclose) read(iu, err = 667, iomsg = errmsg) pl%radius(:)
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 All @@ -1080,7 +1081,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)
if (param%lclose) read(iu, *, err = 667, iomsg = errmsg) self%radius(i)
read(iu, *, err = 667, iomsg = errmsg) self%radius(i)
class is (swiftest_tp)
read(iu, *, err = 667, iomsg = errmsg) self%id(i)
end select
Expand Down Expand Up @@ -1507,7 +1508,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)
if (param%lclose) write(iu, err = 667, iomsg = errmsg) pl%radius(1:n)
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
148 changes: 87 additions & 61 deletions src/kick/kick.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -12,37 +12,8 @@ module subroutine kick_getacch_int_pl(self)
implicit none
! Arguments
class(swiftest_pl), intent(inout) :: self
! Internals
integer(I8B) :: k, nplpl
real(DP) :: rji2
real(DP) :: dx, dy, dz
integer(I4B) :: i, j
real(DP), dimension(:,:), pointer :: ah, xh
real(DP), dimension(NDIM,self%nbody) :: ahi, ahj
integer(I4B), dimension(:,:), pointer :: k_plpl
logical, dimension(:), pointer :: lmask
real(DP), dimension(:), pointer :: Gmass

associate(ah => self%ah, xh => self%xh, k_plpl => self%k_plpl, lmask => self%lmask, Gmass => self%Gmass)
nplpl = self%nplpl
ahi(:,:) = 0.0_DP
ahj(:,:) = 0.0_DP
!$omp parallel do default(shared)&
!$omp private(k, i, j, dx, dy, dz, rji2) &
!$omp reduction(+:ahi) &
!$omp reduction(-:ahj)
do k = 1_I8B, nplpl
i = k_plpl(1,k)
j = k_plpl(2,k)
dx = xh(1, j) - xh(1, i)
dy = xh(2, j) - xh(2, i)
dz = xh(3, j) - xh(3, i)
rji2 = dx**2 + dy**2 + dz**2
if (lmask(i) .and. lmask(j)) call kick_getacch_int_one_pl(rji2, dx, dy, dz, 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
ah(:,1:self%nbody) = ah(:,1:self%nbody) + ahi(:,1:self%nbody) + ahj(:,1:self%nbody)
end associate

call kick_getacch_int_all_pl(self%nbody, self%nplpl, self%k_plpl, self%xh, self%Gmass, self%radius, self%ah)

return
end subroutine kick_getacch_int_pl
Expand All @@ -63,7 +34,6 @@ module pure subroutine kick_getacch_int_tp(self, GMpl, xhp, npl)
integer(I4B), intent(in) :: npl !! Number of active massive bodies
! Internals
integer(I4B) :: i, j
!real(DP), dimension(:,:), allocatable :: aht

if ((self%nbody == 0) .or. (npl == 0)) return

Expand All @@ -73,60 +43,116 @@ module pure subroutine kick_getacch_int_tp(self, GMpl, xhp, npl)
do j = 1, npl
block
real(DP) :: rji2
real(DP) :: dx, dy, dz
dx = tp%xh(1, i) - xhp(1, j)
dy = tp%xh(2, i) - xhp(1, j)
dz = tp%xh(3, i) - xhp(1, j)
rji2 = dx**2 + dy**2 + dz**2
call kick_getacch_int_one_tp(rji2, dx, dy, dz, GMpl(i), tp%ah(1,i), tp%ah(2,i), tp%ah(3,i))
real(DP) :: xr, yr, zr
xr = tp%xh(1, i) - xhp(1, j)
yr = tp%xh(2, i) - xhp(1, j)
zr = tp%xh(3, i) - xhp(1, j)
rji2 = xr**2 + yr**2 + zr**2
call kick_getacch_int_one_tp(rji2, xr, yr, zr, GMpl(i), tp%ah(1,i), tp%ah(2,i), tp%ah(3,i))
end block
end do
end if
end do
!call move_alloc(aht, tp%ah)
end associate

return
end subroutine kick_getacch_int_tp

module pure elemental subroutine kick_getacch_int_one_pl(rji2, dx, dy, dz, Gmi, Gmj, axi, ayi, azi, axj, ayj, azj)

module subroutine kick_getacch_int_all_pl(npl, nplpl, k_plpl, x, Gmass, radius, acc)
!! author: David A. Minton
!!
!! Compute direct cross (third) term heliocentric accelerations for massive bodies, with parallelization
!!
!! 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
real(DP), intent(in) :: rji2
real(DP), intent(in) :: dx, dy, dz
real(DP), intent(in) :: Gmi
real(DP), intent(in) :: Gmj
real(DP), intent(inout) :: axi, ayi, azi
real(DP), intent(inout) :: axj, ayj, azj
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
! Internals
integer(I8B) :: k
real(DP) :: rji2, rlim2
real(DP) :: xr, yr, zr
integer(I4B) :: i, j
real(DP), dimension(NDIM,npl) :: ahi, ahj

ahi(:,:) = 0.0_DP
ahj(:,:) = 0.0_DP
!$omp parallel do default(private)&
!$omp shared(nplpl, k_plpl, x, Gmass, radius) &
!$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
acc(:,1:npl) = acc(:,1:npl) + ahi(:,1:npl) + ahj(:,1:npl)
return
end subroutine kick_getacch_int_all_pl


module pure subroutine kick_getacch_int_one_pl(rji2, xr, yr, zr, Gmi, Gmj, axi, ayi, azi, axj, ayj, azj)
!! author: David A. Minton
!!
!! Compute direct cross (third) term heliocentric accelerations for a single pair of massive bodies
!!
!! 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
real(DP), intent(in) :: rji2 !! Square of distance between the two bodies
real(DP), intent(in) :: xr, yr, zr !! Distances between the two bodies in x, y, and z directions
real(DP), intent(in) :: Gmi !! G*mass of body i
real(DP), intent(in) :: Gmj !! G*mass of body j
real(DP), intent(inout) :: axi, ayi, azi !! Acceleration vector components of body i
real(DP), intent(inout) :: axj, ayj, azj !! Acceleration vector components of body j
! Internals
real(DP) :: faci, facj, irij3

irij3 = 1.0_DP / (rji2 * sqrt(rji2))
faci = Gmi * irij3
facj = Gmj * irij3
axi = axi + facj * dx
ayi = ayi + facj * dy
azi = azi + facj * dz
axj = axj - faci * dx
ayj = ayj - faci * dy
azj = azj - faci * dz
axi = axi + facj * xr
ayi = ayi + facj * yr
azi = azi + facj * zr
axj = axj - faci * xr
ayj = ayj - faci * yr
azj = azj - faci * zr

return
end subroutine kick_getacch_int_one_pl


!module pure elemental subroutine kick_getacch_int_one_tp(rji2, dx, dy, dz, GMpl, ax, ay, az)
module pure subroutine kick_getacch_int_one_tp(rji2, dx, dy, dz, GMpl, ax, ay, az)
module pure subroutine kick_getacch_int_one_tp(rji2, xr, yr, zr, GMpl, ax, ay, az)
!! author: David A. Minton
!!
!! Compute direct cross (third) term heliocentric accelerations of a single test particle massive body pair.
!!
!! Adapted from Hal Levison's Swift routine getacch_ah3_tp.f
!! Adapted from David E. Kaufmann's Swifter routine whm_kick_getacch_ah3.f90 and helio_kick_getacch_int_tp.f90
implicit none
real(DP), intent(in) :: rji2
real(DP), intent(in) :: dx, dy, dz
real(DP), intent(in) :: GMpl
real(DP), intent(inout) :: ax, ay, az
real(DP), intent(in) :: rji2 !! Square of distance between the test particle and massive body
real(DP), intent(in) :: xr, yr, zr !! Distances between the two bodies in x, y, and z directions
real(DP), intent(in) :: Gmpl !! G*mass of massive body
real(DP), intent(inout) :: ax, ay, az !! Acceleration vector components of test particle
! Internals
real(DP) :: fac

fac = GMpl / (rji2 * sqrt(rji2))
ax = ax - fac * dx
ay = ay - fac * dy
az = az - fac * dz
ax = ax - fac * xr
ay = ay - fac * yr
az = az - fac * zr
return
end subroutine kick_getacch_int_one_tp

Expand Down
Loading

0 comments on commit 58c2051

Please sign in to comment.