diff --git a/src/modules/symba_classes.f90 b/src/modules/symba_classes.f90 index 1b3261f62..2aadd86e0 100644 --- a/src/modules/symba_classes.f90 +++ b/src/modules/symba_classes.f90 @@ -13,6 +13,9 @@ module symba_classes !******************************************************************************************************************************** type, public, extends(helio_nbody_system) :: symba_nbody_system contains + private + procedure, public :: step => symba_step_system !! Advance the SyMBA nbody system forward in time by one step + procedure, public :: interp => symba_step_interp_system !! Perform an interpolation step on the SymBA nbody system end type symba_nbody_system !******************************************************************************************************************************** @@ -40,4 +43,24 @@ module symba_classes type, public, extends(helio_tp) :: symba_tp contains end type symba_tp + + interface + module subroutine symba_step_system(self, param, t, dt) + use swiftest_classes, only : swiftest_parameters + implicit none + class(symba_nbody_system), intent(inout) :: self !! RMVS 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 + end subroutine symba_step_system + + module subroutine symba_step_interp_system(self, param, t, dt) + use swiftest_classes, only : swiftest_parameters + implicit none + class(symba_nbody_system), intent(inout) :: self !! RMVS 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 + end subroutine symba_step_interp_system + end interface end module symba_classes \ No newline at end of file diff --git a/src/symba/symba_step.f90 b/src/symba/symba_step.f90 index cb1361833..d3b189bfa 100644 --- a/src/symba/symba_step.f90 +++ b/src/symba/symba_step.f90 @@ -1,127 +1,39 @@ -submodule (symba) s_symba_step +submodule (symba_classes) s_symba_step + use swiftest contains - module procedure symba_step - !! author: David A. Minton - !! - !! Step planets and active test particles ahead in democratic heliocentric coordinates, descending the recursive - !! branch if necessary to handle possible close encounters - !! - !! Adapted from David E. Kaufmann's Swifter routine: symba_step.f90 - !! Adapted from Hal Levison's Swift routine symba5_step_pl.f -use swiftest -implicit none - logical :: lencounter, lvdotr - integer(I4B) :: i, j, irec, nplm - real(DP), dimension(NDIM) :: xr, vr - logical, save :: lfirst = .true. - -! executable code - - do i = 1,npl - symba_plA%nplenc(i) = 0 - symba_plA%ntpenc(i) = 0 - symba_plA%levelg(i) = -1 - symba_plA%levelm(i) = -1 - symba_plA%index_parent(i) = i - symba_plA%index_child(:,i) = 0 - end do - do i =1,ntp - symba_tpA%nplenc(i) = 0 - symba_tpA%levelg(i) = -1 - symba_tpA%levelm(i) = -1 - end do - - - !there should be some parallel bits in here - - nplplenc = 0 - npltpenc = 0 - if (symba_plA%mass(1) < param%mtiny) then - nplm = 0 - else - nplm = 1 - end if - irec = 0 - -! all this needs to be changed to the tree search function for encounters - - do i = 2, npl - if (symba_plA%mass(i) < param%mtiny) exit - nplm = nplm + 1 - do j = i + 1, npl - xr(:) = symba_plA%xh(:,j) - symba_plA%xh(:,i) - vr(:) = symba_plA%vh(:,j) - symba_plA%vh(:,i) - call symba_chk(xr(:), vr(:), symba_plA%rhill(i), & - symba_plA%rhill(j), dt, irec, lencounter, lvdotr) - if (lencounter) then - nplplenc = nplplenc + 1 - if (nplplenc > nenmax) then - write(*, *) "Swiftest Error:" - write(*, *) " pl-pl encounter list is full." - write(*, *) " stopping..." - call util_exit(FAILURE) - end if - plplenc_list%status(nplplenc) = ACTIVE - plplenc_list%lvdotr(nplplenc) = lvdotr - plplenc_list%level(nplplenc) = irec - plplenc_list%index1(nplplenc) = i - plplenc_list%index2(nplplenc) = j - symba_plA%lmerged(i) = .false. - symba_plA%nplenc(i) = symba_plA%nplenc(i) + 1 - symba_plA%levelg(i) = irec - symba_plA%levelm(i) = irec - symba_plA%nchild(i) = 0 - symba_plA%lmerged(j) = .false. - symba_plA%nplenc(j) = symba_plA%nplenc(j) + 1 - symba_plA%levelg(j) = irec - symba_plA%levelm(j) = irec - symba_plA%nchild(j) = 0 - end if - end do - do j = 1, ntp - xr(:) = symba_tpA%xh(:,j) - symba_plA%xh(:,i) - vr(:) = symba_tpA%vh(:,j) - symba_plA%vh(:,i) - call symba_chk(xr(:), vr(:), symba_plA%rhill(i), 0.0_DP, dt, irec, lencounter, lvdotr) - if (lencounter) then - npltpenc = npltpenc + 1 - symba_plA%ntpenc(i) = symba_plA%ntpenc(i) + 1 - symba_plA%levelg(i) = irec - symba_plA%levelm(i) = irec - symba_tpA%nplenc(j) = symba_tpA%nplenc(j) + 1 - symba_tpA%levelg(j) = irec - symba_tpA%levelm(j) = irec - if (npltpenc > nenmax) then - write(*, *) "Swiftest Error:" - write(*, *) " pl-tp encounter list is full." - write(*, *) " stopping..." - call util_exit(FAILURE) + module subroutine symba_step_system(self, param, t, dt) + !! author: David A. Minton + !! + !! Step planets and active test particles ahead in democratic heliocentric coordinates, descending the recursive + !! branch if necessary to handle possible close encounters + !! + !! Adapted from David E. Kaufmann's Swifter routine: symba_step.f90 + !! Adapted from Hal Levison's Swift routine symba5_step_pl.f + implicit none + ! 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 + ! Internals + logical :: lencounter_pl, lencounter_tp, lencounter + + select type(pl => self%pl) + class is (symba_pl) + select type(tp => self%tp) + class is (symba_tp) + lencounter_pl = pl%encounter_check(system, dt) + lencounter_tp = tp%encounter_check(system, dt) + lencounter = lencounter_pl .or. lencounter_tp + if (lencounter) then + call self%interp(param, t, dt) + else + call helio_step_system(self, param, t, dt) end if - pltpenc_list%status(npltpenc) = ACTIVE - pltpenc_list%lvdotr(npltpenc) = lvdotr - pltpenc_list%level(npltpenc) = irec - pltpenc_list%indexpl(npltpenc) = i - pltpenc_list%indextp(npltpenc) = j - end if - end do - end do - -! end of things that need to be changed in the tree - - lencounter = ((nplplenc > 0) .or. (npltpenc > 0)) - if (lencounter) then - call symba_step_interp(param%lextra_force, param%lclose, t, npl, nplm, param%nplmax, & - ntp, param%ntpmax, symba_plA, symba_tpA, param%j2rp2, param%j4rp4, & - dt, eoffset, nplplenc, npltpenc, plplenc_list, pltpenc_list, nmergeadd, & - nmergesub, mergeadd_list, mergesub_list, param%encounter_file, param%out_type, & - fragmax, param) - lfirst = .true. - else - call symba_step_helio(lfirst, param%lextra_force, t, npl, nplm, param%nplmax, ntp,& - param%ntpmax, symba_plA, symba_tpA, & - param%j2rp2, param%j4rp4, dt) - end if + end select + end select - return + return - end procedure symba_step + end subroutine symba_step_system end submodule s_symba_step