From d85ee28f96bf09af373396acc55e51986bef55bc Mon Sep 17 00:00:00 2001 From: David A Minton Date: Mon, 26 Jul 2021 16:01:21 -0400 Subject: [PATCH] Restructured kick methods for entire code with new masked interface. WHM still works, but Helio does not. --- src/helio/helio_drift.f90 | 43 ++++++++++-------- src/helio/helio_kick.f90 | 33 ++++++++++---- src/helio/helio_step.f90 | 31 +++++-------- src/kick/kick.f90 | 26 ----------- src/modules/helio_classes.f90 | 48 +++++++++++++------- src/modules/swiftest_classes.f90 | 48 +++++++++++--------- src/modules/whm_classes.f90 | 58 +++++++++++++++++------- src/symba/symba_step.f90 | 37 ++++++--------- src/whm/whm_kick.f90 | 78 +++++++++++++++++++++++++++++++- src/whm/whm_step.f90 | 22 ++------- 10 files changed, 259 insertions(+), 165 deletions(-) diff --git a/src/helio/helio_drift.f90 b/src/helio/helio_drift.f90 index 0c146eea5..942206945 100644 --- a/src/helio/helio_drift.f90 +++ b/src/helio/helio_drift.f90 @@ -71,7 +71,7 @@ module subroutine helio_drift_tp(self, system, param, dt, mask) return end subroutine helio_drift_tp - module subroutine helio_drift_linear_pl(self, cb, dt, lbeg) + module subroutine helio_drift_linear_pl(self, cb, dt, mask, lbeg) !! author: David A. Minton !! !! Perform linear drift of massive bodies due to barycentric momentum of Sun @@ -80,21 +80,26 @@ module subroutine helio_drift_linear_pl(self, cb, dt, lbeg) !! Adapted from Hal Levison's Swift routine helio_lindrift.f implicit none ! Arguments - class(helio_pl), intent(inout) :: self !! Helio massive body object - class(helio_cb), intent(inout) :: cb !! Helio central bod - real(DP), intent(in) :: dt !! Stepsize - logical, intent(in) :: lbeg !! Argument that determines whether or not this is the beginning or end of the step + class(helio_pl), intent(inout) :: self !! Helio massive body object + class(helio_cb), intent(inout) :: cb !! Helio central body + real(DP), intent(in) :: dt !! Stepsize + logical, dimension(:), intent(in) :: mask !! Mask that determines which bodies to kick + logical, intent(in) :: lbeg !! Argument that determines whether or not this is the beginning or end of the step ! Internals - real(DP), dimension(NDIM) :: pt !! negative barycentric velocity of the central body + real(DP), dimension(NDIM) :: pt !! negative barycentric velocity of the central body + integer(I4B) :: i associate(pl => self, npl => self%nbody) - pt(1) = sum(pl%Gmass(1:npl) * pl%vb(1,1:npl)) - pt(2) = sum(pl%Gmass(1:npl) * pl%vb(2,1:npl)) - pt(3) = sum(pl%Gmass(1:npl) * pl%vb(3,1:npl)) + if (npl == 0) return + pt(1) = sum(pl%Gmass(1:npl) * pl%vb(1,1:npl), mask) + pt(2) = sum(pl%Gmass(1:npl) * pl%vb(2,1:npl), mask) + pt(3) = sum(pl%Gmass(1:npl) * pl%vb(3,1:npl), mask) pt(:) = pt(:) / cb%Gmass - pl%xh(1,1:npl) = pl%xh(1,1:npl) + pt(1) * dt - pl%xh(2,1:npl) = pl%xh(2,1:npl) + pt(2) * dt - pl%xh(3,1:npl) = pl%xh(3,1:npl) + pt(3) * dt + do concurrent(i = 1:npl, mask(i)) + pl%xh(1,1:npl) = pl%xh(1,1:npl) + pt(1) * dt + pl%xh(2,1:npl) = pl%xh(2,1:npl) + pt(2) * dt + pl%xh(3,1:npl) = pl%xh(3,1:npl) + pt(3) * dt + end do if (lbeg) then cb%ptbeg = pt(:) @@ -106,7 +111,7 @@ module subroutine helio_drift_linear_pl(self, cb, dt, lbeg) return end subroutine helio_drift_linear_pl - module subroutine helio_drift_linear_tp(self, cb, dt, lbeg) + module subroutine helio_drift_linear_tp(self, cb, dt, mask, lbeg) !! author: David A. Minton !! !! Perform linear drift of test particles due to barycentric momentum of Sun @@ -116,20 +121,22 @@ module subroutine helio_drift_linear_tp(self, cb, dt, lbeg) !! Adapted from Hal Levison's Swift routine helio_lindrift_tp.f implicit none ! Arguments - class(helio_tp), intent(inout) :: self !! Helio test particleb object - class(helio_cb), intent(in) :: cb !! Helio central body - real(DP), intent(in) :: dt !! Stepsize - logical, intent(in) :: lbeg !! Argument that determines whether or not this is the beginning or end of the step + class(helio_tp), intent(inout) :: self !! Helio test particleb object + class(helio_cb), intent(in) :: cb !! Helio central body + real(DP), intent(in) :: dt !! Stepsize + logical, dimension(:), intent(in) :: mask !! Mask that determines which bodies to kick + logical, intent(in) :: lbeg !! Argument that determines whether or not this is the beginning or end of the step ! Internals real(DP), dimension(NDIM) :: pt !! negative barycentric velocity of the central body associate(tp => self, ntp => self%nbody) + if (ntp == 0) return if (lbeg) then pt(:) = cb%ptbeg else pt(:) = cb%ptend end if - where (tp%status(1:ntp) == ACTIVE) + where (mask(1:ntp)) tp%xh(1, 1:ntp) = tp%xh(1, 1:ntp) + pt(1) * dt tp%xh(2, 1:ntp) = tp%xh(2, 1:ntp) + pt(2) * dt tp%xh(3, 1:ntp) = tp%xh(3, 1:ntp) + pt(3) * dt diff --git a/src/helio/helio_kick.f90 b/src/helio/helio_kick.f90 index a4cd86d1d..1c2fff23e 100644 --- a/src/helio/helio_kick.f90 +++ b/src/helio/helio_kick.f90 @@ -66,7 +66,7 @@ module subroutine helio_kick_getacch_tp(self, system, param, t, lbeg) return end subroutine helio_kick_getacch_tp - module subroutine helio_kick_vb_pl(self, dt) + module subroutine helio_kick_vb_pl(self, system, param, t, dt, mask, lbeg) !! author: David A. Minton !! !! Kick barycentric velocities of bodies @@ -75,14 +75,25 @@ module subroutine helio_kick_vb_pl(self, dt) !! Adapted from David E. Kaufmann's Swifter routine helio_kick_vb.f90 implicit none ! Arguments - class(helio_pl), intent(inout) :: self !! Swiftest generic body object - real(DP), intent(in) :: dt !! Stepsize + class(helio_pl), intent(inout) :: self !! Swiftest generic body object + class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object + class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters + real(DP), intent(in) :: t !! Current time + real(DP), intent(in) :: dt !! Stepsize + logical, dimension(:), intent(in) :: mask !! Mask that determines which bodies to kick + logical, intent(in) :: lbeg !! Logical flag indicating whether this is the beginning of the half step or not. ! Internals integer(I4B) :: i associate(pl => self, npl => self%nbody) if (npl ==0) return - do concurrent(i = 1:npl, pl%status(i) == ACTIVE) + call pl%accel(system, param, t) + if (lbeg) then + call pl%set_beg_end(xbeg = pl%xh) + else + call pl%set_beg_end(xend = pl%xh) + end if + do concurrent(i = 1:npl, mask(i)) pl%vb(:, i) = pl%vb(:, i) + pl%ah(:, i) * dt end do end associate @@ -91,7 +102,7 @@ module subroutine helio_kick_vb_pl(self, dt) end subroutine helio_kick_vb_pl - module subroutine helio_kick_vb_tp(self, dt) + module subroutine helio_kick_vb_tp(self, system, param, t, dt, mask, lbeg) !! author: David A. Minton !! !! Kick barycentric velocities of bodies @@ -100,14 +111,20 @@ module subroutine helio_kick_vb_tp(self, dt) !! Adapted from David E. Kaufmann's Swifter routine helio_kick_vb_tp.f90 implicit none ! Arguments - class(helio_tp), intent(inout) :: self !! Swiftest generic body object - real(DP), intent(in) :: dt !! Stepsize + class(helio_tp), intent(inout) :: self !! Swiftest generic body object + class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object + class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters + real(DP), intent(in) :: t !! Current time + real(DP), intent(in) :: dt !! Stepsize + logical, dimension(:), intent(in) :: mask !! Mask that determines which bodies to kick + logical, intent(in) :: lbeg !! Logical flag indicating whether this is the beginning of the half step or not. ! Internals integer(I4B) :: i associate(tp => self, ntp => self%nbody) if (ntp ==0) return - do concurrent(i = 1:ntp, tp%status(i) == ACTIVE) + call tp%accel(system, param, t, lbeg) + do concurrent(i = 1:ntp, mask(i)) tp%vb(:, i) = tp%vb(:, i) + tp%ah(:, i) * dt end do end associate diff --git a/src/helio/helio_step.f90 b/src/helio/helio_step.f90 index 511ffacb6..d0c4dde83 100644 --- a/src/helio/helio_step.f90 +++ b/src/helio/helio_step.f90 @@ -37,8 +37,7 @@ module subroutine helio_step_pl(self, system, param, t, dt) real(DP), intent(in) :: t !! Current simulation time real(DP), intent(in) :: dt !! Stepsize ! Internals - integer(I4B) :: i - real(DP) :: dth, msys + real(DP) :: dth !! Half step size if (self%nbody == 0) return associate(pl => self) @@ -49,15 +48,11 @@ module subroutine helio_step_pl(self, system, param, t, dt) call pl%vh2vb(cb) pl%lfirst = .false. end if - call pl%lindrift(cb, dth, lbeg=.true.) - call pl%accel(system, param, t) - call pl%kick(dth) - call pl%set_beg_end(xbeg = pl%xh) - call pl%drift(system, param, dt, pl%status(:) == ACTIVE) - call pl%set_beg_end(xend = pl%xh) - call pl%accel(system, param, t + dt) - call pl%kick(dth) - call pl%lindrift(cb, dth, lbeg=.false.) + call pl%lindrift(cb, dth, mask=(pl%status(:) == ACTIVE), lbeg=.true.) + call pl%kick(system, param, t, dth, mask=(pl%status(:) == ACTIVE), lbeg=.true.) + call pl%drift(system, param, dt, mask=(pl%status(:) == ACTIVE)) + call pl%kick(system, param, t + dt, dth, mask=(pl%status(:) == ACTIVE), lbeg=.false.) + call pl%lindrift(cb, dth, mask=(pl%status(:) == ACTIVE), lbeg=.false.) call pl%vb2vh(cb) end select end associate @@ -80,9 +75,9 @@ module subroutine helio_step_tp(self, system, param, t, dt) 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 - real(DP), intent(in) :: dt !! Stepsiz + real(DP), intent(in) :: dt !! Stepsize ! Internals - real(DP) :: dth !! Half step size + real(DP) :: dth !! Half step size if (self%nbody == 0) return @@ -94,13 +89,11 @@ module subroutine helio_step_tp(self, system, param, t, dt) call tp%vh2vb(vbcb = -cb%ptbeg) tp%lfirst = .false. end if - call tp%lindrift(cb, dth, lbeg=.true.) - call tp%accel(system, param, t, lbeg=.true.) - call tp%kick(dth) + call tp%lindrift(cb, dth, mask=(tp%status(:) == ACTIVE), lbeg=.true.) + call tp%kick(system, param, t, dth, mask=(tp%status(:) == ACTIVE), lbeg=.true.) call tp%drift(system, param, dt, tp%status(:) == ACTIVE) - call tp%accel(system, param, t + dt, lbeg=.false.) - call tp%kick(dth) - call tp%lindrift(cb, dth, lbeg=.false.) + call tp%kick(system, param, t + dt, dth, mask=(tp%status(:) == ACTIVE), lbeg=.false.) + call tp%lindrift(cb, dth, mask=(tp%status(:) == ACTIVE), lbeg=.false.) call tp%vb2vh(vbcb = -cb%ptend) end select end associate diff --git a/src/kick/kick.f90 b/src/kick/kick.f90 index f9cd81b35..c10d47dbc 100644 --- a/src/kick/kick.f90 +++ b/src/kick/kick.f90 @@ -64,30 +64,4 @@ module pure subroutine kick_getacch_int_tp(self, GMpl, xhp, npl) return end subroutine kick_getacch_int_tp - module subroutine kick_vh_body(self, dt) - !! author: David A. Minton - !! - !! Kick heliocentric 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 whm_kickvh.f90 and whm_kickvh_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, vh => self%vh, ah => self%ah, status => self%status) - if (n == 0) return - do i = 1, n - if (status(i) == ACTIVE) vh(:, i) = vh(:, i) + ah(:, i) * dt - end do - end associate - - return - end subroutine kick_vh_body - - - end submodule s_kick diff --git a/src/modules/helio_classes.f90 b/src/modules/helio_classes.f90 index 39d1e30e4..c3dc0be62 100644 --- a/src/modules/helio_classes.f90 +++ b/src/modules/helio_classes.f90 @@ -116,20 +116,22 @@ module subroutine helio_drift_tp(self, system, param, dt, mask) logical, dimension(:), intent(in) :: mask !! Logical mask of size self%nbody that determines which bodies to drift end subroutine helio_drift_tp - module subroutine helio_drift_linear_pl(self, cb, dt, lbeg) + module subroutine helio_drift_linear_pl(self, cb, dt, mask, lbeg) implicit none - class(helio_pl), intent(inout) :: self !! Helio massive body object - class(helio_cb), intent(inout) :: cb !! Helio central body object - real(DP), intent(in) :: dt !! Stepsize - logical, intent(in) :: lbeg !! Argument that determines whether or not this is the beginning or end of the step + class(helio_pl), intent(inout) :: self !! Helio massive body object + class(helio_cb), intent(inout) :: cb !! Helio central body + real(DP), intent(in) :: dt !! Stepsize + logical, dimension(:), intent(in) :: mask !! Mask that determines which bodies to kick + logical, intent(in) :: lbeg !! Argument that determines whether or not this is the beginning or end of the step end subroutine helio_drift_linear_pl - module subroutine helio_drift_linear_tp(self, cb, dt, lbeg) + module subroutine helio_drift_linear_tp(self, cb, dt, mask, lbeg) implicit none - class(helio_tp), intent(inout) :: self !! Helio test particle object - class(helio_cb), intent(in) :: cb !! Helio nbody system object - real(DP), intent(in) :: dt !! Stepsize - logical, intent(in) :: lbeg !! Argument that determines whether or not this is the beginning or end of the step + class(helio_tp), intent(inout) :: self !! Helio test particle object + class(helio_cb), intent(in) :: cb !! Helio central body + real(DP), intent(in) :: dt !! Stepsize + logical, dimension(:), intent(in) :: mask !! Mask that determines which bodies to kick + logical, intent(in) :: lbeg !! Argument that determines whether or not this is the beginning or end of the step end subroutine helio_drift_linear_tp module subroutine helio_kick_getacch_pl(self, system, param, t, lbeg) @@ -143,7 +145,7 @@ module subroutine helio_kick_getacch_pl(self, system, param, t, lbeg) end subroutine helio_kick_getacch_pl module subroutine helio_kick_getacch_tp(self, system, param, t, lbeg) - use swiftest_classes, only : swiftest_parameters, swiftest_nbody_system + use swiftest_classes, only : swiftest_nbody_system, swiftest_parameters implicit none class(helio_tp), intent(inout) :: self !! Helio test particle object class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object @@ -152,16 +154,28 @@ module subroutine helio_kick_getacch_tp(self, system, param, t, lbeg) logical, optional, intent(in) :: lbeg !! Optional argument that determines whether or not this is the beginning or end of the step end subroutine helio_kick_getacch_tp - module subroutine helio_kick_vb_pl(self, dt) + module subroutine helio_kick_vb_pl(self, system, param, t, dt, mask, lbeg) + use swiftest_classes, only : swiftest_nbody_system, swiftest_parameters implicit none - class(helio_pl), intent(inout) :: self !! Helio massive body object - real(DP), intent(in) :: dt !! Stepsize + class(helio_pl), intent(inout) :: self !! Helio massive body object + class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object + class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters + real(DP), intent(in) :: t !! Current time + real(DP), intent(in) :: dt !! Stepsize + logical, dimension(:), intent(in) :: mask !! Mask that determines which bodies to kick + logical, intent(in) :: lbeg !! Logical flag indicating whether this is the beginning of the half step or not. end subroutine helio_kick_vb_pl - module subroutine helio_kick_vb_tp(self, dt) + module subroutine helio_kick_vb_tp(self, system, param, t, dt, mask, lbeg) + use swiftest_classes, only : swiftest_nbody_system, swiftest_parameters implicit none - class(helio_tp), intent(inout) :: self !! Helio test particle object - real(DP), intent(in) :: dt !! Stepsize + 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(in) :: param !! Current run configuration parameters + real(DP), intent(in) :: t !! Current time + real(DP), intent(in) :: dt !! Stepsize + logical, dimension(:), intent(in) :: mask !! Mask that determines which bodies to kick + logical, intent(in) :: lbeg !! Logical flag indicating whether this is the beginning of the half step or not. end subroutine helio_kick_vb_tp module subroutine helio_step_pl(self, system, param, t, dt) diff --git a/src/modules/swiftest_classes.f90 b/src/modules/swiftest_classes.f90 index 4092a0a52..7c160b780 100644 --- a/src/modules/swiftest_classes.f90 +++ b/src/modules/swiftest_classes.f90 @@ -13,7 +13,7 @@ module swiftest_classes public :: io_dump_param, io_dump_swiftest, io_dump_system, io_get_args, io_get_token, 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_toupper, io_write_discard, io_write_encounter, io_write_frame_body, io_write_frame_cb, io_write_frame_system - public :: kick_getacch_int_pl, kick_vh_body + public :: kick_getacch_int_pl public :: obl_acc_body, obl_acc_pl, obl_acc_tp public :: orbel_el2xv_vec, orbel_xv2el_vec, orbel_scget, orbel_xv2aeq, orbel_xv2aqt public :: setup_body, setup_construct_system, setup_initialize_system, setup_pl, setup_tp @@ -167,6 +167,7 @@ module swiftest_classes contains private procedure(abstract_discard_body), public, deferred :: discard + procedure(abstract_kick_body), public, deferred :: kick procedure(abstract_set_mu), public, deferred :: set_mu procedure(abstract_step_body), public, deferred :: step procedure(abstract_accel), public, deferred :: accel @@ -177,7 +178,6 @@ 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 :: kick => kick_vh_body !! Kicks the heliocentric velocities procedure, public :: accel_obl => 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 @@ -216,19 +216,19 @@ module swiftest_classes private ! Massive body-specific concrete methods ! These are concrete because they are the same implemenation for all integrators - procedure, public :: discard => discard_pl !! Placeholder method for discarding massive bodies - procedure, public :: eucl_index => eucl_dist_index_plpl !! Sets up the (i, j) -> k indexing used for the single-loop blocking Euclidean distance matrix - procedure, public :: accel_int => kick_getacch_int_pl !! Compute direct cross (third) term heliocentric accelerations of massive bodies - procedure, public :: accel_obl => obl_acc_pl !! Compute the barycentric accelerations of bodies due to the oblateness of the central body - procedure, public :: setup => setup_pl !! A base constructor that sets the number of bodies and allocates and initializes all arrays - procedure, public :: accel_tides => tides_kick_getacch_pl !! Compute the accelerations of bodies due to tidal interactions with the central body - procedure, public :: set_mu => util_set_mu_pl !! Method used to construct the vectorized form of the central body mass - procedure, public :: set_rhill => util_set_rhill !! Calculates the Hill's radii for each body - procedure, public :: h2b => util_coord_h2b_pl !! Convert massive bodies from heliocentric to barycentric coordinates (position and velocity) - procedure, public :: b2h => util_coord_b2h_pl !! Convert massive bodies from barycentric to heliocentric coordinates (position and velocity) - procedure, public :: fill => util_fill_pl !! "Fills" bodies from one object into another depending on the results of a mask (uses the MERGE intrinsic) - procedure, public :: set_beg_end => util_set_beg_end_pl !! Sets the beginning and ending positions and velocities of planets. - procedure, public :: spill => util_spill_pl !! "Spills" bodies from one object to another depending on the results of a mask (uses the PACK intrinsic) + procedure, public :: discard => discard_pl !! Placeholder method for discarding massive bodies + procedure, public :: eucl_index => eucl_dist_index_plpl !! Sets up the (i, j) -> k indexing used for the single-loop blocking Euclidean distance matrix + procedure, public :: accel_int => kick_getacch_int_pl !! Compute direct cross (third) term heliocentric accelerations of massive bodies + procedure, public :: accel_obl => obl_acc_pl !! Compute the barycentric accelerations of bodies due to the oblateness of the central body + procedure, public :: setup => setup_pl !! A base constructor that sets the number of bodies and allocates and initializes all arrays + procedure, public :: accel_tides => tides_kick_getacch_pl !! Compute the accelerations of bodies due to tidal interactions with the central body + procedure, public :: set_mu => util_set_mu_pl !! Method used to construct the vectorized form of the central body mass + procedure, public :: set_rhill => util_set_rhill !! Calculates the Hill's radii for each body + procedure, public :: h2b => util_coord_h2b_pl !! Convert massive bodies from heliocentric to barycentric coordinates (position and velocity) + procedure, public :: b2h => util_coord_b2h_pl !! Convert massive bodies from barycentric to heliocentric coordinates (position and velocity) + procedure, public :: fill => util_fill_pl !! "Fills" bodies from one object into another depending on the results of a mask (uses the MERGE intrinsic) + procedure, public :: set_beg_end => util_set_beg_end_pl !! Sets the beginning and ending positions and velocities of planets. + procedure, public :: spill => util_spill_pl !! "Spills" bodies from one object to another depending on the results of a mask (uses the PACK intrinsic) end type swiftest_pl !******************************************************************************************************************************** @@ -319,6 +319,18 @@ subroutine abstract_initialize(self, param) class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters end subroutine abstract_initialize + subroutine abstract_kick_body(self, system, param, t, dt, mask, lbeg) + import swiftest_body, swiftest_nbody_system, swiftest_parameters, DP + implicit none + class(swiftest_body), intent(inout) :: self !! Swiftest generic body object + class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system objec + class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters + real(DP), intent(in) :: t !! Current time + real(DP), intent(in) :: dt !! Stepsize + logical, dimension(:), intent(in) :: mask !! Mask that determines which bodies to kick + logical, intent(in) :: lbeg !! Logical flag indicating whether this is the beginning of the half step or not. + end subroutine abstract_kick_body + subroutine abstract_read_frame(self, iu, param, form, ierr) import DP, I4B, swiftest_base, swiftest_parameters class(swiftest_base), intent(inout) :: self !! Swiftest base object @@ -617,12 +629,6 @@ module pure subroutine kick_getacch_int_tp(self, GMpl, xhp, npl) integer(I4B), intent(in) :: npl !! Number of active massive bodies end subroutine kick_getacch_int_tp - module subroutine kick_vh_body(self, dt) - implicit none - class(swiftest_body), intent(inout) :: self !! Swiftest body object - real(DP), intent(in) :: dt !! Stepsize - end subroutine kick_vh_body - module subroutine obl_acc_body(self, system) implicit none class(swiftest_body), intent(inout) :: self !! Swiftest body object diff --git a/src/modules/whm_classes.f90 b/src/modules/whm_classes.f90 index f4f98dbf3..e30bd874f 100644 --- a/src/modules/whm_classes.f90 +++ b/src/modules/whm_classes.f90 @@ -30,19 +30,20 @@ module whm_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 whm_setup_pl and whm_util_spill_pl contains - procedure, public :: h2j => whm_coord_h2j_pl !! Convert position and velcoity vectors from heliocentric to Jacobi coordinates - 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 to jacobi coordinates - procedure, public :: fill => whm_util_fill_pl !! "Fills" bodies from one object into another depending on the results of a mask (uses the MERGE intrinsic) - procedure, public :: accel => whm_kick_getacch_pl !! Compute heliocentric accelerations of massive bodies - procedure, public :: accel_gr => whm_gr_kick_getacch_pl !! Acceleration term arising from the post-Newtonian correction - procedure, public :: gr_pos_kick => whm_gr_p4_pl !! Position kick due to p**4 term in the post-Newtonian correction - procedure, public :: setup => whm_setup_pl !! Constructor method - Allocates space for number of particles + procedure, public :: h2j => whm_coord_h2j_pl !! Convert position and velcoity vectors from heliocentric to Jacobi coordinates + 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 to jacobi coordinates + procedure, public :: fill => whm_util_fill_pl !! "Fills" bodies from one object into another depending on the results of a mask (uses the MERGE intrinsic) + procedure, public :: accel => whm_kick_getacch_pl !! Compute heliocentric accelerations of massive bodies + procedure, public :: kick => whm_kick_vh_pl !! Kick heliocentric velocities of massive bodies + procedure, public :: accel_gr => whm_gr_kick_getacch_pl !! Acceleration term arising from the post-Newtonian correction + procedure, public :: gr_pos_kick => whm_gr_p4_pl !! Position kick due to p**4 term in the post-Newtonian correction + procedure, public :: setup => whm_setup_pl !! Constructor method - Allocates space for number of particles procedure, public :: set_mu => whm_util_set_mu_eta_pl !! Sets the Jacobi mass value for all massive bodies. procedure, public :: set_ir3 => whm_setup_set_ir3j !! Sets both the heliocentric and jacobi inverse radius terms (1/rj**3 and 1/rh**3) - procedure, public :: step => whm_step_pl !! Steps the body forward one stepsize - procedure, public :: spill => whm_util_spill_pl !!"Spills" bodies from one object to another depending on the results of a mask (uses the PACK intrinsic) + procedure, public :: step => whm_step_pl !! Steps the body forward one stepsize + procedure, public :: spill => whm_util_spill_pl !!"Spills" bodies from one object to another depending on the results of a mask (uses the PACK intrinsic) end type whm_pl !******************************************************************************************************************************** @@ -55,11 +56,12 @@ module whm_classes !! component list, such as whm_setup_tp and whm_util_spill_tp contains private - procedure, public :: accel => whm_kick_getacch_tp !! Compute heliocentric accelerations of test particles - procedure, public :: accel_gr => whm_gr_kick_getacch_tp !! Acceleration term arising from the post-Newtonian correction - procedure, public :: gr_pos_kick => whm_gr_p4_tp !! Position kick due to p**4 term in the post-Newtonian correction - procedure, public :: setup => whm_setup_tp !! Allocates new components of the whm class and recursively calls parent allocations - procedure, public :: step => whm_step_tp !! Steps the particle forward one stepsize + procedure, public :: accel => whm_kick_getacch_tp !! Compute heliocentric accelerations of test particles + procedure, public :: kick => whm_kick_vh_tp !! Kick heliocentric velocities of test particles + procedure, public :: accel_gr => whm_gr_kick_getacch_tp !! Acceleration term arising from the post-Newtonian correction + procedure, public :: gr_pos_kick => whm_gr_p4_tp !! Position kick due to p**4 term in the post-Newtonian correction + procedure, public :: setup => whm_setup_tp !! Allocates new components of the whm class and recursively calls parent allocations + procedure, public :: step => whm_step_tp !! Steps the particle forward one stepsize end type whm_tp !******************************************************************************************************************************** @@ -136,6 +138,30 @@ module subroutine whm_kick_getacch_tp(self, system, param, t, lbeg) logical, optional, intent(in) :: lbeg !! Optional argument that determines whether or not this is the beginning or end of the step end subroutine whm_kick_getacch_tp + module subroutine whm_kick_vh_pl(self, system, param, t, dt, mask, lbeg) + use swiftest_classes, only : swiftest_nbody_system, swiftest_parameters + implicit none + class(whm_pl), intent(inout) :: self !! WHM massive body object + class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object + class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters + real(DP), intent(in) :: t !! Current time + real(DP), intent(in) :: dt !! Stepsize + logical, dimension(:), intent(in) :: mask !! Mask that determines which bodies to kick + logical, intent(in) :: lbeg !! Logical flag indicating whether this is the beginning of the half step or not. + end subroutine whm_kick_vh_pl + + module subroutine whm_kick_vh_tp(self, system, param, t, dt, mask, lbeg) + use swiftest_classes, only : swiftest_nbody_system, swiftest_parameters + implicit none + class(whm_tp), intent(inout) :: self !! WHM test particle object + class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object + class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters + real(DP), intent(in) :: t !! Current time + real(DP), intent(in) :: dt !! Stepsize + logical, dimension(:), intent(in) :: mask !! Mask that determines which bodies to kick + logical, intent(in) :: lbeg !! Logical flag indicating whether this is the beginning of the half step or not. + end subroutine whm_kick_vh_tp + module subroutine whm_gr_kick_getacch_pl(self, param) use swiftest_classes, only : swiftest_cb, swiftest_parameters implicit none diff --git a/src/symba/symba_step.f90 b/src/symba/symba_step.f90 index 0ef7f8df8..896af6ab4 100644 --- a/src/symba/symba_step.f90 +++ b/src/symba/symba_step.f90 @@ -62,34 +62,27 @@ module subroutine symba_step_interp_system(self, param, t, dt) class is (symba_tp) select type(cb => system%cb) class is (symba_cb) + irec = -1 call pl%vh2vb(cb) - call pl%lindrift(cb, dth, lbeg=.true.) - call tp%vh2vb(vbcb = -cb%ptbeg) - call tp%lindrift(cb, dth, lbeg=.true.) + call pl%lindrift(cb, dth, mask=(pl%status(:) == ACTIVE), lbeg=.true.) + call pl%kick(system, param, t, dth, mask=(pl%status(:) == ACTIVE), lbeg=.true.) + call pl%drift(system, param, dt, mask=(pl%status(:) == ACTIVE .and. pl%levelg(:) == irec)) - call pl%set_beg_end(xbeg = pl%xh) - call pl%accel(system, param, t) - call tp%accel(system, param, t, lbeg=.true.) + call tp%vh2vb(vbcb = -cb%ptbeg) + call tp%lindrift(cb, dth, mask=(tp%status(:) == ACTIVE), lbeg=.true.) + call tp%kick(system, param, t, dth, mask=(tp%status(:) == ACTIVE), lbeg=.true.) + call tp%drift(system, param, dt, mask=(tp%status(:) == ACTIVE .and. tp%levelg(:) == irec)) - call pl%kick(dth) - call tp%kick(dth) - irec = -1 - call pl%drift(system, param, dt, pl%status(:) == ACTIVE .and. pl%levelg(:) == irec) - call tp%drift(system, param, dt, tp%status(:) == ACTIVE .and. tp%levelg(:) == irec) irec = 0 call system%recursive_step(param, irec) - call pl%set_beg_end(xend = pl%xh) - call pl%accel(system, param, t + dt) - call tp%accel(system, param, t + dt, lbeg=.false.) - - call pl%kick(dth) - call tp%kick(dth) - + call pl%kick(system, param, t, dth, mask=(pl%status(:) == ACTIVE), lbeg=.false.) call pl%vb2vh(cb) - call pl%lindrift(cb, dth, lbeg=.false.) + call pl%lindrift(cb, dth, mask=(pl%status(:) == ACTIVE), lbeg=.false.) + + call tp%kick(system, param, t, dth, mask=(tp%status(:) == ACTIVE), lbeg=.true.) call tp%vb2vh(vbcb = -cb%ptend) - call tp%lindrift(cb, dth, lbeg=.false.) + call tp%lindrift(cb, dth, mask=(tp%status(:) == ACTIVE), lbeg=.false.) end select end select end select @@ -146,8 +139,8 @@ module recursive subroutine symba_step_recur_system(self, param, ireci) call pltpenc_list%kick(system, dth, irecp, sgn) end if - call pl%drift(system, param, dtl, pl%status(:) == ACTIVE .and. pl%levelg(:) == ireci) - call tp%drift(system, param, dtl, tp%status(:) == ACTIVE .and. tp%levelg(:) == ireci) + call pl%drift(system, param, dtl, mask=(pl%status(:) == ACTIVE .and. pl%levelg(:) == ireci)) + call tp%drift(system, param, dtl, mask=(tp%status(:) == ACTIVE .and. tp%levelg(:) == ireci)) if (lencounter) call system%recursive_step(param, irecp) diff --git a/src/whm/whm_kick.f90 b/src/whm/whm_kick.f90 index af8805d47..5a0e19fd5 100644 --- a/src/whm/whm_kick.f90 +++ b/src/whm/whm_kick.f90 @@ -23,7 +23,7 @@ module subroutine whm_kick_getacch_pl(self, system, param, t, lbeg) if (npl == 0) return call pl%set_ir3() - ah0 = whm_kick_getacch_ah0(pl%Gmass(2:npl), pl%xh(:,2:npl), npl-1) + ah0(:) = whm_kick_getacch_ah0(pl%Gmass(2:npl), pl%xh(:,2:npl), npl-1) do i = 1, npl pl%ah(:, i) = ah0(:) end do @@ -179,4 +179,80 @@ pure subroutine whm_kick_getacch_ah2(cb, pl) return end subroutine whm_kick_getacch_ah2 + module subroutine whm_kick_vh_pl(self, system, param, t, dt, mask, lbeg) + !! author: David A. Minton + !! + !! Kick heliocentric velocities of massive bodies + !! + !! Adapted from Martin Duncan and Hal Levison's Swift routine kickvh.f + !! Adapted from David E. Kaufmann's Swifter routine whm_kickvh.f90 + implicit none + ! Arguments + class(whm_pl), intent(inout) :: self !! WHM massive body object + class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object + class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters + real(DP), intent(in) :: t !! Current time + real(DP), intent(in) :: dt !! Stepsize + logical, dimension(:), intent(in) :: mask !! Mask that determines which bodies to kick + logical, intent(in) :: lbeg !! Logical flag indicating whether this is the beginning of the half step or not. + ! Internals + integer(I4B) :: i + + associate(pl => self, npl => self%nbody, cb => system%cb) + if (npl == 0) return + if (pl%lfirst) then + call pl%h2j(cb) + call pl%accel(system, param, t) + pl%lfirst = .false. + end if + if (lbeg) then + call pl%set_beg_end(xbeg = pl%xh) + else + call pl%set_beg_end(xend = pl%xh) + call pl%accel(system, param, t) + end if + do concurrent(i = 1:npl, mask(i)) + pl%vh(:, i) = pl%vh(:, i) + pl%ah(:, i) * dt + end do + end associate + + return + end subroutine whm_kick_vh_pl + + module subroutine whm_kick_vh_tp(self, system, param, t, dt, mask, lbeg) + !! author: David A. Minton + !! + !! Kick heliocentric velocities of test particles + !! + !! Adapted from Martin Duncan and Hal Levison's Swift routine kickvh_tp.f + !! Adapted from David E. Kaufmann's Swifter routine whm_kickvh_tp.f90 + implicit none + ! Arguments + class(whm_tp), intent(inout) :: self !! WHM massive body object + class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object + class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters + real(DP), intent(in) :: t !! Current time + real(DP), intent(in) :: dt !! Stepsize + logical, dimension(:), intent(in) :: mask !! Mask that determines which bodies to kick + logical, intent(in) :: lbeg !! Logical flag indicating whether this is the beginning of the half step or not. + ! Internals + integer(I4B) :: i + + associate(tp => self, ntp => self%nbody) + if (ntp == 0) return + if (tp%lfirst) then + call tp%accel(system, param, t, lbeg=.true.) + tp%lfirst = .false. + end if + if (.not.lbeg) call tp%accel(system, param, t, lbeg) + do concurrent(i = 1:ntp, mask(i)) + tp%vh(:, i) = tp%vh(:, i) + tp%ah(:, i) * dt + end do + end associate + + return + end subroutine whm_kick_vh_tp + + + end submodule s_whm_kick diff --git a/src/whm/whm_step.f90 b/src/whm/whm_step.f90 index 64415f15d..ae8722ad9 100644 --- a/src/whm/whm_step.f90 +++ b/src/whm/whm_step.f90 @@ -47,20 +47,13 @@ module subroutine whm_step_pl(self, system, param, t, dt) associate(pl => self, cb => system%cb) dth = 0.5_DP * dt - if (pl%lfirst) then - call pl%h2j(cb) - call pl%accel(system, param, t) - pl%lfirst = .false. - end if - call pl%set_beg_end(xbeg = pl%xh) - call pl%kick(dth) + call pl%kick(system, param, t, dth, mask=(pl%status(:) == ACTIVE), lbeg=.true.) call pl%vh2vj(cb) if (param%lgr) call pl%gr_pos_kick(param, dth) call pl%drift(system, param, dt, pl%status(:) == ACTIVE) if (param%lgr) call pl%gr_pos_kick(param, dth) call pl%j2h(cb) - call pl%accel(system, param, t + dt) - call pl%kick(dth) + call pl%kick(system, param, t + dt, dth, mask=(pl%status(:) == ACTIVE), lbeg=.false.) call pl%set_beg_end(xend = pl%xh) end associate return @@ -89,16 +82,11 @@ module subroutine whm_step_tp(self, system, param, t, dt) class is (whm_nbody_system) associate(tp => self, cb => system%cb, pl => system%pl) dth = 0.5_DP * dt - if (tp%lfirst) then - call tp%accel(system, param, t, lbeg=.true.) - tp%lfirst = .false. - end if - call tp%kick(dth) + call tp%kick(system, param, t, dth, mask=(tp%status(:) == ACTIVE), lbeg=.true.) if (param%lgr) call tp%gr_pos_kick(param, dth) - call tp%drift(system, param, dt, tp%status(:) == ACTIVE) + call tp%drift(system, param, dt, mask=(tp%status(:) == ACTIVE)) if (param%lgr) call tp%gr_pos_kick(param, dth) - call tp%accel(system, param, t + dt, lbeg=.false.) - call tp%kick(dth) + call tp%kick(system, param, t + dt, dth, mask=(tp%status(:) == ACTIVE), lbeg=.false.) end associate end select return