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

Commit

Permalink
Simplified call to recursive step and fleshed out recursive step algo…
Browse files Browse the repository at this point in the history 
…rithm
  • Loading branch information
@daminton
daminton committed Jul 24, 2021
1 parent 43f6b4c commit cceaf46
Show file tree
Hide file tree
Showing 3 changed files with 120 additions and 40 deletions.
50 changes: 34 additions & 16 deletions src/modules/symba_classes.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -122,9 +122,10 @@ module symba_classes
integer(I4B), dimension(:), allocatable :: index2 !! position of the test particle in encounter
contains
procedure, public :: encounter_check => symba_encounter_check_pltpenc !! Checks if massive bodies are going through close encounters with each other
procedure, public :: setup => symba_setup_pltpenc !! A constructor that sets the number of encounters and allocates and initializes all arrays
procedure, public :: copy => symba_util_copy_pltpenc !! Copies all elements of one pltpenc list to another
procedure, public :: resize => symba_util_resize_pltpenc !! Checks the current size of the pltpenc_list against the required size and extends it by a factor of 2 more than requested if it is too small
procedure, public :: kick => symba_kick_pltpenc !! Kick barycentric velocities of active test particles within SyMBA recursion
procedure, public :: setup => symba_setup_pltpenc !! A constructor that sets the number of encounters and allocates and initializes all arrays
procedure, public :: copy => symba_util_copy_pltpenc !! Copies all elements of one pltpenc list to another
procedure, public :: resize => symba_util_resize_pltpenc !! Checks the current size of the pltpenc_list against the required size and extends it by a factor of 2 more than requested if it is too small
end type symba_pltpenc

!********************************************************************************************************************************
Expand All @@ -138,8 +139,9 @@ module symba_classes
real(DP), dimension(:,:), allocatable :: vb2 !! the barycentric velocity of parent 2 in encounter
contains
procedure, public :: encounter_check => symba_encounter_check_plplenc !! Checks if massive bodies are going through close encounters with each other
procedure, public :: setup => symba_setup_plplenc !! A constructor that sets the number of encounters and allocates and initializes all arrays
procedure, public :: copy => symba_util_copy_plplenc !! Copies all elements of one plplenc list to another
procedure, public :: kick => symba_kick_plplenc !! Kick barycentric velocities of massive bodies within SyMBA recursion
procedure, public :: setup => symba_setup_plplenc !! A constructor that sets the number of encounters and allocates and initializes all arrays
procedure, public :: copy => symba_util_copy_plplenc !! Copies all elements of one plplenc list to another
end type symba_plplenc

!********************************************************************************************************************************
Expand Down Expand Up @@ -205,22 +207,40 @@ end function symba_encounter_check_plplenc

module function symba_encounter_check_pltpenc(self, system, dt, irec) result(lany_encounter)
implicit none
class(symba_pltpenc), intent(inout) :: self !! SyMBA pl-pl encounter list object
class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object
real(DP), intent(in) :: dt !! step size
integer(I4B), intent(in) :: irec !! Current recursion level
class(symba_pltpenc), intent(inout) :: self !! SyMBA pl-pl encounter list object
class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object
real(DP), intent(in) :: dt !! step size
integer(I4B), intent(in) :: irec !! Current recursion level
logical :: lany_encounter !! Returns true if there is at least one close encounter
end function symba_encounter_check_pltpenc

module function symba_encounter_check_tp(self, system, dt, irec) result(lany_encounter)
implicit none
class(symba_tp), intent(inout) :: self !! SyMBA test particle object
class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object
real(DP), intent(in) :: dt !! step size
integer(I4B), intent(in) :: irec !! Current recursion level
class(symba_tp), intent(inout) :: self !! SyMBA test particle object
class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object
real(DP), intent(in) :: dt !! step size
integer(I4B), intent(in) :: irec !! Current recursion level
logical :: lany_encounter !! Returns true if there is at least one close encounter
end function symba_encounter_check_tp

module subroutine symba_kick_plplenc(self, system, dt, irec, sgn)
implicit none
class(symba_plplenc), intent(in) :: self !! SyMBA pl-pl encounter list object
class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object
real(DP), intent(in) :: dt !! step size
integer(I4B), intent(in) :: irec !! Current recursion level
integer(I4B), intent(in) :: sgn !! sign to be applied to acceleration
end subroutine symba_kick_plplenc

module subroutine symba_kick_pltpenc(self, system, dt, irec, sgn)
implicit none
class(symba_pltpenc), intent(in) :: self !! SyMBA pl-tp encounter list object
class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object
real(DP), intent(in) :: dt !! step size
integer(I4B), intent(in) :: irec !! Current recursion level
integer(I4B), intent(in) :: sgn !! sign to be applied to acceleration
end subroutine symba_kick_pltpenc

module subroutine symba_io_dump_particle_info(self, param, msg)
use swiftest_classes, only : swiftest_parameters
implicit none
Expand Down Expand Up @@ -325,12 +345,10 @@ module subroutine symba_step_interp_system(self, param, t, dt)
real(DP), intent(in) :: dt !! Current stepsize
end subroutine symba_step_interp_system

module recursive subroutine symba_step_recur_system(self, param, t, dt, ireci)
module recursive subroutine symba_step_recur_system(self, param, ireci)
implicit none
class(symba_nbody_system), intent(inout) :: self !! SyMBA nbody system object
class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters
real(DP), intent(in) :: t !! Simulation time
real(DP), intent(in) :: dt !! Current stepsize
integer(I4B), value, intent(in) :: ireci !! input recursion level
end subroutine symba_step_recur_system

Expand Down
35 changes: 35 additions & 0 deletions src/symba/symba_kick.f90
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,35 @@
submodule(symba_classes) s_symba_kick
use swiftest
contains

