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 
Got rid of intermediate npl argument as it's part of the symba_plA structure.
  • Loading branch information
@daminton
daminton committed May 7, 2021
1 parent 9918a0d commit 49a0b33
Show file tree
Hide file tree
Showing 4 changed files with 121 additions and 122 deletions.
2 changes: 1 addition & 1 deletion src/main/swiftest_symba.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -156,7 +156,7 @@ program swiftest_symba
call symba_discard_pl(t, npl, ntp, symba_plA, symba_tpA, rmin, rmax, rmaxu, qmin, qmin_coord, qmin_alo, qmin_ahi, ldiscard_pl)
if (ldiscard_pl) then
if (param%lenergy) call symba_energy_eucl(npl, symba_plA, j2rp2, j4rp4, ke_orbit_before, ke_spin_before, pe_before, Eorbit_before, Ltot)
call symba_rearray_pl(t, npl, symba_plA, nmergeadd, mergeadd_list, discard_plA, param)
call symba_rearray_pl(t, symba_plA, nmergeadd, mergeadd_list, discard_plA, param)
if (param%lenergy) then
call symba_energy_eucl(npl, symba_plA, j2rp2, j4rp4, ke_orbit_after, ke_spin_after, pe_after, Eorbit_after, Ltot)
Ecollision = Eorbit_before - Eorbit_after ! Energy change resulting in this collisional event Total running energy offset from collision in this step
Expand Down
3 changes: 1 addition & 2 deletions src/modules/module_interfaces.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1075,13 +1075,12 @@ END SUBROUTINE symba_peri
END INTERFACE

