From 5d55ac1b35a108042d8b1bd410dbd4feb281a571 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Fri, 23 Jul 2021 17:26:23 -0400 Subject: [PATCH] Wrote the helio drift subroutine for swiftest_body type and then wrapped it with interfaces so that they can be used as type-bound procedures for helio_pl and helio_tp --- src/helio/helio_drift.f90 | 41 ++++++++++++++++++++++++++++++----- src/helio/helio_getacch.f90 | 5 +++-- src/kick/kick.f90 | 10 +++++---- src/modules/helio_classes.f90 | 23 +++++++++++++++++++- 4 files changed, 67 insertions(+), 12 deletions(-) diff --git a/src/helio/helio_drift.f90 b/src/helio/helio_drift.f90 index 7a6d4db06..38e156e32 100644 --- a/src/helio/helio_drift.f90 +++ b/src/helio/helio_drift.f90 @@ -1,18 +1,17 @@ submodule (helio_classes) s_helio_drift use swiftest contains - module subroutine helio_drift_pl(self, system, param, dt, mask) + module subroutine helio_drift_body(self, system, param, dt, mask) !! author: David A. Minton !! - !! Loop through massive bodies and call Danby drift routine - !! New vectorized version included + !! Loop through bodies and call Danby drift routine on democratic heliocentric coordinates !! !! Adapted from David E. Kaufmann's Swifter routine helio_drift.f90 !! Adapted from Hal Levison's Swift routine drift.f implicit none ! Arguments - class(helio_pl), intent(inout) :: self !! Helio massive body object + class(swiftest_body), intent(inout) :: self !! Swiftest 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) :: dt !! Stepsize @@ -51,7 +50,7 @@ module subroutine helio_drift_pl(self, system, param, dt, mask) if (any(iflag(1:npl) /= 0)) then do i = 1, npl if (iflag(i) /= 0) then - write(*, *) " Planet ", self%id(i), " is lost!!!!!!!!!!" + write(*, *) " Body", self%id(i), " is lost!!!!!!!!!!" write(*, *) pl%xh(:,i) write(*, *) pl%vb(:,i) write(*, *) " stopping " @@ -63,7 +62,39 @@ module subroutine helio_drift_pl(self, system, param, dt, mask) return + end subroutine helio_drift_body + + module subroutine helio_drift_pl(self, system, param, dt, mask) + !! author: David A. Minton + !! + !! Wrapper function used to call the body drift routine from a helio_pl structure + implicit none + ! Arguments + 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) :: dt !! Stepsize + logical, dimension(:), intent(in) :: mask !! Logical mask of size self%nbody that determines which bodies to drift. + + call helio_drift_body(self, system, param, dt, mask) + return end subroutine helio_drift_pl + + module subroutine helio_drift_tp(self, system, param, dt, mask) + !! author: David A. Minton + !! + !! Wrapper function used to call the body drift routine from a helio_pl structure + implicit none + ! Arguments + class(helio_tp), 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) :: dt !! Stepsize + logical, dimension(:), intent(in) :: mask !! Logical mask of size self%nbody that determines which bodies to drift. + + call helio_drift_body(self, system, param, dt, mask) + return + end subroutine helio_drift_tp module subroutine helio_drift_linear_pl(self, cb, dt, lbeg) !! author: David A. Minton diff --git a/src/helio/helio_getacch.f90 b/src/helio/helio_getacch.f90 index 9b00b698b..098549eff 100644 --- a/src/helio/helio_getacch.f90 +++ b/src/helio/helio_getacch.f90 @@ -51,9 +51,10 @@ module subroutine helio_getacch_tp(self, system, param, t, lbeg) real(DP), intent(in) :: t !! Current time logical, optional, intent(in) :: lbeg !! Optional argument that determines whether or not this is the beginning or end of the step - associate(tp => self, ntp => self%nbody, cb => system%cb, pl => system%pl, npl => system%pl%nbody) + associate(tp => self, cb => system%cb, pl => system%pl, npl => system%pl%nbody) + tp%ah(:,:) = 0.0_DP if (present(lbeg)) system%lbeg = lbeg - if (lbeg) then + if (system%lbeg) then call tp%accel_int(pl%Gmass(:), pl%xbeg(:,:), npl) else call tp%accel_int(pl%Gmass(:), pl%xend(:,:), npl) diff --git a/src/kick/kick.f90 b/src/kick/kick.f90 index 2cb6dcde8..61cd560bb 100644 --- a/src/kick/kick.f90 +++ b/src/kick/kick.f90 @@ -48,7 +48,7 @@ module pure subroutine kick_getacch_int_tp(self, GMpl, xhp, npl) integer(I4B), intent(in) :: npl !! Number of active massive bodies ! Internals integer(I4B) :: i, j - real(DP) :: rji2, irij3, fac + real(DP) :: rji2, irij3, fac, r2 real(DP), dimension(NDIM) :: dx, acc associate(tp => self, ntp => self%nbody) @@ -56,9 +56,11 @@ module pure subroutine kick_getacch_int_tp(self, GMpl, xhp, npl) acc(:) = 0.0_DP do j = 1, npl dx(:) = tp%xh(:,i) - xhp(:, j) - rji2 = dot_product(dx(:), dx(:)) - irij3 = 1.0_DP / (rji2 * sqrt(rji2)) - fac = GMpl(j) * irij3 + !rji2 = dot_product(dx(:), dx(:)) + !irij3 = 1.0_DP / (rji2 * sqrt(rji2)) + !fac = GMpl(j) * irij3 + r2 = dot_product(dx(:), dx(:)) + fac = GMpl(j) / (r2 * sqrt(r2)) acc(:) = acc(:) - fac * dx(:) end do tp%ah(:, i) = tp%ah(:, i) + acc(:) diff --git a/src/modules/helio_classes.f90 b/src/modules/helio_classes.f90 index ae2693167..17366e88f 100644 --- a/src/modules/helio_classes.f90 +++ b/src/modules/helio_classes.f90 @@ -53,6 +53,7 @@ module helio_classes 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 :: lindrift => helio_drift_linear_tp !! Method for linear drift of massive bodies due to barycentric momentum of Sun + procedure, public :: drift => helio_drift_tp !! Method for Danby drift in Democratic Heliocentric coordinates procedure, public :: accel => helio_getacch_tp !! Compute heliocentric accelerations of massive bodies procedure, public :: kick => helio_kickvb_tp !! Kicks the barycentric velocities procedure, public :: step => helio_step_tp !! Steps the body forward one stepsize @@ -84,6 +85,16 @@ module subroutine helio_coord_vh2vb_tp(self, vbcb) class(helio_tp), intent(inout) :: self !! Helio massive body object real(DP), dimension(:), intent(in) :: vbcb !! Barycentric velocity of the central body end subroutine helio_coord_vh2vb_tp + + module subroutine helio_drift_body(self, system, param, dt, mask) + use swiftest_classes, only : swiftest_body, swiftest_nbody_system, swiftest_parameters + implicit none + class(swiftest_body), intent(inout) :: self !! Swiftest 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) :: dt !! Stepsize + logical, dimension(:), intent(in) :: mask !! Logical mask of size self%nbody that determines which bodies to drift + end subroutine helio_drift_body module subroutine helio_drift_pl(self, system, param, dt, mask) use swiftest_classes, only : swiftest_nbody_system, swiftest_parameters @@ -92,9 +103,19 @@ module subroutine helio_drift_pl(self, system, param, dt, mask) 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) :: dt !! Stepsize - logical, dimension(:), intent(in) :: mask !! Logical mask of size self%nbody that determines which bodies to drift + logical, dimension(:), intent(in) :: mask !! Logical mask of size self%nbody that determines which bodies to drift end subroutine helio_drift_pl + module subroutine helio_drift_tp(self, system, param, dt, mask) + use swiftest_classes, only : swiftest_nbody_system, swiftest_parameters + implicit none + class(helio_tp), 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) :: dt !! Stepsize + 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) implicit none class(helio_pl), intent(inout) :: self !! Helio massive body object