From 8cfa0995cde1df93f27bd7b5e3b7ddae5321c1db Mon Sep 17 00:00:00 2001 From: David Minton Date: Fri, 16 Apr 2021 16:37:43 -0400 Subject: [PATCH] Major restructuring of RMVS class. Moved encounter objects into container to help streamline data management. Made numerous other changes to clean up the code. --- Makefile | 16 +- Makefile.Defines | 2 + src/driftkick/drift.f90 | 29 +- src/driftkick/kick.f90 | 2 +- src/helio/helio_drift.f90 | 4 +- src/helio/helio_step.f90 | 16 +- src/modules/helio_classes.f90 | 26 +- src/modules/rmvs_classes.f90 | 221 +++-------- src/modules/swiftest_classes.f90 | 26 +- src/modules/whm_classes.f90 | 39 +- src/rmvs/rmvs_encounter_check.f90 | 4 +- src/rmvs/rmvs_getacch.f90 | 98 +++-- src/rmvs/rmvs_interp.f90 | 125 ------ src/rmvs/rmvs_obl.f90 | 31 -- src/rmvs/rmvs_peri.f90 | 87 ----- src/rmvs/rmvs_setup.f90 | 128 +++--- src/rmvs/rmvs_step.f90 | 603 ++++++++++++++++++++--------- src/setup/setup.f90 | 8 +- src/step/step.f90 | 53 --- src/symba/symba_helio_drift.f90 | 2 +- src/symba/symba_helio_drift_tp.f90 | 3 +- src/whm/whm_drift.f90 | 28 +- src/whm/whm_getacch.f90 | 19 +- src/whm/whm_setup.f90 | 6 +- src/whm/whm_step.f90 | 22 +- 25 files changed, 745 insertions(+), 853 deletions(-) delete mode 100644 src/rmvs/rmvs_interp.f90 delete mode 100644 src/rmvs/rmvs_obl.f90 delete mode 100644 src/rmvs/rmvs_peri.f90 delete mode 100644 src/step/step.f90 diff --git a/Makefile b/Makefile index 86173c6bf..9bda6503a 100644 --- a/Makefile +++ b/Makefile @@ -58,14 +58,12 @@ SWIFTEST_MODULES = swiftest_globals.f90 \ include Makefile.Defines MODULES = $(SWIFTEST_MODULES) $(USER_MODULES) -IADVIXE = ${ADVISOR_2019_DIR}/include/intel64 -LADVIXE = ${ADVISOR_2019_DIR}/lib64 .PHONY : all mod lib libdir drivers bin clean force % : %.f90 force - $(FORTRAN) $(FFLAGS) -I$(SWIFTEST_HOME)/include -I$(IADVIXE) $< -o $@ \ - -L$(SWIFTEST_HOME)/lib -L$(LADVIXE) -lswiftest -ladvisor + $(FORTRAN) $(FFLAGS) -I$(SWIFTEST_HOME)/include $< -o $@ \ + -L$(SWIFTEST_HOME)/lib -lswiftest $(INSTALL_PROGRAM) $@ $(SWIFTEST_HOME)/bin rm -f $@ @@ -77,7 +75,7 @@ all: mod: cd $(SWIFTEST_HOME)/src/modules/; \ - $(FORTRAN) $(FFLAGS) -I$(SWIFTEST_HOME)/include -I$(IADVIXE) -c $(MODULES); \ + $(FORTRAN) $(FFLAGS) -I$(SWIFTEST_HOME)/include -c $(MODULES); \ $(AR) rv $(SWIFTEST_HOME)/lib/libswiftest.a *.o; \ $(INSTALL_DATA) *.mod *.smod $(SWIFTEST_HOME)/include; \ rm -f *.o *.mod *.smod @@ -128,11 +126,6 @@ lib: ln -s $(SWIFTEST_HOME)/Makefile.Defines .; \ ln -s $(SWIFTEST_HOME)/Makefile .; \ make libdir - cd $(SWIFTEST_HOME)/src/step; \ - rm -f Makefile.Defines Makefile; \ - ln -s $(SWIFTEST_HOME)/Makefile.Defines .; \ - ln -s $(SWIFTEST_HOME)/Makefile .; \ - make libdir cd $(SWIFTEST_HOME)/src/util; \ rm -f Makefile.Defines Makefile; \ ln -s $(SWIFTEST_HOME)/Makefile.Defines .; \ @@ -160,7 +153,7 @@ lib: make libdir libdir: - $(FORTRAN) $(FFLAGS) -I$(SWIFTEST_HOME)/include -I$(IADVIXE) -c *.f90; \ + $(FORTRAN) $(FFLAGS) -I$(SWIFTEST_HOME)/include -c *.f90; \ $(AR) rv $(SWIFTEST_HOME)/lib/libswiftest.a *.o *.smod; \ $(INSTALL_DATA) *.smod $(SWIFTEST_HOME)/include; \ rm -f *.o *.smod @@ -194,7 +187,6 @@ clean: cd $(SWIFTEST_HOME)/src/operators; rm -f Makefile.Defines Makefile *.gc* cd $(SWIFTEST_HOME)/src/orbel; rm -f Makefile.Defines Makefile *.gc* cd $(SWIFTEST_HOME)/src/setup; rm -f Makefile.Defines Makefile *.gc* - cd $(SWIFTEST_HOME)/src/step; rm -f Makefile.Defines Makefile *.gc* cd $(SWIFTEST_HOME)/src/util; rm -f Makefile.Defines Makefile *.gc* cd $(SWIFTEST_HOME)/src/user; rm -f Makefile.Defines Makefile *.gc* cd $(SWIFTEST_HOME)/src/main; rm -f Makefile.Defines Makefile *.gc* diff --git a/Makefile.Defines b/Makefile.Defines index 4ae84cd22..b9dc143fc 100644 --- a/Makefile.Defines +++ b/Makefile.Defines @@ -71,6 +71,8 @@ GPAR = -fopenmp -ftree-parallelize-loops=4 GMEM = -fsanitize=undefined -fsanitize=address -fsanitize=leak GWARNINGS = -Wall -Warray-bounds -Wimplicit-interface -Wextra -Warray-temporaries +FORTRAN = gfortran +FFLAGS = $(GDEBUG) # DO NOT include in CFLAGS the "-c" option to compile object only diff --git a/src/driftkick/drift.f90 b/src/driftkick/drift.f90 index 5998cec56..68d9d2822 100644 --- a/src/driftkick/drift.f90 +++ b/src/driftkick/drift.f90 @@ -9,7 +9,8 @@ integer(I2B), parameter :: NLAG2 = 40 contains - module pure subroutine drift_one(mu, x, v, dt, iflag) + module pure subroutine drift_one(mu, px, py, pz, vx, vy, vz, dt, iflag) + !$omp declare simd !! author: The Purdue Swiftest Team - David A. Minton, Carlisle A. Wishard, Jennifer L.L. Pouplin, and Jacob R. Elliott !! !! Perform Danby drift for one body, redoing drift with smaller substeps if original accuracy is insufficient @@ -19,12 +20,16 @@ module pure subroutine drift_one(mu, x, v, dt, iflag) implicit none ! Arguments real(DP), intent(in) :: mu !! G * (Mcb + m), G = gravitational constant, Mcb = mass of central body, m = mass of body to drift - real(DP), dimension(:), intent(inout) :: x, v !! Position and velocity of body to drift + real(DP), intent(inout) :: px, py, pz, vx, vy, vz !! Position and velocity of body to drift real(DP), intent(in) :: dt !! Step size integer(I4B), intent(out) :: iflag !! iflag : error status flag for Danby drift (0 = OK, nonzero = ERROR) ! Internals integer(I4B) :: i real(DP) :: dttmp + real(DP), dimension(NDIM) :: x, v + + x = [px, py, pz] + v = [vx, vy, vz] call drift_dan(mu, x(:), v(:), dt, iflag) if (iflag /= 0) then @@ -34,12 +39,14 @@ module pure subroutine drift_one(mu, x, v, dt, iflag) if (iflag /= 0) return end do end if + px = x(1); py = x(2); pz = x(3) + vx = v(1); vy = v(2); vz = v(3) return end subroutine drift_one pure subroutine drift_dan(mu, x0, v0, dt0, iflag) - !$omp declare simd(drift_dan) + !$omp declare simd !! author: David A. Minton !! !! Perform Kepler drift, solving Kepler's equation in appropriate variables @@ -109,7 +116,7 @@ pure subroutine drift_dan(mu, x0, v0, dt0, iflag) end subroutine drift_dan pure subroutine drift_kepmd(dm, es, ec, x, s, c) - !$omp declare simd(drift_kepmd) + !$omp declare simd !! author: David A. Minton !! !! Solve Kepler's equation in difference form for an ellipse for small input dm and eccentricity @@ -148,7 +155,7 @@ pure subroutine drift_kepmd(dm, es, ec, x, s, c) end subroutine drift_kepmd pure subroutine drift_kepu(dt,r0,mu,alpha,u,fp,c1,c2,c3,iflag) - !$omp declare simd(drift_kepu) + !$omp declare simd !! author: David A. Minton !! !! Solve Kepler's equation in universal variables @@ -176,7 +183,7 @@ pure subroutine drift_kepu(dt,r0,mu,alpha,u,fp,c1,c2,c3,iflag) end subroutine drift_kepu pure subroutine drift_kepu_fchk(dt, r0, mu, alpha, u, s, f) - !$omp declare simd(drift_kepu_fchk) + !$omp declare simd !! author: David A. Minton !! !! Computes the value of f, the function whose root we are trying to find in universal variables @@ -199,7 +206,7 @@ pure subroutine drift_kepu_fchk(dt, r0, mu, alpha, u, s, f) end subroutine drift_kepu_fchk pure subroutine drift_kepu_guess(dt, r0, mu, alpha, u, s) - !$omp declare simd(drift_kepu_guess) + !$omp declare simd !! author: David A. Minton !! !! Compute initial guess for solving Kepler's equation using universal variables @@ -237,7 +244,7 @@ pure subroutine drift_kepu_guess(dt, r0, mu, alpha, u, s) end subroutine drift_kepu_guess pure subroutine drift_kepu_lag(s, dt, r0, mu, alpha, u, fp, c1, c2, c3, iflag) - !$omp declare simd(drift_kepu_lag) + !$omp declare simd !! author: David A. Minton !! !! Solve Kepler's equation in universal variables using Laguerre's method @@ -282,7 +289,7 @@ pure subroutine drift_kepu_lag(s, dt, r0, mu, alpha, u, fp, c1, c2, c3, iflag) end subroutine drift_kepu_lag pure subroutine drift_kepu_new(s, dt, r0, mu, alpha, u, fp, c1, c2, c3, iflag) - !$omp declare simd(drift_kepu_new) + !$omp declare simd !! author: David A. Minton !! !! Solve Kepler's equation in universal variables using Newton's method @@ -324,7 +331,7 @@ pure subroutine drift_kepu_new(s, dt, r0, mu, alpha, u, fp, c1, c2, c3, iflag) end subroutine drift_kepu_new pure subroutine drift_kepu_p3solve(dt, r0, mu, alpha, u, s, iflag) - !$omp declare simd(drift_kepu_p3solve) + !$omp declare simd !! author: David A. Minton !! !! Computes real root of cubic involved in setting initial guess for solving Kepler's equation in universal variables @@ -368,7 +375,7 @@ pure subroutine drift_kepu_p3solve(dt, r0, mu, alpha, u, s, iflag) end subroutine drift_kepu_p3solve pure subroutine drift_kepu_stumpff(x, c0, c1, c2, c3) - !$omp declare simd(drift_kepu_stumpff) + !$omp declare simd !! author: David A. Minton !! !! Compute Stumpff functions needed for Kepler drift in universal variables diff --git a/src/driftkick/kick.f90 b/src/driftkick/kick.f90 index 8c428f76c..95e34a1d9 100644 --- a/src/driftkick/kick.f90 +++ b/src/driftkick/kick.f90 @@ -17,7 +17,7 @@ module subroutine kick_vh_body(self, dt) associate(n => self%nbody, vh => self%vh, ah => self%ah, status => self%status) if (n == 0) return - !$omp simd + !!$omp simd do i = 1, n if (status(i) == ACTIVE) vh(:, i) = vh(:, i) + ah(:, i) * dt end do diff --git a/src/helio/helio_drift.f90 b/src/helio/helio_drift.f90 index 244dd5f9e..522b3403b 100644 --- a/src/helio/helio_drift.f90 +++ b/src/helio/helio_drift.f90 @@ -40,7 +40,7 @@ module subroutine helio_drift_pl(self, cb, config, dt) else dtp = dt end if - call drift_one(mu(i), xh(:, i), vb(:, i), dtp, iflag(i)) + call drift_one(mu(i), xh(1, i), xh(2, i), xh(3, i), vb(1, i), vb(2, i), vb(3, i), dtp, iflag(i)) end do if (any(iflag(1:npl) /= 0)) then do i = 1, npl @@ -125,7 +125,7 @@ module subroutine helio_drift_tp(self, cb, config, dt) else dtp = dt end if - call drift_one(mu(i), xh(:, i), vb(:, i), dtp, iflag(i)) + call drift_one(mu(i), xh(1, i), xh(2, i), xh(3, i), vb(1, i), vb(2, i), vb(3, i), dtp, iflag(i)) if (iflag(i) /= 0) status = DISCARDED_DRIFTERR end if end do diff --git a/src/helio/helio_step.f90 b/src/helio/helio_step.f90 index 1f8a6d96a..d31781e45 100644 --- a/src/helio/helio_step.f90 +++ b/src/helio/helio_step.f90 @@ -1,7 +1,7 @@ submodule(helio_classes) s_helio_step use swiftest contains - module subroutine helio_step_system(cb, pl, tp, config) + module subroutine helio_step_system(self, config) !! author: David A. Minton !! !! Step massive bodies and and active test particles ahead in heliocentric coordinates @@ -10,11 +10,15 @@ module subroutine helio_step_system(cb, pl, tp, config) !! Adapted from David E. Kaufmann's Swifter routine helio_step.f90 implicit none ! Arguments - class(helio_cb), intent(inout) :: cb !! Helio central body object - class(helio_pl), intent(inout) :: pl !! Helio central body object - class(helio_tp), intent(inout) :: tp !! Helio central body object + class(helio_nbody_system), intent(inout) :: self !! Helio nbody system object class(swiftest_configuration), intent(in) :: config !! Input collection of on parameters + select type(cb => self%cb) + class is (helio_cb) + select type(pl => self%pl) + class is (helio_pl) + select type(tp => self%tp) + class is (helio_tp) associate(ntp => tp%nbody, npl => pl%nbody, t => config%t, dt => config%dt) call pl%set_rhill(cb) call tp%set_beg_end(xbeg = pl%xh) @@ -24,6 +28,10 @@ module subroutine helio_step_system(cb, pl, tp, config) call tp%step(cb, pl, config, t, dt) end if end associate + end select + end select + end select + return end subroutine helio_step_system module subroutine helio_step_pl(self, cb, config, t, dt) diff --git a/src/modules/helio_classes.f90 b/src/modules/helio_classes.f90 index 3e49e768b..ee8359efc 100644 --- a/src/modules/helio_classes.f90 +++ b/src/modules/helio_classes.f90 @@ -7,6 +7,16 @@ module helio_classes use rmvs_classes, only : rmvs_cb, rmvs_pl, rmvs_tp, rmvs_nbody_system implicit none + + !******************************************************************************************************************************** + ! helio_nbody_system class definitions and method interfaces + !******************************************************************************************************************************** + type, public, extends(rmvs_nbody_system) :: helio_nbody_system + contains + private + procedure, public :: step => helio_step_system + end type helio_nbody_system + !******************************************************************************************************************************** ! helio_cb class definitions and method interfaces !******************************************************************************************************************************* @@ -53,13 +63,6 @@ module helio_classes procedure, public :: step => helio_step_tp !! Steps the body forward one stepsize end type helio_tp - !******************************************************************************************************************************** - ! helio_nbody_system class definitions and method interfaces - !******************************************************************************************************************************** - type, public, extends(rmvs_nbody_system) :: helio_nbody_system - contains - private - end type helio_nbody_system interface @@ -171,19 +174,16 @@ module subroutine helio_setup_tp(self,n) integer, intent(in) :: n !! Number of test particles to allocate end subroutine helio_setup_tp - module subroutine helio_step_system(cb, pl, tp, config) - use swiftest_classes + module subroutine helio_step_system(self, config) + use swiftest_classes, only : swiftest_configuration implicit none - class(helio_cb), intent(inout) :: cb !! Helio central body object - class(helio_pl), intent(inout) :: pl !! Helio massive body object - class(helio_tp), intent(inout) :: tp !! Helio test particle object + class(helio_nbody_system), intent(inout) :: self !! Helio nbody system object class(swiftest_configuration), intent(in) :: config !! Input collection of configuration parameters end subroutine helio_step_system module subroutine helio_step_pl(self, cb, config, t, dt) use swiftest_classes, only : swiftest_cb, swiftest_configuration implicit none - ! Arguments class(helio_pl), intent(inout) :: self !! WHM massive body particle data structure class(swiftest_cb), intent(inout) :: cb !! Swiftest central body particle data structure class(swiftest_configuration), intent(in) :: config !! Input collection of diff --git a/src/modules/rmvs_classes.f90 b/src/modules/rmvs_classes.f90 index fe869be63..9a3f32bd0 100644 --- a/src/modules/rmvs_classes.f90 +++ b/src/modules/rmvs_classes.f90 @@ -16,14 +16,34 @@ module rmvs_classes real(DP), parameter :: RHPSCALE = 1.0_DP real(DP), parameter :: FACQDT = 2.0_DP + type, private :: rmvs_interp + real(DP), dimension(:, :), allocatable :: x !! interpolated heliocentric planet position for outer encounter + real(DP), dimension(:, :), allocatable :: v !! interpolated heliocentric planet velocity for outer encounter + real(DP), dimension(:, :), allocatable :: aobl !! Encountering planet's oblateness acceleration value + end type rmvs_interp + + + !******************************************************************************************************************************** + ! rmvs_nbody_system class definitions and method interfaces + !******************************************************************************************************************************** + type, public, extends(whm_nbody_system) :: rmvs_nbody_system + !> In the RMVS integrator, only test particles are discarded + logical :: lplanetocentric = .false. !! Flag that indicates that the object is a planetocentric set of masive bodies used for close encounter calculations + contains + private + !> Replace the abstract procedures with concrete ones + procedure, public :: initialize => rmvs_setup_system !! Performs RMVS-specific initilization steps, like calculating the Jacobi masses + procedure, public :: step => rmvs_step_system + end type rmvs_nbody_system + !******************************************************************************************************************************** ! rmvs_cb class definitions and method interfaces !******************************************************************************************************************************* - !> RMVS central body particle class + !> RMVS central body particle class type, public, extends(whm_cb) :: rmvs_cb - real(DP), dimension(NDIM) :: xin = (/0.0_DP, 0.0_DP, 0.0_DP/) !! interpolated heliocentric central body position for inner encounter - real(DP), dimension(NDIM) :: vin = (/0.0_DP, 0.0_DP, 0.0_DP/) !! interpolated heliocentric central body velocity for inner encounter - contains + logical :: lplanetocentric !! Flag that indicates that the object is a planetocentric set of test particles used for close encounter calculations + real(DP), dimension(NDIM) :: xin = [0.0_DP, 0.0_DP, 0.0_DP] + real(DP), dimension(NDIM) :: vin = [0.0_DP, 0.0_DP, 0.0_DP] end type rmvs_cb !******************************************************************************************************************************** @@ -35,26 +55,23 @@ module rmvs_classes !! Note to developers: If you add componenets to this class, be sure to update methods and subroutines that traverse the !! component list, such as rmvs_setup_tp and rmvs_spill_tp ! encounter steps) - logical :: lplanetocentric = .false. !! Flag that indicates that the object is a planetocentric set of test particles used for close encounter calculations + logical :: lplanetocentric !! Flag that indicates that the object is a planetocentric set of test particles used for close encounter calculations logical, dimension(:), allocatable :: lperi !! planetocentric pericenter passage flag (persistent for a full rmvs time step) over a full RMVS time step) integer(I4B), dimension(:), allocatable :: plperP !! index of planet associated with pericenter distance peri (persistent over a full RMVS time step) integer(I4B), dimension(:), allocatable :: plencP !! index of planet that test particle is encountering (not persistent for a full RMVS time step) real(DP), dimension(:,:), allocatable :: vbeg !! Planet velocities at beginning ot step ! The following are used to correctly set the oblateness values of the acceleration during an inner encounter with a planet - type(rmvs_cb) :: cb !! Copy of original central body object passed to close encounter (used for oblateness acceleration during planetocentric encoountters) - real(DP), dimension(:,:), allocatable :: xheliocen !! original heliocentric position (used for oblateness calculation during close encounters) - integer(I4B) :: index !! inner substep number within current set - integer(I4B) :: ipleP !! index value of encountering planet - real(DP), dimension(:,:), allocatable :: aoblin_pl !! Encountering planet's oblateness acceleration value - real(DP), dimension(:,:), allocatable :: xh_pl !! Encountering planet's heliocentric position values + type(rmvs_cb) :: cb_heliocentric !! Copy of original central body object passed to close encounter (used for oblateness acceleration during planetocentric encoountters) + real(DP), dimension(:,:), allocatable :: xheliocentric !! original heliocentric position (used for oblateness calculation during close encounters) + integer(I4B) :: index !! inner substep number within current set + integer(I4B) :: ipleP !! index value of encountering planet contains procedure, public :: discard_pl => rmvs_discard_pl_tp procedure, public :: encounter_check => rmvs_encounter_check_tp !! Checks if any test particles are undergoing a close encounter with a massive body procedure, public :: fill => rmvs_fill_tp !! "Fills" bodies from one object into another depending on the results of a mask (uses the MERGE intrinsic) - procedure, public :: getacch => rmvs_getacch_in_tp !! Calculates either the standard or modified version of the oblateness acceleration depending if the + procedure, public :: getacch => rmvs_getacch_tp !! Calculates either the standard or modified version of the acceleration depending if the !! if the test particle is undergoing a close encounter or not - procedure, public :: peri_pass => rmvs_peri_tp !! Determine planetocentric pericenter passages for test particles in close encounters with a planet procedure, public :: set_beg_end => rmvs_setup_set_beg_end !! Sets the beginning and ending values of planet positions. Also adds the end velocity for RMVS procedure, public :: setup => rmvs_setup_tp !! Constructor method - Allocates space for number of particles procedure, public :: spill => rmvs_spill_tp !! "Spills" bodies from one object to another depending on the results of a mask (uses the PACK intrinsic) @@ -65,43 +82,30 @@ module rmvs_classes !******************************************************************************************************************************* !> RMVS massive body particle class - type, public, extends(whm_pl) :: rmvs_pl - integer(I4B), dimension(:), allocatable :: nenc !! number of test particles encountering planet this full rmvs time step - integer(I4B), dimension(:), allocatable :: tpenc1P !! index of first test particle encountering planet - real(DP), dimension(:, :, :), allocatable :: xout !! interpolated heliocentric planet position for outer encounter - real(DP), dimension(:, :, :), allocatable :: vout !! interpolated heliocentric planet velocity for outer encounter - real(DP), dimension(:, :, :), allocatable :: xin !! interpolated heliocentric planet position for inner encounter - real(DP), dimension(:, :, :), allocatable :: vin !! interpolated heliocentric planet velocity for inner encounter - type(rmvs_tp), dimension(:), allocatable :: tpenc !! array of encountering test particles with this planet in planetocentric coordinates - type(whm_pl) , dimension(:), allocatable :: plenc !! array of massive bodies that includes the Sun, but not the encountering planet in planetocentric coordinates - type(rmvs_cb), dimension(:), allocatable :: cbenc !! The planet acting as a central body for close encounters - real(DP), dimension(:, :, :), allocatable :: aoblin !! barycentric acceleration on planets due to central body oblateness during inner encounter - integer(I4B), dimension(:,:), allocatable :: plind + type, private, extends(whm_pl) :: rmvs_base_pl + logical :: lplanetocentric !! Flag that indicates that the object is a planetocentric set of masive bodies used for close encounter calculations + integer(I4B), dimension(:), allocatable :: plind ! Connects the planetocentric indices back to the heliocentric planet list + end type rmvs_base_pl + + type, public :: rmvs_encounter_system + class(rmvs_cb), allocatable :: cb + class(rmvs_tp), allocatable :: tp + class(rmvs_base_pl), allocatable :: pl + end type rmvs_encounter_system + + type, public, extends(rmvs_base_pl) :: rmvs_pl + integer(I4B), dimension(:), allocatable :: nenc !! number of test particles encountering planet this full rmvs time step + integer(I4B), dimension(:), allocatable :: tpenc1P !! index of first test particle encountering planet + type(rmvs_encounter_system), dimension(:), allocatable :: planetocentric + type(rmvs_interp), dimension(:), allocatable :: outer !! interpolated heliocentric central body position for outer encounters + type(rmvs_interp), dimension(:), allocatable :: inner !! interpolated heliocentric central body position for inner encounters !! Note to developers: If you add componenets to this class, be sure to update methods and subroutines that traverse the !! component list, such as rmvs_setup_pl and rmvs_spill_pl contains procedure, public :: fill => rmvs_fill_pl !! "Fills" bodies from one object into another depending on the results of a mask (uses the MERGE intrinsic) - procedure, public :: interp_in => rmvs_interp_in !! Interpolate planet positions between two Keplerian orbits in inner encounter region - procedure, public :: interp_out => rmvs_interp_out !! Interpolate planet positions between two Keplerian orbits in outer encounter region - procedure, public :: obl_acc_in => rmvs_obl_acc_in !! Compute the oblateness acceleration in the inner encounter region with planets procedure, public :: setup => rmvs_setup_pl !! Constructor method - Allocates space for number of particles - procedure, public :: step_in => rmvs_step_in_pl !! Step active test particles ahead in the inner encounter region with planets - procedure, public :: step_out => rmvs_step_out !! Step active test particles ahead in the outer encounter region with planets - procedure, public :: make_planetocentric => rmvs_step_make_planetocentric !! Creates encountering test particle structure for the planets - procedure, public :: end_planetocentric => rmvs_step_end_planetocentric !! Creates encountering test particle structure for the planets procedure, public :: spill => rmvs_spill_pl !! "Spills" bodies from one object to another depending on the results of a mask (uses the PACK intrinsic) end type rmvs_pl - - !******************************************************************************************************************************** - ! rmvs_nbody_system class definitions and method interfaces - !******************************************************************************************************************************** - type, public, extends(whm_nbody_system) :: rmvs_nbody_system - !> In the RMVS integrator, only test particles are discarded - contains - private - !> Replace the abstract procedures with concrete ones - procedure, public :: initialize => rmvs_setup_system !! Performs RMVS-specific initilization steps, like calculating the Jacobi masses - end type rmvs_nbody_system interface module subroutine rmvs_discard_pl_tp(self, pl, t, dt) @@ -113,7 +117,7 @@ module subroutine rmvs_discard_pl_tp(self, pl, t, dt) real(DP), intent(in) :: dt !! Stepsize end subroutine rmvs_discard_pl_tp - module subroutine rmvs_getacch_in_tp(self, cb, pl, config, t, xh) + module subroutine rmvs_getacch_tp(self, cb, pl, config, t, xh) use swiftest_classes, only : swiftest_cb, swiftest_configuration use whm_classes, only : whm_pl implicit none @@ -123,14 +127,7 @@ module subroutine rmvs_getacch_in_tp(self, cb, pl, config, t, xh) class(swiftest_configuration), intent(in) :: config !! Input collection of parameter real(DP), intent(in) :: t !! Current time real(DP), dimension(:,:), intent(in) :: xh !! Heliocentric positions of planets - end subroutine rmvs_getacch_in_tp - - module subroutine rmvs_obl_acc_in(self, cb) - use swiftest_classes, only : swiftest_cb - implicit none - class(rmvs_pl), intent(inout) :: self !! RMVS massive body object - class(swiftest_cb), intent(inout) :: cb !! Swiftest central body object - end subroutine rmvs_obl_acc_in + end subroutine rmvs_getacch_tp module subroutine rmvs_setup_system(self, config) use swiftest_classes, only : swiftest_configuration @@ -141,96 +138,20 @@ end subroutine rmvs_setup_system module subroutine rmvs_setup_pl(self,n) implicit none - class(rmvs_pl), intent(inout) :: self !! RMVS test particle object - integer, intent(in) :: n !! Number of test particles to allocate + class(rmvs_pl), intent(inout) :: self !! RMVS test particle object + integer, intent(in) :: n !! Number of test particles to allocate end subroutine rmvs_setup_pl - module subroutine rmvs_step_make_planetocentric(self, cb, tp, config) - use swiftest_classes, only : swiftest_configuration - implicit none - class(rmvs_pl), intent(inout) :: self !! RMVS test particle object - class(rmvs_cb), intent(inout) :: cb !! RMVS central body particle type - class(rmvs_tp), intent(inout) :: tp !! RMVS test particle object - class(swiftest_configuration), intent(in) :: config !! Input collection of configuration parameters - end subroutine rmvs_step_make_planetocentric - - module subroutine rmvs_step_end_planetocentric(self, cb, tp) - implicit none - class(rmvs_pl), intent(inout) :: self !! RMVS test particle object - class(rmvs_cb), intent(inout) :: cb !! RMVS central body object - class(rmvs_tp), intent(inout) :: tp !! RMVS test particle object - end subroutine rmvs_step_end_planetocentric - module subroutine rmvs_setup_set_beg_end(self, xbeg, xend, vbeg) implicit none - class(rmvs_tp), intent(inout) :: self !! RMVS test particle object - real(DP), dimension(:,:), optional :: xbeg, xend, vbeg + class(rmvs_tp), intent(inout) :: self !! RMVS test particle object + real(DP), dimension(:,:), intent(in), optional :: xbeg, xend, vbeg end subroutine rmvs_setup_set_beg_end - module subroutine rmvs_destruct_planetocentric_encounter(self) - implicit none - class(rmvs_pl), intent(inout) :: self !! RMVS test particle object - end subroutine rmvs_destruct_planetocentric_encounter - - module subroutine rmvs_interp_in(self, cb, dt) - implicit none - class(rmvs_pl), intent(inout) :: self !! RMVS test particle object - class(rmvs_cb), intent(in) :: cb !! RMVS central body particle type - real(DP), intent(in) :: dt !! Step size - end subroutine rmvs_interp_in - - module subroutine rmvs_interp_out(self, cb, dt) - implicit none - class(rmvs_pl), intent(inout) :: self !! RMVS test particle object - class(rmvs_cb), intent(in) :: cb !! RMVS central body particle type - real(DP), intent(in) :: dt !! Step size - end subroutine rmvs_interp_out - - module subroutine rmvs_step_in_pl(self, cb, tp, config, t, dt) - use swiftest_classes, only : swiftest_configuration - implicit none - class(rmvs_pl), intent(inout) :: self !! RMVS massive body object - class(rmvs_cb), intent(inout) :: cb !! RMVS central body object - class(rmvs_tp), intent(inout) :: tp !! RMVS test particle object - class(swiftest_configuration), intent(in) :: config !! Input collection of configuration parameters - real(DP), intent(in) :: t !! Current time - real(DP), intent(in) :: dt !! Step size - end subroutine rmvs_step_in_pl - - module subroutine rmvs_step_out(self, cb, tp, dt, config) + module subroutine rmvs_step_system(self, config) use swiftest_classes, only : swiftest_configuration implicit none - class(rmvs_pl), intent(inout) :: self !! RMVS massive body object - class(rmvs_cb), intent(inout) :: cb !! RMVS central body object - class(rmvs_tp), intent(inout) :: tp !! RMVS test particle object - real(DP), intent(in) :: dt !! Step size - class(swiftest_configuration), intent(in) :: config !! Input collection of configuration parameters - end subroutine rmvs_step_out - - module subroutine rmvs_step_out2(cb, pl, tp, t, dt, index, config) - use swiftest_classes, only : swiftest_configuration - implicit none - class(rmvs_cb), intent(inout) :: cb !! RMVS central body object - class(rmvs_pl), intent(inout) :: pl !! RMVS massive body object - class(rmvs_tp), intent(inout) :: tp !! RMVS test particle object - real(DP), intent(in) :: t !! Simulation time - real(DP), intent(in) :: dt !! Step size - integer(I4B), intent(in) :: index !! outer substep number within current set - class(swiftest_configuration), intent(in) :: config !! Input collection of configuration parameters - end subroutine rmvs_step_out2 - - module elemental function rmvs_chk_ind(r2, v2, vdotr, dt, r2crit) result(lflag) - implicit none - real(DP), intent(in) :: r2, v2, vdotr, dt, r2crit - logical :: lflag - end function rmvs_chk_ind - - module subroutine rmvs_step_system(cb, pl, tp, config) - use swiftest_classes, only : swiftest_configuration - implicit none - class(rmvs_cb), intent(inout) :: cb !! RMVS central body object - class(rmvs_pl), intent(inout) :: pl !! RMVS central body object - class(rmvs_tp), intent(inout) :: tp !! RMVS central body object + class(rmvs_nbody_system), intent(inout) :: self !! RMVS nbody system object class(swiftest_configuration), intent(in) :: config !! Input collection of configuration parameters end subroutine rmvs_step_system @@ -242,40 +163,26 @@ end subroutine rmvs_setup_tp module function rmvs_encounter_check_tp(self, cb, pl, dt, rts) result(lencounter) implicit none - class(rmvs_tp), intent(inout) :: self !! RMVS test particle object - class(rmvs_cb), intent(inout) :: cb !! RMVS central body object - class(rmvs_pl), intent(inout) :: pl !! RMVS massive body object - real(DP), intent(in) :: dt !! step size - real(DP), intent(in) :: rts !! fraction of Hill's sphere radius to use as radius of encounter regio - logical :: lencounter !! Returns true if there is at least one close encounter - + class(rmvs_tp), intent(inout) :: self !! RMVS test particle object + class(rmvs_cb), intent(inout) :: cb !! RMVS central body object + class(rmvs_pl), intent(inout) :: pl !! RMVS massive body object + real(DP), intent(in) :: dt !! step size + real(DP), intent(in) :: rts !! fraction of Hill's sphere radius to use as radius of encounter regio + logical :: lencounter !! Returns true if there is at least one close encounter end function rmvs_encounter_check_tp - module subroutine rmvs_peri_tp(self, cb, pl, t, dt, lfirst, index, ipleP, config) - use swiftest_classes, only : swiftest_configuration - class(rmvs_tp), intent(inout) :: self !! RMVS test particle object - class(rmvs_cb), intent(inout) :: cb !! RMVS central body object - class(rmvs_pl), intent(inout) :: pl !! RMVS massive body object - real(DP), intent(in) :: t !! current time - real(DP), intent(in) :: dt !! step size - logical, intent(in) :: lfirst !! Logical flag indicating whether current invocation is the first - integer(I4B), intent(in) :: index !! outer substep number within current set - integer(I4B), intent(in) :: ipleP !!index of RMVS planet being closely encountered - class(swiftest_configuration), intent(in) :: config !! Input collection of configuration parameters - end subroutine rmvs_peri_tp - module subroutine rmvs_spill_pl(self, discards, lspill_list) use swiftest_classes, only : swiftest_body implicit none - class(rmvs_pl), intent(inout) :: self !! Swiftest massive body body object - class(swiftest_body), intent(inout) :: discards !! Discarded object + class(rmvs_pl), intent(inout) :: self !! RMVS massive body object + class(swiftest_body), intent(inout) :: discards !! Discarded object logical, dimension(:), intent(in) :: lspill_list !! Logical array of bodies to spill into the discards end subroutine rmvs_spill_pl module subroutine rmvs_fill_pl(self, inserts, lfill_list) use swiftest_classes, only : swiftest_body implicit none - class(rmvs_pl), intent(inout) :: self !! RMVS massive body object + class(rmvs_pl), intent(inout) :: self !! RMVS massive body object class(swiftest_body), intent(inout) :: inserts !! Inserted object logical, dimension(:), intent(in) :: lfill_list !! Logical array of bodies to merge into the keeps end subroutine rmvs_fill_pl diff --git a/src/modules/swiftest_classes.f90 b/src/modules/swiftest_classes.f90 index 342e886b7..b0e7c4b67 100644 --- a/src/modules/swiftest_classes.f90 +++ b/src/modules/swiftest_classes.f90 @@ -247,6 +247,7 @@ module swiftest_classes contains private !> Each integrator will have its own version of the step + procedure(abstract_step_system), public, deferred :: step ! Concrete classes that are common to the basic integrator (only test particles considered for discard) procedure, public :: discard => discard_system !! Perform a discard step on the system @@ -256,7 +257,6 @@ module swiftest_classes procedure, public :: set_msys => setup_set_msys !! Sets the value of msys from the masses of system bodies. procedure, public :: write_discard => io_write_discard !! Append a frame of output data to file procedure, public :: write_frame => io_write_frame_system !! Append a frame of output data to file - procedure, public :: step => step_system procedure, public :: copy => util_copy_system !! Copies elements of one object to another. end type swiftest_nbody_system @@ -297,6 +297,13 @@ subroutine abstract_set_mu(self, cb) class(swiftest_cb), intent(inout) :: cb !! Swiftest central body objectt end subroutine abstract_set_mu + subroutine abstract_step_system(self, config) + import swiftest_nbody_system, swiftest_configuration + implicit none + class(swiftest_nbody_system), intent(inout) :: self !! Swiftest system object + class(swiftest_configuration), intent(in) :: config !! Input collection of configuration parameters + end subroutine abstract_step_system + subroutine abstract_write_frame(self, iu, config, t, dt) import DP, I4B, swiftest_base, swiftest_configuration class(swiftest_base), intent(in) :: self !! Swiftest base object @@ -349,13 +356,13 @@ module subroutine discard_system(self, config) class(swiftest_configuration), intent(in) :: config !! Input collection of configuration parameters end subroutine discard_system - module pure subroutine drift_one(mu, x, v, dt, iflag) + module pure subroutine drift_one(mu, px, py, pz, vx, vy, vz, dt, iflag) !$omp declare simd(drift_one) implicit none - real(DP), intent(in) :: mu !! G * (Mcb + m), G = gravitational constant, Mcb = mass of central body, m = mass of body to drift - real(DP), dimension(:), intent(inout) :: x, v !! Position and velocity of body to drift - real(DP), intent(in) :: dt !! Step size - integer(I4B), intent(out) :: iflag !! iflag : error status flag for Danby drift (0 = OK, nonzero = ERROR) + real(DP), intent(in) :: mu !! G * (Mcb + m), G = gravitational constant, Mcb = mass of central body, m = mass of body to drift + real(DP), intent(inout) :: px, py, pz, vx, vy, vz !! Position and velocity of body to drift + real(DP), intent(in) :: dt !! Step size + integer(I4B), intent(out) :: iflag !! iflag : error status flag for Danby drift (0 = OK, nonzero = ERROR) end subroutine drift_one module subroutine eucl_dist_index_plpl(self) @@ -671,13 +678,6 @@ module subroutine setup_tp(self, n) integer, intent(in) :: n !! Number of bodies to allocate space for end subroutine setup_tp - module subroutine step_system(self, config) - implicit none - class(swiftest_nbody_system), intent(inout) :: self !! Swiftest system object - class(swiftest_configuration), intent(in) :: config !! Input collection of configuration parameters - end subroutine step_system - - module subroutine user_getacch_body(self, cb, config, t) implicit none class(swiftest_body), intent(inout) :: self !! Swiftest massive body particle data structure diff --git a/src/modules/whm_classes.f90 b/src/modules/whm_classes.f90 index ed2ce9676..7f1508ae5 100644 --- a/src/modules/whm_classes.f90 +++ b/src/modules/whm_classes.f90 @@ -9,6 +9,20 @@ module whm_classes implicit none public + !******************************************************************************************************************************** + ! whm_nbody_system class definitions and method interfaces + !******************************************************************************************************************************** + !> An abstract class for the WHM integrator nbody system + type, public, extends(swiftest_nbody_system) :: whm_nbody_system + !> In the WHM integrator, only test particles are discarded + class(swiftest_tp), allocatable :: tp_discards + contains + private + !> Replace the abstract procedures with concrete ones + procedure, public :: initialize => whm_setup_system !! Performs WHM-specific initilization steps, like calculating the Jacobi masses + procedure, public :: step => whm_step_system + end type whm_nbody_system + !******************************************************************************************************************************** ! whm_cb class definitions and method interfaces !******************************************************************************************************************************* @@ -73,18 +87,7 @@ module whm_classes procedure, public :: step => whm_step_tp !! Steps the particle forward one stepsize end type whm_tp - !******************************************************************************************************************************** - ! whm_nbody_system class definitions and method interfaces - !******************************************************************************************************************************** - !> An abstract class for the WHM integrator nbody system - type, public, extends(swiftest_nbody_system) :: whm_nbody_system - !> In the WHM integrator, only test particles are discarded - class(swiftest_tp), allocatable :: tp_discards - contains - private - !> Replace the abstract procedures with concrete ones - procedure, public :: initialize => whm_setup_system !! Performs WHM-specific initilization steps, like calculating the Jacobi masses - end type whm_nbody_system + interface !> WHM massive body constructor method @@ -234,9 +237,9 @@ end subroutine whm_step_tp module subroutine whm_setup_set_beg_end(self, xbeg, xend, vbeg) implicit none - class(whm_tp), intent(inout) :: self !! Swiftest test particle object - real(DP), dimension(:,:), optional :: xbeg, xend - real(DP), dimension(:,:), optional :: vbeg ! vbeg is an unused variable to keep this method forward compatible with RMVS + class(whm_tp), intent(inout) :: self !! Swiftest test particle object + real(DP), dimension(:,:), intent(in), optional :: xbeg, xend + real(DP), dimension(:,:), intent(in), optional :: vbeg ! vbeg is an unused variable to keep this method forward compatible with RMVS end subroutine whm_setup_set_beg_end module subroutine whm_gr_getacch_tp(self, cb, config) @@ -293,12 +296,10 @@ module subroutine whm_fill_pl(self, inserts, lfill_list) end subroutine whm_fill_pl !> Steps the Swiftest nbody system forward in time one stepsize - module subroutine whm_step_system(cb, pl, tp, config) + module subroutine whm_step_system(self, config) use swiftest_classes, only : swiftest_configuration implicit none - class(whm_cb), intent(inout) :: cb !! WHM central body object - class(whm_pl), intent(inout) :: pl !! WHM central body object - class(whm_tp), intent(inout) :: tp !! WHM central body object + class(whm_nbody_system), intent(inout) :: self !! WHM system object class(swiftest_configuration), intent(in) :: config !! Input collection of on parameters end subroutine whm_step_system end interface diff --git a/src/rmvs/rmvs_encounter_check.f90 b/src/rmvs/rmvs_encounter_check.f90 index 3eb1f406e..438693c41 100644 --- a/src/rmvs/rmvs_encounter_check.f90 +++ b/src/rmvs/rmvs_encounter_check.f90 @@ -12,7 +12,7 @@ module function rmvs_encounter_check_tp(self, cb, pl, dt, rts) result(lencounter ! Arguments class(rmvs_tp), intent(inout) :: self !! RMVS test particle object class(rmvs_cb), intent(inout) :: cb !! RMVS central body object - class(rmvs_pl), intent(inout) :: pl !! RMVS massive body object + class(rmvs_pl), intent(inout) :: pl !! RMVS massive body object real(DP), intent(in) :: dt !! step size real(DP), intent(in) :: rts !! fraction of Hill's sphere radius to use as radius of encounter regio logical :: lencounter !! Returns true if there is at least one close encounter @@ -46,7 +46,7 @@ module function rmvs_encounter_check_tp(self, cb, pl, dt, rts) result(lencounter return end function rmvs_encounter_check_tp - module elemental function rmvs_chk_ind(r2, v2, vdotr, dt, r2crit) result(lflag) + elemental function rmvs_chk_ind(r2, v2, vdotr, dt, r2crit) result(lflag) !! author: David A. Minton !! !! Determine whether a test particle and planet are having or will have an encounter within the next time step diff --git a/src/rmvs/rmvs_getacch.f90 b/src/rmvs/rmvs_getacch.f90 index 77a902660..9366404be 100644 --- a/src/rmvs/rmvs_getacch.f90 +++ b/src/rmvs/rmvs_getacch.f90 @@ -1,7 +1,7 @@ submodule(rmvs_classes) s_rmvs_getacch use swiftest contains - module subroutine rmvs_getacch_in_tp(self, cb, pl, config, t, xh) + module subroutine rmvs_getacch_tp(self, cb, pl, config, t, xh) !! author: David A. Minton !! !! Compute the oblateness acceleration in the inner encounter region with planets @@ -22,67 +22,65 @@ module subroutine rmvs_getacch_in_tp(self, cb, pl, config, t, xh) integer(I4B) :: i - associate(tp => self, cb_heliocen => self%cb, ipleP => self%ipleP, index => self%index, & - xht => self%xh, vht => self%vh, aht => self%ah) + associate(tp => self, ntp => self%nbody, ipleP => self%ipleP, & + enc_index => self%index, cb_heliocentric => self%cb_heliocentric) + if (tp%lplanetocentric) then ! This is a close encounter step, so any accelerations requiring heliocentric position values ! must be handeled outside the normal WHM method call - allocate(xh_original, source=tp%xh) - config_planetocen = config - ! Temporarily turn off the heliocentric-dependent acceleration terms during an inner encounter - config_planetocen%loblatecb = .false. - config_planetocen%lextra_force = .false. - config_planetocen%lgr = .false. - ! Now compute the heliocentric values of acceleration select type(pl) - class is (rmvs_pl) - if (tp%lfirst) then - do i = 1, NDIM - tp%xheliocen(i,:) = tp%xh(i,:) + pl%xin(i,ipleP,index - 1) + class is (rmvs_pl) + allocate(xh_original, source=tp%xh) + write(*,*) 'ipleP = ',ipleP + write(*,*) 'enc_index = ',enc_index + write(*,*) 'xin-1: ',pl%inner(enc_index - 1)%x(:, ipleP) + write(*,*) 'xin : ',pl%inner(enc_index )%x(:, ipleP) + write(*,*) 'xh :' + do i = 1, pl%nbody + write(*,*) i,pl%xh(:,i) end do - else - do i = 1, NDIM - tp%xheliocen(i,:) = tp%xh(i,:) + pl%xin(i,ipleP,index) - end do - end if - end select - call whm_getacch_tp(tp, cb, pl, config_planetocen, t, xh) - - ! Now compute the heliocentric values of acceleration - select type(pl) - class is (rmvs_pl) - if (tp%lfirst) then - do i = 1, NDIM - tp%xheliocen(i,:) = tp%xh(i,:) + pl%xin(i,ipleP,index - 1) - end do - else - do i = 1, NDIM - tp%xheliocen(i,:) = tp%xh(i,:) + pl%xin(i,ipleP,index) - end do - end if - end select + config_planetocen = config + ! Temporarily turn off the heliocentric-dependent acceleration terms during an inner encounter + config_planetocen%loblatecb = .false. + config_planetocen%lextra_force = .false. + config_planetocen%lgr = .false. + ! Now compute the planetocentric values of acceleration + call whm_getacch_tp(tp, cb, pl, config_planetocen, t, pl%xh) - tp%xh(:,:) = tp%xheliocen(:,:) - if (config%loblatecb) then - ! Put in the current encountering planet's oblateness acceleration as the central body's - if (tp%lfirst) then - cb_heliocen%aobl(:) = tp%aoblin_pl(:, index - 1) - else - cb_heliocen%aobl(:) = tp%aoblin_pl(:, index) - end if - call tp%obl_acc(cb_heliocen) - end if + ! Now compute any heliocentric values of acceleration + if (tp%lfirst) then + do i = 1, ntp + tp%xheliocentric(:,i) = tp%xh(:,i) + pl%inner(enc_index - 1)%x(:,ipleP) + end do + else + do i = 1, ntp + tp%xheliocentric(:,i) = tp%xh(:,i) + pl%inner(enc_index )%x(:,ipleP) + end do + end if + ! Swap the planetocentric and heliocentric position vectors + tp%xh(:,:) = tp%xheliocentric(:,:) + if (config%loblatecb) then + ! Put in the current encountering planet's oblateness acceleration as the central body's + if (tp%lfirst) then + cb_heliocentric%aobl(:) = pl%inner(enc_index - 1)%aobl(:, ipleP) + else + cb_heliocentric%aobl(:) = pl%inner(enc_index )%aobl(:, ipleP) + end if + call tp%obl_acc(cb_heliocentric) + end if - if (config%lextra_force) call tp%user_getacch(cb_heliocen, config, t) - if (config%lgr) call tp%gr_getacch(cb_heliocen, config) - tp%xh(:,:) = xh_original(:,:) + if (config%lextra_force) call tp%user_getacch(cb_heliocentric, config, t) + if (config%lgr) call tp%gr_getacch(cb_heliocentric, config) + + call move_alloc(xh_original, tp%xh) + end select else ! Not a close encounter, so just proceded with the standard WHM method - call whm_getacch_tp(tp, cb, pl, config, t, xh) + call whm_getacch_tp(tp, cb, pl, config, t, pl%xh) end if end associate return - end subroutine rmvs_getacch_in_tp + end subroutine rmvs_getacch_tp end submodule s_rmvs_getacch \ No newline at end of file diff --git a/src/rmvs/rmvs_interp.f90 b/src/rmvs/rmvs_interp.f90 deleted file mode 100644 index a5e239ceb..000000000 --- a/src/rmvs/rmvs_interp.f90 +++ /dev/null @@ -1,125 +0,0 @@ -submodule (rmvs_classes) s_rmvs_interp - use swiftest -contains - module subroutine rmvs_interp_in(self, cb, dt) - !! author: David A. Minton - !! - !! Interpolate planet positions between two Keplerian orbits in inner encounter regio - !! - !! Adapted from David E. Kaufmann's Swifter routine rmvs_interp_in.f90 - !! - !! Adapted from Hal Levison's Swift routine rmvs3_interp.f - implicit none - ! Arguments - class(rmvs_pl), intent(inout) :: self !! RMVS test particle object - class(rmvs_cb), intent(in) :: cb !! RMVS central body particle type - real(DP), intent(in) :: dt !! Step size - ! Internals - integer(I4B) :: i, j, iflag - real(DP) :: dti, frac, dntphenc - real(DP), dimension(NDIM) :: xtmp, vtmp - - dntphenc = real(NTPHENC, kind=DP) - associate (msun => cb%Gmass, npl => self%nbody) - dti = dt / dntphenc - do i = 1, npl - xtmp(:) = self%xin(:, i, 0) - vtmp(:) = self%vin(:, i, 0) - do j = 1, NTPHENC - 1 - call drift_one(msun, xtmp(:), vtmp(:), dti, iflag) - if (iflag /= 0) then - write(*, *) " Planet ", self%name(i), " is lost!!!!!!!!!!" - write(*, *) msun, dti - write(*, *) xtmp(:) - write(*, *) vtmp(:) - write(*, *) " STOPPING " - call util_exit(failure) - end if - frac = 1.0_DP - j / dntphenc - self%xin(:, i, j) = frac * xtmp(:) - self%vin(:, i, j) = frac * vtmp(:) - end do - xtmp(:) = self%xin(:, i, NTPHENC) - vtmp(:) = self%vin(:, i, NTPHENC) - do j = NTPHENC - 1, 1, -1 - call drift_one(msun, xtmp(:), vtmp(:), -dti, iflag) - if (iflag /= 0) then - write(*, *) " Planet ", self%name(i), " is lost!!!!!!!!!!" - write(*, *) msun, -dti - write(*, *) xtmp(:) - write(*, *) vtmp(:) - write(*, *) " STOPPING " - call util_exit(failure) - end if - frac = j / dntphenc - self%xin(:, i, j) = self%xin(:, i, j) + frac * xtmp(:) - self%vin(:, i, j) = self%vin(:, i, j) + frac * vtmp(:) - end do - end do - end associate - return - - end subroutine rmvs_interp_in - - module subroutine rmvs_interp_out(self, cb, dt) - !! author: David A. Minton - !! - !! Interpolate planet positions between two Keplerian orbits in outer encounter region - !! - !! Adapted from David E. Kaufmann's Swifter routine rmvs_interp_out.f90 - !! - !! Adapted from Hal Levison's Swift routine rmvs3_interp.f - implicit none - ! Arguments - class(rmvs_pl), intent(inout) :: self !! RMVS test particle object - class(rmvs_cb), intent(in) :: cb !! RMVS central body particle type - real(DP), intent(in) :: dt !! Step size - ! Internals - integer(I4B) :: i, j, iflag - real(DP) :: dto, frac, dntenc - real(DP), dimension(NDIM) :: xtmp, vtmp - - ! executable code - dntenc = real(NTENC, DP) - associate (msun => cb%Gmass, npl => self%nbody) - dto = dt / dntenc - do i = 1, npl - xtmp(:) = self%xout(:, i, 0) - vtmp(:) = self%vout(:, i, 0) - do j = 1, NTENC - 1 - call drift_one(msun, xtmp(:), vtmp(:), dto, iflag) - if (iflag /= 0) then - write(*, *) " Planet ", self%name(i), " is lost!!!!!!!!!!" - write(*, *) msun, dto - write(*, *) xtmp(:) - write(*, *) vtmp(:) - write(*, *) " STOPPING " - call util_exit(FAILURE) - end if - frac = 1.0_DP - j / dntenc - self%xout(:, i, j) = frac*xtmp(:) - self%vout(:, i, j) = frac*vtmp(:) - end do - xtmp(:) = self%xout(:, i, NTENC) - vtmp(:) = self%vout(:, i, NTENC) - do j = NTENC - 1, 1, -1 - call drift_one(msun, xtmp(:), vtmp(:), -dto, iflag) - if (iflag /= 0) then - write(*, *) " Planet ", self%name(i), " is lost!!!!!!!!!!" - write(*, *) msun, -dto - write(*, *) xtmp(:) - write(*, *) vtmp(:) - write(*, *) " STOPPING " - call util_exit(FAILURE) - end if - frac = j / dntenc - self%xout(:, i, j) = self%xout(:, i, j) + frac * xtmp(:) - self%vout(:, i, j) = self%vout(:, i, j) + frac * vtmp(:) - end do - end do - end associate - - return - - end subroutine rmvs_interp_out -end submodule s_rmvs_interp \ No newline at end of file diff --git a/src/rmvs/rmvs_obl.f90 b/src/rmvs/rmvs_obl.f90 deleted file mode 100644 index 80871fc72..000000000 --- a/src/rmvs/rmvs_obl.f90 +++ /dev/null @@ -1,31 +0,0 @@ -submodule(rmvs_classes) s_rmvs_obl - use swiftest -contains - module subroutine rmvs_obl_acc_in(self, cb) - !! author: David A. Minton - !! - !! Compute the oblateness acceleration in the inner encounter region with planets - !! - implicit none - ! Arguments - class(rmvs_pl), intent(inout) :: self !! RMVS massive body object - class(swiftest_cb), intent(inout) :: cb !! Swiftest central body object - ! Internals - integer(I4B) :: i - real(DP), dimension(:, :), allocatable :: xh_original - - associate(pl => self) - allocate(xh_original, source=pl%xh) - do i = 0, NTPHENC - pl%xh(:,:) = pl%xin(:,:,i) ! Temporarily replace heliocentric position with inner substep values to calculate the oblateness terms - call pl%obl_acc(cb) - pl%aoblin(:,:,i) = pl%aobl(:,:) ! Save the oblateness acceleration on the planet for this substep - end do - ! Put back the original heliocentric position for the planets - pl%xh(:,:) = xh_original(:,:) - deallocate(xh_original) - end associate - - return - end subroutine rmvs_obl_acc_in -end submodule s_rmvs_obl \ No newline at end of file diff --git a/src/rmvs/rmvs_peri.f90 b/src/rmvs/rmvs_peri.f90 deleted file mode 100644 index 4f8cf3612..000000000 --- a/src/rmvs/rmvs_peri.f90 +++ /dev/null @@ -1,87 +0,0 @@ -submodule(rmvs_classes) s_rmvs_peri - use swiftest -contains - module subroutine rmvs_peri_tp(self, cb, pl, t, dt, lfirst, index, ipleP, config) - !! author: David A. Minton - !! - !! Determine planetocentric pericenter passages for test particles in close encounters with a planet - !! - !! Adapted from Hal Levison's Swift routine Adapted from Hal Levison's Swift routine util_peri.f - !! Adapted from David E. Kaufmann's Swifter routine rmvs_peri.f90 - implicit none - ! Arguments - class(rmvs_tp), intent(inout) :: self !! RMVS test particle object - class(rmvs_cb), intent(inout) :: cb !! RMVS central body object - class(rmvs_pl), intent(inout) :: pl !! RMVS massive body object - real(DP), intent(in) :: t !! current time - real(DP), intent(in) :: dt !! step size - logical, intent(in) :: lfirst !! Logical flag indicating whether current invocation is the first - integer(I4B), intent(in) :: index !! outer substep number within current set - integer(I4B), intent(in) :: ipleP !!index of RMVS planet being closely encountered - class(swiftest_configuration), intent(in) :: config !! Input collection of configuration parameters - ! Internals - integer(I4B) :: i, id1, id2 - real(DP) :: r2, mu, rhill2, vdotr, a, peri, capm, tperi, rpl - real(DP), dimension(NDIM) :: xh1, xh2, vh1, vh2 - - rhill2 = pl%Rhill(ipleP)**2 - mu = pl%Gmass(ipleP) - associate(tp => self, nenc => self%nbody) - if (lfirst) then - do i = 1, nenc - if (tp%status(i) == ACTIVE) then - vdotr = dot_product(tp%xh(:, i), tp%vh(:, i)) - if (vdotr > 0.0_DP) then - tp%isperi(i) = 1 - else - tp%isperi(i) = -1 - end if - end if - end do - else - do i = 1, nenc - if (tp%status(i) == ACTIVE) then - vdotr = dot_product(tp%xh(:, i), tp%vh(:, i)) - if (tp%isperi(i) == -1) then - if (vdotr >= 0.0_DP) then - tp%isperi(i) = 0 - call orbel_xv2aqt(mu, tp%xh(:, i), tp%vh(:, i), a, peri, capm, tperi) - r2 = dot_product(tp%xh(:, i), tp%xh(:, i)) - if ((abs(tperi) > FACQDT * dt) .or. (r2 > rhill2)) peri = sqrt(r2) - if (config%encounter_file /= "") then - id1 = pl%name(ipleP) - rpl = pl%radius(ipleP) - xh1(:) = pl%xin(:, ipleP, index) - vh1(:) = pl%vin(:, ipleP, index) - id2 = tp%name(i) - xh2(:) = tp%xh(:, i) + pl%xin(:, ipleP, index) - vh2(:) = tp%vh(:, i) + pl%vin(:, ipleP, index) - call io_write_encounter(t, id1, id2, mu, 0.0_DP, rpl, 0.0_DP, xh1(:), xh2(:), vh1(:), vh2(:), & - config%encounter_file, config%out_type) - end if - if (tp%lperi(i)) then - if (peri < tp%peri(i)) then - tp%peri(i) = peri - tp%plperP(i) = ipleP - end if - else - tp%lperi(i) = .true. - tp%peri(i) = peri - tp%plperP(i) = ipleP - end if - end if - else - if (vdotr > 0.0_DP) then - tp%isperi(i) = 1 - else - tp%isperi(i) = -1 - end if - end if - end if - end do - end if - end associate - return - - end subroutine rmvs_peri_tp -end submodule s_rmvs_peri \ No newline at end of file diff --git a/src/rmvs/rmvs_setup.f90 b/src/rmvs/rmvs_setup.f90 index 030e66496..2f41870a8 100644 --- a/src/rmvs/rmvs_setup.f90 +++ b/src/rmvs/rmvs_setup.f90 @@ -13,30 +13,31 @@ module subroutine rmvs_setup_pl(self,n) integer(I4B), intent(in) :: n !! Number of test particles to allocate ! Internals integer(I4B) :: i,j + !type(swiftest_configuration) :: encounter_config !> Call allocation method for parent class - call whm_setup_pl(self, n) - if (n <= 0) return - - allocate(self%nenc(n)) - allocate(self%xout(NDIM, n, 0:NTENC)) - allocate(self%vout(NDIM, n, 0:NTENC)) - allocate(self%xin(NDIM, n, 0:NTPHENC)) - allocate(self%vin(NDIM, n, 0:NTPHENC)) - allocate(self%aoblin(NDIM, n, 0:NTPHENC)) - allocate(self%plind(n,n)) - allocate(self%tpenc(n)) - allocate(self%plenc(n)) - allocate(self%cbenc(n)) - self%plenc(:)%nbody = n - - self%nenc = 0 - self%xout(:,:,:) = 0.0_DP - self%vout(:,:,:) = 0.0_DP - self%xin(:,:,:) = 0.0_DP - self%vin(:,:,:) = 0.0_DP - self%aoblin(:,:,:) = 0.0_DP - + associate(system => self) + call whm_setup_pl(system, n) + if (n <= 0) return + + if (.not.system%lplanetocentric) then + allocate(self%nenc(n)) + self%nenc(:) = 0 + + ! Set up inner and outer planet interpolation vector storage containers + allocate(self%outer(0:NTENC)) + do i = 0, NTENC + allocate(self%outer(i)%x(NDIM, n)) + allocate(self%outer(i)%v(NDIM, n)) + end do + allocate(self%inner(0:NTPHENC)) + do i = 0, NTPHENC + allocate(self%inner(i)%x(NDIM, n)) + allocate(self%inner(i)%v(NDIM, n)) + allocate(self%inner(i)%aobl(NDIM, n)) + end do + end if + end associate return end subroutine rmvs_setup_pl @@ -58,10 +59,11 @@ module subroutine rmvs_setup_tp(self,n) allocate(self%lperi(n)) allocate(self%plperP(n)) allocate(self%plencP(n)) + if (self%lplanetocentric) then + allocate(self%xheliocentric(NDIM, n)) + end if self%lperi(:) = .false. - self%plperP(:) = 0 - self%plencP(:) = 0 return end subroutine rmvs_setup_tp @@ -69,8 +71,13 @@ end subroutine rmvs_setup_tp module subroutine rmvs_setup_system(self, config) !! author: David A. Minton !! - !! Wrapper method to initialize a basic Swiftest nbody system from files - !! + !! Wrapper method to initialize a basic Swiftest nbody system from files. + !! + !! We currently rearrange the pl order to keep it consistent with the way Swifter does it + !! In Swifter, the central body occupies the first position in the pl list, and during + !! encounters, the encountering planet is skipped in loops. In Swiftest, we instantiate an + !! RMVS nbody system object attached to each pl to store planetocentric versions of the system + !! to use during close encounters. implicit none ! Arguments class(rmvs_nbody_system), intent(inout) :: self !! RMVS system object @@ -80,33 +87,50 @@ module subroutine rmvs_setup_system(self, config) ! Call parent method call whm_setup_system(self, config) - ! Set up the tp-planet encounter structures + ! Set up the pl-tp planetocentric encounter structures for pl and cb. The planetocentric tp structures are + ! generated as necessary during close encounter steps. select type(pl => self%pl) class is(rmvs_pl) - select type(cb => self%cb) - class is (rmvs_cb) - select type (tp => self%tp) - class is (rmvs_tp) - associate(npl => pl%nbody) - tp%cb = cb - do i = 1, npl - allocate(pl%plenc(i)%Gmass(npl)) - allocate(pl%plenc(i)%xh(NDIM,npl)) - allocate(pl%plenc(i)%vh(NDIM,npl)) - pl%plind(i,:) = [(j,j=1,npl)] - pl%plind(i,2:npl) = pack(pl%plind(i,1:npl), pl%plind(i,1:npl) /= i) - ! Shift the planet masses up so that the central body is the first planet - ! During a planetocentric encounter (like the original Swifter) - pl%plenc(i)%Gmass(1) = cb%Gmass - pl%plenc(i)%Gmass(2:npl) = pl%Gmass(pl%plind(i,2:npl)) - pl%cbenc(i) = cb - pl%cbenc(i)%Gmass = pl%Gmass(i) - end do + select type(cb => self%cb) + class is (rmvs_cb) + select type (tp => self%tp) + class is (rmvs_tp) + tp%cb_heliocentric = cb + associate(npl => pl%nbody) + allocate(pl%planetocentric(npl)) + do i = 1, npl + allocate(pl%planetocentric(i)%cb, source=cb) + allocate(rmvs_pl :: pl%planetocentric(i)%pl) + associate(plenci => pl%planetocentric(i)%pl, cbenci => pl%planetocentric(i)%cb) + cbenci%lplanetocentric = .true. + + plenci%lplanetocentric = .true. + call plenci%setup(npl) + plenci%status(:) = ACTIVE + + ! plind stores the heliocentric index value of a planetocentric planet + ! e.g. Consider an encounter with planet 3. + ! Then the following will be the values of plind: + ! pl%planetocentric(3)%pl%plind(1) = 0 (central body - never used) + ! pl%planetocentric(3)%pl%plind(2) = 1 + ! pl%planetocentric(3)%pl%plind(3) = 2 + ! pl%planetocentric(3)%pl%plind(4) = 4 + ! pl%planetocentric(3)%pl%plind(5) = 5 + ! etc. + allocate(plenci%plind(npl)) + plenci%plind(1:npl) = [(j,j=1,npl)] + plenci%plind(2:npl) = pack(plenci%plind(1:npl), plenci%plind(1:npl) /= i) + plenci%plind(1) = 0 + plenci%Gmass(1) = cb%Gmass + plenci%Gmass(2:npl) = pl%Gmass(plenci%plind(2:npl)) + cbenci%Gmass = pl%Gmass(i) end associate - end select - end select + end do + end associate end select - + end select + end select + end subroutine rmvs_setup_system module subroutine rmvs_setup_set_beg_end(self, xbeg, xend, vbeg) @@ -115,8 +139,8 @@ module subroutine rmvs_setup_set_beg_end(self, xbeg, xend, vbeg) !! Sets one or more of the values of xbeg, xend, and vbeg implicit none ! Arguments - class(rmvs_tp), intent(inout) :: self !! RMVS test particle object - real(DP), dimension(:,:), optional :: xbeg, xend, vbeg + class(rmvs_tp), intent(inout) :: self !! RMVS test particle object + real(DP), dimension(:,:), intent(in), optional :: xbeg, xend, vbeg if (present(xbeg)) then if (allocated(self%xbeg)) deallocate(self%xbeg) @@ -135,4 +159,4 @@ module subroutine rmvs_setup_set_beg_end(self, xbeg, xend, vbeg) end subroutine rmvs_setup_set_beg_end -end submodule s_rmvs_setup \ No newline at end of file +end submodule s_rmvs_setup diff --git a/src/rmvs/rmvs_step.f90 b/src/rmvs/rmvs_step.f90 index b406afe9a..0f0b605ed 100644 --- a/src/rmvs/rmvs_step.f90 +++ b/src/rmvs/rmvs_step.f90 @@ -1,7 +1,7 @@ submodule(rmvs_classes) s_rmvs_step use swiftest contains - module subroutine rmvs_step_system(cb, pl, tp, config) + module subroutine rmvs_step_system(self, config) !! author: David A. Minton !! !! Step massive bodies and and active test particles ahead in heliocentric coordinates @@ -10,9 +10,7 @@ module subroutine rmvs_step_system(cb, pl, tp, config) !! Adapted from David E. Kaufmann's Swifter routine rmvs_step.f90 implicit none ! Arguments - class(rmvs_cb), intent(inout) :: cb !! RMVS central body object - class(rmvs_pl), intent(inout) :: pl !! RMVS central body object - class(rmvs_tp), intent(inout) :: tp !! RMVS central body object + class(rmvs_nbody_system), intent(inout) :: self !! RMVS nbody system object class(swiftest_configuration), intent(in) :: config !! Input collection of configuration parameters ! Internals logical :: lencounter, lfirstpl, lfirsttp @@ -20,9 +18,13 @@ module subroutine rmvs_step_system(cb, pl, tp, config) real(DP), dimension(:,:), allocatable :: xbeg, xend, vbeg integer(I4B) :: i - associate(ntp => tp%nbody, npl => pl%nbody, t => config%t, dt => config%dt, & - xh => pl%xh, vh => pl%vh, xj => pl%xj, vj => pl%vj, ah => pl%ah, eta => pl%eta, & ! These two lines of associations aid in debugging with gdb - xht => tp%xh, vht => tp%vh) + select type(cb => self%cb) + class is (rmvs_cb) + select type(pl => self%pl) + class is (rmvs_pl) + select type(tp => self%tp) + class is (rmvs_tp) + associate(ntp => tp%nbody, npl => pl%nbody, t => config%t, dt => config%dt) allocate(xbeg, source=pl%xh) allocate(vbeg, source=pl%vh) call pl%set_rhill(cb) @@ -33,14 +35,14 @@ module subroutine rmvs_step_system(cb, pl, tp, config) if (lencounter) then lfirstpl = pl%lfirst lfirsttp = tp%lfirst - pl%xout(:,:,0) = tp%xbeg(:,:) - pl%vout(:,:,0) = tp%vbeg(:,:) + pl%outer(0)%x(:,:) = xbeg(:,:) + pl%outer(0)%v(:,:) = vbeg(:,:) call pl%step(cb, config, t, dt) - pl%xout(:,:,NTENC) = pl%xh(:,:) - pl%vout(:,:,NTENC) = pl%vh(:,:) + pl%outer(NTENC)%x(:,:) = pl%xh(:,:) + pl%outer(NTENC)%v(:,:) = pl%vh(:,:) call tp%set_beg_end(xend = pl%xh) - call pl%interp_out(cb, dt) - call pl%step_out(cb, tp, dt, config) + call rmvs_interp_out(pl,cb, dt, config) + call rmvs_step_out(pl, cb, tp, dt, config) call tp%reverse_status() call tp%set_beg_end(xbeg = xbeg, xend = xend) tp%lfirst = .true. @@ -49,41 +51,334 @@ module subroutine rmvs_step_system(cb, pl, tp, config) pl%lfirst = lfirstpl tp%lfirst = lfirsttp else - call whm_step_system(cb, pl, tp, config) + call whm_step_system(self, config) end if end associate + end select + end select + end select + return end subroutine rmvs_step_system - module subroutine rmvs_step_out(self, cb, tp, dt, config) + subroutine rmvs_interp_in(pl, cb, dt, enc_index, config) + !! author: David A. Minton + !! + !! Interpolate planet positions between two Keplerian orbits in inner encounter regio + !! + !! Adapted from David E. Kaufmann's Swifter routine rmvs_interp_in.f90 + !! + !! Adapted from Hal Levison's Swift routine rmvs3_interp.f + implicit none + ! Arguments + class(rmvs_pl), intent(inout) :: pl !! RMVS test particle object + class(rmvs_cb), intent(inout) :: cb !! RMVS central body particle type + real(DP), intent(in) :: dt !! Step size + integer(I4B), intent(in) :: enc_index !! Outer substep number within current se + class(swiftest_configuration), intent(in) :: config !! Swiftest configuration file + ! Internals + integer(I4B) :: i, j + real(DP) :: frac, dntphenc + real(DP), dimension(:,:), allocatable :: xtmp, vtmp, xh_original + real(DP), dimension(:), allocatable :: msun, dti + integer(I4B), dimension(:), allocatable :: iflag + + associate (npl => pl%nbody) + dntphenc = real(NTPHENC, kind=DP) + + ! Set the endpoints of the inner region from the outer region values in the current outer step index + pl%inner(0)%x(:,:) = pl%outer(enc_index - 1)%x(:, :) + pl%inner(0)%v(:,:) = pl%outer(enc_index - 1)%v(:, :) + pl%inner(NTPHENC)%x(:,:) = pl%outer(enc_index)%x(:, :) + pl%inner(NTPHENC)%v(:,:) = pl%outer(enc_index)%v(:, :) + + allocate(xtmp,mold=pl%xh) + allocate(vtmp,mold=pl%vh) + allocate(msun(npl)) + allocate(dti(npl)) + allocate(iflag(npl)) + dti(:) = dt / dntphenc + msun(:) = cb%Gmass + xtmp(:, :) = pl%inner(0)%x(:, :) + vtmp(:, :) = pl%inner(0)%v(:, :) + if (config%loblatecb) then + allocate(xh_original,mold=pl%xh) + pl%xh(:, :) = xtmp(:, :) ! Temporarily replace heliocentric position with inner substep values to calculate the oblateness terms + call pl%obl_acc(cb) + pl%inner(0)%aobl(:, :) = pl%aobl(:, :) ! Save the oblateness acceleration on the planet for this substep + end if + + do j = 1, NTPHENC - 1 + !!$omp simd + do i = 1, npl + call drift_one(msun(i), xtmp(1,i), xtmp(2,i), xtmp(3,i), vtmp(1,i), vtmp(2,i), vtmp(3,i), dti(i), iflag(i)) + end do + if (any(iflag(:) /= 0)) then + do i = 1, npl + if (iflag(i) /=0) then + write(*, *) " Planet ", pl%name(i), " is lost!!!!!!!!!!" + write(*, *) msun(i), dti(i) + write(*, *) xtmp(:,i) + write(*, *) vtmp(:,i) + write(*, *) " STOPPING " + call util_exit(failure) + end if + end do + end if + frac = 1.0_DP - j / dntphenc + pl%inner(j)%x(:, :) = frac * xtmp(:,:) + pl%inner(j)%v(:, :) = frac * vtmp(:,:) + end do + + xtmp(:,:) = pl%inner(NTPHENC)%x(:, :) + vtmp(:,:) = pl%inner(NTPHENC)%v(:, :) + + do j = NTPHENC - 1, 1, -1 + !!$omp simd + do i = 1, npl + call drift_one(msun(i), xtmp(1,i), xtmp(2,i), xtmp(3,i), vtmp(1,i), vtmp(2,i), vtmp(3,i), -dti(i), iflag(i)) + end do + if (any(iflag(:) /= 0)) then + do i = 1, npl + if (iflag(i) /=0) then + write(*, *) " Planet ", pl%name(i), " is lost!!!!!!!!!!" + write(*, *) msun(i), -dti(i) + write(*, *) xtmp(:,i) + write(*, *) vtmp(:,i) + write(*, *) " STOPPING " + call util_exit(failure) + end if + end do + end if + frac = j / dntphenc + pl%inner(j)%x(:, :) = pl%inner(j)%x(:, :) + frac * xtmp(:, :) + pl%inner(j)%v(:, :) = pl%inner(j)%v(:, :) + frac * vtmp(:, :) + + if (config%loblatecb) then + pl%xh(:,:) = pl%inner(j)%x(:, :) + call pl%obl_acc(cb) + pl%inner(j)%aobl(:, :) = pl%aobl(:, :) + end if + end do + ! Put the planet positions back into place + if (config%loblatecb) call move_alloc(xh_original, pl%xh) + end associate + return + + end subroutine rmvs_interp_in + + subroutine rmvs_interp_out(pl, cb, dt, config) + !! author: David A. Minton + !! + !! Interpolate planet positions between two Keplerian orbits in outer encounter region + !! + !! Adapted from David E. Kaufmann's Swifter routine rmvs_interp_out.f90 + !! + !! Adapted from Hal Levison's Swift routine rmvs3_interp.f + implicit none + ! Arguments + class(rmvs_pl), intent(inout) :: pl !! RMVS test particle object + class(rmvs_cb), intent(inout) :: cb !! RMVS central body particle type + real(DP), intent(in) :: dt !! Step size + class(swiftest_configuration), intent(in) :: config !! Swiftest configuration file + ! Internals + integer(I4B) :: i, j + real(DP) :: frac, dntenc + real(DP), dimension(:,:), allocatable :: xtmp, vtmp + real(DP), dimension(:), allocatable :: msun, dto + integer(I4B), dimension(:), allocatable :: iflag + + ! executable code + dntenc = real(NTENC, DP) + associate (npl => pl%nbody) + allocate(xtmp, mold = pl%xh) + allocate(vtmp, mold = pl%vh) + allocate(msun(npl)) + allocate(dto(npl)) + allocate(iflag(npl)) + dto(:) = dt / dntenc + msun(:) = cb%Gmass + xtmp(:,:) = pl%outer(0)%x(:, :) + vtmp(:,:) = pl%outer(0)%v(:, :) + do j = 1, NTENC - 1 + !!$omp simd + do i = 1, npl + call drift_one(msun(i), xtmp(1,i), xtmp(2,i), xtmp(3,i), vtmp(1,i), vtmp(2,i), vtmp(3,i), dto(i), iflag(i)) + end do + if (any(iflag(:) /= 0)) then + do i = 1, npl + if (iflag(i) /= 0) then + write(*, *) " Planet ", pl%name(i), " is lost!!!!!!!!!!" + write(*, *) msun(i), dto(i) + write(*, *) xtmp(:,i) + write(*, *) vtmp(:,i) + write(*, *) " STOPPING " + call util_exit(FAILURE) + end if + end do + end if + frac = 1.0_DP - j / dntenc + pl%outer(j)%x(:, :) = frac * xtmp(:,:) + pl%outer(j)%v(:, :) = frac * vtmp(:,:) + end do + xtmp(:,:) = pl%outer(NTENC)%x(:, :) + vtmp(:,:) = pl%outer(NTENC)%v(:, :) + do j = NTENC - 1, 1, -1 + !!$omp simd + do i = 1, npl + call drift_one(msun(i), xtmp(1,i), xtmp(2,i), xtmp(3,i), vtmp(1,i), vtmp(2,i), vtmp(3,i), -dto(i), iflag(i)) + end do + if (any(iflag(:) /= 0)) then + do i = 1, npl + if (iflag(i) /= 0) then + write(*, *) " Planet ", pl%name(i), " is lost!!!!!!!!!!" + write(*, *) msun(i), -dto(i) + write(*, *) xtmp(:,i) + write(*, *) vtmp(:,i) + write(*, *) " STOPPING " + call util_exit(FAILURE) + end if + end do + end if + frac = j / dntenc + pl%outer(j)%x(:, :) = pl%outer(j)%x(:, :) + frac * xtmp(:,:) + pl%outer(j)%v(:, :) = pl%outer(j)%v(:, :) + frac * vtmp(:,:) + end do + end associate + + return + + end subroutine rmvs_interp_out + + subroutine rmvs_peri_tp(tp, pl, t, dt, lfirst, enc_index, ipleP, config) + !! author: David A. Minton + !! + !! Determine planetocentric pericenter passages for test particles in close encounters with a planet + !! + !! Adapted from Hal Levison's Swift routine Adapted from Hal Levison's Swift routine util_peri.f + !! Adapted from David E. Kaufmann's Swifter routine rmvs_peri.f90 + implicit none + ! Arguments + class(rmvs_tp), intent(inout) :: tp !! RMVS test particle object (planetocentric) + class(rmvs_pl), intent(inout) :: pl !! RMVS massive body object (heliocentric) + real(DP), intent(in) :: t !! current time + real(DP), intent(in) :: dt !! step size + logical, intent(in) :: lfirst !! Logical flag indicating whether current invocation is the first + integer(I4B), intent(in) :: enc_index !! Outer substep number within current set + integer(I4B), intent(in) :: ipleP !! index of RMVS planet being closely encountered + class(swiftest_configuration), intent(in) :: config !! Input collection of configuration parameters + ! Internals + integer(I4B) :: i, id1, id2 + real(DP) :: r2, mu, rhill2, vdotr, a, peri, capm, tperi, rpl + real(DP), dimension(NDIM) :: xh1, xh2, vh1, vh2 + + rhill2 = pl%rhill(ipleP)**2 + mu = pl%Gmass(ipleP) + associate( nenc => tp%nbody, xpc => tp%xh, vpc => tp%vh) + if (lfirst) then + do i = 1, nenc + if (tp%status(i) == ACTIVE) then + vdotr = dot_product(xpc(:, i), vpc(:, i)) + if (vdotr > 0.0_DP) then + tp%isperi(i) = 1 + else + tp%isperi(i) = -1 + end if + end if + end do + else + do i = 1, nenc + if (tp%status(i) == ACTIVE) then + vdotr = dot_product(xpc(:, i), vpc(:, i)) + if (tp%isperi(i) == -1) then + if (vdotr >= 0.0_DP) then + tp%isperi(i) = 0 + call orbel_xv2aqt(mu, xpc(:, i), vpc(:, i), a, peri, capm, tperi) + r2 = dot_product(xpc(:, i), xpc(:, i)) + if ((abs(tperi) > FACQDT * dt) .or. (r2 > rhill2)) peri = sqrt(r2) + if (config%encounter_file /= "") then + id1 = pl%name(ipleP) + rpl = pl%radius(ipleP) + xh1(:) = pl%inner(enc_index)%x(:, ipleP) + vh1(:) = pl%inner(enc_index)%v(:, ipleP) + id2 = tp%name(i) + xh2(:) = xpc(:, i) + xh1(:) + vh2(:) = xpc(:, i) + vh1(:) + call io_write_encounter(t, id1, id2, mu, 0.0_DP, rpl, 0.0_DP, xh1(:), xh2(:), vh1(:), vh2(:), & + config%encounter_file, config%out_type) + end if + if (tp%lperi(i)) then + if (peri < tp%peri(i)) then + tp%peri(i) = peri + tp%plperP(i) = ipleP + end if + else + tp%lperi(i) = .true. + tp%peri(i) = peri + tp%plperP(i) = ipleP + end if + end if + else + if (vdotr > 0.0_DP) then + tp%isperi(i) = 1 + else + tp%isperi(i) = -1 + end if + end if + end if + end do + end if + end associate + return + + end subroutine rmvs_peri_tp + + subroutine rmvs_step_out(pl, cb, tp, dt, config) !! author: David A. Minton !! - !! Step ACTIVE test particles ahead in the outer encounter region + !! Step ACTIVE test particles ahead in the outer encounter region, setting up and calling the inner region + !! integration if necessar !! - !! Adapted from Hal Levison's Swift routine rmvs3_step_out.f - !! Adapted from David E. Kaufmann's Swifter routine rmvs_step_out.f90 + !! Adapted from Hal Levison's Swift routines rmvs3_step_out.f and rmvs3_step_out2.f + !! Adapted from David E. Kaufmann's Swifter routines rmvs_step_out.f90 and rmvs_step_out2.f90 implicit none ! Arguments - class(rmvs_pl), intent(inout) :: self !! RMVS massive body object + class(rmvs_pl), intent(inout) :: pl !! RMVS massive body object class(rmvs_cb), intent(inout) :: cb !! RMVS central body object class(rmvs_tp), intent(inout) :: tp !! RMVS test particle object real(DP), intent(in) :: dt !! Step size class(swiftest_configuration), intent(in) :: config !! Input collection of configuration parameters ! Internals - integer(I4B) :: i, j, k - real(DP) :: dto, time + integer(I4B) :: enc_index, j, k + real(DP) :: dto, outer_time, rts + logical :: lencounter, lfirsttp - ! executable code - associate(pl => self, npl => self%nbody, ntp => tp%nbody, t => config%t, xht => tp%xh, vht => tp%vh) + associate(npl => pl%nbody, ntp => tp%nbody, t => config%t) dto = dt / NTENC where(tp%plencP(:) == 0) tp%status(:) = INACTIVE elsewhere tp%lperi(:) = .false. end where - do i = 1, NTENC - time = t + (i - 1) * dto - call rmvs_step_out2(cb, pl, tp, time, dto, i, config) + do enc_index = 1, NTENC + outer_time = t + (enc_index - 1) * dto + call tp%set_beg_end(xbeg = pl%outer(enc_index - 1)%x(:, :), & + vbeg = pl%outer(enc_index - 1)%v(:, :), & + xend = pl%outer(enc_index )%x(:, :)) + rts = RHPSCALE + lencounter = tp%encounter_check(cb, pl, dt, rts) + if (lencounter) then + ! Interpolate planets in inner encounter region + call rmvs_interp_in(pl, cb, dto, enc_index, config) + ! Step through the inner region + call rmvs_step_in(pl, cb, tp, config, outer_time, dto) + lfirsttp = tp%lfirst + tp%lfirst = .true. + call tp%step(cb, pl, config, outer_time, dto) + tp%lfirst = lfirsttp + else + call tp%step(cb, pl, config, outer_time, dto) + end if do j = 1, npl if (pl%nenc(j) == 0) cycle where((tp%plencP(:) == j) .and. (tp%status(:) == INACTIVE)) tp%status(:) = ACTIVE @@ -94,57 +389,7 @@ module subroutine rmvs_step_out(self, cb, tp, dt, config) end subroutine rmvs_step_out - module subroutine rmvs_step_out2(cb, pl, tp, t, dt, index, config) - !! author: David A. Minton - !! - !! Step active test particles ahead in the outer encounter region, setting up and calling the inner region - !! integration if necessary - !! - !! Adapted from Hal Levison's Swift routine rmvs3_step_out.f - !! Adapted from David E. Kaufmann's Swifter routine rmvs_step_out.f90 - implicit none - ! Arguments - class(rmvs_cb), intent(inout) :: cb !! RMVS central body object - class(rmvs_pl), intent(inout) :: pl !! RMVS massive body object - class(rmvs_tp), intent(inout) :: tp !! RMVS test particle object - real(DP), intent(in) :: t !! Simulation time - real(DP), intent(in) :: dt !! Step size - integer(I4B), intent(in) :: index !! outer substep number within current set - class(swiftest_configuration), intent(in) :: config !! Input collection of configuration parameters - ! Internals - logical :: lfirsttp - real(DP) :: rts - logical :: lencounter - - associate(xht => tp%xh, vht => tp%vh, xbeg => tp%xbeg, xend => tp%xend, aht => tp%ah,& - xout1 => pl%xout(:, :, index-1), vout1 => pl%vout(:,:,index-1), xout2 => pl%xout(:,:,index)) - - call tp%set_beg_end(xbeg = pl%xout(:, :, index-1), & - vbeg = pl%vout(:, :, index-1), & - xend = pl%xout(:, :, index)) - rts = RHPSCALE - lencounter = tp%encounter_check(cb, pl, dt, rts) - if (lencounter) then - ! Interpolate planets in inner encounter region - pl%xin(:,:,0) = tp%xbeg(:,:) - pl%vin(:,:,0) = tp%vbeg(:,:) - pl%xin(:,:,NTPHENC) = tp%xend(:, :) - pl%vin(:,:,NTPHENC) = pl%vout(:, :, index) - call pl%interp_in(cb, dt) - call pl%step_in(cb, tp, config, t, dt) - lfirsttp = tp%lfirst - tp%lfirst = .true. - call tp%step(cb, pl, config, t, dt) - tp%lfirst = lfirsttp - else - call tp%step(cb, pl, config, t, dt) - end if - end associate - return - - end subroutine rmvs_step_out2 - - module subroutine rmvs_step_in_pl(self, cb, tp, config, t, dt) + subroutine rmvs_step_in(pl, cb, tp, config, outer_time, dto) !! author: David A. Minton !! !! Step active test particles ahead in the inner encounter region @@ -153,63 +398,81 @@ module subroutine rmvs_step_in_pl(self, cb, tp, config, t, dt) !! Adapted from David E. Kaufmann's Swifter routine rmvs_step_in.f90 implicit none ! Arguments - class(rmvs_pl), intent(inout) :: self !! RMVS massive body object + class(rmvs_pl), intent(inout) :: pl !! RMVS massive body object class(rmvs_cb), intent(inout) :: cb !! RMVS central body object class(rmvs_tp), intent(inout) :: tp !! RMVS test particle object class(swiftest_configuration), intent(in) :: config !! Input collection of configuration parameters - real(DP), intent(in) :: t !! Current time - real(DP), intent(in) :: dt !! Step size + real(DP), intent(in) :: outer_time !! Current time + real(DP), intent(in) :: dto !! Step size ! Internals logical :: lfirsttp - integer(I4B) :: i, j - real(DP) :: dti, time - real(DP), dimension(NDIM, self%nbody) :: xbeg, xend, vbeg + integer(I4B) :: i, j, ipleP + real(DP) :: dti, inner_time + real(DP), dimension(:, :), allocatable :: xbeg, xend, vbeg - dti = dt / NTPHENC - associate(pl => self, npl => self%nbody, xht => tp%xh, vht => tp%vh, plind => self%plind, nenc => self%nenc) - if (config%loblatecb) call pl%obl_acc_in(cb) - call pl%make_planetocentric(cb, tp, config) - do i = 1, npl - if (nenc(i) == 0) cycle - ! There are inner encounters with this planet...switch to planetocentric coordinates to proceed - time = t - call pl%tpenc(i)%peri_pass(cb, pl, time, dti, .true., 0, i, config) - ! now step the encountering test particles fully through the inner encounter - lfirsttp = .true. - associate(index => pl%tpenc(i)%index, & - xpc => self%tpenc(i)%xh, vpc => self%tpenc(i)%vh, apc => self%tpenc(i)%ah) - pl%tpenc(i)%lfirst = .true. - do index = 1, NTPHENC ! Integrate over the encounter region, using the "substitute" planetocentric systems at each level - xbeg(:,1) = cb%xin(:) - pl%xin(:, i, index - 1) - xend(:,1) = cb%xin(:) - pl%xin(:, i, index) - vbeg(:,1) = cb%vin(:) - pl%vin(:, i, index - 1) - do j = 1, NDIM - xbeg(j,2:npl) = pl%xin(j, plind(i,2:npl), index - 1) - pl%xin(j, i, index - 1) - xend(j,2:npl) = pl%xin(j, plind(i,2:npl), index) - pl%xin(j, i, index) - vbeg(j,2:npl) = pl%vin(j, plind(i,2:npl), index - 1) - pl%vin(j, i, index - 1) - end do - pl%plenc(i)%xh(:,:) = xbeg(:,:) - pl%plenc(i)%vh(:,:) = vbeg(:,:) - call pl%tpenc(i)%set_beg_end(xbeg = xbeg, xend = xend) - call pl%tpenc(i)%step(pl%cbenc(i), pl%plenc(i), config, time, dti) - do j = 1, NDIM - pl%tpenc(i)%xheliocen(j, :) = pl%tpenc(i)%xh(j, :) + pl%xin(j, i, index) + associate(npl => pl%nbody, nenc => pl%nenc, & + cbenci => pl%planetocentric(i)%cb, & + plenci => pl%planetocentric(i)%pl, & + tpenci => pl%planetocentric(i)%tp) + associate(enc_index => tpenci%index, plind => plenci%plind) + + dti = dto / NTPHENC + allocate(xbeg, mold=pl%xh) + allocate(xend, mold=pl%xh) + allocate(vbeg, mold=pl%vh) + if (config%loblatecb) call pl%obl_acc(cb) + call rmvs_make_planetocentric(pl, cb, tp, config) + do i = 1, npl + if (nenc(i) == 0) cycle + ! There are inner encounters with this planet...switch to planetocentric coordinates to proceed + tpenci%lfirst = .true. + inner_time = outer_time + call rmvs_peri_tp(tpenci, pl, inner_time, dti, .true., 0, i, config) + ! now step the encountering test particles fully through the inner encounter + lfirsttp = .true. + do enc_index = 1, NTPHENC ! Integrate over the encounter region, using the "substitute" planetocentric systems at each level + write(*,*) 'enc_index: ',enc_index + write(*,*) 'xin-1 ',pl%inner(enc_index - 1)%x(:, i) + write(*,*) 'xin ',pl%inner(enc_index )%x(:, i) + xbeg(:,1) = cb%xin(:) - pl%inner(enc_index - 1)%x(:, i) + xend(:,1) = cb%xin(:) - pl%inner(enc_index )%x(:, i) + vbeg(:,1) = cb%vin(:) - pl%inner(enc_index - 1)%v(:, i) + do j = 2, npl + ipleP = plind(j) + write(*,*) 'sub planet ', j + write(*,*) 'plind(i,j) ',ipleP + write(*,*) 'xin-1 ',pl%inner(enc_index - 1)%x(:,ipleP) + write(*,*) 'xin ',pl%inner(enc_index )%x(:,ipleP) + xbeg(:,j) = pl%inner(enc_index - 1)%x(:,ipleP) - pl%inner(enc_index - 1)%x(:,i) + xend(:,j) = pl%inner(enc_index )%x(:,ipleP) - pl%inner(enc_index )%x(:,i) + vbeg(:,j) = pl%inner(enc_index - 1)%v(:,ipleP) - pl%inner(enc_index - 1)%v(:,i) + end do + plenci%xh(:,:) = xbeg(:,:) + plenci%vh(:,:) = vbeg(:,:) + call tpenci%set_beg_end(xbeg = xbeg, xend = xend) + write(*,*) 'xend check' + do j = 1, npl + write(*,*) 'xend: ',j,xend(:,j) + write(*,*) 'xenc: ',j,tpenci%xend(:,j) + end do + call tpenci%step(cbenci, plenci, config, inner_time, dti) + do j = 1, nenc(i) + tpenci%xheliocentric(:, j) = tpenci%xh(:, j) + pl%inner(enc_index )%x(:,i) + end do + inner_time = outer_time + j * dti + call rmvs_peri_tp(tpenci, pl, inner_time, dti, .false., enc_index, i, config) end do - time = config%t + j * dti - call pl%tpenc(i)%peri_pass(cb, pl, time, dti, .false., index, i, config) - end do - where(pl%tpenc(i)%status(:) == ACTIVE) pl%tpenc(i)%status(:) = INACTIVE - end associate - end do - call pl%end_planetocentric(cb,tp) + where(tpenci%status(:) == ACTIVE) tpenci%status(:) = INACTIVE + end do + call rmvs_end_planetocentric(pl, cb, tp) + end associate end associate - return - end subroutine rmvs_step_in_pl + end subroutine rmvs_step_in - module subroutine rmvs_step_make_planetocentric(self, cb, tp, config) + subroutine rmvs_make_planetocentric(pl, cb, tp, config) !! author: David A. Minton !! !! When encounters are detected, this method will call the interpolation methods for the planets and @@ -218,7 +481,7 @@ module subroutine rmvs_step_make_planetocentric(self, cb, tp, config) !! implicit none ! Arguments - class(rmvs_pl), intent(inout) :: self !! RMVS test particle object + class(rmvs_pl), intent(inout) :: pl !! RMVS test particle object class(rmvs_cb), intent(inout) :: cb !! RMVS central body particle type class(rmvs_tp), intent(inout) :: tp !! RMVS test particle object class(swiftest_configuration), intent(in) :: config !! Input collection of configuration parameters @@ -226,9 +489,8 @@ module subroutine rmvs_step_make_planetocentric(self, cb, tp, config) integer(I4B) :: i, j logical, dimension(:), allocatable :: encmask - associate(pl => self, npl => self%nbody, ntp => tp%nbody, nenc => self%nenc, tpenc => self%tpenc, cbenc => self%cbenc, & - plenc => self%plenc, GMpl => self%Gmass) - + associate(npl => pl%nbody, ntp => tp%nbody, GMpl => pl%Gmass, nenc => pl%nenc, & + cbenci => pl%planetocentric(i)%cb) do i = 1, npl if (nenc(i) == 0) cycle ! There are inner encounters with this planet @@ -236,58 +498,39 @@ module subroutine rmvs_step_make_planetocentric(self, cb, tp, config) allocate(encmask(ntp)) encmask(:) = tp%plencP(:) == i - ! Save the index value of the planet corresponding to this encounter - tpenc(i)%ipleP = i - tpenc(i)%lplanetocentric = .true. - tpenc(i)%nbody = nenc(i) - !call tpenc(i)%setup(nenc(i)) - ! Create space for the heliocentric position values for those acceleration calculations that need them - allocate(tpenc(i)%xheliocen(NDIM, nenc(i))) - allocate(tpenc(i)%xh(NDIM, nenc(i))) - allocate(tpenc(i)%vh(NDIM, nenc(i))) - allocate(tpenc(i)%ah(NDIM, nenc(i))) - allocate(tpenc(i)%status(nenc(i))) - allocate(tpenc(i)%mu(nenc(i))) - allocate(tpenc(i)%isperi(nenc(i))) - allocate(tpenc(i)%name(nenc(i))) - allocate(tpenc(i)%plperP(nenc(i))) - allocate(tpenc(i)%lperi(nenc(i))) - allocate(tpenc(i)%peri(nenc(i))) - tpenc(i)%status(:) = ACTIVE - tpenc(i)%lperi(:) = .false. - tpenc(i)%plperP(:) = 0 - tpenc(i)%name(:) = pack(tp%name(:), encmask(:)) - !call tp%spill(tpenc(i), pl%encmask(:,i)) - ! Grab all the encountering test particles and convert them to a planetocentric frame - do j = 1, NDIM - tpenc(i)%xheliocen(j, :) = pack(tp%xh(j,:), encmask(:)) - tpenc(i)%xh(j, :) = tpenc(i)%xheliocen(j, :) - pl%xin(j, i, 0) - tpenc(i)%vh(j, :) = pack(tp%vh(j,:), encmask(:)) - pl%vin(j, i, 0) - end do - - ! Make sure that the test particles get the planetocentric value of mu - call tpenc(i)%set_mu(pl%cbenc(i)) + ! Create encountering test particle structure + allocate(rmvs_tp :: pl%planetocentric(i)%tp) + associate(tpenci => pl%planetocentric(i)%tp) + tpenci%lplanetocentric = .true. + call tpenci%setup(nenc(i)) + tpenci%ipleP = i + tpenci%status(:) = ACTIVE + tpenci%name(:) = pack(tp%name(:), encmask(:)) + ! Grab all the encountering test particles and convert them to a planetocentric frame + do j = 1, NDIM + tpenci%xheliocentric(j, :) = pack(tp%xh(j,:), encmask(:)) + tpenci%xh(j, :) = tpenci%xheliocentric(j, :) - pl%inner(0)%x(j, i) + tpenci%vh(j, :) = pack(tp%vh(j,:), encmask(:)) - pl%inner(0)%v(j, i) + end do - ! Save the heliocentric position values of the encountering planet - !allocate(tpenc(i)%xh_pl(NDIM,0:NTPHENC)) - ! Save the encountering planet's values of oblateness acceleration - if (config%loblatecb) then - allocate(tpenc(i)%aoblin_pl(NDIM,0:NTPHENC)) - tpenc(i)%aoblin_pl(:,:) = pl%aoblin(:,i,0:NTPHENC) - end if + ! Make sure that the test particles get the planetocentric value of mu + write(*,*) i,'nenc: ',nenc(i) + write(*,*) cbenci%Gmass + call tpenci%set_mu(cbenci) + end associate end do end associate return - end subroutine rmvs_step_make_planetocentric + end subroutine rmvs_make_planetocentric - module subroutine rmvs_step_end_planetocentric(self, cb, tp) + subroutine rmvs_end_planetocentric(pl, cb, tp) !! author: David A. Minton !! !! Deallocates all of the encountering particle data structures for next time !! implicit none ! Arguments - class(rmvs_pl), intent(inout) :: self !! RMVS test particle object + class(rmvs_pl), intent(inout) :: pl !! RMVS test particle object class(rmvs_cb), intent(inout) :: cb !! RMVS central body object class(rmvs_tp), intent(inout) :: tp !! RMVS test particle object ! Internals @@ -295,8 +538,7 @@ module subroutine rmvs_step_end_planetocentric(self, cb, tp) integer(I4B), dimension(:), allocatable :: tpind logical, dimension(:), allocatable :: encmask - associate(pl => self, nenc => self%nenc, npl => self%nbody, ntp => tp%nbody, name => tp%name, & - tpenc => self%tpenc, plenc => self%plenc, cbenc => self%cbenc) + associate(nenc => pl%nenc, npl => pl%nbody, ntp => tp%nbody, tpenci => pl%planetocentric(i)%tp) do i = 1, npl if (nenc(i) == 0) cycle @@ -306,29 +548,18 @@ module subroutine rmvs_step_end_planetocentric(self, cb, tp) allocate(encmask(ntp)) encmask(:) = tp%plencP(:) == i tpind(:) = pack([(j,j=1,ntp)], encmask(:)) - + ! Copy the results of the integration back over and shift back to heliocentric reference - tp%status(tpind(1:nenc(i))) = tpenc(i)%status(1:nenc(i)) + tp%status(tpind(1:nenc(i))) = tpenci%status(1:nenc(i)) do j = 1, NDIM - tp%xh(j, tpind(1:nenc(i))) = tpenc(i)%xh(j,1:nenc(i)) + pl%xin(j, i, NTPHENC) - tp%vh(j, tpind(1:nenc(i))) = tpenc(i)%vh(j,1:nenc(i)) + pl%vin(j, i, NTPHENC) + tp%xh(j, tpind(1:nenc(i))) = tpenci%xh(j,1:nenc(i)) + pl%inner(NTPHENC)%x(j, i) + tp%vh(j, tpind(1:nenc(i))) = tpenci%vh(j,1:nenc(i)) + pl%inner(NTPHENC)%v(j, i) end do - deallocate(tpenc(i)%xheliocen) - deallocate(tpenc(i)%xh) - deallocate(tpenc(i)%vh) - deallocate(tpenc(i)%ah) - deallocate(tpenc(i)%status) - deallocate(tpenc(i)%mu) - deallocate(tpenc(i)%isperi) - deallocate(tpenc(i)%name) - deallocate(tpenc(i)%plperP) - deallocate(tpenc(i)%lperi) - deallocate(tpenc(i)%peri) - deallocate(tpind) + deallocate(pl%planetocentric(i)%tp) end do end associate return - end subroutine rmvs_step_end_planetocentric + end subroutine rmvs_end_planetocentric end submodule s_rmvs_step diff --git a/src/setup/setup.f90 b/src/setup/setup.f90 index 7f2f26d7c..f7f3ba9fd 100644 --- a/src/setup/setup.f90 +++ b/src/setup/setup.f90 @@ -1,12 +1,16 @@ submodule (swiftest_classes) s_setup use swiftest contains - module procedure setup_construct_system + module subroutine setup_construct_system(system, config) !! author: David A. Minton !! !! Constructor for a Swiftest nbody system. Creates the nbody system object based on the user-input integrator !! implicit none + ! Arguments + class(swiftest_nbody_system), allocatable, intent(inout) :: system !! Swiftest system object + type(swiftest_configuration), intent(in) :: config !! Swiftest configuration parameters + select case(config%integrator) case (BS) write(*,*) 'Bulirsch-Stoer integrator not yet enabled' @@ -52,7 +56,7 @@ end select return - end procedure setup_construct_system + end subroutine setup_construct_system module procedure setup_body !! author: David A. Minton diff --git a/src/step/step.f90 b/src/step/step.f90 deleted file mode 100644 index 7fddc8225..000000000 --- a/src/step/step.f90 +++ /dev/null @@ -1,53 +0,0 @@ -submodule (swiftest_classes) s_step - use swiftest -contains - module procedure step_system - !! author: David A. Minton - !! - !! Wrapper function that selects the correct step function for the integrator type. - !! The nested select statements serve to make this method a "bridge" between the polymorphic swiftest_nbody_system class - !! in which the cb, pl, and tp components are allocatable abstract classes and the actual integrator-specific methods that are - !! called internally. Before this point, the actual types of cb, pl, and tp are ambiguous. The select type constructs remove the - !! ambiguity. - !! - implicit none - - select type(system => self) - class is (whm_nbody_system) - select type(cb => self%cb) - class is (whm_cb) - select type(pl => self%pl) - class is (whm_pl) - select type(tp => self%tp) - class is (whm_tp) - call whm_step_system(cb, pl, tp, config) - end select - end select - end select - class is (rmvs_nbody_system) - select type(cb => self%cb) - class is (rmvs_cb) - select type(pl => self%pl) - class is (rmvs_pl) - select type(tp => self%tp) - class is (rmvs_tp) - call rmvs_step_system(cb, pl, tp, config) - end select - end select - end select - class is (helio_nbody_system) - select type(cb => self%cb) - class is (helio_cb) - select type(pl => self%pl) - class is (helio_pl) - select type(tp => self%tp) - class is (helio_tp) - call helio_step_system(cb, pl, tp, config) - end select - end select - end select - end select - return - end procedure step_system - -end submodule s_step diff --git a/src/symba/symba_helio_drift.f90 b/src/symba/symba_helio_drift.f90 index 87855b2f0..07b8d97a5 100644 --- a/src/symba/symba_helio_drift.f90 +++ b/src/symba/symba_helio_drift.f90 @@ -19,7 +19,7 @@ !$omp private (i, iflag) do i = 2, npl if ((symba_plA%levelg(i) == irec) .and. (symba_plA%status(i) == ACTIVE)) then - call drift_one(mu, symba_plA%xh(:,i), symba_plA%vb(:,i), dt, iflag) + call drift_one(mu, symba_plA%xh(1,i), symba_plA%xh(2,i), symba_plA%xh(3,i), symba_plA%vb(1,i), symba_plA%vb(2,i), symba_plA%vb(3,i), dt, iflag) if (iflag /= 0) then write(*, *) " massive body ", symba_plA%name(i), " is lost!!!!!!!!!!" write(*, *) mu, dt diff --git a/src/symba/symba_helio_drift_tp.f90 b/src/symba/symba_helio_drift_tp.f90 index 9eb3c615d..05d248360 100644 --- a/src/symba/symba_helio_drift_tp.f90 +++ b/src/symba/symba_helio_drift_tp.f90 @@ -11,9 +11,10 @@ implicit none integer(I4B) :: i, iflag + associate(symba_tpA%xh => xh, symba_tpA%vb => vb) do i = 1, ntp if ((symba_tpA%levelg(i) == irec) .and. (symba_tpA%status(i) == ACTIVE)) then - call drift_one(mu, symba_tpA%xh(:,i), symba_tpA%vb(:,i), dt, iflag) + call drift_one(mu, xh(1,i), xh(2,i), xh(3,i), vb(1,i), vb(2,i), vb(3,i), dt, iflag) if (iflag /= 0) then symba_tpA%status(i) = DISCARDED_DRIFTERR write(*, *) "particle ", symba_tpA%name(i), " lost due to error in danby drift" diff --git a/src/whm/whm_drift.f90 b/src/whm/whm_drift.f90 index 85997fa0c..082953963 100644 --- a/src/whm/whm_drift.f90 +++ b/src/whm/whm_drift.f90 @@ -43,15 +43,15 @@ module subroutine whm_drift_pl(self, cb, config, dt) dtp(:) = dt end if - !$omp simd + !!$omp simd do i = 1, npl - call drift_one(mu(i), xj(:, i), vj(:, i), dtp(i), iflag(i)) + call drift_one(mu(i), xj(1, i), xj(2, i), xj(3, i), vj(1, i), vj(2, i), vj(3, i), dtp(i), iflag(i)) end do - if (any(iflag(1:npl) /= 0)) then + if (any(iflag(:) /= 0)) then do i = 1, npl write(*, *) " Planet ", self%name(i), " is lost!!!!!!!!!!" - write(*, *) xj - write(*, *) vj + write(*, *) xj(:,i) + write(*, *) vj(:,i) write(*, *) " stopping " call util_exit(FAILURE) end do @@ -82,11 +82,13 @@ module subroutine whm_drift_tp(self, cb, config, dt) integer(I4B), dimension(:), allocatable :: iflag real(DP), dimension(:), allocatable :: dtp + associate(ntp => self%nbody, & - xh => self%xh, & - vh => self%vh, & - status => self%status,& - mu => self%mu) + xh => self%xh, & + vh => self%vh, & + status => self%status,& + mu => self%mu) + if (ntp == 0) return allocate(iflag(ntp)) iflag(:) = 0 @@ -101,11 +103,13 @@ module subroutine whm_drift_tp(self, cb, config, dt) else dtp(:) = dt end if - !$omp simd + !!$omp simd do i = 1,ntp - if (status(i) == ACTIVE) call drift_one(mu(i), xh(:, i), vh(:, i), dtp(i), iflag(i)) + if (status(i) == ACTIVE) call drift_one(mu(i), xh(1, i), xh(2, i), xh(3, i), & + vh(1, i), vh(2, i), vh(3, i), dtp(i), & + iflag(i)) end do - if (any(iflag(1:ntp) /= 0)) then + if (any(iflag(:) /= 0)) then where(iflag(:) /= 0) status(:) = DISCARDED_DRIFTERR do i = 1, ntp if (iflag(i) /= 0) write(*, *) "Particle ", self%name(i), " lost due to error in Danby drift" diff --git a/src/whm/whm_getacch.f90 b/src/whm/whm_getacch.f90 index 65497c402..14f7bfabf 100644 --- a/src/whm/whm_getacch.f90 +++ b/src/whm/whm_getacch.f90 @@ -23,9 +23,9 @@ module subroutine whm_getacch_pl(self, cb, config, t) if (npl == 0) return call pl%set_ir3() - ah0 = whm_getacch_ah0(pl%Gmass(2:npl), pl%xh(:,2:npl)) - do i = 1, NDIM - pl%ah(i, 1:npl) = ah0(i) + ah0 = whm_getacch_ah0(pl%Gmass(2:npl), pl%xh(:,2:npl), npl-1) + do i = 1, npl + pl%ah(:, i) = ah0(:) end do call whm_getacch_ah1(cb, pl) call whm_getacch_ah2(cb, pl) @@ -62,7 +62,7 @@ module subroutine whm_getacch_tp(self, cb, pl, config, t, xh) ir3h => pl%ir3h, GMpl => pl%Gmass) if (ntp == 0 .or. npl == 0) return - ah0 = whm_getacch_ah0(pl%Gmass(:), xh(:,:)) + ah0 = whm_getacch_ah0(pl%Gmass(:), xh(:,:), npl) do i = 1, ntp tp%ah(:, i) = ah0(:) end do @@ -74,21 +74,20 @@ module subroutine whm_getacch_tp(self, cb, pl, config, t, xh) return end subroutine whm_getacch_tp - pure function whm_getacch_ah0(mu, xh) result(ah0) + function whm_getacch_ah0(mu, xh, n) result(ah0) !! author: David A. Minton !! !! Compute zeroth term heliocentric accelerations of planets implicit none ! Arguments - real(DP), dimension(:), intent(in) :: mu + real(DP), dimension(:), intent(in) :: mu real(DP), dimension(:,:), intent(in) :: xh + integer(I4B), intent(in) :: n ! Result real(DP), dimension(NDIM) :: ah0 ! Internals real(DP) :: fac, r2, ir3h - integer(I4B) :: i, n - - n = size(mu) + integer(I4B) :: i ah0(:) = 0.0_DP do i = 1, n @@ -221,7 +220,7 @@ pure subroutine whm_getacch_ah3_tp(cb, pl, tp, xh) if (ntp == 0) return do i = 1, ntp acc(:) = 0.0_DP - !$omp simd private(dx,rji2,irij3,fac) reduction(-:acc) + !!$omp simd private(dx,rji2,irij3,fac) reduction(-:acc) do j = 1, npl dx(:) = xht(:, i) - xh(:, j) rji2 = dot_product(dx(:), dx(:)) diff --git a/src/whm/whm_setup.f90 b/src/whm/whm_setup.f90 index da3258df5..822108527 100644 --- a/src/whm/whm_setup.f90 +++ b/src/whm/whm_setup.f90 @@ -133,9 +133,9 @@ module subroutine whm_setup_set_beg_end(self, xbeg, xend, vbeg) !! Sets one or more of the values of xbeg and xend implicit none ! Arguments - class(whm_tp), intent(inout) :: self !! Swiftest test particle object - real(DP), dimension(:,:), optional :: xbeg, xend - real(DP), dimension(:,:), optional :: vbeg ! vbeg is an unused variable to keep this method forward compatible with RMVS + class(whm_tp), intent(inout) :: self !! Swiftest test particle object + real(DP), dimension(:,:), intent(in), optional :: xbeg, xend + real(DP), dimension(:,:), intent(in), optional :: vbeg ! vbeg is an unused variable to keep this method forward compatible with RMVS if (present(xbeg)) then if (allocated(self%xbeg)) deallocate(self%xbeg) diff --git a/src/whm/whm_step.f90 b/src/whm/whm_step.f90 index 120a7194b..d7c620786 100644 --- a/src/whm/whm_step.f90 +++ b/src/whm/whm_step.f90 @@ -1,7 +1,9 @@ submodule(whm_classes) s_whm_step use swiftest contains - module subroutine whm_step_system(cb, pl, tp, config) + + + module subroutine whm_step_system(self, config) !! author: David A. Minton !! !! Step massive bodies and and active test particles ahead in heliocentric coordinates @@ -10,11 +12,15 @@ module subroutine whm_step_system(cb, pl, tp, config) !! Adapted from David E. Kaufmann's Swifter routine whm_step.f90 implicit none ! Arguments - class(whm_cb), intent(inout) :: cb !! Swiftest central body object - class(whm_pl), intent(inout) :: pl !! Swiftest central body object - class(whm_tp), intent(inout) :: tp !! Swiftest central body object + class(whm_nbody_system), intent(inout) :: self !! WHM nbody system object class(swiftest_configuration), intent(in) :: config !! Input collection of on parameters + select type(cb => self%cb) + class is (whm_cb) + select type(pl => self%pl) + class is (whm_pl) + select type(tp => self%tp) + class is (whm_tp) associate(ntp => tp%nbody, npl => pl%nbody, t => config%t, dt => config%dt) call pl%set_rhill(cb) call tp%set_beg_end(xbeg = pl%xh) @@ -24,6 +30,10 @@ module subroutine whm_step_system(cb, pl, tp, config) call tp%step(cb, pl, config, t, dt) end if end associate + end select + end select + end select + return end subroutine whm_step_system module subroutine whm_step_pl(self, cb, config, t, dt) @@ -88,7 +98,7 @@ module subroutine whm_step_tp(self, cb, pl, config, t, dt) xbeg => self%xbeg, xend => self%xend) dth = 0.5_DP * dt if (tp%lfirst) then - call tp%getacch(cb, pl, config, t, tp%xbeg) + call tp%getacch(cb, pl, config, t, xbeg) tp%lfirst = .false. end if call tp%kickvh(dth) @@ -96,7 +106,7 @@ module subroutine whm_step_tp(self, cb, pl, config, t, dt) if (config%lgr) call tp%gr_p4(config, dth) call tp%drift(cb, config, dt) if (config%lgr) call tp%gr_p4(config, dth) - call tp%getacch(cb, pl, config, t + dt, tp%xend) + call tp%getacch(cb, pl, config, t + dt, xend) call tp%kickvh(dth) end associate return