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 
Fixed a bunch of mismatched interfaces and trying to troubleshoot energy growth during discards with the central body
  • Loading branch information
@daminton
daminton committed May 20, 2021
1 parent 314a5ff commit 5c2b8e4
Show file tree
Hide file tree
Showing 8 changed files with 149 additions and 154 deletions.
18 changes: 16 additions & 2 deletions src/io/io_conservation_report.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -15,8 +15,8 @@ module subroutine io_conservation_report(t, symba_plA, npl, j2rp2, j4rp4, param,
type(user_input_parameters), intent(in) :: param !! Input colleciton of user-defined parameters
logical, intent(in) :: lterminal !! Indicates whether to output information to the terminal screen
! Internals
real(DP), dimension(NDIM), save :: Ltot_orig
real(DP), save :: Eorbit_orig, Mtot_orig, Lmag_orig
real(DP), dimension(NDIM), save :: Ltot_orig, Ltot_last
real(DP), save :: Eorbit_orig, Mtot_orig, Lmag_orig, ke_orb_last, ke_spin_last, pe_last, Eorbit_last
real(DP) :: ke_orbit, ke_spin, pe, Eorbit
real(DP), dimension(NDIM) :: Ltot_now
real(DP) :: Eorbit_error, Etotal_error, Ecoll_error
Expand Down Expand Up @@ -62,7 +62,21 @@ module subroutine io_conservation_report(t, symba_plA, npl, j2rp2, j4rp4, param,
Etotal_error = (Eorbit - (Eorbit_orig - Ecollisions)) / abs(Eorbit_orig)
Merror = (Mtot_now - Mtot_orig) / Mtot_orig
write(*, egytermfmt) Lerror, Ecoll_error, Etotal_error, Merror
if (Ecoll_error > 0.0_DP) then
write(*,*) 'Something has gone wrong! Collisional energy should not be positive!'
write(*,*) 'Comparisons with last time: '
write(*,*) 'ke_orb : ',ke_orbit - ke_orb_last, (ke_orbit - ke_orb_last) / abs(Eorbit_last)
write(*,*) 'ke_spin: ',ke_spin - ke_spin_last, (ke_spin - ke_spin_last) / abs(Eorbit_last)
write(*,*) 'pe : ',pe - pe_last, (pe - pe_last) / abs(Eorbit_last)
write(*,*) 'Etot : ',Eorbit - Eorbit_last, (Eorbit - Eorbit_last) / abs(Eorbit_last)
write(*,*)
end if
end if
ke_orb_last = ke_orbit
ke_spin_last = ke_spin
pe_last = pe
Ltot_last(:) = Ltot_now(:)
Eorbit_last = Eorbit
end associate
return

Expand Down
11 changes: 8 additions & 3 deletions src/main/swiftest_symba.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -42,6 +42,8 @@ program swiftest_symba
character(len=*), parameter :: endwallfmt = '("Wall time to complete run (s): ", es12.5)'
integer(I4B), allocatable :: seed(:)
integer(I4B) :: nseeds
integer(I4B), dimension(:), allocatable :: discard_stat_list
logical, dimension(:), allocatable :: discard_l_pl

! Executable code
! temporary until the conversion to the derived type argument list is complete
Expand Down Expand Up @@ -160,10 +162,13 @@ program swiftest_symba
call symba_step_eucl(t, dt, param,npl,ntp,symba_plA, symba_tpA, nplplenc, npltpenc,&
plplenc_list, pltpenc_list, nmergeadd, nmergesub, mergeadd_list, mergesub_list)

ldiscard_pl = .false.
if (allocated(discard_l_pl)) deallocate(discard_l_pl)
allocate(discard_l_pl(npl))
discard_l_pl(:) = .false.
ldiscard_tp = .false.
lfrag_add = .false.
call symba_discard_pl(t, npl, ntp, symba_plA, symba_tpA, rmin, rmax, rmaxu, qmin, qmin_coord, qmin_alo, qmin_ahi, ldiscard_pl)
call symba_discard_pl(t, npl, ntp, symba_plA, symba_tpA, rmin, rmax, rmaxu, qmin, qmin_coord, qmin_alo, qmin_ahi, discard_l_pl, discard_stat_list)
ldiscard_pl = any(discard_l_pl(:))
call symba_discard_tp(t, npl, ntp, symba_plA, symba_tpA, dt, rmin, rmax, rmaxu, qmin, qmin_coord, &
qmin_alo, qmin_ahi, param%lrhill_present, ldiscard_tp)
call symba_collision(t, symba_plA, nplplenc, plplenc_list, lfrag_add, mergeadd_list, nmergeadd, param)
Expand All @@ -173,7 +178,7 @@ program swiftest_symba
if (ldiscard_pl .or. ldiscard_tp) 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(t, npl, nplm, ntp, nsppl, nsptp, symba_plA, symba_tpA, nmergeadd, mergeadd_list, discard_plA, &
discard_tpA, ldiscard_pl, ldiscard_tp, mtiny, param)
discard_tpA, ldiscard_pl, ldiscard_tp, mtiny, param, discard_l_pl, discard_stat_list)
call io_discard_write_symba(t, mtiny, npl, nsppl, nsptp, nmergesub, symba_plA, &
discard_plA%helio%swiftest, discard_tpA%helio%swiftest, mergeadd_list, mergesub_list, discard_file, param%lbig_discard)
nmergeadd = 0
Expand Down
69 changes: 34 additions & 35 deletions src/modules/module_interfaces.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -757,8 +757,8 @@ function symba_casemerge (symba_plA, family, nmergeadd, mergeadd_list, x, v, mas
integer(I4B), dimension(:), intent(in) :: family
integer(I4B), intent(inout) :: nmergeadd
type(symba_merger), intent(inout) :: mergeadd_list
real(DP), dimension(:,:), intent(inout) :: x, v, lspin, Ip
real(DP), dimension(:), intent(inout) :: mass, radius
real(DP), dimension(:,:), intent(in) :: x, v, lspin, Ip
real(DP), dimension(:), intent(in) :: mass, radius
type(user_input_parameters),intent(inout) :: param
integer(I4B) :: status
end function symba_casemerge
Expand All @@ -773,8 +773,8 @@ function symba_casesupercatastrophic (symba_plA, family, nmergeadd, mergeadd_lis
integer(I4B), dimension(:), intent(in) :: family
integer(I4B), intent(inout) :: nmergeadd
type(symba_merger), intent(inout) :: mergeadd_list
real(DP), dimension(:,:), intent(in) :: x, v, lspin, Ip
real(DP), dimension(:), intent(in) :: mass, radius, mass_res
real(DP), dimension(:,:), intent(inout) :: x, v, lspin, Ip
real(DP), dimension(:), intent(inout) :: mass, radius, mass_res
type(user_input_parameters),intent(inout) :: param
real(DP), intent(inout) :: Qloss
integer(I4B) :: status
Expand Down Expand Up @@ -846,20 +846,21 @@ END SUBROUTINE symba_discard_peri_pl
END INTERFACE

INTERFACE
SUBROUTINE symba_discard_pl(t, npl, ntp, symba_plA, symba_tpA, rmin, rmax, rmaxu, qmin, qmin_coord, &
qmin_alo, qmin_ahi, ldiscard)
USE swiftest_globals
USE swiftest_data_structures
USE module_helio
USE module_symba
IMPLICIT NONE
INTEGER(I4B), INTENT(INOUT) :: npl, ntp
REAL(DP), INTENT(IN) :: t, rmin, rmax, rmaxu, qmin, qmin_alo, qmin_ahi
CHARACTER(*), INTENT(IN) :: qmin_coord
TYPE(symba_pl), INTENT(INOUT) :: symba_plA
TYPE(symba_tp), INTENT(INOUT) :: symba_tpA
LOGICAL(LGT), INTENT(INOUT) :: ldiscard
END SUBROUTINE symba_discard_pl
subroutine symba_discard_pl(t, npl, ntp, symba_pla, symba_tpa, rmin, rmax, rmaxu, qmin, qmin_coord, &
qmin_alo, qmin_ahi, discard_l_pl, discard_stat_list)
use swiftest_globals
use swiftest_data_structures
use module_helio
use module_symba
implicit none
integer(I4B), intent(inout) :: npl, ntp
real(DP), intent(in) :: t, rmin, rmax, rmaxu, qmin, qmin_alo, qmin_ahi
character(*), intent(in) :: qmin_coord
type(symba_pl), intent(inout) :: symba_pla
type(symba_tp), intent(inout) :: symba_tpa
logical, dimension(:), intent(out) :: discard_l_pl
integer(I4B), dimension(:), allocatable, intent(out) :: discard_stat_list
end subroutine symba_discard_pl
END INTERFACE

INTERFACE
Expand Down Expand Up @@ -934,7 +935,7 @@ subroutine symba_frag_pos(param, symba_plA, family, x, v, L_spin, Ip, mass, radi
integer(I4B), intent(inout) :: nfrag
real(DP), dimension(:), intent(inout) :: m_frag, rad_frag
real(DP), dimension(:,:), intent(inout) :: Ip_frag
real(DP), dimension(:,:), intent(out) :: xb_frag, vb_frag, rot_frag
real(DP), dimension(:,:), intent(inout) :: xb_frag, vb_frag, rot_frag
real(DP), intent(inout) :: Qloss
logical, intent(out) :: lmerge
end subroutine symba_frag_pos
Expand All @@ -952,7 +953,7 @@ SUBROUTINE symba_getacch(lextra_force, t, npl, nplm, symba_plA, j2rp2, j4rp4, np
INTEGER(I4B), INTENT(IN) :: npl, nplm, nplplenc
REAL(DP), INTENT(IN) :: t, j2rp2, j4rp4
TYPE(symba_pl), INTENT(INOUT) :: symba_plA
TYPE(symba_plplenc), INTENT(IN) :: plplenc_list
TYPE(symba_plplenc), INTENT(INOUT) :: plplenc_list
END SUBROUTINE symba_getacch
END INTERFACE

Expand Down Expand Up @@ -1077,7 +1078,7 @@ END SUBROUTINE symba_peri

INTERFACE
subroutine symba_rearray(t, npl, nplm, ntp, nsppl, nsptp, symba_plA, symba_tpA, nmergeadd, mergeadd_list, discard_plA, &
discard_tpA, ldiscard_pl, ldiscard_tp, mtiny, param)
discard_tpA, ldiscard_pl, ldiscard_tp, mtiny, param, discard_l_pl, discard_stat_list)
use swiftest_globals
use swiftest_data_structures
use module_symba
Expand All @@ -1093,6 +1094,8 @@ subroutine symba_rearray(t, npl, nplm, ntp, nsppl, nsptp, symba_plA, symba_tpA,
logical(LGT), intent(in) :: ldiscard_pl, ldiscard_tp
real(DP), intent(in) :: mtiny
type(user_input_parameters), intent(in) :: param
logical, dimension(:), allocatable, intent(inout) :: discard_l_pl
integer(I4B), dimension(:), allocatable, intent(inout) :: discard_stat_list
end subroutine symba_rearray

END INTERFACE
Expand Down Expand Up @@ -1308,7 +1311,7 @@ SUBROUTINE symba_getacch_eucl(lextra_force, t, npl, symba_plA, j2rp2, j4rp4, npl
INTEGER(I4B), INTENT(IN) :: npl, nplplenc
REAL(DP), INTENT(IN) :: t, j2rp2, j4rp4
TYPE(symba_pl), INTENT(INOUT) :: symba_plA
TYPE(symba_plplenc), INTENT(IN) :: plplenc_list
TYPE(symba_plplenc), INTENT(INOUT) :: plplenc_list
END SUBROUTINE symba_getacch_eucl
END INTERFACE

Expand Down Expand Up @@ -1649,19 +1652,15 @@ FUNCTION collresolve_resolve(model,m1,m2,r1,r2,p1,p2,v1,v2,n,mres,rres,pres,vres
END INTERFACE

INTERFACE
SUBROUTINE symba_regime(Mcb, m1, m2, rad1, rad2, xh1, xh2, vb1, vb2, den1, den2, regime, Mlr, Mslr, mtiny, Qloss)
USE swiftest_globals
USE module_symba
USE swiftest_data_structures
USE module_helio
USE module_nrutil
USE module_swiftestalloc
IMPLICIT NONE
INTEGER(I4B), INTENT(OUT) :: regime
REAL(DP), INTENT(INOUT) :: Mcb, Mlr, Mslr, m1, m2, rad1, rad2, den1, den2, mtiny
REAL(DP), DIMENSION(:), INTENT(IN) :: xh1, xh2, vb1, vb2
real(DP), intent(out) :: Qloss
END SUBROUTINE symba_regime
subroutine symba_regime(mcb, m1, m2, rad1, rad2, xh1, xh2, vb1, vb2, den1, den2, regime, mlr, mslr, mtiny, qloss)
use swiftest_globals
implicit none
integer(I4B), intent(out) :: regime
real(DP), intent(out) :: Mlr, Mslr
real(DP), intent(in) :: Mcb, m1, m2, rad1, rad2, den1, den2, mtiny
real(DP), dimension(:), intent(in) :: xh1, xh2, vb1, vb2
real(DP), intent(out) :: Qloss
end subroutine symba_regime
END INTERFACE

END MODULE module_interfaces
Expand Down
2 changes: 1 addition & 1 deletion src/symba/symba_discard_conserve_mtm.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -16,7 +16,7 @@ subroutine symba_discard_conserve_mtm(swiftest_plA, ipl, lescape)
associate(npl => swiftest_plA%nbody, xb => swiftest_plA%xb, vb => swiftest_plA%vb, &
rot => swiftest_plA%rot, Ip => swiftest_plA%Ip, rad => swiftest_plA%radius, mass => swiftest_plA%mass, &
statipl => swiftest_plA%status(ipl), xbpl => swiftest_plA%xb(:,ipl), xbcb => swiftest_plA%xb(:,1))

xcom(:) = (mass(ipl) * xb(:, ipl) + mass(1) * xb(:,1)) / (mass(1) + mass(ipl))
vcom(:) = (mass(ipl) * vb(:, ipl) + mass(1) * vb(:,1)) / (mass(1) + mass(ipl))

Expand Down
130 changes: 48 additions & 82 deletions src/symba/symba_discard_pl.f90
View file
Original file line number Diff line number Diff line change
@@ -1,87 +1,53 @@
!**********************************************************************************************************************************
!
! Unit Name : symba_discard_pl
! Unit Type : subroutine
! Project : Swiftest
! Package : symba
! Language : Fortran 90/95
!
! Description : Check to see if planets should be discarded based on their positions or because they are unbound
!
! Input
! Arguments : t : time
! npl : number of planets
! symba_pl1P : pointer to head of SyMBA planet structure linked-list
! symba_pld1P : pointer to head of discard SyMBA planet structure linked-list
! rmin : minimum allowed heliocentric radius
! rmax : maximum allowed heliocentric radius
! rmaxu : maximum allowed heliocentric radius for unbound planets
! qmin : minimum allowed pericenter distance
! qmin_coord : coordinate frame for qmin
! qmin_alo : minimum semimajor axis for qmin
! qmin_ahi : maximum semimajor axis for qmin
! ldiscard : logical flag indicating that at least one body needs to be discarded this step
! Terminal : none
! File : none
!
! Output
! Arguments : npl : number of planets
! symba_pl1P : pointer to head of SyMBA planet structure linked-list
! symba_pld1P : pointer to head of discard SyMBA planet structure linked-list
! eoffset : energy offset
! Terminal : none
! File : none
!
! Invocation : CALL symba_discard_pl(t, npl, symba_pl1P, symba_pld1P, rmin, rmax, rmaxu, qmin, qmin_coord,
! qmin_alo, qmin_ahi)
!
! Notes : Adapted from Hal Levison's Swift routine discard_massive5.f
!
!**********************************************************************************************************************************
SUBROUTINE symba_discard_pl(t, npl, ntp, symba_plA, symba_tpA, rmin, rmax, rmaxu, qmin, qmin_coord, qmin_alo, qmin_ahi, ldiscard)
subroutine symba_discard_pl(t, npl, ntp, symba_plA, symba_tpA, rmin, rmax, rmaxu, &
qmin, qmin_coord, qmin_alo, qmin_ahi, discard_l_pl, discard_stat_list)
!! author: David A. Minton
!!
!! Check to see if planets should be discarded based on their positions or because they are unbound
!!
!! Adapted from David E. Kaufmann Swifter routine symba_discard_pl.f90
!! Adapted from Martin Duncan and Hal Levison's Swift routine discard_massive5.f
use swiftest
use module_helio
use module_symba
use module_interfaces, except_this_one => symba_discard_pl
implicit none

! Modules
USE swiftest
USE module_helio
USE module_symba
USE module_interfaces, EXCEPT_THIS_ONE => symba_discard_pl
IMPLICIT NONE
! Arguments
integer(I4B), intent(inout) :: npl, ntp
real(DP), intent(in) :: t, rmin, rmax, rmaxu, qmin, qmin_alo, qmin_ahi
character(*), intent(in) :: qmin_coord
type(symba_pl), intent(inout) :: symba_plA
type(symba_tp), intent(inout) :: symba_tpa
logical, dimension(:), intent(out) :: discard_l_pl
integer(I4B), dimension(:), allocatable, intent(out) :: discard_stat_list

! Arguments
INTEGER(I4B), INTENT(INOUT) :: npl, ntp
REAL(DP), INTENT(IN) :: t, rmin, rmax, rmaxu, qmin, qmin_alo, qmin_ahi
CHARACTER(*), INTENT(IN) :: qmin_coord
TYPE(symba_pl) :: symba_plA
TYPE(symba_tp) :: symba_tpA
LOGICAL(LGT), INTENT(INOUT) :: ldiscard
! Internals
real(DP) :: msys
integer(I4B) :: i, ndiscard
logical :: ldiscard

! Internals
REAL(DP) :: msys
! Executable code
if ((rmin >= 0.0_DP) .or. (rmax >= 0.0_DP) .or. (rmaxu >= 0.0_DP) .or. ((qmin >= 0.0_DP) .and. (qmin_coord == "bary"))) &
call coord_h2b(npl, symba_plA%helio%swiftest, msys)
if ((rmin >= 0.0_DP) .or. (rmax >= 0.0_DP) .or. (rmaxu >= 0.0_DP)) &
call symba_discard_sun_pl(t, npl, ntp, msys, symba_plA%helio%swiftest, symba_tpa%helio%swiftest, rmin, rmax, rmaxu, ldiscard)
!if (qmin >= 0.0_DP) call symba_discard_peri_pl(t, npl, symba_plA, msys, qmin, qmin_alo, qmin_ahi, qmin_coord, ldiscard)
if (.not.ldiscard) return

! Executable code
IF ((rmin >= 0.0_DP) .OR. (rmax >= 0.0_DP) .OR. (rmaxu >= 0.0_DP) .OR. ((qmin >= 0.0_DP) .AND. (qmin_coord == "BARY"))) &
CALL coord_h2b(npl, symba_plA%helio%swiftest, msys)
IF ((rmin >= 0.0_DP) .OR. (rmax >= 0.0_DP) .OR. (rmaxu >= 0.0_DP)) &
CALL symba_discard_sun_pl(t, npl, ntp, msys, symba_plA%helio%swiftest, symba_tpA%helio%swiftest, rmin, rmax, rmaxu, ldiscard)
!IF (qmin >= 0.0_DP) CALL symba_discard_peri_pl(t, npl, symba_plA, msys, qmin, qmin_alo, qmin_ahi, qmin_coord, ldiscard)
RETURN
! We need to keep track of the bodies that are discarded via collision or escape from the central body separately from those that collide with each other.
associate(status => symba_plA%helio%swiftest%status)
discard_l_pl(1:npl) = (status(1:npl) == DISCARDED_RMIN) .or. (status(1:npl) == DISCARDED_PERI) .or. (status(1:npl) == DISCARDED_RMAX) .or. (status(1:npl) == DISCARDED_RMAXU)
ndiscard = count(discard_l_pl(1:npl))
if (ndiscard > 0) then
if (allocated(discard_stat_list)) deallocate(discard_stat_list)
allocate(discard_stat_list(ndiscard))
discard_stat_list = pack(status(1:npl), discard_l_pl)
where(discard_l_pl(1:npl))
status(:) = ACTIVE
end where
end if
end associate

return

END SUBROUTINE symba_discard_pl
!**********************************************************************************************************************************
!
! Author(s) : David E. Kaufmann
!
! Revision Control System (RCS) Information
!
! Source File : $RCSfile$
! Full Path : $Source$
! Revision : $Revision$
! Date : $Date$
! Programmer : $Author$
! Locked By : $Locker$
! State : $State$
!
! Modification History:
!
! $Log$
!**********************************************************************************************************************************
end subroutine symba_discard_pl
18 changes: 9 additions & 9 deletions src/symba/symba_energy_eucl.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -12,17 +12,17 @@ subroutine symba_energy_eucl(npl, symba_plA, j2rp2, j4rp4, ke_orbit, ke_spin, pe
implicit none

! arguments
integer(I4B), intent(in) :: npl
real(DP), intent(in) :: j2rp2, j4rp4
real(DP), intent(out) :: ke_orbit, ke_spin, pe, te
real(DP), dimension(:), intent(out) :: Ltot
type(symba_pl), intent(inout) :: symba_plA
integer(I4B), intent(in) :: npl
real(DP), intent(in) :: j2rp2, j4rp4
real(DP), intent(out) :: ke_orbit, ke_spin, pe, te
real(DP), dimension(:), intent(out) :: Ltot
type(symba_pl), intent(inout) :: symba_plA

! internals
integer(I4B) :: i, j
integer(I8B) :: k
real(DP) :: rmag, v2, rot2, oblpot, hx, hy, hz, hsx, hsy, hsz
real(DP), dimension(npl) :: irh, kepl, kespinpl, pecb
integer(I4B) :: i, j
integer(I8B) :: k
real(DP) :: rmag, v2, rot2, oblpot, hx, hy, hz, hsx, hsy, hsz
real(DP), dimension(npl) :: irh, kepl, kespinpl, pecb
real(DP), dimension(npl) :: Lplx, Lply, Lplz
real(DP), dimension(symba_plA%helio%swiftest%num_plpl_comparisons) :: pepl
logical, dimension(symba_plA%helio%swiftest%num_plpl_comparisons) :: lstatpl
Expand Down
Loading

0 comments on commit 5c2b8e4

Please sign in to comment.