diff --git a/src/helio/helio_drift.f90 b/src/helio/helio_drift.f90 index 34255f38c..5ce7438ba 100644 --- a/src/helio/helio_drift.f90 +++ b/src/helio/helio_drift.f90 @@ -56,7 +56,7 @@ module subroutine helio_drift_pl(self, system, param, dt) return end subroutine helio_drift_pl - module subroutine helio_drift_linear_pl(self, cb, dt, pt) + module subroutine helio_drift_linear_pl(self, system, dt, pt) !! author: David A. Minton !! !! Perform linear drift of massive bodies due to barycentric momentum of Sun @@ -65,23 +65,23 @@ module subroutine helio_drift_linear_pl(self, cb, dt, pt) !! Adapted from Hal Levison's Swift routine helio_lindrift.f implicit none ! Arguments - class(helio_pl), intent(inout) :: self !! Helio massive body 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 + class(helio_pl), intent(inout) :: self !! Helio massive body object + class(helio_nbody_system), intent(in) :: system !! Swiftest nbody system object + real(DP), intent(in) :: dt !! Stepsize + real(DP), dimension(:), intent(out) :: pt !! negative barycentric velocity of the central body ! Internals integer(I4B) :: i real(DP),dimension(NDIM) :: pttmp !intent(out) variables don't play nicely !with openmp's reduction for some reason - associate(npl => self%nbody, xh => self%xh, vb => self%vb, GMpl => self%Gmass, GMcb => cb%Gmass) + associate(pl => self, npl => self%nbody, cb => system%cb) pttmp(:) = 0.0_DP do i = 2, npl - pttmp(:) = pttmp(:) + GMpl(i) * vb(:,i) + pttmp(:) = pttmp(:) + pl%Gmass(i) * pl%vb(:,i) end do - pttmp(:) = pttmp(:) / GMcb + pttmp(:) = pttmp(:) / cb%Gmass do i = 2, npl - xh(:,i) = xh(:,i) + pttmp(:) * dt + pl%xh(:,i) = pl%xh(:,i) + pttmp(:) * dt end do pt(:) = pttmp(:) end associate @@ -100,13 +100,13 @@ module subroutine helio_drift_tp(self, system, param, dt) implicit none ! Arguments class(helio_tp), intent(inout) :: self !! Helio test particle object - class(swiftest_nbody_system), intent(inout) :: system !! WHM nbody system object + class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters of - real(DP), intent(in) :: dt !! Stepsiz + real(DP), intent(in) :: dt !! Stepsize ! Internals integer(I4B) :: i !! Loop counter real(DP) :: rmag, vmag2, energy - real(DP), dimension(:), allocatable :: dtp + real(DP), dimension(:), allocatable :: dtp integer(I4B), dimension(:),allocatable :: iflag !! Vectorized error code flag associate(tp => self, ntp => self%nbody) @@ -139,7 +139,7 @@ module subroutine helio_drift_tp(self, system, param, dt) return end subroutine helio_drift_tp - module subroutine helio_drift_linear_tp(self, dt, pt) + module subroutine helio_drift_linear_tp(self, system, dt, pt) !! author: David A. Minton !! !! Perform linear drift of test particles due to barycentric momentum of Sun @@ -149,15 +149,16 @@ module subroutine helio_drift_linear_tp(self, dt, pt) !! Adapted from Hal Levison's Swift routine helio_lindrift_tp.f implicit none ! Arguments - class(helio_tp), intent(inout) :: self !! Helio test particle data structure - real(DP), intent(in) :: dt !! Stepsize - real(DP), dimension(:), intent(in) :: pt !! negative barycentric velocity of the Sun + class(helio_tp), intent(inout) :: self !! Helio test particleb object + class(helio_nbody_system), intent(in) :: system !! Swiftest nbody system object + real(DP), intent(in) :: dt !! Stepsize + real(DP), dimension(:), intent(in) :: pt !! negative barycentric velocity of the central body - associate(ntp => self%nbody, xh => self%xh, status => self%status) - where (status(1:ntp) == ACTIVE) - xh(1, 1:ntp) = xh(1, 1:ntp) + pt(1) * dt - xh(2, 1:ntp) = xh(2, 1:ntp) + pt(2) * dt - xh(3, 1:ntp) = xh(3, 1:ntp) + pt(3) * dt + associate(tp => self, ntp => self%nbody) + where (tp%status(1:ntp) == ACTIVE) + 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 end where end associate diff --git a/src/helio/helio_step.f90 b/src/helio/helio_step.f90 index c639b4bee..9882be084 100644 --- a/src/helio/helio_step.f90 +++ b/src/helio/helio_step.f90 @@ -50,26 +50,29 @@ module subroutine helio_step_pl(self, system, param, t, dt) ! Internals integer(I4B) :: i real(DP) :: dth, msys - real(DP), dimension(NDIM) :: ptbeg, ptend !! TODO: Incorporate these into the tp structure + !real(DP), dimension(NDIM) :: ptbeg, ptend !! TODO: Incorporate these into the tp structure logical, save :: lfirst = .true. - - associate(pl => self, cb => system%cb) - dth = 0.5_DP * dt - if (lfirst) then - call pl%vh2vb(cb) - lfirst = .false. - end if - call pl%lindrift(cb, dth, ptbeg) - call pl%get_accel(system, param, t) - call pl%kick(dth) - call system%set_beg_end(xbeg = pl%xh) - call pl%drift(system, param, dt) - call system%set_beg_end(xend = pl%xh) - call pl%get_accel(system, param, t + dt) - call pl%kick(dth) - call pl%lindrift(cb, dth, ptend) - call pl%vb2vh(cb) - end associate + + select type(system) + class is (helio_nbody_system) + associate(pl => self, cb => system%cb, ptb => system%ptb, pte => system%pte) + dth = 0.5_DP * dt + if (lfirst) then + call pl%vh2vb(cb) + lfirst = .false. + end if + call pl%lindrift(system, dth, ptb) + call pl%get_accel(system, param, t) + call pl%kick(dth) + call system%set_beg_end(xbeg = pl%xh) + call pl%drift(system, param, dt) + call system%set_beg_end(xend = pl%xh) + call pl%get_accel(system, param, t + dt) + call pl%kick(dth) + call pl%lindrift(system, dth, pte) + call pl%vb2vh(cb) + end associate + end select return @@ -96,20 +99,20 @@ module subroutine helio_step_tp(self, system, param, t, dt) select type(system) class is (helio_nbody_system) - associate(tp => self, cb => system%cb, pl => system%pl, xbeg => system%xbeg, xend => system%xend) + associate(tp => self, cb => system%cb, pl => system%pl, xbeg => system%xbeg, xend => system%xend, ptb => system%ptb, pte => system%pte) dth = 0.5_DP * dt if (lfirst) then - call tp%vh2vb(vbcb = -tp%ptbeg) + call tp%vh2vb(vbcb = -ptb) lfirst = .false. end if - call tp%lindrift(dth, tp%ptbeg) + call tp%lindrift(system, dth, ptb) call tp%get_accel(system, param, t, xbeg) call tp%kick(dth) call tp%drift(system, param, dt) 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) + call tp%lindrift(system, dth, pte) + call tp%vb2vh(vbcb = -pte) end associate end select diff --git a/src/modules/helio_classes.f90 b/src/modules/helio_classes.f90 index 65320346c..7fc1749c8 100644 --- a/src/modules/helio_classes.f90 +++ b/src/modules/helio_classes.f90 @@ -13,6 +13,8 @@ module helio_classes ! helio_nbody_system class definitions and method interfaces !******************************************************************************************************************************** type, public, extends(whm_nbody_system) :: helio_nbody_system + real(DP), dimension(NDIM) :: ptb !! negative barycentric velocity of the central body at the beginning of time step + real(DP), dimension(NDIM) :: pte !! negative barycentric velocity of the central body at the end of time step contains private procedure, public :: initialize => helio_setup_system !! Performs Helio-specific initilization steps, @@ -110,20 +112,20 @@ module subroutine helio_drift_tp(self, system, param, dt) real(DP), intent(in) :: dt !! Stepsize end subroutine helio_drift_tp - module subroutine helio_drift_linear_pl(self, cb, dt, pt) - use swiftest_classes, only : swiftest_cb + module subroutine helio_drift_linear_pl(self, system, dt, pt) implicit none - 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 + class(helio_pl), intent(inout) :: self !! Helio massive body object + class(helio_nbody_system), intent(in) :: system !! Helio nbody system 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) + module subroutine helio_drift_linear_tp(self, system, dt, pt) implicit none - 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 + class(helio_tp), intent(inout) :: self !! Helio test particle object + class(helio_nbody_system), intent(in) :: system !! Helio nbody system 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 module subroutine helio_getacch_pl(self, system, param, t)