From e8bd4306302d388dbeb69ae4c73fa7aca148259a Mon Sep 17 00:00:00 2001 From: David A Minton Date: Wed, 7 Jul 2021 08:09:59 -0400 Subject: [PATCH] Restructured kick methods and refactored getacch to get_accel method. Created basic helio setup method --- src/drift/drift.f90 | 8 ++--- src/gr/gr.f90 | 8 ++--- src/helio/helio_getacch.f90 | 4 +-- src/helio/helio_kick.f90 | 53 +++++++++++++++++++++++++++++ src/helio/helio_setup.f90 | 14 ++++++++ src/helio/helio_step.f90 | 16 ++++----- src/kick/kick.f90 | 28 ++------------- src/modules/helio_classes.f90 | 57 +++++++++++++++++++++---------- src/modules/rmvs_classes.f90 | 8 ++--- src/modules/swiftest_classes.f90 | 23 +++++-------- src/modules/whm_classes.f90 | 17 ++++------ src/rmvs/rmvs_getacch.f90 | 2 +- src/rmvs/rmvs_setup.f90 | 58 +++++++++++++++++--------------- src/whm/whm_getacch.f90 | 4 +-- src/whm/whm_setup.f90 | 2 +- src/whm/whm_step.f90 | 16 ++++----- 16 files changed, 188 insertions(+), 130 deletions(-) create mode 100644 src/helio/helio_kick.f90 diff --git a/src/drift/drift.f90 b/src/drift/drift.f90 index 31afa8c08..8bba1a273 100644 --- a/src/drift/drift.f90 +++ b/src/drift/drift.f90 @@ -18,10 +18,10 @@ module pure elemental subroutine drift_one(mu, px, py, pz, vx, vy, vz, dt, iflag !! Adapted from Hal Levison and Martin Duncan's Swift routine drift_one.f 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), 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) + 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) ! Internals integer(I4B) :: i real(DP) :: dttmp diff --git a/src/gr/gr.f90 b/src/gr/gr.f90 index 6afc9f5ed..a778e3db2 100644 --- a/src/gr/gr.f90 +++ b/src/gr/gr.f90 @@ -11,11 +11,11 @@ subroutine gr_getaccb_ns_body(self, cb, param, agr, agr0) !! Adapted from David A. Minton's Swifter routine routine gr_getaccb_ns.f90 implicit none ! Arguments - class(swiftest_body), intent(inout) :: self - class(swiftest_cb), intent(inout) :: cb + class(swiftest_body), intent(inout) :: self + class(swiftest_cb), intent(inout) :: cb class(swiftest_parameters), intent(in) :: param - real(DP), dimension(:, :), intent(inout) :: agr - real(DP), dimension(NDIM), intent(out) :: agr0 + real(DP), dimension(:, :), intent(inout) :: agr + real(DP), dimension(NDIM), intent(out) :: agr0 ! Internals real(DP), dimension(NDIM) :: xh, vh real(DP) :: rmag, rdotv, vmag2 diff --git a/src/helio/helio_getacch.f90 b/src/helio/helio_getacch.f90 index af6ab9e4d..dd63ec392 100644 --- a/src/helio/helio_getacch.f90 +++ b/src/helio/helio_getacch.f90 @@ -27,7 +27,7 @@ module subroutine helio_getacch_pl(self, system, param, t) pl%ah(:,:) = pl%ahi(:,:) if (param%loblatecb) call pl%obl_acc(cb) if (param%lextra_force) call pl%user_getacch(system, param, t) - if (param%lgr) call pl%gr_getacch(param) + !if (param%lgr) call pl%gr_get_accel(param) end associate return @@ -61,7 +61,7 @@ module subroutine helio_getacch_tp(self, system, param, t, xhp) tp%ah(:,:) = tp%ahi(:,:) if (param%loblatecb) call tp%obl_acc(cb) if (param%lextra_force) call tp%user_getacch(system, param, t) - if (param%lgr) call tp%gr_getacch(param) + !if (param%lgr) call tp%gr_get_accel(param) end select end associate return diff --git a/src/helio/helio_kick.f90 b/src/helio/helio_kick.f90 new file mode 100644 index 000000000..9d5cea3a6 --- /dev/null +++ b/src/helio/helio_kick.f90 @@ -0,0 +1,53 @@ +submodule(helio_classes) s_helio_kick + use swiftest +contains + module subroutine helio_kickvb_pl(self, dt) + !! author: David A. Minton + !! + !! Kick barycentric velocities of bodies + !! + !! Adapted from Martin Duncan and Hal Levison's Swift routine kickvh.f + !! Adapted from David E. Kaufmann's Swifter routine helio_kickvb.f90 + implicit none + ! Arguments + class(helio_pl), intent(inout) :: self !! Swiftest generic body object + real(DP), intent(in) :: dt !! Stepsize + ! Internals + integer(I4B) :: i + + associate(pl => self, npl => self%nbody) + if (npl ==0) return + do concurrent(i = 1:npl, pl%status(i) == ACTIVE) + pl%vb(:, i) = pl%vb(:, i) + pl%ah(:, i) * dt + end do + end associate + + return + + end subroutine helio_kickvb_pl + + module subroutine helio_kickvb_tp(self, dt) + !! author: David A. Minton + !! + !! Kick barycentric velocities of bodies + !! + !! Adapted from Martin Duncan and Hal Levison's Swift routine kickvh_tp.f + !! Adapted from David E. Kaufmann's Swifter routine helio_kickvb_tp.f90 + implicit none + ! Arguments + class(helio_tp), intent(inout) :: self !! Swiftest generic body object + real(DP), intent(in) :: dt !! Stepsize + ! Internals + integer(I4B) :: i + + associate(tp => self, ntp => self%nbody) + if (ntp ==0) return + do concurrent(i = 1:ntp, tp%status(i) == ACTIVE) + tp%vb(:, i) = tp%vb(:, i) + tp%ah(:, i) * dt + end do + end associate + + return + + end subroutine helio_kickvb_tp +end submodule s_helio_kick \ No newline at end of file diff --git a/src/helio/helio_setup.f90 b/src/helio/helio_setup.f90 index b31b43b26..6a1f99f2a 100644 --- a/src/helio/helio_setup.f90 +++ b/src/helio/helio_setup.f90 @@ -1,6 +1,20 @@ submodule(helio_classes) s_helio_setup use swiftest contains + module subroutine helio_setup_system(self, param) + !! author: David A. Minton + !! + !! Initialize a Helio nbody system from files + implicit none + ! Arguments + class(helio_nbody_system), intent(inout) :: self !! Helio system object + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters + + call io_read_initialize_system(self, param) + ! Make sure that the discard list gets allocated initially + call self%tp_discards%setup(self%tp%nbody) + end subroutine helio_setup_system + module procedure helio_setup_pl !! author: David A. Minton & Carlisle A. Wishard !! diff --git a/src/helio/helio_step.f90 b/src/helio/helio_step.f90 index 53052fc9e..be1b2b71c 100644 --- a/src/helio/helio_step.f90 +++ b/src/helio/helio_step.f90 @@ -64,12 +64,12 @@ module subroutine helio_step_pl(self, system, param, t, dt) lfirst = .false. end if call pl%lindrift(cb, dth, ptbeg) - call pl%getacch(system, param, t) - call pl%kickvb(dth) + call pl%get_accel(system, param, t) + call pl%kick(dth) call pl%drift(system, param, dt) - call pl%getacch(system, param, t + dt) - call pl%kickvb(dth) + call pl%get_accel(system, param, t + dt) + call pl%kick(dth) call pl%lindrift(cb, dth, ptend) call pl%vb2vh(cb) end associate @@ -106,11 +106,11 @@ module subroutine helio_step_tp(self, system, param, t, dt) lfirst = .false. end if call tp%lindrift(dth, tp%ptbeg) - call tp%getacch(system, param, t, xbeg) - call tp%kickvb(dth) + call tp%get_accel(system, param, t, xbeg) + call tp%kick(dth) call tp%drift(system, param, dt) - call tp%getacch(system, param, t + dt, xend) - call tp%kickvb(dth) + call tp%get_accel(system, param, t + dt, xend) + call tp%kick(dth) call tp%lindrift(dth, tp%ptend) call tp%vb2vh(vbcb = -tp%ptend) end associate diff --git a/src/kick/kick.f90 b/src/kick/kick.f90 index 97a2b2bab..a54677dcd 100644 --- a/src/kick/kick.f90 +++ b/src/kick/kick.f90 @@ -1,7 +1,7 @@ submodule(swiftest_classes) s_kick use swiftest contains - module subroutine kick_vh_body(self, dt) + module subroutine kickvh_body(self, dt) !! author: David A. Minton !! !! Kick heliocentric velocities of bodies @@ -23,30 +23,6 @@ module subroutine kick_vh_body(self, dt) end associate return - end subroutine kick_vh_body + end subroutine kickvh_body - module subroutine kick_vb_body(self, dt) - !! author: David A. Minton - !! - !! Kick barycentric velocities of bodies - !! - !! Adapted from Martin Duncan and Hal Levison's Swift routine kickvh.f and kickvh_tp.f - !! Adapted from David E. Kaufmann's Swifter routine helio_kickvb.f90 and helio_kickvb_tp.f90 - implicit none - ! Arguments - class(swiftest_body), intent(inout) :: self !! Swiftest generic body object - real(DP), intent(in) :: dt !! Stepsize - ! Internals - integer(I4B) :: i - - associate(n => self%nbody, vb => self%vb, ah => self%ah, status => self%status) - if (n ==0) return - do concurrent(i = 1:n, status(i) == ACTIVE) - vb(:, i) = vb(:, i) + ah(:, i) * dt - end do - end associate - - return - - end subroutine kick_vb_body end submodule s_kick diff --git a/src/modules/helio_classes.f90 b/src/modules/helio_classes.f90 index ad10018a2..65320346c 100644 --- a/src/modules/helio_classes.f90 +++ b/src/modules/helio_classes.f90 @@ -4,16 +4,18 @@ module helio_classes !! Definition of classes and methods specific to the Democratic Heliocentric Method !! Adapted from David E. Kaufmann's Swifter routine: helio.f90 use swiftest_globals - use rmvs_classes, only : rmvs_cb, rmvs_pl, rmvs_tp, rmvs_nbody_system + use swiftest_classes, only : swiftest_cb, swiftest_pl, swiftest_tp + use whm_classes, only : whm_nbody_system implicit none !******************************************************************************************************************************** ! helio_nbody_system class definitions and method interfaces !******************************************************************************************************************************** - type, public, extends(rmvs_nbody_system) :: helio_nbody_system + type, public, extends(whm_nbody_system) :: helio_nbody_system contains private + procedure, public :: initialize => helio_setup_system !! Performs Helio-specific initilization steps, procedure, public :: step => helio_step_system end type helio_nbody_system @@ -21,7 +23,7 @@ module helio_classes ! helio_cb class definitions and method interfaces !******************************************************************************************************************************* !> Helio central body particle class - type, public, extends(rmvs_cb) :: helio_cb + type, public, extends(swiftest_cb) :: helio_cb contains end type helio_cb @@ -30,14 +32,15 @@ module helio_classes !******************************************************************************************************************************* !! Helio massive body particle class - type, public, extends(rmvs_pl) :: helio_pl + type, public, extends(swiftest_pl) :: helio_pl real(DP), dimension(:,:), allocatable :: ahi !! heliocentric acceleration due to interactions contains procedure, public :: vh2vb => helio_coord_vh2vb_pl !! Convert massive bodies from heliocentric to barycentric coordinates (velocity only) procedure, public :: vb2vh => helio_coord_vb2vh_pl !! Convert massive bodies from barycentric to heliocentric coordinates (velocity only) procedure, public :: drift => helio_drift_pl !! Method for Danby drift in Democratic Heliocentric coordinates procedure, public :: lindrift => helio_drift_linear_pl !! Method for linear drift of massive bodies due to barycentric momentum of Sun - procedure, public :: getacch => helio_getacch_pl !! Compute heliocentric accelerations of massive bodies + procedure, public :: get_accel => helio_getacch_pl !! Compute heliocentric accelerations of massive bodies + procedure, public :: kick => helio_kickvb_pl !! Kicks the barycentric velocities procedure, public :: setup => helio_setup_pl !! Constructor method - Allocates space for number of particles procedure, public :: step => helio_step_pl !! Steps the body forward one stepsize end type helio_pl @@ -47,18 +50,19 @@ module helio_classes !******************************************************************************************************************************* !! Helio test particle class - type, public, extends(rmvs_tp) :: helio_tp + type, public, extends(swiftest_tp) :: helio_tp real(DP), dimension(:,:), allocatable :: ahi !! heliocentric acceleration due to interactions real(DP), dimension(NDIM) :: ptbeg !! negative barycentric velocity of the Sun at beginning of time step real(DP), dimension(NDIM) :: ptend !! negative barycentric velocity of the Sun at beginning of time step contains - procedure, public :: vh2vb => helio_coord_vh2vb_tp !! Convert test particles from heliocentric to barycentric coordinates (velocity only) - procedure, public :: vb2vh => helio_coord_vb2vh_tp !! Convert test particles from barycentric to heliocentric coordinates (velocity only) + procedure, public :: vh2vb => helio_coord_vh2vb_tp !! Convert test particles from heliocentric to barycentric coordinates (velocity only) + procedure, public :: vb2vh => helio_coord_vb2vh_tp !! Convert test particles from barycentric to heliocentric coordinates (velocity only) procedure, public :: drift => helio_drift_tp !! Method for Danby drift in Democratic Heliocentric coordinates procedure, public :: lindrift => helio_drift_linear_tp !! Method for linear drift of massive bodies due to barycentric momentum of Sun - procedure, public :: getacch => helio_getacch_tp !! Compute heliocentric accelerations of massive bodies + procedure, public :: get_accel => helio_getacch_tp !! Compute heliocentric accelerations of massive bodies + procedure, public :: kick => helio_kickvb_tp !! Kicks the barycentric velocities procedure, public :: setup => helio_setup_tp !! Constructor method - Allocates space for number of particles - procedure, public :: step => helio_step_tp !! Steps the body forward one stepsize + procedure, public :: step => helio_step_tp !! Steps the body forward one stepsize end type helio_tp interface @@ -109,15 +113,15 @@ end subroutine helio_drift_tp module subroutine helio_drift_linear_pl(self, cb, dt, pt) use swiftest_classes, only : swiftest_cb implicit none - class(helio_pl), intent(inout) :: self !! Helio test particle data structure - class(swiftest_cb), intent(in) :: cb !! Helio central body data structure + class(helio_pl), intent(inout) :: self !! Helio test particle object + class(swiftest_cb), intent(in) :: cb !! Helio central body object real(DP), intent(in) :: dt !! Stepsize real(DP), dimension(:), intent(out) :: pt !! negative barycentric velocity of the central body end subroutine helio_drift_linear_pl module subroutine helio_drift_linear_tp(self, dt, pt) implicit none - class(helio_tp), intent(inout) :: self !! Helio test particle data structure + class(helio_tp), intent(inout) :: self !! Helio test particle object real(DP), intent(in) :: dt !! Stepsize real(DP), dimension(:), intent(in) :: pt !! negative barycentric velocity of the Sun end subroutine helio_drift_linear_tp @@ -125,7 +129,7 @@ end subroutine helio_drift_linear_tp module subroutine helio_getacch_pl(self, system, param, t) use swiftest_classes, only : swiftest_parameters, swiftest_nbody_system implicit none - class(helio_pl), intent(inout) :: self !! Helio massive body particle data structure + class(helio_pl), intent(inout) :: self !! Helio massive body object class(swiftest_nbody_system), intent(inout) :: system !! WHM nbody system object class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters of real(DP), intent(in) :: t !! Current simulation time @@ -134,19 +138,38 @@ end subroutine helio_getacch_pl module subroutine helio_getacch_tp(self, system, param, t, xhp) use swiftest_classes, only : swiftest_parameters, swiftest_nbody_system implicit none - class(helio_tp), intent(inout) :: self !! Helio test particle data structure + class(helio_tp), intent(inout) :: self !! Helio test particle object class(swiftest_nbody_system), intent(inout) :: system !! WHM nbody system object class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters real(DP), intent(in) :: t !! Current time real(DP), dimension(:,:), intent(in) :: xhp !! Heliocentric positions of planets end subroutine helio_getacch_tp + module subroutine helio_kickvb_pl(self, dt) + implicit none + class(helio_pl), intent(inout) :: self !! Helio massive body object + real(DP), intent(in) :: dt !! Stepsize + end subroutine helio_kickvb_pl + + module subroutine helio_kickvb_tp(self, dt) + implicit none + class(helio_tp), intent(inout) :: self !! Helio test particle object + real(DP), intent(in) :: dt !! Stepsize + end subroutine helio_kickvb_tp + module subroutine helio_setup_pl(self, n) implicit none class(helio_pl), intent(inout) :: self !! Helio massive body object integer, intent(in) :: n !! Number of test particles to allocate end subroutine helio_setup_pl + module subroutine helio_setup_system(self, param) + use swiftest_classes, only : swiftest_parameters + implicit none + class(helio_nbody_system), intent(inout) :: self !! Helio system object + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters + end subroutine helio_setup_system + module subroutine helio_setup_tp(self,n) implicit none class(helio_tp), intent(inout) :: self !! Helio test particle object @@ -165,7 +188,7 @@ end subroutine helio_step_system module subroutine helio_step_pl(self, system, param, t, dt) use swiftest_classes, only : swiftest_nbody_system, swiftest_parameters implicit none - class(helio_pl), intent(inout) :: self !! Helio massive body particle data structure + class(helio_pl), intent(inout) :: self !! Helio massive body particle object class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nboody system class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters real(DP), intent(in) :: t !! Current simulation time @@ -175,7 +198,7 @@ end subroutine helio_step_pl module subroutine helio_step_tp(self, system, param, t, dt) use swiftest_classes, only : swiftest_cb, swiftest_parameters implicit none - class(helio_tp), intent(inout) :: self !! Helio test particle data structure + class(helio_tp), intent(inout) :: self !! Helio test particle object class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters real(DP), intent(in) :: t !! Current simulation time diff --git a/src/modules/rmvs_classes.f90 b/src/modules/rmvs_classes.f90 index e595af5fe..136bd3f37 100644 --- a/src/modules/rmvs_classes.f90 +++ b/src/modules/rmvs_classes.f90 @@ -27,9 +27,9 @@ module rmvs_classes 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 - 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 :: initialize => rmvs_setup_system !! Performs RMVS-specific initilization steps, including generating the close encounter planetocentric structures + procedure, public :: step => rmvs_step_system !! Advance the RMVS nbody system forward in time by one step + 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. end type rmvs_nbody_system type, private :: rmvs_interp @@ -71,7 +71,7 @@ module rmvs_classes procedure, public :: discard => rmvs_discard_tp !! Check to see if test particles should be discarded based on pericenter passage distances with respect to planets encountered 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_tp !! Calculates either the standard or modified version of the acceleration depending if the + procedure, public :: get_accel => 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 :: 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) diff --git a/src/modules/swiftest_classes.f90 b/src/modules/swiftest_classes.f90 index ad0a17598..be9c83038 100644 --- a/src/modules/swiftest_classes.f90 +++ b/src/modules/swiftest_classes.f90 @@ -9,10 +9,10 @@ module swiftest_classes public :: discard_pl, discard_system, discard_tp public :: drift_one public :: eucl_dist_index_plpl, eucl_dist_index_pltp, eucl_irij3_plpl - public :: kick_vb_body, kick_vh_body public :: io_dump_param, io_dump_swiftest, io_dump_system, io_get_args, io_param_reader, io_param_writer, io_read_body_in, & io_read_cb_in, io_read_param_in, io_read_frame_body, io_read_frame_cb, io_read_frame_system, io_read_initialize_system, & io_write_discard, io_write_encounter, io_write_frame_body, io_write_frame_cb, io_write_frame_system + public :: kickvh_body public :: obl_acc_body public :: orbel_el2xv_vec, orbel_xv2el_vec, orbel_scget, orbel_xv2aeq, orbel_xv2aqt public :: setup_body, setup_construct_system, setup_pl, setup_set_ir3h, setup_set_msys, setup_set_mu_pl, setup_set_mu_tp, & @@ -164,8 +164,7 @@ module swiftest_classes procedure, public :: initialize => io_read_body_in !! Read in body initial conditions from a file procedure, public :: read_frame => io_read_frame_body !! I/O routine for writing out a single frame of time-series data for the central body procedure, public :: write_frame => io_write_frame_body !! I/O routine for writing out a single frame of time-series data for the central body - procedure, public :: kickvb => kick_vb_body !! Kicks the barycentric velocities - procedure, public :: kickvh => kick_vh_body !! Kicks the heliocentric velocities + procedure, public :: kick => kickvh_body !! Kicks the heliocentric velocities procedure, public :: obl_acc => obl_acc_body !! Compute the barycentric accelerations of bodies due to the oblateness of the central body procedure, public :: el2xv => orbel_el2xv_vec !! Convert orbital elements to position and velocity vectors procedure, public :: xv2el => orbel_xv2el_vec !! Convert position and velocity vectors to orbital elements @@ -383,18 +382,6 @@ module subroutine eucl_irij3_plpl(self) class(swiftest_pl), intent(inout) :: self !! Swiftest massive body object end subroutine eucl_irij3_plpl - module subroutine kick_vb_body(self, dt) - implicit none - class(swiftest_body), intent(inout) :: self !! Swiftest generic body object - real(DP), intent(in) :: dt !! Stepsize - end subroutine kick_vb_body - - module subroutine kick_vh_body(self, dt) - implicit none - class(swiftest_body), intent(inout) :: self !! Swiftest generic body object - real(DP), intent(in) :: dt !! Stepsize - end subroutine kick_vh_body - module subroutine io_dump_param(self, param_file_name) implicit none class(swiftest_parameters),intent(in) :: self !! Output collection of parameters @@ -532,6 +519,12 @@ module subroutine io_write_frame_system(self, iu, param) class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters end subroutine io_write_frame_system + module subroutine kickvh_body(self, dt) + implicit none + class(swiftest_body), intent(inout) :: self !! Swiftest generic body object + real(DP), intent(in) :: dt !! Stepsize + end subroutine kickvh_body + module subroutine obl_acc_body(self, cb) implicit none class(swiftest_body), intent(inout) :: self !! Swiftest generic body object diff --git a/src/modules/whm_classes.f90 b/src/modules/whm_classes.f90 index 1707ef875..52f87249b 100644 --- a/src/modules/whm_classes.f90 +++ b/src/modules/whm_classes.f90 @@ -5,11 +5,9 @@ module whm_classes !! Partially adapted from David E. Kaufmann's Swifter module: module_whm.f90 use swiftest_globals use swiftest_classes, only : swiftest_cb, swiftest_pl, swiftest_tp, swiftest_nbody_system - implicit none public - !******************************************************************************************************************************** ! whm_cb class definitions and method interfaces !******************************************************************************************************************************* @@ -37,9 +35,9 @@ module whm_classes procedure, public :: j2h => whm_coord_j2h_pl !! Convert position and velcoity vectors from Jacobi to helliocentric coordinates procedure, public :: vh2vj => whm_coord_vh2vj_pl !! Convert velocity vectors from heliocentric to Jacobi coordinates procedure, public :: drift => whm_drift_pl !! Loop through massive bodies and call Danby drift routine - procedure, public :: getacch => whm_getacch_pl !! Compute heliocentric accelerations of massive bodies procedure, public :: fill => whm_fill_pl !! "Fills" bodies from one object into another depending on the results of a mask (uses the MERGE intrinsic) - procedure, public :: gr_getacch => whm_gr_getacch_pl !! Acceleration term arising from the post-Newtonian correction + procedure, public :: get_accel => whm_getacch_pl !! Compute heliocentric accelerations of massive bodies + procedure, public :: gr_get_accel => whm_gr_getacch_pl !! Acceleration term arising from the post-Newtonian correction procedure, public :: gr_p4 => whm_gr_p4_pl !! Position kick due to p**4 term in the post-Newtonian correction procedure, public :: gr_vh2pv => whm_gr_vh2pv_pl !! Converts from heliocentric velocity to psudeovelocity for GR calculations procedure, public :: gr_pv2vh => whm_gr_pv2vh_pl !! Converts from psudeovelocity to heliocentric velocity for GR calculations @@ -62,8 +60,8 @@ module whm_classes contains private procedure, public :: drift => whm_drift_tp !! Loop through test particles and call Danby drift routine - procedure, public :: getacch => whm_getacch_tp !! Compute heliocentric accelerations of test particles - procedure, public :: gr_getacch => whm_gr_getacch_tp !! Acceleration term arising from the post-Newtonian correction + procedure, public :: get_accel => whm_getacch_tp !! Compute heliocentric accelerations of test particles + procedure, public :: gr_get_accel => whm_gr_getacch_tp !! Acceleration term arising from the post-Newtonian correction procedure, public :: gr_p4 => whm_gr_p4_tp !! Position kick due to p**4 term in the post-Newtonian correction procedure, public :: gr_vh2pv => whm_gr_vh2pv_tp !! Converts from heliocentric velocity to psudeovelocity for GR calculations procedure, public :: gr_pv2vh => whm_gr_pv2vh_tp !! Converts from psudeovelocity to heliocentric velocity for GR calculations @@ -81,9 +79,9 @@ module whm_classes 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 - procedure, public :: set_beg_end => whm_setup_set_beg_end !! Sets the beginning and ending positions of planets. + procedure, public :: initialize => whm_setup_system !! Performs WHM-specific initilization steps, like calculating the Jacobi masses + procedure, public :: step => whm_step_system !! Advance the WHM nbody system forward in time by one step + procedure, public :: set_beg_end => whm_setup_set_beg_end !! Sets the beginning and ending positions of planets. end type whm_nbody_system interface @@ -213,7 +211,6 @@ module pure subroutine whm_gr_vh2pv_tp(self, param) class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters of on parameters end subroutine whm_gr_vh2pv_tp - !> Reads WHM massive body object in from file module subroutine whm_setup_pl(self,n) implicit none diff --git a/src/rmvs/rmvs_getacch.f90 b/src/rmvs/rmvs_getacch.f90 index dcbf5aff9..58cccbf3e 100644 --- a/src/rmvs/rmvs_getacch.f90 +++ b/src/rmvs/rmvs_getacch.f90 @@ -64,7 +64,7 @@ module subroutine rmvs_getacch_tp(self, system, param, t, xhp) end if if (param%lextra_force) call tp%user_getacch(system, param, t) - if (param%lgr) call tp%gr_getacch(param) + if (param%lgr) call tp%gr_get_accel(param) tp%xh(:,:) = xh_original(:,:) end associate diff --git a/src/rmvs/rmvs_setup.f90 b/src/rmvs/rmvs_setup.f90 index 9d5916f56..c5e126a74 100644 --- a/src/rmvs/rmvs_setup.f90 +++ b/src/rmvs/rmvs_setup.f90 @@ -46,37 +46,10 @@ module subroutine rmvs_setup_pl(self,n) return end subroutine rmvs_setup_pl - module subroutine rmvs_setup_tp(self,n) - !! author: David A. Minton - !! - !! Allocate WHM test particle structure - !! - !! Equivalent in functionality to David E. Kaufmann's Swifter routine whm_setup.f90 - implicit none - ! Arguments - class(rmvs_tp), intent(inout) :: self !! RMVS test particle object - integer, intent(in) :: n !! Number of test particles to allocate - - !> Call allocation method for parent class - call whm_setup_tp(self, n) - if (n <= 0) return - - 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. - - return - end subroutine rmvs_setup_tp - module subroutine rmvs_setup_system(self, param) !! author: David A. Minton !! - !! Wrapper method to initialize a basic Swiftest nbody system from files. + !! nitialize an RMVS nbody system from files and sets up the planetocentric structures. !! !! 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 @@ -144,6 +117,35 @@ module subroutine rmvs_setup_system(self, param) end select end subroutine rmvs_setup_system + + module subroutine rmvs_setup_tp(self,n) + !! author: David A. Minton + !! + !! Allocate WHM test particle structure + !! + !! Equivalent in functionality to David E. Kaufmann's Swifter routine whm_setup.f90 + implicit none + ! Arguments + class(rmvs_tp), intent(inout) :: self !! RMVS test particle object + integer, intent(in) :: n !! Number of test particles to allocate + + !> Call allocation method for parent class + call whm_setup_tp(self, n) + if (n <= 0) return + + 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. + + return + end subroutine rmvs_setup_tp + + module subroutine rmvs_setup_set_beg_end(self, xbeg, xend, vbeg) !! author: David A. Minton diff --git a/src/whm/whm_getacch.f90 b/src/whm/whm_getacch.f90 index 67ecc7487..595c4546b 100644 --- a/src/whm/whm_getacch.f90 +++ b/src/whm/whm_getacch.f90 @@ -33,7 +33,7 @@ module subroutine whm_getacch_pl(self, system, param, t) if (param%loblatecb) call pl%obl_acc(cb) if (param%lextra_force) call pl%user_getacch(system, param, t) - if (param%lgr) call pl%gr_getacch(param) + if (param%lgr) call pl%gr_get_accel(param) end associate return @@ -67,7 +67,7 @@ module subroutine whm_getacch_tp(self, system, param, t, xhp) call whm_getacch_ah3_tp(system, xhp) if (param%loblatecb) call tp%obl_acc(cb) if (param%lextra_force) call tp%user_getacch(system, param, t) - if (param%lgr) call tp%gr_getacch(param) + if (param%lgr) call tp%gr_get_accel(param) end associate return end subroutine whm_getacch_tp diff --git a/src/whm/whm_setup.f90 b/src/whm/whm_setup.f90 index e0812d00d..80c0e81c9 100644 --- a/src/whm/whm_setup.f90 +++ b/src/whm/whm_setup.f90 @@ -75,7 +75,7 @@ end subroutine whm_setup_set_mu_eta_pl module subroutine whm_setup_system(self, param) !! author: David A. Minton !! - !! Wrapper method to initialize a basic Swiftest nbody system from files + !! Initialize a WHM nbody system from files !! implicit none ! Arguments diff --git a/src/whm/whm_step.f90 b/src/whm/whm_step.f90 index 8e87796ea..411e8b59c 100644 --- a/src/whm/whm_step.f90 +++ b/src/whm/whm_step.f90 @@ -50,18 +50,18 @@ module subroutine whm_step_pl(self, system, param, t, dt) dth = 0.5_DP * dt if (pl%lfirst) then call pl%h2j(cb) - call pl%getacch(system, param, t) + call pl%get_accel(system, param, t) pl%lfirst = .false. end if - call pl%kickvh(dth) + call pl%kick(dth) call pl%vh2vj(cb) !If GR enabled, calculate the p4 term before and after each drift if (param%lgr) call pl%gr_p4(param, dth) call pl%drift(system, param, dt) if (param%lgr) call pl%gr_p4(param, dth) call pl%j2h(cb) - call pl%getacch(system, param, t + dt) - call pl%kickvh(dth) + call pl%get_accel(system, param, t + dt) + call pl%kick(dth) end associate return end subroutine whm_step_pl @@ -88,16 +88,16 @@ module subroutine whm_step_tp(self, system, param, t, dt) associate(tp => self, cb => system%cb, pl => system%pl, xbeg => system%xbeg, xend => system%xend) dth = 0.5_DP * dt if (tp%lfirst) then - call tp%getacch(system, param, t, xbeg) + call tp%get_accel(system, param, t, xbeg) tp%lfirst = .false. end if - call tp%kickvh(dth) + call tp%kick(dth) !If GR enabled, calculate the p4 term before and after each drift if (param%lgr) call tp%gr_p4(param, dth) call tp%drift(system, param, dt) if (param%lgr) call tp%gr_p4(param, dth) - call tp%getacch(system, param, t + dt, xend) - call tp%kickvh(dth) + call tp%get_accel(system, param, t + dt, xend) + call tp%kick(dth) end associate end select return