INTERFACE
subroutine symba_rearray_pl(t, npl, symba_plA, nmergeadd, mergeadd_list, discard_plA, param)
subroutine symba_rearray_pl(t, symba_plA, nmergeadd, mergeadd_list, discard_plA, param)
use swiftest_globals
use swiftest_data_structures
use module_symba
implicit none
real(DP), intent(in) :: t
integer(I4B), intent(inout) :: npl
integer(I4B), intent(in) :: nmergeadd
type(symba_pl), intent(inout) :: symba_plA
type(symba_pl), intent(inout) :: discard_plA
Expand Down
2 changes: 1 addition & 1 deletion src/symba/symba_collision.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -304,7 +304,7 @@ subroutine symba_collision(t, symba_plA, nplplenc, plplenc_list, lfrag_add, nmer
end do
deallocate(family)

call symba_rearray_pl(t, npl, symba_plA, nmergeadd, mergeadd_list, discard_plA, param)
call symba_rearray_pl(t, symba_plA, nmergeadd, mergeadd_list, discard_plA, param)

end do
end associate
Expand Down
236 changes: 118 additions & 118 deletions src/symba/symba_rearray.f90
View file
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
subroutine symba_rearray_pl(t, npl, symba_plA, nmergeadd, mergeadd_list, discard_plA, param)
subroutine symba_rearray_pl(t, symba_plA, nmergeadd, mergeadd_list, discard_plA, param)
!! Author: the Purdue Swiftest Team - David A. Minton, Carlisle A. Wishard, Jennifer L.L. Pouplin, and Jacob R. Elliott
!!
!! Clean up tp and pl arrays to remove discarded bodies and add new bodies
Expand All @@ -12,7 +12,6 @@ subroutine symba_rearray_pl(t, npl, symba_plA, nmergeadd, mergeadd_list, discard

! Arguments
real(DP), intent(in) :: t
integer(I4B), intent(inout) :: npl
integer(I4B), intent(in) :: nmergeadd
type(symba_pl), intent(inout) :: symba_plA
type(symba_pl), intent(inout) :: discard_plA
Expand All @@ -24,132 +23,133 @@ subroutine symba_rearray_pl(t, npl, symba_plA, nmergeadd, mergeadd_list, discard
logical, dimension(:), allocatable :: discard_l_pl, add_l_pl
logical :: lescape, ldiscard_pl

ldiscard_pl = any(symba_plA%helio%swiftest%status(1:npl) /= ACTIVE)
if (ldiscard_pl) then

! Deal with the central body/system discards if there are any
do i = 1, npl
if ((symba_plA%helio%swiftest%status(i) == DISCARDED_RMIN) .or. (symba_plA%helio%swiftest%status(i) == DISCARDED_PERI)) then
lescape = .false.
else if ((symba_plA%helio%swiftest%status(i) == DISCARDED_RMAX) .or. (symba_plA%helio%swiftest%status(i) == DISCARDED_RMAXU)) then
lescape = .true.
else
cycle
associate(npl => symba_plA%helio%swiftest%nbody)
ldiscard_pl = any(symba_plA%helio%swiftest%status(1:npl) /= ACTIVE)
if (ldiscard_pl) then

! Deal with the central body/system discards if there are any
do i = 1, npl
if ((symba_plA%helio%swiftest%status(i) == DISCARDED_RMIN) .or. (symba_plA%helio%swiftest%status(i) == DISCARDED_PERI)) then
lescape = .false.
else if ((symba_plA%helio%swiftest%status(i) == DISCARDED_RMAX) .or. (symba_plA%helio%swiftest%status(i) == DISCARDED_RMAXU)) then
lescape = .true.
else
cycle
end if
call symba_discard_conserve_mtm(symba_plA%helio%swiftest, i, lescape)
end do

allocate(discard_l_pl(npl))
nkpl = count(symba_plA%helio%swiftest%status(:) == ACTIVE)
nsppl = npl - nkpl
dlo = 1
if (allocated(discard_plA%helio%swiftest%id)) then ! We alredy made a discard list in this step, so we need to append to it
nsppl = nsppl + discard_plA%helio%swiftest%nbody
dlo = dlo + discard_plA%helio%swiftest%nbody
call util_resize_pl(discard_plA, nsppl)
else
call symba_pl_allocate(discard_plA, nsppl)
end if
call symba_discard_conserve_mtm(symba_plA%helio%swiftest, i, lescape)
end do

allocate(discard_l_pl(npl))
nkpl = count(symba_plA%helio%swiftest%status(:) == ACTIVE)
nsppl = npl - nkpl
dlo = 1
if (allocated(discard_plA%helio%swiftest%id)) then ! We alredy made a discard list in this step, so we need to append to it
nsppl = nsppl + discard_plA%helio%swiftest%nbody
dlo = dlo + discard_plA%helio%swiftest%nbody
call util_resize_pl(discard_plA, nsppl)
else
call symba_pl_allocate(discard_plA, nsppl)
end if

discard_l_pl(:) = (symba_plA%helio%swiftest%status(1:npl) /= ACTIVE)

! Spill discarded bodies into discard list
discard_plA%helio%swiftest%id(dlo:nsppl) = pack(symba_plA%helio%swiftest%id(1:npl), discard_l_pl)
discard_plA%helio%swiftest%status(dlo:nsppl) = pack(symba_plA%helio%swiftest%status(1:npl), discard_l_pl)
discard_plA%helio%swiftest%mass(dlo:nsppl) = pack(symba_plA%helio%swiftest%mass(1:npl), discard_l_pl)
discard_plA%helio%swiftest%radius(dlo:nsppl) = pack(symba_plA%helio%swiftest%radius(1:npl), discard_l_pl)
discard_plA%helio%swiftest%rhill(dlo:nsppl) = pack(symba_plA%helio%swiftest%rhill(1:npl), discard_l_pl)
do i = 1, NDIM
discard_plA%helio%swiftest%xh(i,dlo:nsppl) = pack(symba_plA%helio%swiftest%xh(i,1:npl), discard_l_pl)
discard_plA%helio%swiftest%vh(i,dlo:nsppl) = pack(symba_plA%helio%swiftest%vh(i,1:npl), discard_l_pl)
discard_plA%helio%swiftest%xb(i,dlo:nsppl) = pack(symba_plA%helio%swiftest%xb(i,1:npl), discard_l_pl)
discard_plA%helio%swiftest%vb(i,dlo:nsppl) = pack(symba_plA%helio%swiftest%vb(i,1:npl), discard_l_pl)
discard_plA%helio%swiftest%Ip(i,dlo:nsppl) = pack(symba_plA%helio%swiftest%Ip(i,1:npl), discard_l_pl)
discard_plA%helio%swiftest%rot(i,dlo:nsppl) = pack(symba_plA%helio%swiftest%rot(i,1:npl), discard_l_pl)
end do

if (nkpl > 0) then
! Pack kept bodies down
symba_plA%helio%swiftest%id(1:nkpl) = pack(symba_plA%helio%swiftest%id(1:npl), .not. discard_l_pl)
symba_plA%helio%swiftest%status(1:nkpl) = pack(symba_plA%helio%swiftest%status(1:npl), .not. discard_l_pl)
symba_plA%helio%swiftest%mass(1:nkpl) = pack(symba_plA%helio%swiftest%mass(1:npl), .not. discard_l_pl)
symba_plA%helio%swiftest%radius(1:nkpl) = pack(symba_plA%helio%swiftest%radius(1:npl), .not. discard_l_pl)
symba_plA%helio%swiftest%rhill(1:nkpl) = pack(symba_plA%helio%swiftest%rhill(1:npl), .not. discard_l_pl)
discard_l_pl(:) = (symba_plA%helio%swiftest%status(1:npl) /= ACTIVE)

! Spill discarded bodies into discard list
discard_plA%helio%swiftest%id(dlo:nsppl) = pack(symba_plA%helio%swiftest%id(1:npl), discard_l_pl)
discard_plA%helio%swiftest%status(dlo:nsppl) = pack(symba_plA%helio%swiftest%status(1:npl), discard_l_pl)
discard_plA%helio%swiftest%mass(dlo:nsppl) = pack(symba_plA%helio%swiftest%mass(1:npl), discard_l_pl)
discard_plA%helio%swiftest%radius(dlo:nsppl) = pack(symba_plA%helio%swiftest%radius(1:npl), discard_l_pl)
discard_plA%helio%swiftest%rhill(dlo:nsppl) = pack(symba_plA%helio%swiftest%rhill(1:npl), discard_l_pl)
do i = 1, NDIM
symba_plA%helio%swiftest%xh(i,1:nkpl) = pack(symba_plA%helio%swiftest%xh(i,1:npl), .not. discard_l_pl)
symba_plA%helio%swiftest%vh(i,1:nkpl) = pack(symba_plA%helio%swiftest%vh(i,1:npl), .not. discard_l_pl)
symba_plA%helio%swiftest%xb(i,1:nkpl) = pack(symba_plA%helio%swiftest%xb(i,1:npl), .not. discard_l_pl)
symba_plA%helio%swiftest%vb(i,1:nkpl) = pack(symba_plA%helio%swiftest%vb(i,1:npl), .not. discard_l_pl)
symba_plA%helio%swiftest%Ip(i,1:nkpl) = pack(symba_plA%helio%swiftest%Ip(i,1:npl), .not. discard_l_pl)
symba_plA%helio%swiftest%rot(i,1:nkpl) = pack(symba_plA%helio%swiftest%rot(i,1:npl), .not. discard_l_pl)
discard_plA%helio%swiftest%xh(i,dlo:nsppl) = pack(symba_plA%helio%swiftest%xh(i,1:npl), discard_l_pl)
discard_plA%helio%swiftest%vh(i,dlo:nsppl) = pack(symba_plA%helio%swiftest%vh(i,1:npl), discard_l_pl)
discard_plA%helio%swiftest%xb(i,dlo:nsppl) = pack(symba_plA%helio%swiftest%xb(i,1:npl), discard_l_pl)
discard_plA%helio%swiftest%vb(i,dlo:nsppl) = pack(symba_plA%helio%swiftest%vb(i,1:npl), discard_l_pl)
discard_plA%helio%swiftest%Ip(i,dlo:nsppl) = pack(symba_plA%helio%swiftest%Ip(i,1:npl), discard_l_pl)
discard_plA%helio%swiftest%rot(i,dlo:nsppl) = pack(symba_plA%helio%swiftest%rot(i,1:npl), discard_l_pl)
end do

if (nkpl > 0) then
! Pack kept bodies down
symba_plA%helio%swiftest%id(1:nkpl) = pack(symba_plA%helio%swiftest%id(1:npl), .not. discard_l_pl)
symba_plA%helio%swiftest%status(1:nkpl) = pack(symba_plA%helio%swiftest%status(1:npl), .not. discard_l_pl)
symba_plA%helio%swiftest%mass(1:nkpl) = pack(symba_plA%helio%swiftest%mass(1:npl), .not. discard_l_pl)
symba_plA%helio%swiftest%radius(1:nkpl) = pack(symba_plA%helio%swiftest%radius(1:npl), .not. discard_l_pl)
symba_plA%helio%swiftest%rhill(1:nkpl) = pack(symba_plA%helio%swiftest%rhill(1:npl), .not. discard_l_pl)
do i = 1, NDIM
symba_plA%helio%swiftest%xh(i,1:nkpl) = pack(symba_plA%helio%swiftest%xh(i,1:npl), .not. discard_l_pl)
symba_plA%helio%swiftest%vh(i,1:nkpl) = pack(symba_plA%helio%swiftest%vh(i,1:npl), .not. discard_l_pl)
symba_plA%helio%swiftest%xb(i,1:nkpl) = pack(symba_plA%helio%swiftest%xb(i,1:npl), .not. discard_l_pl)
symba_plA%helio%swiftest%vb(i,1:nkpl) = pack(symba_plA%helio%swiftest%vb(i,1:npl), .not. discard_l_pl)
symba_plA%helio%swiftest%Ip(i,1:nkpl) = pack(symba_plA%helio%swiftest%Ip(i,1:npl), .not. discard_l_pl)
symba_plA%helio%swiftest%rot(i,1:nkpl) = pack(symba_plA%helio%swiftest%rot(i,1:npl), .not. discard_l_pl)
end do
end if

npl = nkpl
else
nkpl = npl
end if

npl = nkpl
else
nkpl = npl
end if

if (nmergeadd > 0) then
call util_resize_pl(symba_plA, nkpl+nmergeadd)
npl = nkpl + nmergeadd
symba_plA%helio%swiftest%nbody = npl

!add merge products to the end of the planet list
symba_plA%helio%swiftest%status(nkpl+1:npl) = mergeadd_list%status(1:nmergeadd)
symba_plA%helio%swiftest%id(nkpl+1:npl) = mergeadd_list%id(1:nmergeadd)
symba_plA%helio%swiftest%mass(nkpl+1:npl) = mergeadd_list%mass(1:nmergeadd)
symba_plA%helio%swiftest%radius(nkpl+1:npl) = mergeadd_list%radius(1:nmergeadd)
do i = 1, NDIM
symba_plA%helio%swiftest%xb(i,nkpl+1:npl) = mergeadd_list%xb(i,1:nmergeadd)
symba_plA%helio%swiftest%vb(i,nkpl+1:npl) = mergeadd_list%vb(i,1:nmergeadd)
symba_plA%helio%swiftest%Ip(i,nkpl+1:npl) = mergeadd_list%Ip(i,1:nmergeadd)
symba_plA%helio%swiftest%rot(i,nkpl+1:npl) = mergeadd_list%rot(i,1:nmergeadd)
end do
symba_plA%helio%swiftest%info(nkpl+1:npl) = mergeadd_list%info(1:nmergeadd)
end if

if (ldiscard_pl) then
call coord_b2h(npl, symba_plA%helio%swiftest)
! Create the particle information and set the status flags of all new particles
allocate(add_l_pl(npl))

where ((symba_plA%helio%swiftest%status(:) == DISRUPTION) .or. &
(symba_plA%helio%swiftest%status(:) == SUPERCATASTROPHIC) .or. &
(symba_plA%helio%swiftest%status(:) == HIT_AND_RUN))
symba_plA%helio%swiftest%info(:)%origin_time = t
symba_plA%helio%swiftest%info(:)%origin_xh(1) = symba_plA%helio%swiftest%xh(1,:)
symba_plA%helio%swiftest%info(:)%origin_xh(2) = symba_plA%helio%swiftest%xh(2,:)
symba_plA%helio%swiftest%info(:)%origin_xh(3) = symba_plA%helio%swiftest%xh(3,:)
symba_plA%helio%swiftest%info(:)%origin_vh(1) = symba_plA%helio%swiftest%vh(1,:)
symba_plA%helio%swiftest%info(:)%origin_vh(2) = symba_plA%helio%swiftest%vh(2,:)
symba_plA%helio%swiftest%info(:)%origin_vh(3) = symba_plA%helio%swiftest%vh(3,:)
add_l_pl(:) = .true.
elsewhere
add_l_pl(:) = .false.
end where

! check for duplicate names and fix if ncessary
do j = 1, npl
do i = j + 1, npl
if (symba_plA%helio%swiftest%id(i) == symba_plA%helio%swiftest%id(j)) then
symba_plA%helio%swiftest%id(i) = maxval(symba_plA%helio%swiftest%id(:)) + 1
end if
end do
end do
symba_plA%helio%swiftest%maxid = max(symba_plA%helio%swiftest%maxid, maxval(symba_plA%helio%swiftest%id))
if (nmergeadd > 0) then
call util_resize_pl(symba_plA, nkpl+nmergeadd)
npl = nkpl + nmergeadd
symba_plA%helio%swiftest%nbody = npl

call io_write_particle_pl(symba_plA%helio%swiftest, pack([(i, i=1,npl)], add_l_pl(:)), param)
!add merge products to the end of the planet list
symba_plA%helio%swiftest%status(nkpl+1:npl) = mergeadd_list%status(1:nmergeadd)
symba_plA%helio%swiftest%id(nkpl+1:npl) = mergeadd_list%id(1:nmergeadd)
symba_plA%helio%swiftest%mass(nkpl+1:npl) = mergeadd_list%mass(1:nmergeadd)
symba_plA%helio%swiftest%radius(nkpl+1:npl) = mergeadd_list%radius(1:nmergeadd)
do i = 1, NDIM
symba_plA%helio%swiftest%xb(i,nkpl+1:npl) = mergeadd_list%xb(i,1:nmergeadd)
symba_plA%helio%swiftest%vb(i,nkpl+1:npl) = mergeadd_list%vb(i,1:nmergeadd)
symba_plA%helio%swiftest%Ip(i,nkpl+1:npl) = mergeadd_list%Ip(i,1:nmergeadd)
symba_plA%helio%swiftest%rot(i,nkpl+1:npl) = mergeadd_list%rot(i,1:nmergeadd)
end do
symba_plA%helio%swiftest%info(nkpl+1:npl) = mergeadd_list%info(1:nmergeadd)
end if

if (ldiscard_pl) then
call coord_b2h(npl, symba_plA%helio%swiftest)
! Create the particle information and set the status flags of all new particles
allocate(add_l_pl(npl))

where ((symba_plA%helio%swiftest%status(:) == DISRUPTION) .or. &
(symba_plA%helio%swiftest%status(:) == SUPERCATASTROPHIC) .or. &
(symba_plA%helio%swiftest%status(:) == HIT_AND_RUN))
symba_plA%helio%swiftest%info(:)%origin_time = t
symba_plA%helio%swiftest%info(:)%origin_xh(1) = symba_plA%helio%swiftest%xh(1,:)
symba_plA%helio%swiftest%info(:)%origin_xh(2) = symba_plA%helio%swiftest%xh(2,:)
symba_plA%helio%swiftest%info(:)%origin_xh(3) = symba_plA%helio%swiftest%xh(3,:)
symba_plA%helio%swiftest%info(:)%origin_vh(1) = symba_plA%helio%swiftest%vh(1,:)
symba_plA%helio%swiftest%info(:)%origin_vh(2) = symba_plA%helio%swiftest%vh(2,:)
symba_plA%helio%swiftest%info(:)%origin_vh(3) = symba_plA%helio%swiftest%vh(3,:)
add_l_pl(:) = .true.
elsewhere
add_l_pl(:) = .false.
end where

! check for duplicate names and fix if ncessary
do j = 1, npl
do i = j + 1, npl
if (symba_plA%helio%swiftest%id(i) == symba_plA%helio%swiftest%id(j)) then
symba_plA%helio%swiftest%id(i) = maxval(symba_plA%helio%swiftest%id(:)) + 1
end if
end do
end do
symba_plA%helio%swiftest%maxid = max(symba_plA%helio%swiftest%maxid, maxval(symba_plA%helio%swiftest%id))

symba_plA%helio%swiftest%status(1:npl) = ACTIVE
call symba_reorder_pl(npl, symba_plA)

call util_hills(npl, symba_plA%helio%swiftest)
call io_write_particle_pl(symba_plA%helio%swiftest, pack([(i, i=1,npl)], add_l_pl(:)), param)

nplm = count(symba_plA%helio%swiftest%mass > param%mtiny)
CALL util_dist_index_plpl(npl, nplm, symba_plA)
symba_plA%helio%swiftest%status(1:npl) = ACTIVE
call symba_reorder_pl(npl, symba_plA)

call util_hills(npl, symba_plA%helio%swiftest)

end if
nplm = count(symba_plA%helio%swiftest%mass > param%mtiny)
CALL util_dist_index_plpl(npl, nplm, symba_plA)
end if
end associate


end subroutine symba_rearray_pl

0 comments on commit 49a0b33

Please sign in to comment.