From 294cdbc7c91ffdfa614a7ee6d335e3a8e91e4626 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Fri, 2 Jul 2021 17:51:16 -0400 Subject: [PATCH] Fixed Helio interfaces --- src/helio/helio_step.f90 | 105 +++++++++++++++++----------------- src/modules/helio_classes.f90 | 66 ++++++++++----------- 2 files changed, 87 insertions(+), 84 deletions(-) diff --git a/src/helio/helio_step.f90 b/src/helio/helio_step.f90 index ca6df15fa..480de821d 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(self, param) + module subroutine helio_step_system(self, param, t, dt) !! author: David A. Minton !! !! Step massive bodies and and active test particles ahead in heliocentric coordinates @@ -10,31 +10,34 @@ module subroutine helio_step_system(self, param) !! Adapted from David E. Kaufmann's Swifter routine helio_step.f90 implicit none ! Arguments - class(helio_nbody_system), intent(inout) :: self !! Helio nbody system object - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters of on parameters + class(helio_nbody_system), intent(inout) :: self !! Helio nbody system object + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters + real(DP), intent(in) :: t !! Simulation time + real(DP), intent(in) :: dt !! Current stepsize + select type(system => self) + 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) - associate(ntp => tp%nbody, npl => pl%nbody, t => param%t, dt => param%dt) call pl%set_rhill(cb) call tp%set_beg_end(xbeg = pl%xh) - call pl%step(cb, param, t, dt) + call pl%step(system, param, t, dt) if (ntp > 0) then call tp%set_beg_end(xend = pl%xh) - call tp%step(cb, pl, param, t, dt) + call tp%step(system, param, t, dt) end if - end associate + end select end select end select end select return end subroutine helio_step_system - module subroutine helio_step_pl(self, cb, param, t, dt) + module subroutine helio_step_pl(self, system, param, t, dt) !! author: David A. Minton !! !! Step massive bodies ahead Democratic Heliocentric method @@ -43,38 +46,39 @@ module subroutine helio_step_pl(self, cb, param, t, dt) !! Adapted from Hal Levison's Swift routine helio_step_pl.f implicit none ! Arguments - class(helio_pl), intent(inout) :: self !! WHM massive body particle data structure - class(swiftest_cb), intent(inout) :: cb !! Helio central body particle data structure - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters of - real(DP), intent(in) :: t !! Current time - real(DP), intent(in) :: dt !! Stepsize + class(helio_pl), intent(inout) :: self !! Helio massive body particle data structure + 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 !! Stepsize ! Internals - integer(I4B) :: i,npl + integer(I4B) :: i real(DP) :: dth, msys real(DP), dimension(NDIM) :: ptbeg, ptend !! TODO: Incorporate these into the tp structure logical, save :: lfirst = .true. - - npl = self%nbody - dth = 0.5_DP * dt - if (lfirst) then - call self%vh2vb(cb) - lfirst = .false. - end if - call self%lindrift(cb, dth, ptbeg) - call self%getacch(cb, param, t) - call self%kickvb(dth) + + associate(cb => system%cb) + dth = 0.5_DP * dt + if (lfirst) then + call self%vh2vb(cb) + lfirst = .false. + end if + call self%lindrift(cb, dth, ptbeg) + call self%getacch(system, param, t) + call self%kickvb(dth) - call self%drift(cb, param, dt) - call self%getacch(cb, param, t + dt) - call self%kickvb(dth) - call self%lindrift(cb, dth, ptend) - call self%vb2vh(cb) + call self%drift(system, param, dt) + call self%getacch(system, param, t + dt) + call self%kickvb(dth) + call self%lindrift(cb, dth, ptend) + call self%vb2vh(cb) + end associate return end subroutine helio_step_pl - module subroutine helio_step_tp(self, cb, pl, param, t, dt) + module subroutine helio_step_tp(self, system, param, t, dt) !! author: David A. Minton !! !! Step active test particles ahead using Democratic Heliocentric method @@ -83,32 +87,31 @@ module subroutine helio_step_tp(self, cb, pl, param, t, dt) !! Adapted from Hal Levison's Swift routine helio_step_tp.f implicit none ! Arguments - class(helio_tp), intent(inout) :: self !! Helio test particle data structure - class(swiftest_cb), intent(inout) :: cb !! Swiftest central body particle data structure - class(whm_pl), intent(inout) :: pl !! WHM massive body data structure - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters of - real(DP), intent(in) :: t !! Current time - real(DP), intent(in) :: dt !! Stepsize + class(helio_tp), intent(inout) :: self !! Helio test particle data structure + 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 ! Internals logical, save :: lfirst = .true. !! Flag to indicate that this is the first call real(DP) :: dth !! Half step size - real(DP) :: mu !! Central mass term ! executable code - dth = 0.5_DP * dt - mu = cb%Gmass - if (lfirst) then - call self%vh2vb(vbcb = -self%ptbeg) - lfirst = .false. - end if - call self%lindrift(dth, self%ptbeg) - call self%getacch(cb, pl, param, t, self%xbeg) - call self%kickvb(dth) - call self%drift(cb, param, dt) - call self%getacch(cb, pl, param, t + dt, self%xend) - call self%kickvb(dth) - call self%lindrift(dth, self%ptend) - call self%vb2vh(vbcb = -self%ptend) + associate(cb => system%cb, pl => system%pl) + dth = 0.5_DP * dt + if (lfirst) then + call self%vh2vb(vbcb = -self%ptbeg) + lfirst = .false. + end if + call self%lindrift(dth, self%ptbeg) + call self%getacch(system, param, t, self%xbeg) + call self%kickvb(dth) + call self%drift(system, param, dt) + call self%getacch(system, param, t + dt, self%xend) + call self%kickvb(dth) + call self%lindrift(dth, self%ptend) + call self%vb2vh(vbcb = -self%ptend) + end associate return diff --git a/src/modules/helio_classes.f90 b/src/modules/helio_classes.f90 index 751e94fed..15c61343b 100644 --- a/src/modules/helio_classes.f90 +++ b/src/modules/helio_classes.f90 @@ -63,9 +63,7 @@ module helio_classes procedure, public :: step => helio_step_tp !! Steps the body forward one stepsize end type helio_tp - interface - module subroutine helio_coord_vb2vh_pl(self, cb) use swiftest_classes, only : swiftest_cb implicit none @@ -92,22 +90,24 @@ module subroutine helio_coord_vh2vb_tp(self, vbcb) real(DP), dimension(:), intent(in) :: vbcb !! Barycentric velocity of the central body end subroutine helio_coord_vh2vb_tp - module subroutine helio_drift_pl(self, cb, param, dt) - use swiftest_classes, only : swiftest_cb, swiftest_parameters + module subroutine helio_drift_pl(self, system, param, dt) + use swiftest_classes, only : swiftest_parameters + use whm_classes, only : whm_nbody_system implicit none - class(helio_pl), intent(inout) :: self !! Helio massive body object - class(swiftest_cb), intent(inout) :: cb !! Helio central body object - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters - real(DP), intent(in) :: dt !! Stepsize + class(helio_pl), intent(inout) :: self !! Helio massive body object + class(whm_nbody_system), intent(inout) :: system !! WHM nbody system object + class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters of + real(DP), intent(in) :: dt !! Stepsize end subroutine helio_drift_pl - module subroutine helio_drift_tp(self, cb, param, dt) - use swiftest_classes, only : swiftest_cb, swiftest_parameters + module subroutine helio_drift_tp(self, system, param, dt) + use swiftest_classes, only : swiftest_parameters + use whm_classes, only : whm_nbody_system implicit none - class(helio_tp), intent(inout) :: self !! Helio test particle object - class(swiftest_cb), intent(inout) :: cb !! Helio central body object - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters - real(DP), intent(in) :: dt !! Stepsize + class(helio_tp), intent(inout) :: self !! Helio test particle object + class(whm_nbody_system), intent(inout) :: system !! WHM nbody system object + class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters of + real(DP), intent(in) :: dt !! Stepsiz end subroutine helio_drift_tp module subroutine helio_drift_linear_pl(self, cb, dt, pt) @@ -126,25 +126,25 @@ module subroutine helio_drift_linear_tp(self, dt, pt) real(DP), dimension(:), intent(in) :: pt !! negative barycentric velocity of the Sun end subroutine helio_drift_linear_tp - module subroutine helio_getacch_pl(self, cb, param, t) - use swiftest_classes, only : swiftest_cb, swiftest_parameters + module subroutine helio_getacch_pl(self, system, param, t) + use swiftest_classes, only : swiftest_parameters + use whm_classes, only : whm_nbody_system implicit none - class(helio_pl), intent(inout) :: self !! Helio massive body particle data structure - class(swiftest_cb), intent(inout) :: cb !! Swiftest central body particle data structure - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters - real(DP), intent(in) :: t !! Current time + class(helio_pl), intent(inout) :: self !! Helio massive body particle data structure + class(whm_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 end subroutine helio_getacch_pl - module subroutine helio_getacch_tp(self, cb, pl, param, t, xh) - use swiftest_classes, only : swiftest_cb, swiftest_parameters - use whm_classes, only : whm_pl + module subroutine helio_getacch_tp(self, system, param, t, xh) + use swiftest_classes, only : swiftest_parameters + use whm_classes, only : whm_nbody_system implicit none - class(helio_tp), intent(inout) :: self !! Helio test particle data structure - class(swiftest_cb), intent(inout) :: cb !! Swiftest central body particle data structuree - class(whm_pl), intent(inout) :: pl !! Swiftest massive body particle data structure. - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters - real(DP), intent(in) :: t !! Current time - real(DP), dimension(:,:), intent(in) :: xh !! Heliocentric positions of planets + class(helio_tp), intent(inout) :: self !! Helio test particle data structure + class(whm_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) :: xh !! Heliocentric positions of planets end subroutine helio_getacch_tp module subroutine helio_getacch_int_pl(self, t) @@ -156,9 +156,9 @@ end subroutine helio_getacch_int_pl module subroutine helio_getacch_int_tp(self, pl, t, xh) use whm_classes, only : whm_pl implicit none - class(helio_tp), intent(inout) :: self !! Helio test particle data structure + class(helio_tp), intent(inout) :: self !! Helio test particle data structure class(whm_pl), intent(inout) :: pl !! WhM massive body particle data structure - real(DP), intent(in) :: t !! Current time + real(DP), intent(in) :: t !! Current time real(DP), dimension(:,:), intent(in) :: xh !! Heliocentric positions of planet end subroutine helio_getacch_int_tp @@ -174,7 +174,7 @@ 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(self, param) + module subroutine helio_step_system(self, param, t, dt) use swiftest_classes, only : swiftest_parameters implicit none class(helio_nbody_system), intent(inout) :: self !! Helio nbody system object @@ -186,7 +186,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 !! WHM massive body particle data structure + class(helio_pl), intent(inout) :: self !! Helio massive body particle data structure 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