From cceaf460bcb700188d5a2f47ef33c5fb0bad3cc4 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Sat, 24 Jul 2021 17:24:28 -0400 Subject: [PATCH] Simplified call to recursive step and fleshed out recursive step algorithm --- src/modules/symba_classes.f90 | 50 +++++++++++++++-------- src/symba/symba_kick.f90 | 35 ++++++++++++++++ src/symba/symba_step.f90 | 75 ++++++++++++++++++++++++----------- 3 files changed, 120 insertions(+), 40 deletions(-) create mode 100644 src/symba/symba_kick.f90 diff --git a/src/modules/symba_classes.f90 b/src/modules/symba_classes.f90 index 877f5b585..345b1b31c 100644 --- a/src/modules/symba_classes.f90 +++ b/src/modules/symba_classes.f90 @@ -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 !******************************************************************************************************************************** @@ -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 !******************************************************************************************************************************** @@ -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 @@ -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 diff --git a/src/symba/symba_kick.f90 b/src/symba/symba_kick.f90 new file mode 100644 index 000000000..6d438a950 --- /dev/null +++ b/src/symba/symba_kick.f90 @@ -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 \ No newline at end of file diff --git a/src/symba/symba_step.f90 b/src/symba/symba_step.f90 index cc7090065..b56db66d5 100644 --- a/src/symba/symba_step.f90 +++ b/src/symba/symba_step.f90 @@ -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) @@ -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 @@ -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