module subroutine symba_kick_plplenc(self, system, dt, irec, sgn)
!! author: David A. Minton
!!
!! !! Kick barycentric velocities of massive bodies within SyMBA recursion.
!!
!! Adapted from David E. Kaufmann's Swifter routine: symba_kick.f90
!! Adapted from Hal Levison's Swift routine symba5_kick.f
implicit none
class(symba_plplenc), intent(in) :: self !! SyMBA pl-pl encounter list object
class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object
real(DP), intent(in) :: dt !! step size
integer(I4B), intent(in) :: irec !! Current recursion level
integer(I4B), intent(in) :: sgn !! sign to be applied to acceleration
end subroutine symba_kick_plplenc

module subroutine symba_kick_pltpenc(self, system, dt, irec, sgn)
!! author: David A. Minton
!!
!! !! Kick barycentric velocities of active test particles within SyMBA recursion.
!!
!! Adapted from David E. Kaufmann's Swifter routine: symba_kick.f90
!! Adapted from Hal Levison's Swift routine symba5_kick.f
implicit none
class(symba_pltpenc), intent(in) :: self !! SyMBA pl-tp encounter list object
class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object
real(DP), intent(in) :: dt !! step size
integer(I4B), intent(in) :: irec !! Current recursion level
integer(I4B), intent(in) :: sgn !! sign to be applied to acceleration
end subroutine symba_kick_pltpenc

end submodule s_symba_kick
75 changes: 51 additions & 24 deletions src/symba/symba_step.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -76,7 +76,7 @@ module subroutine symba_step_interp_system(self, param, t, dt)
call pl%drift(system, param, dt, pl%status(:) == ACTIVE)
call tp%drift(system, param, dt, tp%status(:) == ACTIVE)

call system%recursive_step(param, t, dt, 0)
call system%recursive_step(param, 0)

call pl%set_beg_end(xend = pl%xh)
call pl%accel(system, param, t + dt)
Expand All @@ -96,7 +96,7 @@ module subroutine symba_step_interp_system(self, param, t, dt)
return
end subroutine symba_step_interp_system

module recursive subroutine symba_step_recur_system(self, param, t, dt, ireci)
module recursive subroutine symba_step_recur_system(self, param, ireci)
!! author: David A. Minton
!!
!! Step interacting planets and active test particles ahead in democratic heliocentric coordinates at the current
Expand All @@ -108,33 +108,60 @@ module recursive subroutine symba_step_recur_system(self, param, t, dt, ireci)
! Arguments
class(symba_nbody_system), intent(inout) :: self !! SyMBA nbody system object
class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters
real(DP), intent(in) :: t !! Simulation time
real(DP), intent(in) :: dt !! Current stepsize
integer(I4B), value, intent(in) :: ireci !! input recursion level
! Internals
integer(I4B) :: i, j, irecp, nloops
real(DP) :: dtl, dth, sgn
integer(I4B) :: i, j, irecp, nloops, sgn
real(DP) :: dtl, dth
real(DP), dimension(NDIM) :: xr, vr
logical :: lencounter

associate(system => self, pl => self%pl, tp => self%tp, plplenc_list => self%plplenc_list, pltpenc_list => self%pltpenc_list)
dtl = param%dt / (NTENC**ireci)
dth = 0.5_DP * dtl
IF (dtl / param%dt < VSMALL) THEN
write(*, *) "SWIFTEST Warning:"
write(*, *) " In symba_step_recur_system, local time step is too small"
write(*, *) " Roundoff error will be important!"
call util_exit(FAILURE)
END IF
irecp = ireci + 1
if (ireci == 0) then
nloops = 1
else
nloops = NTENC
end if
do j = 1, nloops
lencounter = plplenc_list%encounter_check(system, dtl, irecp) .or. pltpenc_list%encounter_check(system, dtl, irecp)
end do
associate(system => self, plplenc_list => self%plplenc_list, pltpenc_list => self%pltpenc_list)
select type(pl => self%pl)
class is (symba_pl)
select type(tp => self%tp)
class is (symba_tp)
dtl = param%dt / (NTENC**ireci)
dth = 0.5_DP * dtl
IF (dtl / param%dt < VSMALL) THEN
write(*, *) "SWIFTEST Warning:"
write(*, *) " In symba_step_recur_system, local time step is too small"
write(*, *) " Roundoff error will be important!"
call util_exit(FAILURE)
END IF
irecp = ireci + 1
if (ireci == 0) then
nloops = 1
else
nloops = NTENC
end if
do j = 1, nloops
lencounter = plplenc_list%encounter_check(system, dtl, irecp) .or. pltpenc_list%encounter_check(system, dtl, irecp)
sgn = 1
call plplenc_list%kick(system, dth, irecp, sgn)
call pltpenc_list%kick(system, dth, irecp, sgn)
if (ireci /= 0) then
sgn = -1
call plplenc_list%kick(system, dth, irecp, sgn)
call pltpenc_list%kick(system, dth, irecp, sgn)
end if

call pl%drift(system, param, dtl, pl%status(:) == ACTIVE .and. pl%levelg(:) == ireci)
call tp%drift(system, param, dtl, tp%status(:) == ACTIVE .and. tp%levelg(:) == ireci)

if (lencounter) call system%recursive_step(param, irecp)

sgn = 1
call plplenc_list%kick(system, dth, irecp, sgn)
call pltpenc_list%kick(system, dth, irecp, sgn)
if (ireci /= 0) then
sgn = -1
call plplenc_list%kick(system, dth, irecp, sgn)
call pltpenc_list%kick(system, dth, irecp, sgn)
end if

end do
end select
end select
end associate

end subroutine symba_step_recur_system
Expand Down

0 comments on commit cceaf46

Please sign in to comment.