diff --git a/src/discard/discard.f90 b/src/discard/discard.f90 index 1e69639be..daee4b15a 100644 --- a/src/discard/discard.f90 +++ b/src/discard/discard.f90 @@ -9,37 +9,35 @@ module subroutine discard_system(self, param) !! Adapted from David E. Kaufmann's Swifter routine: discard.f90 !! Adapted from Hal Levison's Swift routine discard.f implicit none - class(swiftest_nbody_system), intent(inout) :: self !! Swiftest system object - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters + class(swiftest_nbody_system), intent(inout) :: self !! Swiftest system object + class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters if (self%tp%nbody == 0) return - select type(self) + select type(system => self) class is (whm_nbody_system) - associate(cb => self%cb, pl => self%pl, tp => self%tp, t => param%t, dt => param%dt, & - msys => self%msys, discards => self%tp_discards, & - ntp => self%tp%nbody, ldiscard => self%tp%ldiscard) + associate(cb => system%cb, pl => system%pl, npl => system%pl%nbody, tp => system%tp, ntp => system%tp%nbody) if ((param%rmin >= 0.0_DP) .or. (param%rmax >= 0.0_DP) .or. & (param%rmaxu >= 0.0_DP) .or. ((param%qmin >= 0.0_DP) .and. (param%qmin_coord == "BARY"))) then - call pl%h2b(cb) + if (npl > 0) call pl%h2b(cb) if (ntp > 0) call tp%h2b(cb) end if if ((param%rmin >= 0.0_DP) .or. (param%rmax >= 0.0_DP) .or. (param%rmaxu >= 0.0_DP)) then - if (ntp > 0) call tp%discard_sun(cb, param, t, msys) + if (ntp > 0) call tp%discard_sun(system, param) end if - if (param%qmin >= 0.0_DP .and. ntp > 0) call tp%discard_peri(cb, pl, param, t, msys) - if (param%lclose .and. ntp > 0) call tp%discard_pl(pl, t, dt) + if (param%qmin >= 0.0_DP .and. ntp > 0) call tp%discard_peri(system, param) + if (param%lclose .and. ntp > 0) call tp%discard_pl(system, param) if (any(tp%ldiscard(1:ntp))) then ! Spill the discards to the spill list - call tp%spill(discards, ldiscard) - call self%write_discard(param, discards) + call tp%spill(system%tp_discards, tp%ldiscard) + call self%write_discard(param, system%tp_discards) end if end associate end select return end subroutine discard_system - module subroutine discard_sun_tp(self, cb, param, t, msys) + module subroutine discard_sun_tp(self, system, param) !! author: David A. Minton !! !! Check to see if test particles should be discarded based on their positions relative to the Sun @@ -49,16 +47,14 @@ module subroutine discard_sun_tp(self, cb, param, t, msys) !! Adapted from Hal Levison's Swift routine discard_sun.f implicit none ! Arguments - class(swiftest_tp), intent(inout) :: self !! Swiftest test particle object - class(swiftest_cb), intent(inout) :: cb !! Swiftest central body object - class(swiftest_parameters), intent(in) :: param !! parameters - real(DP), intent(in) :: t !! Current simulation tim - real(DP), intent(in) :: msys !! Total system mass + class(swiftest_tp), intent(inout) :: self !! Swiftest test particle object + class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object + class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters ! Internals integer(I4B) :: i real(DP) :: energy, vb2, rb2, rh2, rmin2, rmax2, rmaxu2 - associate(tp => self, ntp => self%nbody) + associate(tp => self, ntp => self%nbody, cb => system%cb, t => param%t, msys => system%msys) rmin2 = max(param%rmin * param%rmin, cb%radius * cb%radius) rmax2 = param%rmax**2 rmaxu2 = param%rmaxu**2 @@ -90,7 +86,7 @@ module subroutine discard_sun_tp(self, cb, param, t, msys) return end subroutine discard_sun_tp - module subroutine discard_peri_tp(self, cb, pl, param, t, msys) + module subroutine discard_peri_tp(self, system, param) !! author: David A. Minton !! !! Check to see if a test particle should be discarded because its perihelion distance becomes too small @@ -99,19 +95,16 @@ module subroutine discard_peri_tp(self, cb, pl, param, t, msys) !! Adapted from Hal Levison's Swift routine discard_peri.f implicit none ! Arguments - class(swiftest_tp), intent(inout) :: self !! Swiftest test particle object - class(swiftest_cb), intent(inout) :: cb !! Swiftest central body object - class(swiftest_pl), intent(inout) :: pl !! Swiftest massive body object - class(swiftest_parameters), intent(in) :: param !! parameters - real(DP), intent(in) :: t !! Current simulation tim - real(DP), intent(in) :: msys !! Total system mass + class(swiftest_tp), intent(inout) :: self !! Swiftest test particle object + class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object + class(swiftest_parameters), intent(in) :: param !! Current run configuration parameterss ! Internals logical, save :: lfirst = .true. integer(I4B) :: i, j, ih real(DP) :: r2 real(DP), dimension(NDIM) :: dx - associate(tp => self, ntp => self%nbody, npl => pl%nbody, qmin_coord => param%qmin_coord) + associate(cb => system%cb, tp => self, ntp => self%nbody, pl => system%pl, npl => system%pl%nbody, qmin_coord => param%qmin_coord, t => param%t, msys => system%msys) if (lfirst) then call util_hills(npl, pl) call util_peri(lfirst, ntp, tp, cb%Gmass, msys, param%qmin_coord) @@ -145,7 +138,7 @@ module subroutine discard_peri_tp(self, cb, pl, param, t, msys) end subroutine discard_peri_tp - module subroutine discard_pl_tp(self, pl, t, dt) + module subroutine discard_pl_tp(self, system, param) !! author: David A. Minton !! !! Check to see if test particles should be discarded based on their positions relative to the massive bodies @@ -154,16 +147,15 @@ module subroutine discard_pl_tp(self, pl, t, dt) !! Adapted from Hal Levison's Swift routine discard_pl.f implicit none ! Arguments - class(swiftest_tp), intent(inout) :: self !! Swiftest test particle object - class(swiftest_pl), intent(inout) :: pl !! Swiftest massive body object - real(DP), intent(in) :: t !! Current simulation tim - real(DP), intent(in) :: dt !! Stepsize + class(swiftest_tp), intent(inout) :: self !! Swiftest test particle object + class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object + class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters ! Internals integer(I4B) :: i, j, isp real(DP) :: r2min, radius real(DP), dimension(NDIM) :: dx, dv - associate(tp => self, ntp => self%nbody, npl => pl%nbody) + associate(tp => self, ntp => self%nbody, pl => system%pl, npl => system%pl%nbody, t => param%t, dt => param%dt) do i = 1, ntp if (tp%status(i) == ACTIVE) then do j = 1, npl @@ -187,7 +179,7 @@ module subroutine discard_pl_tp(self, pl, t, dt) end subroutine discard_pl_tp - module subroutine discard_pl_close(dx, dv, dt, r2crit, iflag, r2min) + subroutine discard_pl_close(dx, dv, dt, r2crit, iflag, r2min) !! author: David A. Minton !! !! Check to see if a test particle and massive body are having, or will have within the next time step, an encounter such diff --git a/src/gr/gr.f90 b/src/gr/gr.f90 index d50968ad1..6afc9f5ed 100644 --- a/src/gr/gr.f90 +++ b/src/gr/gr.f90 @@ -1,7 +1,8 @@ submodule(swiftest_classes) s_gr use swiftest contains - module procedure gr_getaccb_ns_body + subroutine gr_getaccb_ns_body(self, cb, param, agr, agr0) + !! author: David A. Minton !! !! Add relativistic correction acceleration for non-symplectic integrators @@ -9,7 +10,13 @@ !! !! 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_parameters), intent(in) :: param + 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 integer(I4B) :: i @@ -29,7 +36,6 @@ agr0 = 0.0_DP select type(self) class is (swiftest_pl) - !do concurrent(i = 1:NDIM) do i = 1, NDIM agr0(i) = -sum(self%Gmass(1:n) * agr(1:n, i) / msun) end do @@ -39,6 +45,6 @@ return - end procedure gr_getaccb_ns_body + end subroutine gr_getaccb_ns_body end submodule s_gr \ No newline at end of file diff --git a/src/helio/helio_drift.f90 b/src/helio/helio_drift.f90 index 1244533a0..59f828425 100644 --- a/src/helio/helio_drift.f90 +++ b/src/helio/helio_drift.f90 @@ -1,7 +1,8 @@ submodule (helio_classes) s_helio_drift use swiftest contains - module subroutine helio_drift_pl(self, cb, param, dt) + module subroutine helio_drift_pl(self, system, param, dt) + !! author: David A. Minton !! !! Loop through massive bodies and call Danby drift routine @@ -11,22 +12,17 @@ module subroutine helio_drift_pl(self, cb, param, dt) !! Adapted from Hal Levison's Swift routine drift.f implicit none ! Arguments - class(helio_pl), intent(inout) :: self !! Helio test particle data structure - class(swiftest_cb), intent(inout) :: cb !! Helio central body particle data structuree - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters of - real(DP), intent(in) :: dt !! Stepsize + 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) :: dt !! Stepsize) ! Internals integer(I4B) :: i !! Loop counter real(DP) :: rmag, vmag2, energy integer(I4B), dimension(:),allocatable :: iflag !! Vectorized error code flag real(DP), dimension(:), allocatable :: dtp - associate(npl => self%nbody, & - xh => self%xh, & - vb => self%vb, & - status => self%status,& - mu => self%mu) - + associate(pl => self, npl => self%nbody) if (npl == 0) return allocate(iflag(npl)) @@ -35,23 +31,23 @@ module subroutine helio_drift_pl(self, cb, param, dt) if (param%lgr) then do i = 1,npl - rmag = norm2(xh(:, i)) - vmag2 = dot_product(vb(:, i), vb(:, i)) - energy = 0.5_DP * vmag2 - mu(i) / rmag + rmag = norm2(pl%xh(:, i)) + vmag2 = dot_product(pl%vb(:, i), pl%vb(:, i)) + energy = 0.5_DP * vmag2 - pl%mu(i) / rmag dtp(i) = dt * (1.0_DP + 3 * param%inv_c2 * energy) end do else dtp(:) = dt end if - call drift_one(mu(1:npl), xh(1,1:npl), xh(2,1:npl), xh(3,1:npl), & - vb(1,1:npl), vb(2,1:npl), vb(3,1:npl), & - dtp(1:npl), iflag(1:npl)) + call drift_one(pl%mu(1:npl), pl%xh(1,1:npl), pl%xh(2,1:npl), pl%xh(3,1:npl), & + pl%vb(1,1:npl), pl%vb(2,1:npl), pl%vb(3,1:npl), & + dtp(1:npl), iflag(1:npl)) if (any(iflag(1:npl) /= 0)) then do i = 1, npl - write(*, *) " Planet ", self%name(i), " is lost!!!!!!!!!!" - write(*, *) xh(:,i) - write(*, *) vb(:,i) + write(*, *) " Planet ", pl%name(i), " is lost!!!!!!!!!!" + write(*, *) pl%xh(:,i) + write(*, *) pl%vb(:,i) write(*, *) " stopping " call util_exit(FAILURE) end do @@ -94,7 +90,7 @@ module subroutine helio_drift_linear_pl(self, cb, dt, pt) end subroutine helio_drift_linear_pl - module subroutine helio_drift_tp(self, cb, param, dt) + module subroutine helio_drift_tp(self, system, param, dt) !! author: David A. Minton !! !! Loop through test particles and call Danby drift routine @@ -103,21 +99,17 @@ module subroutine helio_drift_tp(self, cb, param, dt) !! Adapted from Hal Levison's Swift routine drift_tp.f implicit none ! Arguments - class(helio_tp), intent(inout) :: self !! Helio test particle data structure - class(swiftest_cb), intent(inout) :: cb !! Helio central body particle data structuree - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters of - real(DP), intent(in) :: dt !! Stepsize + 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 of + real(DP), intent(in) :: dt !! Stepsiz ! Internals integer(I4B) :: i !! Loop counter real(DP) :: rmag, vmag2, energy real(DP), dimension(:), allocatable :: dtp integer(I4B), dimension(:),allocatable :: iflag !! Vectorized error code flag - associate(ntp => self%nbody, & - xh => self%xh, & - vh => self%vh, & - status => self%status,& - mu => self%mu) + associate(tp => self, ntp => self%nbody) if (ntp == 0) return allocate(iflag(ntp)) allocate(dtp(ntp)) @@ -128,21 +120,21 @@ module subroutine helio_drift_tp(self, cb, param, dt) if (param%lgr) then do i = 1,ntp - rmag = norm2(xh(:, i)) - vmag2 = dot_product(vh(:, i), vh(:, i)) - energy = 0.5_DP * vmag2 - mu(i) / rmag + rmag = norm2(tp%xh(:, i)) + vmag2 = dot_product(tp%vh(:, i), tp%vh(:, i)) + energy = 0.5_DP * vmag2 - tp%mu(i) / rmag dtp(i) = dt * (1.0_DP + 3 * param%inv_c2 * energy) end do else dtp(:) = dt end if - call drift_one(mu(1:ntp), xh(1,1:ntp), xh(2,1:ntp), xh(3,1:ntp), & - vh(1,1:ntp), vh(2,1:ntp), vh(3,1:ntp), & - dtp(1:ntp), iflag(1:ntp)) + call drift_one(tp%mu(1:ntp), tp%xh(1,1:ntp), tp%xh(2,1:ntp), tp%xh(3,1:ntp), & + tp%vh(1,1:ntp), tp%vh(2,1:ntp), tp%vh(3,1:ntp), & + dtp(1:ntp), iflag(1:ntp)) if (any(iflag(1:ntp) /= 0)) then - status = DISCARDED_DRIFTERR + tp%status = DISCARDED_DRIFTERR do i = 1, ntp - if (iflag(i) /= 0) write(*, *) "Particle ", self%name(i), " lost due to error in Danby drift" + if (iflag(i) /= 0) write(*, *) "Particle ", tp%name(i), " lost due to error in Danby drift" end do end if end associate diff --git a/src/helio/helio_getacch.f90 b/src/helio/helio_getacch.f90 index 3dbcc4e45..806047cc3 100644 --- a/src/helio/helio_getacch.f90 +++ b/src/helio/helio_getacch.f90 @@ -1,7 +1,7 @@ submodule (helio_classes) s_helio_getacch use swiftest contains - module subroutine helio_getacch_pl(self, cb, param, t) + module subroutine helio_getacch_pl(self, system, param, t) !! author: David A. Minton !! !! Compute heliocentric accelerations of massive bodies @@ -10,10 +10,10 @@ module subroutine helio_getacch_pl(self, cb, param, t) !! Adapted from Hal Levison's Swift routine helio_getacch.f implicit none ! Arguments - 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 of - real(DP), intent(in) :: t !! Current time + class(helio_pl), intent(inout) :: self !! Helio massive body particle data structure + 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 ! Internals logical, save :: lmalloc = .true. integer(I4B) :: i @@ -21,22 +21,22 @@ module subroutine helio_getacch_pl(self, cb, param, t) real(DP), dimension(:), allocatable, save :: irh real(DP), dimension(:, :), allocatable, save :: xh_loc, aobl - associate(npl => self%nbody) + associate(pl => self, npl => self%nbody) !if (lflag) then - self%ahi(:,2:npl) = 0.0_DP - call helio_getacch_int_pl(self, t) + pl%ahi(:,2:npl) = 0.0_DP + call helio_getacch_int_pl(pl, t) !end if !if (param%loblatecb) call self%obl_acc(cb) TODO: Fix this !else - self%ah(:,:) = self%ahi(:,:) + pl%ah(:,:) = pl%ahi(:,:) !end if - if (param%lextra_force) call self%user_getacch(cb, param, t) + if (param%lextra_force) call pl%user_getacch(system, param, t) end associate return end subroutine helio_getacch_pl - module subroutine helio_getacch_tp(self, cb, pl, param, t, xh) + module subroutine helio_getacch_tp(self, system, param, t, xhp) !! author: David A. Minton !! !! Compute heliocentric accelerations of test particles @@ -45,34 +45,35 @@ module subroutine helio_getacch_tp(self, cb, pl, param, t, xh) !! Adapted from Hal Levison's Swift routine helio_getacch_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 structuree - class(whm_pl), intent(inout) :: pl !! WHM massive body particle data structure. - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters of - real(DP), intent(in) :: t !! Current time - real(DP), dimension(:,:), intent(in) :: xh !! Heliocentric positions of planets + class(helio_tp), intent(inout) :: self !! WHM test particle data structure + 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 time + real(DP), dimension(:,:), intent(in) :: xhp !! Heliocentric positions of planets at the current substep ! Internals logical, save :: lmalloc = .true. integer(I4B) :: i real(DP) :: r2, mu real(DP), dimension(:), allocatable, save :: irh, irht - real(DP), dimension(:, :), allocatable, save :: aobl, xht, aoblt ! executable code - associate(ntp => self%nbody, npl => pl%nbody) - !if (lflag) then - self%ahi(:,:) = 0.0_DP - call helio_getacch_int_tp(self, pl, t, xh) - !end if - !if (param%loblatecb) call self%obl_acc(cb) TODO: Fix this - self%ah(:,:) = self%ahi(:,:) - if (param%lextra_force) call self%user_getacch(cb, param, t) - if (param%lgr) call self%gr_getacch(cb, param) + associate(tp => self, ntp => self%nbody, npl => system%pl%nbody) + select type(pl => system%pl) + class is (rmvs_pl) + !if (lflag) then + self%ahi(:,:) = 0.0_DP + call helio_getacch_int_tp(tp, pl, t, xhp) + !end if + !if (param%loblatecb) call self%obl_acc(cb) TODO: Fix this + tp%ah(:,:) = tp%ahi(:,:) + if (param%lextra_force) call tp%user_getacch(system, param, t) + if (param%lgr) call tp%gr_getacch(system, param) + end select end associate return end subroutine helio_getacch_tp - module subroutine helio_getacch_int_pl(self, t) + subroutine helio_getacch_int_pl(pl, t) !! author: David A. Minton !! !! Compute direct cross term heliocentric accelerations of massive bodiese @@ -81,17 +82,17 @@ module subroutine helio_getacch_int_pl(self, t) !! Adapted from Hal Levison's Swift routine getacch_ah3.f implicit none ! Arguments - class(helio_pl), intent(inout) :: self !! Helio massive body particle data structure - real(DP), intent(in) :: t !! Current time + class(helio_pl), intent(inout) :: pl !! Helio massive body particle data structure + real(DP), intent(in) :: t !! Current time ! Internals integer(I4B) :: i, j real(DP) :: rji2, irij3, faci, facj real(DP), dimension(NDIM) :: dx - associate(npl => self%nbody) + associate(npl => pl%nbody) do i = 2, npl - 1 do j = i + 1, npl - dx(:) = self%xh(:,j) - self%xh(:,i) + dx(:) = pl%xh(:,j) - pl%xh(:,i) rji2 = dot_product(dx(:), dx(:)) irij3 = 1.0_DP / (rji2 * sqrt(rji2)) faci = self%Gmass(i) * irij3 @@ -105,7 +106,7 @@ module subroutine helio_getacch_int_pl(self, t) return end subroutine helio_getacch_int_pl - module subroutine helio_getacch_int_tp(self, pl, t, xh) + subroutine helio_getacch_int_tp(tp, pl, t, xhp) !! author: David A. Minton !! !! Compute direct cross term heliocentric accelerations of test particles @@ -114,23 +115,23 @@ module subroutine helio_getacch_int_tp(self, pl, t, xh) !! Adapted from Hal Levison's Swift routine getacch_ah3_tp.f implicit none ! Arguments - class(helio_tp), intent(inout) :: self !! Helio test particle data structure - class(whm_pl), intent(inout) :: pl !! Helio massive body particle data structure + class(helio_tp), intent(inout) :: tp !! Helio test particle data structure + class(helio_pl), intent(inout) :: pl !! Helio massive body particle data structure real(DP), intent(in) :: t !! Current time - real(DP), dimension(:,:), intent(in) :: xh !! Heliocentric positions of planets + real(DP), dimension(:,:), intent(in) :: xhp !! Heliocentric positions of planets ! Internals integer(I4B) :: i, j real(DP) :: r2, fac real(DP), dimension(NDIM) :: dx - associate(ntp => self%nbody, npl => pl%nbody) + associate(ntp => tp%nbody, npl => pl%nbody) do i = 1, ntp - if (self%status(i) == ACTIVE) then + if (tp%status(i) == ACTIVE) then do j = 2, npl - dx(:) = self%xh(:,i) - xh(:,j) + dx(:) = tp%xh(:,i) - xhp(:,j) r2 = dot_product(dx(:), dx(:)) fac = pl%Gmass(j) / (r2 * sqrt(r2)) - self%ahi(:,i) = self%ahi(:,i) - fac * dx(:) + tp%ahi(:,i) = tp%ahi(:,i) - fac * dx(:) end do end if end do diff --git a/src/io/io.f90 b/src/io/io.f90 index 6f930cdf2..629a4952f 100644 --- a/src/io/io.f90 +++ b/src/io/io.f90 @@ -392,11 +392,13 @@ module subroutine io_param_writer(self, unit, iotype, v_list, iostat, iomsg) write(unit, Lfmt) "ENERGY", self%lenergy !write(unit, Lfmt) "YARKOVSKY", self%lyarkovsky !write(unit, Lfmt) "YORP", self%lyorp + iostat = 0 + iomsg = "UDIO not implemented" return end subroutine io_param_writer - module subroutine io_dump_param(self, param_file_name, t, dt) + module subroutine io_dump_param(self, param_file_name) !! author: David A. Minton !! !! Dump integration parameters to file @@ -406,13 +408,11 @@ module subroutine io_dump_param(self, param_file_name, t, dt) implicit none ! Arguments class(swiftest_parameters),intent(in) :: self !! Output collection of parameters - character(len=*), intent(in) :: param_file_name !! Parameter input file name (i.e. param.in) - real(DP),intent(in) :: t !! Current simulation time - real(DP),intent(in) :: dt !! Step size + character(len=*), intent(in) :: param_file_name !! Parameter input file name (i.e. param.in) ! Internals - integer(I4B), parameter :: LUN = 7 !! Unit number of output file - integer(I4B) :: ierr !! Error code - character(STRMAX) :: error_message !! Error message in UDIO procedure + integer(I4B), parameter :: LUN = 7 !! Unit number of output file + integer(I4B) :: ierr !! Error code + character(STRMAX) :: error_message !! Error message in UDIO procedure open(unit = LUN, file = param_file_name, status='replace', form = 'FORMATTED', iostat =ierr) if (ierr /=0) then @@ -434,7 +434,7 @@ module subroutine io_dump_param(self, param_file_name, t, dt) return end subroutine io_dump_param - module subroutine io_dump_swiftest(self, param, t, dt, msg) + module subroutine io_dump_swiftest(self, param, msg) !! author: David A. Minton !! !! Dump massive body data to files @@ -443,11 +443,9 @@ module subroutine io_dump_swiftest(self, param, t, dt, msg) !! Adapted from Hal Levison's Swift routine io_dump_pl.f and io_dump_tp.f implicit none ! Arguments - class(swiftest_base), intent(inout) :: self !! Swiftest base object + class(swiftest_base), intent(inout) :: self !! Swiftest base object class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters - real(DP), intent(in) :: t !! Current simulation time - real(DP), intent(in) :: dt !! Stepsize - character(*), optional, intent(in) :: msg !! Message to display with dump operation + character(*), optional, intent(in) :: msg !! Message to display with dump operation ! Internals integer(I4B) :: ierr !! Error code integer(I4B),parameter :: LUN = 7 !! Unit number for dump file @@ -468,13 +466,13 @@ module subroutine io_dump_swiftest(self, param, t, dt, msg) write(*, *) " Unable to open binary dump file " // dump_file_name call util_exit(FAILURE) end if - call self%write_frame(iu, param, t, dt) + call self%write_frame(iu, param) close(LUN) return end subroutine io_dump_swiftest - module subroutine io_dump_system(self, param, t, dt, msg) + module subroutine io_dump_system(self, param, msg) !! author: David A. Minton !! !! Dumps the state of the system to files in case the simulation is interrupted. @@ -482,11 +480,9 @@ module subroutine io_dump_system(self, param, t, dt, msg) !! so that if a dump file gets corrupted during writing, the user can restart from the older one. implicit none ! Arguments - class(swiftest_nbody_system), intent(inout) :: self !! Swiftest system object - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters - real(DP), intent(in) :: t !! Current simulation time - real(DP), intent(in) :: dt !! Stepsize - character(*), optional, intent(in) :: msg !! Message to display with dump operation + class(swiftest_nbody_system), intent(inout) :: self !! Swiftest system object + class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters + character(*), optional, intent(in) :: msg !! Message to display with dump operation ! Internals class(swiftest_parameters), allocatable :: dump_param !! Local parameters variable used to parameters change input file names !! to dump file-specific values without changing the user-defined values @@ -495,26 +491,25 @@ module subroutine io_dump_system(self, param, t, dt, msg) character(len=:), allocatable :: param_file_name real(DP) :: tfrac - allocate(dump_param, source=param) - param_file_name = trim(adjustl(DUMP_PARAM_FILE(idx))) + param_file_name = trim(adjustl(DUMP_PARAM_FILE(idx))) dump_param%incbfile = trim(adjustl(DUMP_CB_FILE(idx))) dump_param%inplfile = trim(adjustl(DUMP_PL_FILE(idx))) dump_param%intpfile = trim(adjustl(DUMP_TP_FILE(idx))) dump_param%out_form = XV dump_param%out_stat = 'APPEND' - call dump_param%dump(param_file_name,t,dt) + call dump_param%dump(param_file_name) - call self%cb%dump(dump_param, t, dt) - if (self%pl%nbody > 0) call self%pl%dump(dump_param, t, dt) - if (self%tp%nbody > 0) call self%tp%dump(dump_param, t, dt) + call self%cb%dump(dump_param) + if (self%pl%nbody > 0) call self%pl%dump(dump_param) + if (self%tp%nbody > 0) call self%tp%dump(dump_param) idx = idx + 1 if (idx > NDUMPFILES) idx = 1 ! Print the status message (format code passed in from main driver) - tfrac = (t - param%t0) / (param%tstop - param%t0) - write(*,msg) t, tfrac, self%pl%nbody, self%tp%nbody + tfrac = (param%t - param%t0) / (param%tstop - param%t0) + write(*,msg) param%t, tfrac, self%pl%nbody, self%tp%nbody return end subroutine io_dump_system @@ -588,12 +583,12 @@ module function io_get_token(buffer, ifirst, ilast, ierr) result(token) !! Adapted from David E. Kaufmann's Swifter routine io_get_token.f90 implicit none ! Arguments - character(len=*), intent(in) :: buffer !! Input string buffer - integer(I4B), intent(inout) :: ifirst !! Index of the buffer at which to start the search for a token - integer(I4B), intent(out) :: ilast !! Index of the buffer at the end of the returned token - integer(I4B), intent(out) :: ierr !! Error code + character(len=*), intent(in) :: buffer !! Input string buffer + integer(I4B), intent(inout) :: ifirst !! Index of the buffer at which to start the search for a token + integer(I4B), intent(out) :: ilast !! Index of the buffer at the end of the returned token + integer(I4B), intent(out) :: ierr !! Error code ! Result - character(len=:),allocatable :: token !! Returned token string + character(len=:), allocatable :: token !! Returned token string ! Internals integer(I4B) :: i,ilength @@ -635,7 +630,7 @@ module subroutine io_read_body_in(self, param) !! Adapted from Martin Duncan's Swift routine swiftest_init_pl.f and swiftest_init_tp.f implicit none ! Arguments - class(swiftest_body), intent(inout) :: self !! Swiftest particle object + class(swiftest_body), intent(inout) :: self !! Swiftest particle object class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters ! Internals integer(I4B), parameter :: LUN = 7 !! Unit number of input file @@ -699,7 +694,7 @@ module subroutine io_read_body_in(self, param) read(iu, iostat = ierr) nbody call self%setup(nbody) if (nbody > 0) then - call self%read_frame(iu, param, XV, t, ierr) + call self%read_frame(iu, param, XV, ierr) self%status(:) = ACTIVE end if case default @@ -724,7 +719,7 @@ module subroutine io_read_cb_in(self, param) !! Adapted from Martin Duncan's Swift routine swiftest_init_pl.f implicit none ! Arguments - class(swiftest_cb), intent(inout) :: self + class(swiftest_cb), intent(inout) :: self class(swiftest_parameters), intent(inout) :: param ! Internals integer(I4B), parameter :: LUN = 7 !! Unit number of input file @@ -755,7 +750,7 @@ module subroutine io_read_cb_in(self, param) else open(unit = iu, file = param%incbfile, status = 'old', form = 'UNFORMATTED', iostat = ierr) - call self%read_frame(iu, param, XV, t, ierr) + call self%read_frame(iu, param, XV, ierr) end if close(iu) if (ierr /= 0) then @@ -808,7 +803,7 @@ module subroutine io_read_param_in(self, param_file_name) return end subroutine io_read_param_in - module function io_read_encounter(t, name1, name2, mass1, mass2, radius1, radius2, & + function io_read_encounter(t, name1, name2, mass1, mass2, radius1, radius2, & xh1, xh2, vh1, vh2, encounter_file, out_type) result(ierr) !! author: David A. Minton !! @@ -818,10 +813,10 @@ module function io_read_encounter(t, name1, name2, mass1, mass2, radius1, radius !! Adapted from David E. Kaufmann's Swifter routine: io_read_encounter.f90 implicit none ! Arguments - integer(I4B), intent(out) :: name1, name2 - real(DP), intent(out) :: t, mass1, mass2, radius1, radius2 - real(DP), dimension(NDIM), intent(out) :: xh1, xh2, vh1, vh2 - character(*), intent(in) :: encounter_file, out_type + integer(I4B), intent(out) :: name1, name2 + real(DP), intent(out) :: t, mass1, mass2, radius1, radius2 + real(DP), dimension(:), intent(out) :: xh1, xh2, vh1, vh2 + character(*), intent(in) :: encounter_file, out_type ! Result integer(I4B) :: ierr ! Internals @@ -858,7 +853,7 @@ module function io_read_encounter(t, name1, name2, mass1, mass2, radius1, radius return end function io_read_encounter - module subroutine io_read_frame_body(self, iu, param, form, t, ierr) + module subroutine io_read_frame_body(self, iu, param, form, ierr) !! author: David A. Minton !! !! Reads a frame of output of either test particle or massive body data to the binary output file @@ -868,12 +863,11 @@ module subroutine io_read_frame_body(self, iu, param, form, t, ierr) !! Adapted from Hal Levison's Swift routine io_read_frame.F implicit none ! Arguments - class(swiftest_body), intent(inout) :: self !! Swiftest particle object - integer(I4B), intent(inout) :: iu !! Unit number for the output file to write frame to + class(swiftest_body), intent(inout) :: self !! Swiftest particle object + integer(I4B), intent(inout) :: iu !! Unit number for the output file to write frame to class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters - character(*), intent(in) :: form !! Input format code ("XV" or "EL") - real(DP), intent(out) :: t !! Simulation time - integer(I4B), intent(out) :: ierr !! Error code + character(*), intent(in) :: form !! Input format code ("XV" or "EL") + integer(I4B), intent(out) :: ierr !! Error code associate(n => self%nbody) read(iu, iostat = ierr) self%name(1:n) @@ -921,7 +915,7 @@ module subroutine io_read_frame_body(self, iu, param, form, t, ierr) return end subroutine io_read_frame_body - module subroutine io_read_frame_cb(self, iu, param, form, t, ierr) + module subroutine io_read_frame_cb(self, iu, param, form, ierr) !! author: David A. Minton !! !! Reads a frame of output of central body data to the binary output file @@ -930,12 +924,11 @@ module subroutine io_read_frame_cb(self, iu, param, form, t, ierr) !! Adapted from Hal Levison's Swift routine io_read_frame.F implicit none ! Arguments - class(swiftest_cb), intent(inout) :: self !! Swiftest central body object - integer(I4B), intent(inout) :: iu !! Unit number for the output file to write frame to + class(swiftest_cb), intent(inout) :: self !! Swiftest central body object + integer(I4B), intent(inout) :: iu !! Unit number for the output file to write frame to class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters - character(*), intent(in) :: form !! Input format code ("XV" or "EL") - real(DP), intent(out) :: t !! Simulation time - integer(I4B), intent(out) :: ierr !! Error cod + character(*), intent(in) :: form !! Input format code ("XV" or "EL") + integer(I4B), intent(out) :: ierr !! Error cod read(iu, iostat = ierr) self%Gmass self%mass = self%Gmass / param%GU @@ -958,18 +951,17 @@ module subroutine io_read_frame_cb(self, iu, param, form, t, ierr) return end subroutine io_read_frame_cb - module subroutine io_read_frame_system(self, iu, param, form, t, ierr) + module subroutine io_read_frame_system(self, iu, param, form, ierr) !! author: The Purdue Swiftest Team - David A. Minton, Carlisle A. Wishard, Jennifer L.L. Pouplin, and Jacob R. Elliott !! !! Read a frame (header plus records for each massive body and active test particle) from a output binary file implicit none ! Arguments - class(swiftest_nbody_system), intent(inout) :: self !! Swiftest system object - integer(I4B), intent(inout) :: iu !! Unit number for the output file to write frame to - class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters - character(*), intent(in) :: form !! Input format code ("XV" or "EL") - real(DP), intent(out) :: t !! Current simulation time - integer(I4B), intent(out) :: ierr !! Error code + class(swiftest_nbody_system), intent(inout) :: self !! Swiftest system object + integer(I4B), intent(inout) :: iu !! Unit number for the output file to write frame to + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters + character(*), intent(in) :: form !! Input format code ("XV" or "EL") + integer(I4B), intent(out) :: ierr !! Error code ! Internals logical, save :: lfirst = .true. @@ -983,21 +975,21 @@ module subroutine io_read_frame_system(self, iu, param, form, t, ierr) return end if end if - ierr = io_read_hdr(iu, t, self%pl%nbody, self%tp%nbody, param%out_form, param%out_type) + ierr = io_read_hdr(iu, param%t, self%pl%nbody, self%tp%nbody, param%out_form, param%out_type) if (ierr /= 0) then write(*, *) "Swiftest error:" write(*, *) " Binary output file already exists or cannot be accessed" return end if - call self%cb%read_frame(iu, param, form, t, ierr) + call self%cb%read_frame(iu, param, form, ierr) if (ierr /= 0) return - call self%pl%read_frame(iu, param, form, t, ierr) + call self%pl%read_frame(iu, param, form, ierr) if (ierr /= 0) return - call self%tp%read_frame(iu, param, form, t, ierr) + call self%tp%read_frame(iu, param, form, ierr) return end subroutine io_read_frame_system - module function io_read_hdr(iu, t, npl, ntp, out_form, out_type) result(ierr) + function io_read_hdr(iu, t, npl, ntp, out_form, out_type) result(ierr) !! author: David A. Minton !! !! Read frame header from input binary files @@ -1009,7 +1001,7 @@ module function io_read_hdr(iu, t, npl, ntp, out_form, out_type) result(ierr) integer(I4B), intent(in) :: iu integer(I4B), intent(out) :: npl, ntp character(*), intent(out) :: out_form - real(DP), intent(out) :: t + real(DP), intent(out) :: t character(*), intent(in) :: out_type ! Result integer(I4B) :: ierr @@ -1041,8 +1033,8 @@ module subroutine io_read_initialize_system(self, param) !! implicit none ! Arguments - class(swiftest_nbody_system), intent(inout) :: self !! Swiftest system object - class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters + class(swiftest_nbody_system), intent(inout) :: self !! Swiftest system object + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters call self%cb%initialize(param) call self%pl%initialize(param) @@ -1062,9 +1054,9 @@ module subroutine io_write_discard(self, param, discards) !! Adapted from Hal Levison's Swift routine io_discard_write.f implicit none ! Arguments - class(swiftest_nbody_system), intent(inout) :: self !! Swiftest system object - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters - class(swiftest_body), intent(inout) :: discards !! Swiftest discard object + class(swiftest_nbody_system), intent(inout) :: self !! Swiftest system object + class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters + class(swiftest_body), intent(inout) :: discards !! Swiftest discard object ! Internals integer(I4B), parameter :: LUN = 40 integer(I4B) :: i, ierr @@ -1134,8 +1126,8 @@ module subroutine io_write_discard(self, param, discards) end subroutine io_write_discard - module subroutine io_write_encounter(t, name1, name2, mass1, mass2, radius1, radius2, & - xh1, xh2, vh1, vh2, encounter_file, out_type) + subroutine io_write_encounter(t, name1, name2, mass1, mass2, radius1, radius2, & + xh1, xh2, vh1, vh2, encounter_file, out_type) !! author: David A. Minton !! !! Write close encounter data to output binary files @@ -1145,10 +1137,10 @@ module subroutine io_write_encounter(t, name1, name2, mass1, mass2, radius1, rad !! Adapted from Hal Levison's Swift routine io_write_encounter.f implicit none ! Arguments - integer(I4B), intent(in) :: name1, name2 - real(DP), intent(in) :: t, mass1, mass2, radius1, radius2 - real(DP), dimension(NDIM), intent(in) :: xh1, xh2, vh1, vh2 - character(*), intent(in) :: encounter_file, out_type + integer(I4B), intent(in) :: name1, name2 + real(DP), intent(in) :: t, mass1, mass2, radius1, radius2 + real(DP), dimension(:), intent(in) :: xh1, xh2, vh1, vh2 + character(*), intent(in) :: encounter_file, out_type ! Internals logical , save :: lfirst = .true. integer(I4B), parameter :: lun = 30 @@ -1184,7 +1176,7 @@ module subroutine io_write_encounter(t, name1, name2, mass1, mass2, radius1, rad end subroutine io_write_encounter - module subroutine io_write_frame_body(self, iu, param, t, dt) + module subroutine io_write_frame_body(self, iu, param) !! author: David A. Minton !! !! Write a frame of output of either test particle or massive body data to the binary output file @@ -1194,11 +1186,9 @@ module subroutine io_write_frame_body(self, iu, param, t, dt) !! Adapted from Hal Levison's Swift routine io_write_frame.F implicit none ! Arguments - class(swiftest_body), intent(in) :: self !! Swiftest particle object - integer(I4B), intent(inout) :: iu !! Unit number for the output file to write frame to + class(swiftest_body), intent(in) :: self !! Swiftest particle object + integer(I4B), intent(inout) :: iu !! Unit number for the output file to write frame to class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters - real(DP), intent(in) :: t !! Current simulation time - real(DP), intent(in) :: dt !! Step size associate(n => self%nbody) if (n == 0) return @@ -1241,7 +1231,7 @@ module subroutine io_write_frame_body(self, iu, param, t, dt) return end subroutine io_write_frame_body - module subroutine io_write_frame_cb(self, iu, param, t, dt) + module subroutine io_write_frame_cb(self, iu, param) !! author: David A. Minton !! !! Write a frame of output of central body data to the binary output file @@ -1250,11 +1240,9 @@ module subroutine io_write_frame_cb(self, iu, param, t, dt) !! Adapted from Hal Levison's Swift routine io_write_frame.F implicit none ! Arguments - class(swiftest_cb), intent(in) :: self !! Swiftest central body object - integer(I4B), intent(inout) :: iu !! Unit number for the output file to write frame to + class(swiftest_cb), intent(in) :: self !! Swiftest central body object + integer(I4B), intent(inout) :: iu !! Unit number for the output file to write frame to class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters - real(DP), intent(in) :: t !! Current simulation time - real(DP), intent(in) :: dt !! Step size write(iu) self%Gmass write(iu) self%radius @@ -1276,7 +1264,7 @@ module subroutine io_write_frame_cb(self, iu, param, t, dt) return end subroutine io_write_frame_cb - module subroutine io_write_frame_system(self, iu, param, t, dt) + module subroutine io_write_frame_system(self, iu, param) !! author: The Purdue Swiftest Team - David A. Minton, Carlisle A. Wishard, Jennifer L.L. Pouplin, and Jacob R. Elliott !! !! Write a frame (header plus records for each massive body and active test particle) to output binary file @@ -1286,18 +1274,16 @@ module subroutine io_write_frame_system(self, iu, param, t, dt) !! Adapted from Hal Levison's Swift routine io_write_frame.F implicit none ! Arguments - class(swiftest_nbody_system), intent(in) :: self !! Swiftest system object - integer(I4B), intent(inout) :: iu !! Unit number for the output file to write frame to - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters - real(DP), intent(in) :: t !! Current simulation time - real(DP), intent(in) :: dt !! Step size + class(swiftest_nbody_system), intent(in) :: self !! Swiftest system object + integer(I4B), intent(inout) :: iu !! Unit number for the output file to write frame to + class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters ! Internals - logical, save :: lfirst = .true. !! Flag to determine if this is the first call of this method - integer(I4B) :: ierr !! I/O error code + logical, save :: lfirst = .true. !! Flag to determine if this is the first call of this method + integer(I4B) :: ierr !! I/O error code - class(swiftest_cb), allocatable :: cb !! Temporary local version of pl structure used for non-destructive conversions - class(swiftest_pl), allocatable :: pl !! Temporary local version of pl structure used for non-destructive conversions - class(swiftest_tp), allocatable :: tp !! Temporary local version of pl structure used for non-destructive conversions + class(swiftest_cb), allocatable :: cb !! Temporary local version of pl structure used for non-destructive conversions + class(swiftest_pl), allocatable :: pl !! Temporary local version of pl structure used for non-destructive conversions + class(swiftest_tp), allocatable :: tp !! Temporary local version of pl structure used for non-destructive conversions allocate(cb, source = self%cb) allocate(pl, source = self%pl) @@ -1329,7 +1315,7 @@ module subroutine io_write_frame_system(self, iu, param, t, dt) call util_exit(FAILURE) end if end if - call io_write_hdr(iu, t, pl%nbody, tp%nbody, param%out_form, param%out_type) + call io_write_hdr(iu, param%t, pl%nbody, tp%nbody, param%out_form, param%out_type) if (param%lgr) then associate(vh => pl%vh, vht => tp%vh) @@ -1350,9 +1336,9 @@ module subroutine io_write_frame_system(self, iu, param, t, dt) end if ! Write out each data type frame - call cb%write_frame(iu, param, t, dt) - call pl%write_frame(iu, param, t, dt) - call tp%write_frame(iu, param, t, dt) + call cb%write_frame(iu, param) + call pl%write_frame(iu, param) + call tp%write_frame(iu, param) deallocate(cb, pl, tp) @@ -1361,7 +1347,7 @@ module subroutine io_write_frame_system(self, iu, param, t, dt) return end subroutine io_write_frame_system - module subroutine io_write_hdr(iu, t, npl, ntp, out_form, out_type) + subroutine io_write_hdr(iu, t, npl, ntp, out_form, out_type) !! author: The Purdue Swiftest Team - David A. Minton, Carlisle A. Wishard, Jennifer L.L. Pouplin, and Jacob R. Elliott !! !! Write frame header to output binary file diff --git a/src/modules/helio_classes.f90 b/src/modules/helio_classes.f90 index 15c61343b..ad10018a2 100644 --- a/src/modules/helio_classes.f90 +++ b/src/modules/helio_classes.f90 @@ -38,7 +38,6 @@ module helio_classes 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 :: getacch_int => helio_getacch_int_pl !! Compute direct cross term heliocentric accelerations of planets 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 @@ -58,7 +57,6 @@ module helio_classes 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 :: getacch_int => helio_getacch_int_tp !! Compute direct cross term heliocentric accelerations of test particles 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 end type helio_tp @@ -91,23 +89,21 @@ module subroutine helio_coord_vh2vb_tp(self, vbcb) end subroutine helio_coord_vh2vb_tp module subroutine helio_drift_pl(self, system, param, dt) - use swiftest_classes, only : swiftest_parameters - use whm_classes, only : whm_nbody_system + use swiftest_classes, only : swiftest_nbody_system, swiftest_parameters implicit none - 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 + 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) :: dt !! Stepsize end subroutine helio_drift_pl module subroutine helio_drift_tp(self, system, param, dt) - use swiftest_classes, only : swiftest_parameters - use whm_classes, only : whm_nbody_system + use swiftest_classes, only : swiftest_nbody_system, swiftest_parameters implicit none - 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 + 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 of + real(DP), intent(in) :: dt !! Stepsize end subroutine helio_drift_tp module subroutine helio_drift_linear_pl(self, cb, dt, pt) @@ -127,41 +123,24 @@ module subroutine helio_drift_linear_tp(self, dt, pt) end subroutine helio_drift_linear_tp module subroutine helio_getacch_pl(self, system, param, t) - use swiftest_classes, only : swiftest_parameters - use whm_classes, only : whm_nbody_system + use swiftest_classes, only : swiftest_parameters, swiftest_nbody_system implicit none - 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 + class(helio_pl), intent(inout) :: self !! Helio massive body particle data structure + 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 end subroutine helio_getacch_pl - module subroutine helio_getacch_tp(self, system, param, t, xh) - use swiftest_classes, only : swiftest_parameters - use whm_classes, only : whm_nbody_system + 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(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 + class(helio_tp), intent(inout) :: self !! Helio test particle data structure + 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_getacch_int_pl(self, t) - implicit none - class(helio_pl), intent(inout) :: self !! Helio massive body particle data structure - real(DP), intent(in) :: t !! Current time - 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(whm_pl), intent(inout) :: pl !! WhM massive body particle data structure - real(DP), intent(in) :: t !! Current time - real(DP), dimension(:,:), intent(in) :: xh !! Heliocentric positions of planet - end subroutine helio_getacch_int_tp - module subroutine helio_setup_pl(self, n) implicit none class(helio_pl), intent(inout) :: self !! Helio massive body object diff --git a/src/modules/rmvs_classes.f90 b/src/modules/rmvs_classes.f90 index 8a1f94a8e..8fc57cadf 100644 --- a/src/modules/rmvs_classes.f90 +++ b/src/modules/rmvs_classes.f90 @@ -21,12 +21,15 @@ module rmvs_classes !******************************************************************************************************************************** type, public, extends(whm_nbody_system) :: rmvs_nbody_system !> In the RMVS integrator, only test particles are discarded - logical :: lplanetocentric = .false. !! Flag that indicates that the object is a planetocentric set of masive bodies used for close encounter calculations + logical :: lplanetocentric = .false. !! Flag that indicates that the object is a planetocentric set of masive bodies used for close encounter calculations + real(DP) :: rts !! fraction of Hill's sphere radius to use as radius of encounter region + real(DP), dimension(:,:), allocatable :: vbeg !! Planet velocities at beginning ot step 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 end type rmvs_nbody_system type, private :: rmvs_interp @@ -40,9 +43,9 @@ module rmvs_classes !******************************************************************************************************************************* !> RMVS central body particle class type, public, extends(whm_cb) :: rmvs_cb - logical :: lplanetocentric = .false. !! Flag that indicates that the object is a planetocentric set of test particles used for close encounter calculations - type(rmvs_interp), dimension(:), allocatable :: outer !! interpolated heliocentric central body position for outer encounters - type(rmvs_interp), dimension(:), allocatable :: inner !! interpolated heliocentric central body position for inner encounters + type(rmvs_interp), dimension(:), allocatable :: outer !! interpolated heliocentric central body position for outer encounters + type(rmvs_interp), dimension(:), allocatable :: inner !! interpolated heliocentric central body position for inner encounters + logical :: lplanetocentric = .false. !! Flag that indicates that the object is a planetocentric set of masive bodies used for close encounter calculations end type rmvs_cb !******************************************************************************************************************************** @@ -54,24 +57,22 @@ module rmvs_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 rmvs_setup_tp and rmvs_spill_tp ! encounter steps) - logical :: lplanetocentric = .false. !! Flag that indicates that the object is a planetocentric set of test particles used for close encounter calculations logical, dimension(:), allocatable :: lperi !! planetocentric pericenter passage flag (persistent for a full rmvs time step) over a full RMVS time step) integer(I4B), dimension(:), allocatable :: plperP !! index of planet associated with pericenter distance peri (persistent over a full RMVS time step) integer(I4B), dimension(:), allocatable :: plencP !! index of planet that test particle is encountering (not persistent for a full RMVS time step) - real(DP), dimension(:,:), allocatable :: vbeg !! Planet velocities at beginning ot step ! The following are used to correctly set the oblateness values of the acceleration during an inner encounter with a planet - type(rmvs_cb) :: cb_heliocentric !! Copy of original central body object passed to close encounter (used for oblateness acceleration during planetocentric encoountters) - real(DP), dimension(:,:), allocatable :: xheliocentric !! original heliocentric position (used for oblateness calculation during close encounters) - integer(I4B) :: index !! inner substep number within current set - integer(I4B) :: ipleP !! index value of encountering planet + type(rmvs_cb) :: cb_heliocentric !! Copy of original central body object passed to close encounter (used for oblateness acceleration during planetocentric encoountters) + real(DP), dimension(:,:), allocatable :: xheliocentric !! original heliocentric position (used for oblateness calculation during close encounters) + integer(I4B) :: index !! inner substep number within current set + integer(I4B) :: ipleP !! index value of encountering planet + logical :: lplanetocentric = .false. !! Flag that indicates that the object is a planetocentric set of masive bodies used for close encounter calculations contains procedure, public :: discard_pl => rmvs_discard_pl_tp 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 !! if the test particle is undergoing a close encounter or not - 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 :: 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) end type rmvs_tp @@ -80,81 +81,82 @@ module rmvs_classes ! rmvs_pl class definitions and method interfaces !******************************************************************************************************************************* - !> RMVS massive body particle class - type, private, extends(whm_pl) :: rmvs_base_pl - logical :: lplanetocentric = .false. !! Flag that indicates that the object is a planetocentric set of masive bodies used for close encounter calculations - integer(I4B), dimension(:), allocatable :: nenc !! number of test particles encountering planet this full rmvs time step - integer(I4B), dimension(:), allocatable :: tpenc1P !! index of first test particle encountering planet - integer(I4B), dimension(:), allocatable :: plind ! Connects the planetocentric indices back to the heliocentric planet list - type(rmvs_interp), dimension(:), allocatable :: outer !! interpolated heliocentric central body position for outer encounters - type(rmvs_interp), dimension(:), allocatable :: inner !! interpolated heliocentric central body position for inner encounters + type, private, extends(whm_pl) :: rmvs_pl + integer(I4B), dimension(:), allocatable :: nenc !! number of test particles encountering planet this full rmvs time step + integer(I4B), dimension(:), allocatable :: tpenc1P !! index of first test particle encountering planet + integer(I4B), dimension(:), allocatable :: plind ! Connects the planetocentric indices back to the heliocentric planet list + type(rmvs_interp), dimension(:), allocatable :: outer !! interpolated heliocentric central body position for outer encounters + type(rmvs_interp), dimension(:), allocatable :: inner !! interpolated heliocentric central body position for inner encounters + class(rmvs_nbody_system), dimension(:), allocatable :: planetocentric + logical :: lplanetocentric = .false. !! Flag that indicates that the object is a planetocentric set of masive bodies used for close encounter calculations contains procedure, public :: fill => rmvs_fill_pl !! "Fills" bodies from one object into another depending on the results of a mask (uses the MERGE intrinsic) procedure, public :: setup => rmvs_setup_pl !! Constructor method - Allocates space for number of particles procedure, public :: spill => rmvs_spill_pl !! "Spills" bodies from one object to another depending on the results of a mask (uses the PACK intrinsic) - end type rmvs_base_pl - - type, public :: rmvs_encounter_system - class(rmvs_cb), allocatable :: cb - class(rmvs_tp), allocatable :: tp - class(rmvs_base_pl), allocatable :: pl - end type rmvs_encounter_system - - type, public, extends(rmvs_base_pl) :: rmvs_pl - type(rmvs_encounter_system), dimension(:), allocatable :: planetocentric - !! Note to developers: If you add componenets to this class, be sure to update methods and subroutines that traverse the - !! component list, such as rmvs_setup_pl and rmvs_spill_pl end type rmvs_pl - + interface - module subroutine rmvs_discard_pl_tp(self, pl, t, dt) - use swiftest_classes, only : swiftest_cb, swiftest_pl, swiftest_parameters + module subroutine rmvs_discard_pl_tp(self, system, param) + use swiftest_classes, only : swiftest_nbody_system, swiftest_parameters implicit none - class(rmvs_tp), intent(inout) :: self - class(swiftest_pl), intent(inout) :: pl !! WHM massive body particle data structure. - real(DP), intent(in) :: t !! Current simulation time - real(DP), intent(in) :: dt !! Stepsize + class(rmvs_tp), intent(inout) :: self !! RMVS test particle object + class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object + class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters end subroutine rmvs_discard_pl_tp - module subroutine rmvs_getacch_tp(self, system, param, t, xh) - use swiftest_classes, only : swiftest_cb, swiftest_parameters - use whm_classes, only : whm_nbody_system + module function rmvs_encounter_check_tp(self, system, dt) result(lencounter) implicit none - class(rmvs_tp), intent(inout) :: self !! RMVS test particle data structure - class(whm_nbody_system), intent(inout) :: system !! Swiftest central body particle data structuree - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters of parameter - real(DP), intent(in) :: t !! Current time - real(DP), dimension(:,:), intent(in) :: xh !! Heliocentric positions of planets - end subroutine rmvs_getacch_tp + class(rmvs_tp), intent(inout) :: self !! RMVS test particle object + class(rmvs_nbody_system), intent(inout) :: system !! RMVS nbody system object + real(DP), intent(in) :: dt !! step size + logical :: lencounter !! Returns true if there is at least one close encounter + end function rmvs_encounter_check_tp - module subroutine rmvs_setup_system(self, param) - use swiftest_classes, only : swiftest_parameters + module subroutine rmvs_fill_pl(self, inserts, lfill_list) + use swiftest_classes, only : swiftest_body implicit none - class(rmvs_nbody_system), intent(inout) :: self !! RMVS system object - class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters - end subroutine rmvs_setup_system + class(rmvs_pl), intent(inout) :: self !! RMVS massive body object + class(swiftest_body), intent(inout) :: inserts !! Inserted object + logical, dimension(:), intent(in) :: lfill_list !! Logical array of bodies to merge into the keeps + end subroutine rmvs_fill_pl + + module subroutine rmvs_fill_tp(self, inserts, lfill_list) + use swiftest_classes, only : swiftest_body + implicit none + class(rmvs_tp), intent(inout) :: self !! RMVS massive body object + class(swiftest_body), intent(inout) :: inserts !! Inserted object + logical, dimension(:), intent(in) :: lfill_list !! Logical array of bodies to merge into the keeps + end subroutine rmvs_fill_tp + + module subroutine rmvs_getacch_tp(self, system, param, t, xhp) + use swiftest_classes, only : swiftest_nbody_system, swiftest_parameters + implicit none + class(rmvs_tp), intent(inout) :: self !! RMVS test particle data structure + class(swiftest_nbody_system), intent(inout) :: system !! Swiftest central body particle data structuree + 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 at current substep + end subroutine rmvs_getacch_tp module subroutine rmvs_setup_pl(self,n) implicit none - class(rmvs_base_pl), intent(inout) :: self !! RMVS test particle object + class(rmvs_pl), intent(inout) :: self !! RMVS test particle object integer, intent(in) :: n !! Number of test particles to allocate end subroutine rmvs_setup_pl module subroutine rmvs_setup_set_beg_end(self, xbeg, xend, vbeg) implicit none - class(rmvs_tp), intent(inout) :: self !! RMVS test particle object + class(rmvs_nbody_system), intent(inout) :: self !! RMVS nbody system object real(DP), dimension(:,:), intent(in), optional :: xbeg, xend, vbeg end subroutine rmvs_setup_set_beg_end - module subroutine rmvs_step_system(self, param, t, dt) + module subroutine rmvs_setup_system(self, param) use swiftest_classes, only : swiftest_parameters implicit none - class(rmvs_nbody_system), intent(inout) :: self !! RMVS nbody system object + class(rmvs_nbody_system), intent(inout) :: self !! RMVS 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 - end subroutine rmvs_step_system + end subroutine rmvs_setup_system module subroutine rmvs_setup_tp(self,n) implicit none @@ -162,32 +164,14 @@ module subroutine rmvs_setup_tp(self,n) integer, intent(in) :: n !! Number of test particles to allocate end subroutine rmvs_setup_tp - module function rmvs_encounter_check_tp(self, cb, pl, dt, rts) result(lencounter) - implicit none - class(rmvs_tp), intent(inout) :: self !! RMVS test particle object - class(rmvs_cb), intent(inout) :: cb !! RMVS central body object - class(rmvs_pl), intent(inout) :: pl !! RMVS massive body object - real(DP), intent(in) :: dt !! step size - real(DP), intent(in) :: rts !! fraction of Hill's sphere radius to use as radius of encounter regio - logical :: lencounter !! Returns true if there is at least one close encounter - end function rmvs_encounter_check_tp - module subroutine rmvs_spill_pl(self, discards, lspill_list) use swiftest_classes, only : swiftest_body implicit none - class(rmvs_base_pl), intent(inout) :: self !! RMVS massive body object + class(rmvs_pl), intent(inout) :: self !! RMVS massive body object class(swiftest_body), intent(inout) :: discards !! Discarded object - logical, dimension(:), intent(in) :: lspill_list !! Logical array of bodies to spill into the discards + logical, dimension(:), intent(in) :: lspill_list !! Logical array of bodies to spill into the discards end subroutine rmvs_spill_pl - module subroutine rmvs_fill_pl(self, inserts, lfill_list) - use swiftest_classes, only : swiftest_body - implicit none - class(rmvs_base_pl), intent(inout) :: self !! RMVS massive body object - class(swiftest_body), intent(inout) :: inserts !! Inserted object - logical, dimension(:), intent(in) :: lfill_list !! Logical array of bodies to merge into the keeps - end subroutine rmvs_fill_pl - module subroutine rmvs_spill_tp(self, discards, lspill_list) use swiftest_classes, only : swiftest_body implicit none @@ -196,13 +180,15 @@ module subroutine rmvs_spill_tp(self, discards, lspill_list) logical, dimension(:), intent(in) :: lspill_list !! Logical array of bodies to spill into the discards end subroutine rmvs_spill_tp - module subroutine rmvs_fill_tp(self, inserts, lfill_list) - use swiftest_classes, only : swiftest_body + module subroutine rmvs_step_system(self, param, t, dt) + use swiftest_classes, only : swiftest_parameters implicit none - class(rmvs_tp), intent(inout) :: self !! RMVS massive body object - class(swiftest_body), intent(inout) :: inserts !! Inserted object - logical, dimension(:), intent(in) :: lfill_list !! Logical array of bodies to merge into the keeps - end subroutine rmvs_fill_tp + class(rmvs_nbody_system), intent(inout) :: self !! RMVS 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 + end subroutine rmvs_step_system + end interface end module rmvs_classes diff --git a/src/modules/swiftest_classes.f90 b/src/modules/swiftest_classes.f90 index 905044f17..150408537 100644 --- a/src/modules/swiftest_classes.f90 +++ b/src/modules/swiftest_classes.f90 @@ -273,13 +273,12 @@ subroutine abstract_initialize(self, param) class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters end subroutine abstract_initialize - subroutine abstract_read_frame(self, iu, param, form, t, ierr) + 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 integer(I4B), intent(inout) :: iu !! Unit number for the output file to write frame to class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters character(*), intent(in) :: form !! Input format code ("XV" or "EL") - real(DP), intent(out) :: t !! Simulation time integer(I4B), intent(out) :: ierr !! Error code end subroutine abstract_read_frame @@ -317,38 +316,25 @@ end subroutine abstract_write_frame end interface interface - module subroutine discard_peri_tp(self, cb, pl, param, t, msys) - implicit none - class(swiftest_tp), intent(inout) :: self !! Swiftest test particle object - class(swiftest_cb), intent(inout) :: cb !! Swiftest central body object - class(swiftest_pl), intent(inout) :: pl !! Swiftest massive body object - class(swiftest_parameters), intent(in) :: param !! parameters - real(DP), intent(in) :: t !! Current simulation tim - real(DP), intent(in) :: msys !! Total system mass - end subroutine discard_peri_tp - - module subroutine discard_pl_close(dx, dv, dt, r2crit, iflag, r2min) + module subroutine discard_peri_tp(self, system, param) implicit none - real(DP), dimension(:), intent(in) :: dx, dv - real(DP), intent(in) :: dt, r2crit - integer(I4B), intent(out) :: iflag - real(DP), intent(out) :: r2min - end subroutine discard_pl_close + class(swiftest_tp), intent(inout) :: self !! Swiftest test particle object + class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object + class(swiftest_parameters), intent(in) :: param !! Current run configuration parameter + end subroutine discard_peri_tp - module subroutine discard_pl_tp(self, pl, t, dt) + module subroutine discard_pl_tp(self, system, param) implicit none - class(swiftest_tp), intent(inout) :: self !! Swiftest test particle object - class(swiftest_pl), intent(inout) :: pl !! Swiftest massive body object - real(DP), intent(in) :: t !! Current simulation tim - real(DP), intent(in) :: dt !! Stepsize + class(swiftest_tp), intent(inout) :: self !! Swiftest test particle object + class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object + class(swiftest_parameters), intent(in) :: param !! Current run configuration parameter end subroutine discard_pl_tp - module subroutine discard_sun_tp(self, cb, param, msys) + module subroutine discard_sun_tp(self, system, param) implicit none - class(swiftest_tp), intent(inout) :: self !! Swiftest test particle object - class(swiftest_cb), intent(inout) :: cb !! Swiftest central body object - class(swiftest_parameters), intent(in) :: param !! parameters - real(DP), intent(in) :: msys !! Total system mass + class(swiftest_tp), intent(inout) :: self !! Swiftest test particle object + class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object + class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters end subroutine discard_sun_tp module subroutine discard_system(self, param) @@ -468,56 +454,33 @@ module subroutine io_read_param_in(self, param_file_name) character(len=*), intent(in) :: param_file_name !! Parameter input file name (i.e. param.in) end subroutine io_read_param_in - module function io_read_encounter(t, name1, name2, mass1, mass2, radius1, radius2, & - xh1, xh2, vh1, vh2, encounter_file, out_type) result(ierr) - implicit none - integer(I4B) :: ierr - integer(I4B), intent(out) :: name1, name2 - real(DP), intent(out) :: t, mass1, mass2, radius1, radius2 - real(DP), dimension(NDIM), intent(out) :: xh1, xh2, vh1, vh2 - character(*), intent(in) :: encounter_file,out_type - end function io_read_encounter - - module subroutine io_read_frame_body(self, iu, param, form, t, ierr) + module subroutine io_read_frame_body(self, iu, param, form, ierr) implicit none class(swiftest_body), intent(inout) :: self !! Swiftest particle object integer(I4B), intent(inout) :: iu !! Unit number for the output file to write frame to class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters character(*), intent(in) :: form !! Input format code ("XV" or "EL") - real(DP), intent(out) :: t !! Simulation time integer(I4B), intent(out) :: ierr !! Error code end subroutine io_read_frame_body - module subroutine io_read_frame_cb(self, iu, param, form, t, ierr) + module subroutine io_read_frame_cb(self, iu, param, form, ierr) implicit none class(swiftest_cb), intent(inout) :: self !! Swiftest central body object integer(I4B), intent(inout) :: iu !! Unit number for the output file to write frame to class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters character(*), intent(in) :: form !! Input format code ("XV" or "EL") - real(DP), intent(out) :: t !! Simulation time integer(I4B), intent(out) :: ierr !! Error code end subroutine io_read_frame_cb - module subroutine io_read_frame_system(self, iu, param, form, t, ierr) + module subroutine io_read_frame_system(self, iu, param, form, ierr) implicit none class(swiftest_nbody_system),intent(inout) :: self !! Swiftest system object integer(I4B), intent(inout) :: iu !! Unit number for the output file to write frame to class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters character(*), intent(in) :: form !! Input format code ("XV" or "EL") - real(DP), intent(out) :: t !! Current simulation time integer(I4B), intent(out) :: ierr !! Error code end subroutine io_read_frame_system - module function io_read_hdr(iu, t, npl, ntp, out_form, out_type) result(ierr) - implicit none - integer(I4B) :: ierr - integer(I4B), intent(in) :: iu - integer(I4B), intent(out) :: npl, ntp - character(*), intent(out) :: out_form - real(DP), intent(out) :: t - character(*), intent(in) :: out_type - end function io_read_hdr - module subroutine io_read_initialize_system(self, param) implicit none class(swiftest_nbody_system), intent(inout) :: self !! Swiftest system object @@ -531,26 +494,17 @@ module subroutine io_write_discard(self, param, discards) class(swiftest_body), intent(inout) :: discards !! Swiftest discard object end subroutine io_write_discard - module subroutine io_write_encounter(t, name1, name2, mass1, mass2, radius1, radius2, & - xh1, xh2, vh1, vh2, encounter_file, out_type) - implicit none - integer(I4B), intent(in) :: name1, name2 - real(DP), intent(in) :: t, mass1, mass2, radius1, radius2 - real(DP), dimension(NDIM), intent(in) :: xh1, xh2, vh1, vh2 - character(*), intent(in) :: encounter_file, out_type - end subroutine io_write_encounter - module subroutine io_write_frame_body(self, iu, param) implicit none - class(swiftest_body), intent(in) :: self !! Swiftest particle object - integer(I4B), intent(inout) :: iu !! Unit number for the output file to write frame to + class(swiftest_body), intent(in) :: self !! Swiftest particle object + integer(I4B), intent(inout) :: iu !! Unit number for the output file to write frame to class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters end subroutine io_write_frame_body module subroutine io_write_frame_cb(self, iu, param) implicit none - class(swiftest_cb), intent(in) :: self !! Swiftest central body object - integer(I4B), intent(inout) :: iu !! Unit number for the output file to write frame to + class(swiftest_cb), intent(in) :: self !! Swiftest central body object + integer(I4B), intent(inout) :: iu !! Unit number for the output file to write frame to class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters end subroutine io_write_frame_cb @@ -561,15 +515,6 @@ 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 io_write_hdr(iu, t, npl, ntp, out_form, out_type) - integer(I4B), intent(in) :: iu !! Output file unit number - real(DP), intent(in) :: t !! Current time of simulation - integer(I4B), intent(in) :: npl !! Number of massive bodies - integer(I4B), intent(in) :: ntp !! Number of test particles - character(*), intent(in) :: out_form !! Output format type ("EL" or "XV") - character(*), intent(in) :: out_type !! Output file format type (REAL4, REAL8 - see swiftest module for symbolic name definitions) - end subroutine io_write_hdr - module subroutine obl_acc_body(self, cb) implicit none class(swiftest_body), intent(inout) :: self !! Swiftest generic body object @@ -660,12 +605,12 @@ module subroutine setup_tp(self, n) integer, intent(in) :: n !! Number of bodies to allocate space for end subroutine setup_tp - module subroutine user_getacch_body(self, cb, param, t) + module subroutine user_getacch_body(self, system, param, t) implicit none - class(swiftest_body), intent(inout) :: self !! Swiftest massive body particle data structure - class(swiftest_cb), intent(inout) :: cb !! Swiftest central body particle data structuree - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters of - real(DP), intent(in) :: t !! Current time + class(swiftest_body), intent(inout) :: self !! Swiftest massive body particle data structure + class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nobody system object + class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters of + real(DP), intent(in) :: t !! Current time end subroutine user_getacch_body module subroutine util_coord_b2h_pl(self, cb) diff --git a/src/modules/whm_classes.f90 b/src/modules/whm_classes.f90 index e411201ae..031318161 100644 --- a/src/modules/whm_classes.f90 +++ b/src/modules/whm_classes.f90 @@ -9,19 +9,6 @@ module whm_classes implicit none public - !******************************************************************************************************************************** - ! whm_nbody_system class definitions and method interfaces - !******************************************************************************************************************************** - !> An abstract class for the WHM integrator nbody system - type, public, extends(swiftest_nbody_system) :: whm_nbody_system - !> In the WHM integrator, only test particles are discarded - class(swiftest_tp), allocatable :: tp_discards - 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 - end type whm_nbody_system !******************************************************************************************************************************** ! whm_cb class definitions and method interfaces @@ -71,8 +58,6 @@ module whm_classes type, public, extends(swiftest_tp) :: whm_tp !! 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_tp and whm_spill_tp - real(DP), dimension(:,:), allocatable :: xbeg, xend - logical :: beg logical :: lfirst = .true. contains private @@ -82,67 +67,27 @@ module whm_classes 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 - procedure, public :: set_beg_end => whm_setup_set_beg_end !! Sets the beginning and ending positions of planets. 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 - interface - !> WHM massive body constructor method - module subroutine whm_setup_pl(self,n) - implicit none - class(whm_pl), intent(inout) :: self !! Swiftest test particle object - integer(I4B), intent(in) :: n !! Number of test particles to allocate - end subroutine whm_setup_pl - - module subroutine whm_setup_set_mu_eta_pl(self, cb) - use swiftest_classes, only : swiftest_cb - implicit none - class(whm_pl), intent(inout) :: self !! WHM massive body object - class(swiftest_cb), intent(inout) :: cb !! Swiftest central body object - end subroutine whm_setup_set_mu_eta_pl - - module subroutine whm_setup_set_ir3j(self) - implicit none - class(whm_pl), intent(inout) :: self !! WHM massive body object - end subroutine whm_setup_set_ir3j - - module subroutine whm_step_pl(self, system, param, t, dt) - 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 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 - end subroutine whm_step_pl - - !> Get heliocentric accelration of massive bodies - module subroutine whm_getacch_pl(self, system, param, t) - use swiftest_classes, only : swiftest_cb, swiftest_parameters - implicit none - class(whm_pl), intent(inout) :: self !! WHM 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 whm_getacch_pl - - module subroutine whm_drift_pl(self, system, param, dt) - use swiftest_classes, only : swiftest_cb, swiftest_parameters - implicit none - class(whm_pl), intent(inout) :: self !! WHM 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) :: dt !! Stepsize - end subroutine whm_drift_pl - - module subroutine whm_getacch_int_pl(self, cb) - use swiftest_classes, only : swiftest_cb - implicit none - class(whm_pl), intent(inout) :: self !! WHM massive body particle data structure - class(swiftest_cb), intent(inout) :: cb !! Swiftest central body particle data structure - end subroutine whm_getacch_int_pl + !******************************************************************************************************************************** + ! whm_nbody_system class definitions and method interfaces + !******************************************************************************************************************************** + !> An abstract class for the WHM integrator nbody system + type, public, extends(swiftest_nbody_system) :: whm_nbody_system + !> In the WHM integrator, only test particles are discarded + class(whm_tp), allocatable :: tp_discards !! WHM test particle object that + real(DP), dimension(:,:), allocatable :: xbeg, xend !! Positions of massive bodies at beginning and end of a step. Required in order to separate the test particle step from the massive body step + 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. + end type whm_nbody_system + interface module subroutine whm_coord_h2j_pl(self, cb) use swiftest_classes, only : swiftest_cb implicit none @@ -164,27 +109,82 @@ module subroutine whm_coord_vh2vj_pl(self, cb) class(swiftest_cb), intent(inout) :: cb !! Swiftest central body particle data structuree end subroutine whm_coord_vh2vj_pl + module subroutine whm_drift_pl(self, system, param, dt) + use swiftest_classes, only : swiftest_nbody_system, swiftest_parameters + implicit none + class(whm_pl), intent(inout) :: self !! WHM massive body particle data structure + 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) :: dt !! Stepsize + end subroutine whm_drift_pl + + module subroutine whm_drift_tp(self, system, param, dt) + use swiftest_classes, only : swiftest_nbody_system, swiftest_parameters + implicit none + class(whm_tp), intent(inout) :: self !! WHM test particle data structure + 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) :: dt !! Stepsize + end subroutine whm_drift_tp + + module subroutine whm_fill_pl(self, inserts, lfill_list) + use swiftest_classes, only : swiftest_body + implicit none + class(whm_pl), intent(inout) :: self !! WHM massive body object + class(swiftest_body), intent(inout) :: inserts !! inserted object + logical, dimension(:), intent(in) :: lfill_list !! Logical array of bodies to merge into the keeps + end subroutine whm_fill_pl + + !> Get heliocentric accelration of massive bodies + module subroutine whm_getacch_pl(self, system, param, t) + use swiftest_classes, only : swiftest_cb, swiftest_parameters + implicit none + class(whm_pl), intent(inout) :: self !! WHM massive body particle data structure + 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 + end subroutine whm_getacch_pl + + !> Get heliocentric accelration of the test particle + module subroutine whm_getacch_tp(self, system, param, t, xhp) + use swiftest_classes, only : swiftest_cb, swiftest_parameters + implicit none + class(whm_tp), intent(inout) :: self !! WHM test particle data structure + 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 time + real(DP), dimension(:,:), intent(in) :: xhp !! Heliocentric positions of planets at the current substep + end subroutine whm_getacch_tp + module subroutine whm_gr_getacch_pl(self, param) use swiftest_classes, only : swiftest_cb, swiftest_parameters implicit none - class(whm_pl), intent(inout) :: self !! WHM massive body particle data structure + class(whm_pl), intent(inout) :: self !! WHM massive body particle data structure class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters of end subroutine whm_gr_getacch_pl + module subroutine whm_gr_getacch_tp(self, param) + use swiftest_classes, only : swiftest_cb, swiftest_parameters + implicit none + class(whm_tp), intent(inout) :: self !! WHM test particle data structure + class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters + end subroutine whm_gr_getacch_tp + module pure subroutine whm_gr_p4_pl(self, param, dt) use swiftest_classes, only : swiftest_parameters implicit none - class(whm_pl), intent(inout) :: self !! Swiftest particle object + class(whm_pl), intent(inout) :: self !! Swiftest particle object class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters of on parameters - real(DP), intent(in) :: dt !! Step size + real(DP), intent(in) :: dt !! Step size end subroutine whm_gr_p4_pl - module pure subroutine whm_gr_vh2pv_pl(self, param) + module pure subroutine whm_gr_p4_tp(self, param, dt) use swiftest_classes, only : swiftest_parameters implicit none - class(whm_pl), intent(inout) :: self !! Swiftest particle object + class(whm_tp), intent(inout) :: self !! WHM test particle object class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters of on parameters - end subroutine whm_gr_vh2pv_pl + real(DP), intent(in) :: dt !! Step size + end subroutine whm_gr_p4_tp module pure subroutine whm_gr_pv2vh_pl(self, param) use swiftest_classes, only : swiftest_parameters @@ -193,78 +193,53 @@ module pure subroutine whm_gr_pv2vh_pl(self, param) class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters of on parameters end subroutine whm_gr_pv2vh_pl - !> WHM test particle constructor - module subroutine whm_setup_tp(self,n) + module pure subroutine whm_gr_pv2vh_tp(self, param) + use swiftest_classes, only : swiftest_parameters implicit none - class(whm_tp), intent(inout) :: self !! WHM test particle data structure - integer, intent(in) :: n !! Number of test particles to allocate - end subroutine whm_setup_tp + class(whm_tp), intent(inout) :: self !! WHM test particle object + class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters of on parameters + end subroutine whm_gr_pv2vh_tp - module subroutine whm_drift_tp(self, system, param, dt) - use swiftest_classes, only : swiftest_cb, swiftest_parameters + module pure subroutine whm_gr_vh2pv_pl(self, param) + use swiftest_classes, only : swiftest_parameters implicit none - class(whm_tp), intent(inout) :: self !! WHM 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 of - real(DP), intent(in) :: dt !! Stepsize - end subroutine whm_drift_tp + class(whm_pl), intent(inout) :: self !! Swiftest particle object + class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters of on parameters + end subroutine whm_gr_vh2pv_pl - !> Get heliocentric accelration of the test particle - module subroutine whm_getacch_tp(self, system, param, t, xh) - use swiftest_classes, only : swiftest_cb, swiftest_parameters + module pure subroutine whm_gr_vh2pv_tp(self, param) + use swiftest_classes, only : swiftest_parameters implicit none - class(whm_tp), intent(inout) :: self !! WHM 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 of - real(DP), intent(in) :: t !! Current time - real(DP), dimension(:,:), intent(in) :: xh !! Heliocentric positions of planets - end subroutine whm_getacch_tp + class(whm_tp), intent(inout) :: self !! WHM test particle object + class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters of on parameters + end subroutine whm_gr_vh2pv_tp - module subroutine whm_step_tp(self, system, param, t, dt) - use swiftest_classes, only : swiftest_nbody_system, swiftest_parameters + + !> Reads WHM massive body object in from file + module subroutine whm_setup_pl(self,n) implicit none - class(whm_tp), intent(inout) :: self !! WHM test particle data structure - 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 - real(DP), intent(in) :: dt !! Stepsize - end subroutine whm_step_tp + class(whm_pl), intent(inout) :: self !! Swiftest test particle object + integer(I4B), intent(in) :: n !! Number of test particles to allocate + end subroutine whm_setup_pl module subroutine whm_setup_set_beg_end(self, xbeg, xend, vbeg) implicit none - class(whm_tp), intent(inout) :: self !! Swiftest test particle object + class(whm_nbody_system), intent(inout) :: self !! WHM nbody system object real(DP), dimension(:,:), intent(in), optional :: xbeg, xend real(DP), dimension(:,:), intent(in), optional :: vbeg ! vbeg is an unused variable to keep this method forward compatible with RMVS end subroutine whm_setup_set_beg_end - module subroutine whm_gr_getacch_tp(self, param) - use swiftest_classes, only : swiftest_cb, swiftest_parameters - implicit none - class(whm_tp), intent(inout) :: self !! WHM test particle data structure - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters - end subroutine whm_gr_getacch_tp - - module pure subroutine whm_gr_p4_tp(self, param, dt) - use swiftest_classes, only : swiftest_parameters - implicit none - class(whm_tp), intent(inout) :: self !! WHM test particle object - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters of on parameters - real(DP), intent(in) :: dt !! Step size - end subroutine whm_gr_p4_tp - - module pure subroutine whm_gr_vh2pv_tp(self, param) - use swiftest_classes, only : swiftest_parameters + module subroutine whm_setup_set_ir3j(self) implicit none - class(whm_tp), intent(inout) :: self !! WHM test particle object - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters of on parameters - end subroutine whm_gr_vh2pv_tp + class(whm_pl), intent(inout) :: self !! WHM massive body object + end subroutine whm_setup_set_ir3j - module pure subroutine whm_gr_pv2vh_tp(self, param) - use swiftest_classes, only : swiftest_parameters + module subroutine whm_setup_set_mu_eta_pl(self, cb) + use swiftest_classes, only : swiftest_cb implicit none - class(whm_tp), intent(inout) :: self !! WHM test particle object - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters of on parameters - end subroutine whm_gr_pv2vh_tp + class(whm_pl), intent(inout) :: self !! WHM massive body object + class(swiftest_cb), intent(inout) :: cb !! Swiftest central body object + end subroutine whm_setup_set_mu_eta_pl module subroutine whm_setup_system(self, param) use swiftest_classes, only : swiftest_parameters @@ -273,6 +248,33 @@ module subroutine whm_setup_system(self, param) class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters of on parameters end subroutine whm_setup_system + !> Reads WHM test particle object in from file + module subroutine whm_setup_tp(self,n) + implicit none + class(whm_tp), intent(inout) :: self !! WHM test particle data structure + integer, intent(in) :: n !! Number of test particles to allocate + end subroutine whm_setup_tp + + module subroutine whm_step_pl(self, system, param, t, dt) + 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 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 + end subroutine whm_step_pl + + module subroutine whm_step_tp(self, system, param, t, dt) + use swiftest_classes, only : swiftest_nbody_system, swiftest_parameters + implicit none + class(whm_tp), intent(inout) :: self !! WHM test particle data structure + 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 + real(DP), intent(in) :: dt !! Stepsize + end subroutine whm_step_tp + module subroutine whm_spill_pl(self, discards, lspill_list) use swiftest_classes, only : swiftest_body implicit none @@ -281,14 +283,6 @@ module subroutine whm_spill_pl(self, discards, lspill_list) logical, dimension(:), intent(in) :: lspill_list !! Logical array of bodies to spill into the discards end subroutine whm_spill_pl - module subroutine whm_fill_pl(self, inserts, lfill_list) - use swiftest_classes, only : swiftest_body - implicit none - class(whm_pl), intent(inout) :: self !! WHM massive body object - class(swiftest_body), intent(inout) :: inserts !! inserted object - logical, dimension(:), intent(in) :: lfill_list !! Logical array of bodies to merge into the keeps - end subroutine whm_fill_pl - !> Steps the Swiftest nbody system forward in time one stepsize module subroutine whm_step_system(self, param, t, dt) use swiftest_classes, only : swiftest_parameters diff --git a/src/rmvs/rmvs_discard.f90 b/src/rmvs/rmvs_discard.f90 index 4f75cd645..dc8b4b9bb 100644 --- a/src/rmvs/rmvs_discard.f90 +++ b/src/rmvs/rmvs_discard.f90 @@ -1,7 +1,7 @@ submodule(rmvs_classes) s_rmvs_discard use swiftest contains - module subroutine rmvs_discard_pl_tp(self, pl, t, dt) + module subroutine rmvs_discard_pl_tp(self, system, param) !! author: David A. Minton !! !! Check to see if test particles should be discarded based on pericenter passage distances with respect to @@ -11,14 +11,13 @@ module subroutine rmvs_discard_pl_tp(self, pl, t, dt) !! Adapted from Hal Levison's Swift routine rmvs_discard_pl.f90 implicit none ! Arguments - class(rmvs_tp), intent(inout) :: self - class(swiftest_pl), intent(inout) :: pl !! WHM massive body particle data structure. - real(DP), intent(in) :: t !! Current simulation time - real(DP), intent(in) :: dt !! Stepsize + class(rmvs_tp), intent(inout) :: self !! RMVS test particle object + class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object + class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters ! Internals integer(I4B) :: i - associate(tp => self, ntp => self%nbody) + associate(tp => self, ntp => self%nbody, pl => system%pl, t => param%t) do i = 1, ntp associate(iplperP => tp%plperP(i)) if ((tp%status(i) == ACTIVE) .and. (tp%lperi(i))) then @@ -30,7 +29,7 @@ module subroutine rmvs_discard_pl_tp(self, pl, t, dt) end if end associate end do - call discard_pl_tp(tp, pl, t, dt) + call discard_pl_tp(tp, system, param) end associate end subroutine rmvs_discard_pl_tp diff --git a/src/rmvs/rmvs_encounter_check.f90 b/src/rmvs/rmvs_encounter_check.f90 index 3d5543166..9719567a9 100644 --- a/src/rmvs/rmvs_encounter_check.f90 +++ b/src/rmvs/rmvs_encounter_check.f90 @@ -1,7 +1,7 @@ submodule (rmvs_classes) s_rmvs_chk use swiftest contains - module function rmvs_encounter_check_tp(self, cb, pl, dt, rts) result(lencounter) + module function rmvs_encounter_check_tp(self, system, dt) result(lencounter) !! author: David A. Minton !! !! Determine whether a test particle and planet are having or will have an encounter within the next time step @@ -10,38 +10,39 @@ module function rmvs_encounter_check_tp(self, cb, pl, dt, rts) result(lencounter !! Adapted from Hal Levison's Swift routine rmvs3_chk.f implicit none ! Arguments - class(rmvs_tp), intent(inout) :: self !! RMVS test particle object - class(rmvs_cb), intent(inout) :: cb !! RMVS central body object - class(rmvs_pl), intent(inout) :: pl !! RMVS massive body object - real(DP), intent(in) :: dt !! step size - real(DP), intent(in) :: rts !! fraction of Hill's sphere radius to use as radius of encounter regio - logical :: lencounter !! Returns true if there is at least one close encounter + class(rmvs_tp), intent(inout) :: self !! RMVS test particle object + class(rmvs_nbody_system), intent(inout) :: system !! RMVS nbody system object + real(DP), intent(in) :: dt !! step size + ! Result + logical :: lencounter !! Returns true if there is at least one close encounter ! Internals - integer(I4B) :: i, j - real(DP) :: r2, v2, vdotr - real(DP), dimension(NDIM) :: xr, vr - real(DP), dimension(pl%nbody) :: r2crit - logical :: lflag + integer(I4B) :: i, j + real(DP) :: r2, v2, vdotr + real(DP), dimension(NDIM) :: xr, vr + real(DP), dimension(system%pl%nbody) :: r2crit + logical :: lflag - associate(tp => self, ntp => self%nbody, npl => pl%nbody, rhill => pl%rhill, xht => self%xh, vht => self%vh, & - xbeg => self%xbeg, vbeg => self%vbeg, status => self%status, plencP => self%plencP, nenc => pl%nenc) - r2crit(:) = (rts * rhill(:))**2 - plencP(:) = 0 - do j = 1, npl - do i = 1, ntp - if ((status(i) /= ACTIVE).or.(plencP(i) /= 0)) cycle - xr(:) = xht(:, i) - xbeg(:, j) - vr(:) = vht(:, i) - vbeg(:, j) - r2 = dot_product(xr(:), xr(:)) - v2 = dot_product(vr(:), vr(:)) - vdotr = dot_product(vr(:), xr(:)) - lflag = rmvs_chk_ind(r2, v2, vdotr, dt, r2crit(j)) - if (lflag) plencP(i) = j + select type(pl => system%pl) + class is (rmvs_pl) + associate(tp => self, ntp => self%nbody, npl => pl%nbody, xbeg => system%xbeg, vbeg => system%vbeg, rts => system%rts) + r2crit(:) = (rts * pl%rhill(:))**2 + tp%plencP(:) = 0 + do j = 1, npl + do i = 1, ntp + if ((tp%status(i) /= ACTIVE).or.(tp%plencP(i) /= 0)) cycle + xr(:) = tp%xh(:, i) - xbeg(:, j) + vr(:) = tp%vh(:, i) - vbeg(:, j) + r2 = dot_product(xr(:), xr(:)) + v2 = dot_product(vr(:), vr(:)) + vdotr = dot_product(vr(:), xr(:)) + lflag = rmvs_chk_ind(r2, v2, vdotr, dt, r2crit(j)) + if (lflag) tp%plencP(i) = j + end do + pl%nenc(j) = count(tp%plencP(:) == j) end do - nenc(j) = count(plencP(:) == j) - end do - lencounter = any(nenc(:) > 0) - end associate + lencounter = any(pl%nenc(:) > 0) + end associate + end select return end function rmvs_encounter_check_tp diff --git a/src/rmvs/rmvs_getacch.f90 b/src/rmvs/rmvs_getacch.f90 index 15fbc3b7a..8329b580a 100644 --- a/src/rmvs/rmvs_getacch.f90 +++ b/src/rmvs/rmvs_getacch.f90 @@ -1,7 +1,8 @@ submodule(rmvs_classes) s_rmvs_getacch use swiftest contains - module subroutine rmvs_getacch_tp(self, cb, pl, param, t, xh) + module subroutine rmvs_getacch_tp(self, system, param, t, xhp) + !! author: David A. Minton !! !! Compute the oblateness acceleration in the inner encounter region with planets @@ -10,69 +11,69 @@ module subroutine rmvs_getacch_tp(self, cb, pl, param, t, xh) !! uses object polymorphism, and so is not directly adapted. implicit none ! Arguments - class(rmvs_tp), intent(inout) :: self !! RMVS test particle data structure - class(swiftest_cb), intent(inout) :: cb !! Swiftest central body particle data structuree - class(whm_pl), intent(inout) :: pl !! WHM massive body particle data structure. - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters of parameter - real(DP), intent(in) :: t !! Current time - real(DP), dimension(:,:), intent(in) :: xh !! Heliocentric positions of planets + class(rmvs_tp), intent(inout) :: self !! RMVS test particle data structure + class(swiftest_nbody_system), intent(inout) :: system !! Swiftest central body particle data structuree + 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 at current substep ! Internals type(swiftest_parameters) :: param_planetocen real(DP), dimension(:, :), allocatable :: xh_original integer(I4B) :: i - associate(tp => self, ntp => self%nbody, ipleP => self%ipleP, & - inner_index => self%index, cb_heliocentric => self%cb_heliocentric) - - if (tp%lplanetocentric) then ! This is a close encounter step, so any accelerations requiring heliocentric position values - ! must be handeled outside the normal WHM method call - select type(pl) - class is (rmvs_pl) - select type (cb) - class is (rmvs_cb) - associate(xpc => pl%xh, xpct => self%xh) - allocate(xh_original, source=tp%xh) - param_planetocen = param - ! Temporarily turn off the heliocentric-dependent acceleration terms during an inner encounter - param_planetocen%loblatecb = .false. - param_planetocen%lextra_force = .false. - param_planetocen%lgr = .false. - ! Now compute the planetocentric values of acceleration - call whm_getacch_tp(tp, cb, pl, param_planetocen, t, xh) + associate(tp => self, ntp => self%nbody, ipleP => self%ipleP, inner_index => self%index, cb_heliocentric => self%cb_heliocentric) + select type(system) + class is (rmvs_nbody_system) + if (system%lplanetocentric) then ! This is a close encounter step, so any accelerations requiring heliocentric position values + ! must be handeled outside the normal WHM method call + select type(pl => system%pl) + class is (rmvs_pl) + select type (cb => system%cb) + class is (rmvs_cb) + associate(xpc => pl%xh, xpct => self%xh) + allocate(xh_original, source=tp%xh) + param_planetocen = param + ! Temporarily turn off the heliocentric-dependent acceleration terms during an inner encounter + param_planetocen%loblatecb = .false. + param_planetocen%lextra_force = .false. + param_planetocen%lgr = .false. + ! Now compute the planetocentric values of acceleration + call whm_getacch_tp(tp, system, param, t, xhp) - ! Now compute any heliocentric values of acceleration - if (tp%lfirst) then - do i = 1, ntp - tp%xheliocentric(:,i) = tp%xh(:,i) + cb%inner(inner_index - 1)%x(:,1) - end do - else - do i = 1, ntp - tp%xheliocentric(:,i) = tp%xh(:,i) + cb%inner(inner_index )%x(:,1) - end do - end if - ! Swap the planetocentric and heliocentric position vectors - tp%xh(:,:) = tp%xheliocentric(:,:) - if (param%loblatecb) then - ! Put in the current encountering planet's oblateness acceleration as the central body's + ! Now compute any heliocentric values of acceleration if (tp%lfirst) then - cb_heliocentric%aobl(:) = cb%inner(inner_index - 1)%aobl(:,1) + do i = 1, ntp + tp%xheliocentric(:,i) = tp%xh(:,i) + cb%inner(inner_index - 1)%x(:,1) + end do else - cb_heliocentric%aobl(:) = cb%inner(inner_index )%aobl(:,1) + do i = 1, ntp + tp%xheliocentric(:,i) = tp%xh(:,i) + cb%inner(inner_index )%x(:,1) + end do + end if + ! Swap the planetocentric and heliocentric position vectors + tp%xh(:,:) = tp%xheliocentric(:,:) + if (param%loblatecb) then + ! Put in the current encountering planet's oblateness acceleration as the central body's + if (tp%lfirst) then + cb_heliocentric%aobl(:) = cb%inner(inner_index - 1)%aobl(:,1) + else + cb_heliocentric%aobl(:) = cb%inner(inner_index )%aobl(:,1) + end if + call tp%obl_acc(cb_heliocentric) end if - call tp%obl_acc(cb_heliocentric) - end if - if (param%lextra_force) call tp%user_getacch(cb_heliocentric, param, t) - if (param%lgr) call tp%gr_getacch(cb_heliocentric, param) - - tp%xh(:,:) = xh_original(:,:) - end associate - end select - end select - else ! Not a close encounter, so just proceded with the standard WHM method - call whm_getacch_tp(tp, cb, pl, param, t, xh) - end if + if (param%lextra_force) call tp%user_getacch(system, param, t) + if (param%lgr) call tp%gr_getacch(param) + + tp%xh(:,:) = xh_original(:,:) + end associate + end select + end select + else ! Not a close encounter, so just proceded with the standard WHM method + call whm_getacch_tp(tp, system, param, t, xhp) + end if + end select end associate diff --git a/src/rmvs/rmvs_setup.f90 b/src/rmvs/rmvs_setup.f90 index b3bcf8f28..2835aea6c 100644 --- a/src/rmvs/rmvs_setup.f90 +++ b/src/rmvs/rmvs_setup.f90 @@ -9,7 +9,7 @@ module subroutine rmvs_setup_pl(self,n) !! Equivalent in functionality to David E. Kaufmann's Swifter routine rmvs_setup.f90 implicit none ! Arguments - class(rmvs_base_pl), intent(inout) :: self !! RMVS test particle object + class(rmvs_pl), intent(inout) :: self !! RMVS test particle object integer(I4B), intent(in) :: n !! Number of test particles to allocate ! Internals integer(I4B) :: i,j @@ -97,46 +97,49 @@ module subroutine rmvs_setup_system(self, param) ! generated as necessary during close encounter steps. select type(pl => self%pl) class is(rmvs_pl) - select type(cb => self%cb) - class is (rmvs_cb) - select type (tp => self%tp) - class is (rmvs_tp) - tp%cb_heliocentric = cb - pl%lplanetocentric = .false. - tp%lplanetocentric = .false. - cb%lplanetocentric = .false. - associate(npl => pl%nbody) - allocate(pl%planetocentric(npl)) - do i = 1, npl - allocate(pl%planetocentric(i)%cb, source=cb) - allocate(rmvs_pl :: pl%planetocentric(i)%pl) - associate(plenci => pl%planetocentric(i)%pl, cbenci => pl%planetocentric(i)%cb) - cbenci%lplanetocentric = .true. - plenci%lplanetocentric = .true. - call plenci%setup(npl) - plenci%status(:) = ACTIVE - - ! plind stores the heliocentric index value of a planetocentric planet - ! e.g. Consider an encounter with planet 3. - ! Then the following will be the values of plind: - ! pl%planetocentric(3)%pl%plind(1) = 0 (central body - never used) - ! pl%planetocentric(3)%pl%plind(2) = 1 - ! pl%planetocentric(3)%pl%plind(3) = 2 - ! pl%planetocentric(3)%pl%plind(4) = 4 - ! pl%planetocentric(3)%pl%plind(5) = 5 - ! etc. - allocate(plenci%plind(npl)) - plenci%plind(1:npl) = [(j,j=1,npl)] - plenci%plind(2:npl) = pack(plenci%plind(1:npl), plenci%plind(1:npl) /= i) - plenci%plind(1) = 0 - plenci%Gmass(1) = cb%Gmass - plenci%Gmass(2:npl) = pl%Gmass(plenci%plind(2:npl)) - cbenci%Gmass = pl%Gmass(i) + select type(cb => self%cb) + class is (rmvs_cb) + select type (tp => self%tp) + class is (rmvs_tp) + tp%cb_heliocentric = cb + pl%lplanetocentric = .false. + tp%lplanetocentric = .false. + cb%lplanetocentric = .false. + associate(npl => pl%nbody) + allocate(pl%planetocentric(npl)) + do i = 1, npl + pl%planetocentric(i)%lplanetocentric = .true. + allocate(pl%planetocentric(i)%cb, source=cb) + allocate(rmvs_cb :: pl%planetocentric(i)%cb) + allocate(rmvs_pl :: pl%planetocentric(i)%pl) + select type(cbenci => pl%planetocentric(i)%cb) + class is (rmvs_cb) + select type(plenci => pl%planetocentric(i)%pl) + class is (rmvs_pl) + call plenci%setup(npl) + plenci%status(:) = ACTIVE + ! plind stores the heliocentric index value of a planetocentric planet + ! e.g. Consider an encounter with planet 3. + ! Then the following will be the values of plind: + ! pl%planetocentric(3)%pl%plind(1) = 0 (central body - never used) + ! pl%planetocentric(3)%pl%plind(2) = 1 + ! pl%planetocentric(3)%pl%plind(3) = 2 + ! pl%planetocentric(3)%pl%plind(4) = 4 + ! pl%planetocentric(3)%pl%plind(5) = 5 + ! etc. + allocate(plenci%plind(npl)) + plenci%plind(1:npl) = [(j,j=1,npl)] + plenci%plind(2:npl) = pack(plenci%plind(1:npl), plenci%plind(1:npl) /= i) + plenci%plind(1) = 0 + plenci%Gmass(1) = cb%Gmass + plenci%Gmass(2:npl) = pl%Gmass(plenci%plind(2:npl)) + cbenci%Gmass = pl%Gmass(i) + end select + end select + end do end associate - end do - end associate - end select - end select + end select + end select end select end subroutine rmvs_setup_system @@ -147,7 +150,7 @@ module subroutine rmvs_setup_set_beg_end(self, xbeg, xend, vbeg) !! Sets one or more of the values of xbeg, xend, and vbeg implicit none ! Arguments - class(rmvs_tp), intent(inout) :: self !! RMVS test particle object + class(rmvs_nbody_system), intent(inout) :: self !! RMVS test particle object real(DP), dimension(:,:), intent(in), optional :: xbeg, xend, vbeg if (present(xbeg)) then diff --git a/src/rmvs/rmvs_spill_and_fill.f90 b/src/rmvs/rmvs_spill_and_fill.f90 index 4eb1e81dc..ae0ff563b 100644 --- a/src/rmvs/rmvs_spill_and_fill.f90 +++ b/src/rmvs/rmvs_spill_and_fill.f90 @@ -9,9 +9,9 @@ module subroutine rmvs_spill_pl(self, discards, lspill_list) !! Adapted from David E. Kaufmann's Swifter routine discard_discard_spill.f90 implicit none ! Arguments - class(rmvs_base_pl), intent(inout) :: self !! Swiftest massive body body object - class(swiftest_body), intent(inout) :: discards !! Discarded object - logical, dimension(:), intent(in) :: lspill_list !! Logical array of bodies to spill into the discards + class(rmvs_pl), intent(inout) :: self !! Swiftest massive body body object + class(swiftest_body), intent(inout) :: discards !! Discarded object + logical, dimension(:), intent(in) :: lspill_list !! Logical array of bodies to spill into the discards ! Internals integer(I4B) :: i @@ -41,7 +41,7 @@ module subroutine rmvs_fill_pl(self, inserts, lfill_list) !! implicit none ! Arguments - class(rmvs_base_pl), intent(inout) :: self !! RMVS massive body object + class(rmvs_pl), intent(inout) :: self !! RMVS massive body object class(swiftest_body), intent(inout) :: inserts !! Inserted object logical, dimension(:), intent(in) :: lfill_list !! Logical array of bodies to merge into the keeps ! Internals diff --git a/src/rmvs/rmvs_step.f90 b/src/rmvs/rmvs_step.f90 index 7995efdd8..3336d3e98 100644 --- a/src/rmvs/rmvs_step.f90 +++ b/src/rmvs/rmvs_step.f90 @@ -1,7 +1,7 @@ submodule(rmvs_classes) s_rmvs_step use swiftest contains - module subroutine rmvs_step_system(self, param, dt) + module subroutine rmvs_step_system(self, param, t, dt) !! author: David A. Minton !! !! Step massive bodies and and active test particles ahead in heliocentric coordinates @@ -12,60 +12,58 @@ module subroutine rmvs_step_system(self, param, dt) ! Arguments class(rmvs_nbody_system), intent(inout) :: self !! RMVS nbody system object class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters - integer(I4B), intent(in) :: dt !! Current stepsize + real(DP), intent(in) :: t !! Current simulation time + real(DP), intent(in) :: dt !! Current stepsiz ! Internals logical :: lencounter, lfirstpl, lfirsttp real(DP) :: rts real(DP), dimension(:,:), allocatable :: xbeg, xend, vbeg integer(I4B) :: i - select type(system => self) - class is (rmvs_nbody_system) select type(cb => self%cb) class is (rmvs_cb) - select type(pl => self%pl) - class is (rmvs_pl) - select type(tp => self%tp) - class is (rmvs_tp) - associate(ntp => tp%nbody, npl => pl%nbody, t => param%t) - allocate(xbeg, source=pl%xh) - allocate(vbeg, source=pl%vh) - call pl%set_rhill(cb) - call tp%set_beg_end(xbeg = xbeg, vbeg = vbeg) - ! ****** Check for close encounters ***** ! - rts = RHSCALE - lencounter = tp%encounter_check(cb, pl, dt, rts) - if (lencounter) then - lfirstpl = pl%lfirst - lfirsttp = tp%lfirst - pl%outer(0)%x(:,:) = xbeg(:,:) - pl%outer(0)%v(:,:) = vbeg(:,:) - call pl%step(system, param, dt) - pl%outer(NTENC)%x(:,:) = pl%xh(:,:) - pl%outer(NTENC)%v(:,:) = pl%vh(:,:) - call tp%set_beg_end(xend = pl%xh) - call rmvs_interp_out(system, param, dt) - call rmvs_step_out(system, param, dt) - call tp%reverse_status() - call tp%set_beg_end(xbeg = xbeg, xend = xend) - tp%lfirst = .true. - call tp%step(system, param, dt) - where (tp%status(:) == INACTIVE) tp%status(:) = ACTIVE - pl%lfirst = lfirstpl - tp%lfirst = lfirsttp - else - call whm_step_system(self, param) - end if - end associate - end select - end select - end select + select type(pl => self%pl) + class is (rmvs_pl) + select type(tp => self%tp) + class is (rmvs_tp) + associate(system => self, ntp => tp%nbody, npl => pl%nbody) + allocate(xbeg, source=pl%xh) + allocate(vbeg, source=pl%vh) + call pl%set_rhill(cb) + call system%set_beg_end(xbeg = xbeg, vbeg = vbeg) + ! ****** Check for close encounters ***** ! + system%rts = RHSCALE + lencounter = tp%encounter_check(system, dt) + if (lencounter) then + lfirstpl = pl%lfirst + lfirsttp = tp%lfirst + pl%outer(0)%x(:,:) = xbeg(:,:) + pl%outer(0)%v(:,:) = vbeg(:,:) + call pl%step(system, param, t, dt) + pl%outer(NTENC)%x(:,:) = pl%xh(:,:) + pl%outer(NTENC)%v(:,:) = pl%vh(:,:) + call system%set_beg_end(xend = pl%xh) + call rmvs_interp_out(system, param, dt) + call rmvs_step_out(system, param, t, dt) + call tp%reverse_status() + call system%set_beg_end(xbeg = xbeg, xend = xend) + tp%lfirst = .true. + call tp%step(system, param, t, dt) + where (tp%status(:) == INACTIVE) tp%status(:) = ACTIVE + pl%lfirst = lfirstpl + tp%lfirst = lfirsttp + else + call whm_step_system(self, param, t, dt) + end if + end associate + end select + end select end select return end subroutine rmvs_step_system - subroutine rmvs_step_out(system, param, dt) + subroutine rmvs_step_out(system, param, t, dt) !! author: David A. Minton !! !! Step ACTIVE test particles ahead in the outer encounter region, setting up and calling the inner region @@ -77,49 +75,56 @@ subroutine rmvs_step_out(system, param, dt) ! Arguments class(rmvs_nbody_system), intent(inout) :: system !! Swiftest system object class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters - integer(I4B), intent(in) :: dt !! Current stepsize + real(DP), intent(in) :: t !! Current simulation time + real(DP), intent(in) :: dt !! Current stepsiz ! Internals integer(I4B) :: outer_index, j, k real(DP) :: dto, outer_time, rts logical :: lencounter, lfirsttp - associate(cb => system%cb, pl => system%pl, npl => system%pl%nbody, tp => system%tp, ntp => system%tp%nbody, t => param%t) - dto = dt / NTENC - where(tp%plencP(:) == 0) - tp%status(:) = INACTIVE - elsewhere - tp%lperi(:) = .false. - end where - do outer_index = 1, NTENC - outer_time = t + (outer_index - 1) * dto - call tp%set_beg_end(xbeg = pl%outer(outer_index - 1)%x(:, :), & - vbeg = pl%outer(outer_index - 1)%v(:, :), & - xend = pl%outer(outer_index )%x(:, :)) - rts = RHPSCALE - lencounter = tp%encounter_check(cb, pl, dt, rts) - if (lencounter) then - ! Interpolate planets in inner encounter region - call rmvs_interp_in(system, param, dto, outer_index) - ! Step through the inner region - call rmvs_step_in(system, param, outer_time, dto) - lfirsttp = tp%lfirst - tp%lfirst = .true. - call tp%step(cb, pl, param, outer_time, dto) - tp%lfirst = lfirsttp - else - call tp%step(cb, pl, param, outer_time, dto) - end if - do j = 1, npl - if (pl%nenc(j) == 0) cycle - where((tp%plencP(:) == j) .and. (tp%status(:) == INACTIVE)) tp%status(:) = ACTIVE - end do - end do - end associate + select type(pl => system%pl) + class is (rmvs_pl) + select type(tp => system%tp) + class is (rmvs_tp) + associate(cb => system%cb, npl => pl%nbody, ntp => tp%nbody) + dto = dt / NTENC + where(tp%plencP(:) == 0) + tp%status(:) = INACTIVE + elsewhere + tp%lperi(:) = .false. + end where + do outer_index = 1, NTENC + outer_time = t + (outer_index - 1) * dto + call system%set_beg_end(xbeg = pl%outer(outer_index - 1)%x(:, :), & + vbeg = pl%outer(outer_index - 1)%v(:, :), & + xend = pl%outer(outer_index )%x(:, :)) + system%rts = RHPSCALE + lencounter = tp%encounter_check(system, dt) + if (lencounter) then + ! Interpolate planets in inner encounter region + call rmvs_interp_in(system, param, dto, outer_index) + ! Step through the inner region + call rmvs_step_in(system, param, outer_time, dto) + lfirsttp = tp%lfirst + tp%lfirst = .true. + call tp%step(system, param, outer_time, dto) + tp%lfirst = lfirsttp + else + call tp%step(system, param, outer_time, dto) + end if + do j = 1, npl + if (pl%nenc(j) == 0) cycle + where((tp%plencP(:) == j) .and. (tp%status(:) == INACTIVE)) tp%status(:) = ACTIVE + end do + end do + end associate + end select + end select return end subroutine rmvs_step_out - subroutine rmvs_step_in(pl, cb, tp, param, outer_time, dto) + subroutine rmvs_step_in(system, param, outer_time, dto) !! author: David A. Minton !! !! Step active test particles ahead in the inner encounter region @@ -128,57 +133,63 @@ subroutine rmvs_step_in(pl, cb, tp, param, outer_time, dto) !! Adapted from David E. Kaufmann's Swifter routine rmvs_step_in.f90 implicit none ! Arguments - class(rmvs_pl), intent(inout) :: pl !! RMVS massive body object - class(rmvs_cb), intent(inout) :: cb !! RMVS central body object - class(rmvs_tp), intent(inout) :: tp !! RMVS test particle object - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters - real(DP), intent(in) :: outer_time !! Current time - real(DP), intent(in) :: dto !! Step size + class(rmvs_nbody_system), intent(inout) :: system !! RMVS nbody system object + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters + real(DP), intent(in) :: outer_time !! Current time + real(DP), intent(in) :: dto !! Outer step size ! Internals - logical :: lfirsttp - integer(I4B) :: i, j, ipleP - real(DP) :: dti, inner_time - real(DP), dimension(:, :), allocatable :: xbeg, xend, vbeg + logical :: lfirsttp + integer(I4B) :: i, j, ipleP + real(DP) :: dti, inner_time - associate(npl => pl%nbody, nenc => pl%nenc) - dti = dto / NTPHENC - allocate(xbeg, mold=pl%xh) - allocate(xend, mold=pl%xh) - allocate(vbeg, mold=pl%vh) - if (param%loblatecb) call pl%obl_acc(cb) - call rmvs_make_planetocentric(pl, cb, tp, param) - do i = 1, npl - if (nenc(i) == 0) cycle - associate(cbenci => pl%planetocentric(i)%cb, plenci => pl%planetocentric(i)%pl, & - tpenci => pl%planetocentric(i)%tp) - associate(inner_index => tpenci%index) - ! There are inner encounters with this planet...switch to planetocentric coordinates to proceed - tpenci%lfirst = .true. - inner_time = outer_time - call rmvs_peri_tp(tpenci, pl, inner_time, dti, .true., 0, i, param) - ! now step the encountering test particles fully through the inner encounter - lfirsttp = .true. - do inner_index = 1, NTPHENC ! Integrate over the encounter region, using the "substitute" planetocentric systems at each level - plenci%xh(:,:) = plenci%inner(inner_index - 1)%x(:,:) - call tpenci%set_beg_end(xbeg = plenci%inner(inner_index - 1)%x, & - xend = plenci%inner(inner_index)%x) - call tpenci%step(cbenci, plenci, param, inner_time, dti) - do j = 1, nenc(i) - tpenci%xheliocentric(:, j) = tpenci%xh(:, j) + pl%inner(inner_index)%x(:,i) - end do - inner_time = outer_time + j * dti - call rmvs_peri_tp(tpenci, pl, inner_time, dti, .false., inner_index, i, param) - end do - where(tpenci%status(:) == ACTIVE) tpenci%status(:) = INACTIVE - end associate + select type(pl => system%pl) + class is (rmvs_pl) + select type (tp => system%tp) + class is (rmvs_tp) + associate(npl => pl%nbody, cb => system%cb) + dti = dto / NTPHENC + if (param%loblatecb) call pl%obl_acc(cb) + call rmvs_make_planetocentric(system, param) + do i = 1, npl + if (pl%nenc(i) == 0) cycle + select type(planetocen_system => pl%planetocentric(i)) + class is (rmvs_nbody_system) + select type(plenci => planetocen_system%pl) + class is (rmvs_pl) + select type(tpenci => planetocen_system%tp) + class is (rmvs_tp) + associate(inner_index => tpenci%index) + ! There are inner encounters with this planet...switch to planetocentric coordinates to proceed + tpenci%lfirst = .true. + inner_time = outer_time + call rmvs_peri_tp(tpenci, pl, inner_time, dti, .true., 0, i, param) + ! now step the encountering test particles fully through the inner encounter + lfirsttp = .true. + do inner_index = 1, NTPHENC ! Integrate over the encounter region, using the "substitute" planetocentric systems at each level + plenci%xh(:,:) = plenci%inner(inner_index - 1)%x(:,:) + call planetocen_system%set_beg_end(xbeg = plenci%inner(inner_index - 1)%x, & + xend = plenci%inner(inner_index)%x) + call tpenci%step(planetocen_system, param, inner_time, dti) + do j = 1, pl%nenc(i) + tpenci%xheliocentric(:, j) = tpenci%xh(:, j) + pl%inner(inner_index)%x(:,i) + end do + inner_time = outer_time + j * dti + call rmvs_peri_tp(tpenci, pl, inner_time, dti, .false., inner_index, i, param) + end do + where(tpenci%status(:) == ACTIVE) tpenci%status(:) = INACTIVE + end associate + end select + end select + end select + end do + call rmvs_end_planetocentric(system) end associate - end do - call rmvs_end_planetocentric(pl, cb, tp) - end associate + end select + end select return end subroutine rmvs_step_in - subroutine rmvs_interp_in(pl, cb, dt, outer_index, param) + subroutine rmvs_interp_in(system, param, dt, outer_index) !! author: David A. Minton !! !! Interpolate planet positions between two Keplerian orbits in inner encounter regio @@ -188,107 +199,109 @@ subroutine rmvs_interp_in(pl, cb, dt, outer_index, param) !! Adapted from Hal Levison's Swift routine rmvs3_interp.f implicit none ! Arguments - class(rmvs_pl), intent(inout) :: pl !! RMVS test particle object - class(rmvs_cb), intent(inout) :: cb !! RMVS central body particle type - real(DP), intent(in) :: dt !! Step size - integer(I4B), intent(in) :: outer_index !! Outer substep number within current se - class(swiftest_parameters), intent(in) :: param !! Swiftest parameters file + class(rmvs_nbody_system), intent(inout) :: system !! RMVS test particle object + class(swiftest_parameters), intent(in) :: param !! Swiftest parameters file + real(DP), intent(in) :: dt !! Step size + integer(I4B), intent(in) :: outer_index !! Outer substep number within current set ! Internals - integer(I4B) :: i, inner_index - real(DP) :: frac, dntphenc - real(DP), dimension(:,:), allocatable :: xtmp, vtmp, xh_original - real(DP), dimension(:), allocatable :: msun, dti - integer(I4B), dimension(:), allocatable :: iflag - - associate (npl => pl%nbody) - dntphenc = real(NTPHENC, kind=DP) + integer(I4B) :: i, inner_index + real(DP) :: frac, dntphenc + real(DP), dimension(:,:), allocatable :: xtmp, vtmp, xh_original + real(DP), dimension(:), allocatable :: GMcb, dti + integer(I4B), dimension(:), allocatable :: iflag - ! Set the endpoints of the inner region from the outer region values in the current outer step index - pl%inner(0)%x(:,:) = pl%outer(outer_index - 1)%x(:, :) - pl%inner(0)%v(:,:) = pl%outer(outer_index - 1)%v(:, :) - pl%inner(NTPHENC)%x(:,:) = pl%outer(outer_index)%x(:, :) - pl%inner(NTPHENC)%v(:,:) = pl%outer(outer_index)%v(:, :) - - allocate(xtmp,mold=pl%xh) - allocate(vtmp,mold=pl%vh) - allocate(msun(npl)) - allocate(dti(npl)) - allocate(iflag(npl)) - dti(:) = dt / dntphenc - msun(:) = cb%Gmass - xtmp(:, :) = pl%inner(0)%x(:, :) - vtmp(:, :) = pl%inner(0)%v(:, :) - if (param%loblatecb) then - allocate(xh_original,source=pl%xh) - pl%xh(:, :) = xtmp(:, :) ! Temporarily replace heliocentric position with inner substep values to calculate the oblateness terms - call pl%obl_acc(cb) - pl%inner(0)%aobl(:, :) = pl%aobl(:, :) ! Save the oblateness acceleration on the planet for this substep - end if - - do inner_index = 1, NTPHENC - 1 - call drift_one(msun(1:npl), xtmp(1,1:npl), xtmp(2,1:npl), xtmp(3,1:npl), & - vtmp(1,1:npl), vtmp(2,1:npl), vtmp(3,1:npl), & - dti(1:npl), iflag(1:npl)) - if (any(iflag(1:npl) /= 0)) then - do i = 1, npl - if (iflag(i) /=0) then - write(*, *) " Planet ", pl%name(i), " is lost!!!!!!!!!!" - write(*, *) msun(i), dti(i) - write(*, *) xtmp(:,i) - write(*, *) vtmp(:,i) - write(*, *) " STOPPING " - call util_exit(failure) - end if - end do - end if - frac = 1.0_DP - inner_index / dntphenc - pl%inner(inner_index)%x(:, :) = frac * xtmp(:,:) - pl%inner(inner_index)%v(:, :) = frac * vtmp(:,:) - end do + associate (cb => system%cb, npl => system%pl%nbody) + select type(pl => system%pl) + class is (rmvs_pl) + dntphenc = real(NTPHENC, kind=DP) - xtmp(:,:) = pl%inner(NTPHENC)%x(:, :) - vtmp(:,:) = pl%inner(NTPHENC)%v(:, :) - - do inner_index = NTPHENC - 1, 1, -1 - call drift_one(msun(1:npl), xtmp(1,1:npl), xtmp(2,1:npl), xtmp(3,1:npl), & - vtmp(1,1:npl), vtmp(2,1:npl), vtmp(3,1:npl), & - -dti(1:npl), iflag(1:npl)) - if (any(iflag(1:npl) /= 0)) then - do i = 1, npl - if (iflag(i) /=0) then - write(*, *) " Planet ", pl%name(i), " is lost!!!!!!!!!!" - write(*, *) msun(i), -dti(i) - write(*, *) xtmp(:,i) - write(*, *) vtmp(:,i) - write(*, *) " STOPPING " - call util_exit(failure) - end if - end do + ! Set the endpoints of the inner region from the outer region values in the current outer step index + pl%inner(0)%x(:,:) = pl%outer(outer_index - 1)%x(:, :) + pl%inner(0)%v(:,:) = pl%outer(outer_index - 1)%v(:, :) + pl%inner(NTPHENC)%x(:,:) = pl%outer(outer_index)%x(:, :) + pl%inner(NTPHENC)%v(:,:) = pl%outer(outer_index)%v(:, :) + + allocate(xtmp,mold=pl%xh) + allocate(vtmp,mold=pl%vh) + allocate(GMcb(npl)) + allocate(dti(npl)) + allocate(iflag(npl)) + dti(:) = dt / dntphenc + GMcb(:) = cb%Gmass + xtmp(:, :) = pl%inner(0)%x(:, :) + vtmp(:, :) = pl%inner(0)%v(:, :) + if (param%loblatecb) then + allocate(xh_original,source=pl%xh) + pl%xh(:, :) = xtmp(:, :) ! Temporarily replace heliocentric position with inner substep values to calculate the oblateness terms + call pl%obl_acc(cb) + pl%inner(0)%aobl(:, :) = pl%aobl(:, :) ! Save the oblateness acceleration on the planet for this substep end if - frac = inner_index / dntphenc - pl%inner(inner_index)%x(:, :) = pl%inner(inner_index)%x(:, :) + frac * xtmp(:, :) - pl%inner(inner_index)%v(:, :) = pl%inner(inner_index)%v(:, :) + frac * vtmp(:, :) - + + do inner_index = 1, NTPHENC - 1 + call drift_one(GMcb(1:npl), xtmp(1,1:npl), xtmp(2,1:npl), xtmp(3,1:npl), & + vtmp(1,1:npl), vtmp(2,1:npl), vtmp(3,1:npl), & + dti(1:npl), iflag(1:npl)) + if (any(iflag(1:npl) /= 0)) then + do i = 1, npl + if (iflag(i) /=0) then + write(*, *) " Planet ", pl%name(i), " is lost!!!!!!!!!!" + write(*, *) GMcb(i), dti(i) + write(*, *) xtmp(:,i) + write(*, *) vtmp(:,i) + write(*, *) " STOPPING " + call util_exit(failure) + end if + end do + end if + frac = 1.0_DP - inner_index / dntphenc + pl%inner(inner_index)%x(:, :) = frac * xtmp(:,:) + pl%inner(inner_index)%v(:, :) = frac * vtmp(:,:) + end do + + xtmp(:,:) = pl%inner(NTPHENC)%x(:, :) + vtmp(:,:) = pl%inner(NTPHENC)%v(:, :) + + do inner_index = NTPHENC - 1, 1, -1 + call drift_one(GMcb(1:npl), xtmp(1,1:npl), xtmp(2,1:npl), xtmp(3,1:npl), & + vtmp(1,1:npl), vtmp(2,1:npl), vtmp(3,1:npl), & + -dti(1:npl), iflag(1:npl)) + if (any(iflag(1:npl) /= 0)) then + do i = 1, npl + if (iflag(i) /=0) then + write(*, *) " Planet ", pl%name(i), " is lost!!!!!!!!!!" + write(*, *) GMcb(i), -dti(i) + write(*, *) xtmp(:,i) + write(*, *) vtmp(:,i) + write(*, *) " STOPPING " + call util_exit(failure) + end if + end do + end if + frac = inner_index / dntphenc + pl%inner(inner_index)%x(:, :) = pl%inner(inner_index)%x(:, :) + frac * xtmp(:, :) + pl%inner(inner_index)%v(:, :) = pl%inner(inner_index)%v(:, :) + frac * vtmp(:, :) + + if (param%loblatecb) then + pl%xh(:,:) = pl%inner(inner_index)%x(:, :) + call pl%obl_acc(cb) + pl%inner(inner_index)%aobl(:, :) = pl%aobl(:, :) + end if + end do if (param%loblatecb) then - pl%xh(:,:) = pl%inner(inner_index)%x(:, :) + ! Calculate the final value of oblateness accelerations at the final inner substep + pl%xh(:,:) = pl%inner(NTPHENC)%x(:, :) call pl%obl_acc(cb) - pl%inner(inner_index)%aobl(:, :) = pl%aobl(:, :) + pl%inner(NTPHENC)%aobl(:, :) = pl%aobl(:, :) + ! Put the planet positions back into place + call move_alloc(xh_original, pl%xh) end if - end do - if (param%loblatecb) then - ! Calculate the final value of oblateness accelerations at the final inner substep - pl%xh(:,:) = pl%inner(NTPHENC)%x(:, :) - call pl%obl_acc(cb) - pl%inner(NTPHENC)%aobl(:, :) = pl%aobl(:, :) - ! Put the planet positions back into place - call move_alloc(xh_original, pl%xh) - end if + end select end associate return end subroutine rmvs_interp_in - subroutine rmvs_interp_out(pl, cb, dt, param) + subroutine rmvs_interp_out(system, param, dt) !! author: David A. Minton !! !! Interpolate planet positions between two Keplerian orbits in outer encounter region @@ -298,72 +311,73 @@ subroutine rmvs_interp_out(pl, cb, dt, param) !! Adapted from Hal Levison's Swift routine rmvs3_interp.f implicit none ! Arguments - class(rmvs_pl), intent(inout) :: pl !! RMVS test particle object - class(rmvs_cb), intent(inout) :: cb !! RMVS central body particle type - real(DP), intent(in) :: dt !! Step size + class(rmvs_nbody_system), intent(inout) :: system !! RMVS nbody system object class(swiftest_parameters), intent(in) :: param !! Swiftest parameters file + real(DP), intent(in) :: dt !! Step size ! Internals - integer(I4B) :: i, outer_index - real(DP) :: frac, dntenc - real(DP), dimension(:,:), allocatable :: xtmp, vtmp - real(DP), dimension(:), allocatable :: msun, dto - integer(I4B), dimension(:), allocatable :: iflag + integer(I4B) :: i, outer_index + real(DP) :: frac, dntenc + real(DP), dimension(:,:), allocatable :: xtmp, vtmp + real(DP), dimension(:), allocatable :: GMcb, dto + integer(I4B), dimension(:), allocatable :: iflag - ! executable code - dntenc = real(NTENC, DP) - associate (npl => pl%nbody) - allocate(xtmp, mold = pl%xh) - allocate(vtmp, mold = pl%vh) - allocate(msun(npl)) - allocate(dto(npl)) - allocate(iflag(npl)) - dto(:) = dt / dntenc - msun(:) = cb%Gmass - xtmp(:,:) = pl%outer(0)%x(:, :) - vtmp(:,:) = pl%outer(0)%v(:, :) - do outer_index = 1, NTENC - 1 - call drift_one(msun(1:npl), xtmp(1,1:npl), xtmp(2,1:npl), xtmp(3,1:npl), & - vtmp(1,1:npl), vtmp(2,1:npl), vtmp(3,1:npl), & - dto(1:npl), iflag(1:npl)) + dntenc = real(NTENC, kind=DP) + associate (cb => system%cb, pl => system%pl, npl => system%pl%nbody) + select type(pl => system%pl) + class is (rmvs_pl) + allocate(xtmp, mold = pl%xh) + allocate(vtmp, mold = pl%vh) + allocate(GMcb(npl)) + allocate(dto(npl)) + allocate(iflag(npl)) + dto(:) = dt / dntenc + GMcb(:) = cb%Gmass + xtmp(:,:) = pl%outer(0)%x(:, :) + vtmp(:,:) = pl%outer(0)%v(:, :) + do outer_index = 1, NTENC - 1 + call drift_one(GMcb(1:npl), xtmp(1,1:npl), xtmp(2,1:npl), xtmp(3,1:npl), & + vtmp(1,1:npl), vtmp(2,1:npl), vtmp(3,1:npl), & + dto(1:npl), iflag(1:npl)) - if (any(iflag(1:npl) /= 0)) then - do i = 1, npl - if (iflag(i) /= 0) then - write(*, *) " Planet ", pl%name(i), " is lost!!!!!!!!!!" - write(*, *) msun(i), dto(i) - write(*, *) xtmp(:,i) - write(*, *) vtmp(:,i) - write(*, *) " STOPPING " - call util_exit(FAILURE) - end if - end do - end if - frac = 1.0_DP - outer_index / dntenc - pl%outer(outer_index)%x(:, :) = frac * xtmp(:,:) - pl%outer(outer_index)%v(:, :) = frac * vtmp(:,:) - end do - xtmp(:,:) = pl%outer(NTENC)%x(:, :) - vtmp(:,:) = pl%outer(NTENC)%v(:, :) - do outer_index = NTENC - 1, 1, -1 - call drift_one(msun(1:npl), xtmp(1,1:npl), xtmp(2,1:npl), xtmp(3,1:npl), & - vtmp(1,1:npl), vtmp(2,1:npl), vtmp(3,1:npl), & - -dto(1:npl), iflag(1:npl)) - if (any(iflag(1:npl) /= 0)) then - do i = 1, npl - if (iflag(i) /= 0) then - write(*, *) " Planet ", pl%name(i), " is lost!!!!!!!!!!" - write(*, *) msun(i), -dto(i) - write(*, *) xtmp(:,i) - write(*, *) vtmp(:,i) - write(*, *) " STOPPING " - call util_exit(FAILURE) - end if - end do - end if - frac = outer_index / dntenc - pl%outer(outer_index)%x(:, :) = pl%outer(outer_index)%x(:, :) + frac * xtmp(:,:) - pl%outer(outer_index)%v(:, :) = pl%outer(outer_index)%v(:, :) + frac * vtmp(:,:) - end do + if (any(iflag(1:npl) /= 0)) then + do i = 1, npl + if (iflag(i) /= 0) then + write(*, *) " Planet ", pl%name(i), " is lost!!!!!!!!!!" + write(*, *) GMcb(i), dto(i) + write(*, *) xtmp(:,i) + write(*, *) vtmp(:,i) + write(*, *) " STOPPING " + call util_exit(FAILURE) + end if + end do + end if + frac = 1.0_DP - outer_index / dntenc + pl%outer(outer_index)%x(:, :) = frac * xtmp(:,:) + pl%outer(outer_index)%v(:, :) = frac * vtmp(:,:) + end do + xtmp(:,:) = pl%outer(NTENC)%x(:, :) + vtmp(:,:) = pl%outer(NTENC)%v(:, :) + do outer_index = NTENC - 1, 1, -1 + call drift_one(GMcb(1:npl), xtmp(1,1:npl), xtmp(2,1:npl), xtmp(3,1:npl), & + vtmp(1,1:npl), vtmp(2,1:npl), vtmp(3,1:npl), & + -dto(1:npl), iflag(1:npl)) + if (any(iflag(1:npl) /= 0)) then + do i = 1, npl + if (iflag(i) /= 0) then + write(*, *) " Planet ", pl%name(i), " is lost!!!!!!!!!!" + write(*, *) GMcb(i), -dto(i) + write(*, *) xtmp(:,i) + write(*, *) vtmp(:,i) + write(*, *) " STOPPING " + call util_exit(FAILURE) + end if + end do + end if + frac = outer_index / dntenc + pl%outer(outer_index)%x(:, :) = pl%outer(outer_index)%x(:, :) + frac * xtmp(:,:) + pl%outer(outer_index)%v(:, :) = pl%outer(outer_index)%v(:, :) + frac * vtmp(:,:) + end do + end select end associate return @@ -453,7 +467,7 @@ subroutine rmvs_peri_tp(tp, pl, t, dt, lfirst, inner_index, ipleP, param) end subroutine rmvs_peri_tp - subroutine rmvs_make_planetocentric(pl, cb, tp, param) + subroutine rmvs_make_planetocentric(system, param) !! author: David A. Minton !! !! When encounters are detected, this method will call the interpolation methods for the planets and @@ -462,112 +476,135 @@ subroutine rmvs_make_planetocentric(pl, cb, tp, param) !! implicit none ! Arguments - class(rmvs_pl), intent(inout) :: pl !! RMVS test particle object - class(rmvs_cb), intent(inout) :: cb !! RMVS central body particle type - class(rmvs_tp), intent(inout) :: tp !! RMVS test particle object - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters + class(rmvs_nbody_system), intent(inout) :: system !! RMVS nbody system object + class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters + ! Internals integer(I4B) :: i, j, inner_index, ipc2hc logical, dimension(:), allocatable :: encmask - associate(npl => pl%nbody, ntp => tp%nbody, GMpl => pl%Gmass, nenc => pl%nenc) - do i = 1, npl - if (nenc(i) == 0) cycle - ! There are inner encounters with this planet - if (allocated(encmask)) deallocate(encmask) - allocate(encmask(ntp)) - encmask(:) = tp%plencP(:) == i - - ! Create encountering test particle structure - allocate(rmvs_tp :: pl%planetocentric(i)%tp) - associate(cbenci => pl%planetocentric(i)%cb, & - plenci => pl%planetocentric(i)%pl, & - tpenci => pl%planetocentric(i)%tp) - tpenci%lplanetocentric = .true. - call tpenci%setup(nenc(i)) - tpenci%cb_heliocentric = cb - tpenci%ipleP = i - tpenci%status(:) = ACTIVE - tpenci%name(:) = pack(tp%name(:), encmask(:)) - ! Grab all the encountering test particles and convert them to a planetocentric frame - do j = 1, NDIM - tpenci%xheliocentric(j, :) = pack(tp%xh(j,:), encmask(:)) - tpenci%xh(j, :) = tpenci%xheliocentric(j, :) - pl%inner(0)%x(j, i) - tpenci%vh(j, :) = pack(tp%vh(j,:), encmask(:)) - pl%inner(0)%v(j, i) - end do - ! Make sure that the test particles get the planetocentric value of mu - allocate(cbenci%inner(0:NTPHENC)) - do inner_index = 0, NTPHENC - allocate(plenci%inner(inner_index)%x, mold=pl%inner(inner_index)%x) - allocate(plenci%inner(inner_index)%v, mold=pl%inner(inner_index)%x) - allocate(plenci%inner(inner_index)%aobl, mold=pl%inner(inner_index)%aobl) - allocate(cbenci%inner(inner_index)%x(NDIM,1)) - allocate(cbenci%inner(inner_index)%v(NDIM,1)) - allocate(cbenci%inner(inner_index)%aobl(NDIM,1)) - cbenci%inner(inner_index)%x(:,1) = pl%inner(inner_index)%x(:, i) - cbenci%inner(inner_index)%v(:,1) = pl%inner(inner_index)%v(:, i) - cbenci%inner(inner_index)%aobl(:,1) = pl%inner(inner_index)%aobl(:, i) - - plenci%inner(inner_index)%x(:,1) = -cbenci%inner(inner_index)%x(:,1) - plenci%inner(inner_index)%v(:,1) = -cbenci%inner(inner_index)%v(:,1) - do j = 2, npl - ipc2hc = plenci%plind(j) - plenci%inner(inner_index)%x(:,j) = pl%inner(inner_index)%x(:, ipc2hc) - cbenci%inner(inner_index)%x(:,1) - plenci%inner(inner_index)%v(:,j) = pl%inner(inner_index)%v(:, ipc2hc) - cbenci%inner(inner_index)%v(:,1) + select type(cb => system%cb) + class is (rmvs_cb) + select type(pl => system%pl) + class is (rmvs_pl) + select type (tp => system%tp) + class is (rmvs_tp) + associate (npl => pl%nbody, ntp => tp%nbody) + do i = 1, npl + if (pl%nenc(i) == 0) cycle + ! There are inner encounters with this planet + if (allocated(encmask)) deallocate(encmask) + allocate(encmask(ntp)) + encmask(:) = tp%plencP(:) == i + allocate(rmvs_tp :: pl%planetocentric(i)%tp) + ! Create encountering test particle structure + select type(cbenci => pl%planetocentric(i)%cb) + class is (rmvs_cb) + select type(plenci => pl%planetocentric(i)%pl) + class is (rmvs_pl) + select type(tpenci => pl%planetocentric(i)%tp) + class is (rmvs_tp) + tpenci%lplanetocentric = .true. + call tpenci%setup(pl%nenc(i)) + tpenci%cb_heliocentric = cb + tpenci%ipleP = i + tpenci%status(:) = ACTIVE + tpenci%name(:) = pack(tp%name(:), encmask(:)) + ! Grab all the encountering test particles and convert them to a planetocentric frame + do j = 1, NDIM + tpenci%xheliocentric(j, :) = pack(tp%xh(j,:), encmask(:)) + tpenci%xh(j, :) = tpenci%xheliocentric(j, :) - pl%inner(0)%x(j, i) + tpenci%vh(j, :) = pack(tp%vh(j,:), encmask(:)) - pl%inner(0)%v(j, i) + end do + ! Make sure that the test particles get the planetocentric value of mu + allocate(cbenci%inner(0:NTPHENC)) + do inner_index = 0, NTPHENC + allocate(plenci%inner(inner_index)%x, mold=pl%inner(inner_index)%x) + allocate(plenci%inner(inner_index)%v, mold=pl%inner(inner_index)%x) + allocate(plenci%inner(inner_index)%aobl, mold=pl%inner(inner_index)%aobl) + allocate(cbenci%inner(inner_index)%x(NDIM,1)) + allocate(cbenci%inner(inner_index)%v(NDIM,1)) + allocate(cbenci%inner(inner_index)%aobl(NDIM,1)) + cbenci%inner(inner_index)%x(:,1) = pl%inner(inner_index)%x(:, i) + cbenci%inner(inner_index)%v(:,1) = pl%inner(inner_index)%v(:, i) + cbenci%inner(inner_index)%aobl(:,1) = pl%inner(inner_index)%aobl(:, i) + plenci%inner(inner_index)%x(:,1) = -cbenci%inner(inner_index)%x(:,1) + plenci%inner(inner_index)%v(:,1) = -cbenci%inner(inner_index)%v(:,1) + do j = 2, npl + ipc2hc = plenci%plind(j) + plenci%inner(inner_index)%x(:,j) = pl%inner(inner_index)%x(:, ipc2hc) - cbenci%inner(inner_index)%x(:,1) + plenci%inner(inner_index)%v(:,j) = pl%inner(inner_index)%v(:, ipc2hc) - cbenci%inner(inner_index)%v(:,1) + end do + end do + call tpenci%set_mu(cbenci) + end select + end select + end select end do - end do - call tpenci%set_mu(cbenci) - end associate - end do - end associate + end associate + end select + end select + end select return - end subroutine rmvs_make_planetocentric + end subroutine rmvs_make_planetocentric - subroutine rmvs_end_planetocentric(pl, cb, tp) + subroutine rmvs_end_planetocentric(system) !! author: David A. Minton !! !! Deallocates all of the encountering particle data structures for next time !! implicit none ! Arguments - class(rmvs_pl), intent(inout) :: pl !! RMVS test particle object - class(rmvs_cb), intent(inout) :: cb !! RMVS central body object - class(rmvs_tp), intent(inout) :: tp !! RMVS test particle object + class(rmvs_nbody_system), intent(inout) :: system !! RMVS nbody system object ! Internals integer(I4B) :: i, j, inner_index integer(I4B), dimension(:), allocatable :: tpind logical, dimension(:), allocatable :: encmask - associate(nenc => pl%nenc, npl => pl%nbody, ntp => tp%nbody) - do i = 1, npl - if (nenc(i) == 0) cycle - associate(cbenci => pl%planetocentric(i)%cb, & - tpenci => pl%planetocentric(i)%tp, & - plenci => pl%planetocentric(i)%pl) - if (allocated(tpind)) deallocate(tpind) - allocate(tpind(nenc(i))) - ! Index array of encountering test particles - if (allocated(encmask)) deallocate(encmask) - allocate(encmask(ntp)) - encmask(:) = tp%plencP(:) == i - tpind(:) = pack([(j,j=1,ntp)], encmask(:)) - - ! Copy the results of the integration back over and shift back to heliocentric reference - tp%status(tpind(1:nenc(i))) = tpenci%status(1:nenc(i)) - do j = 1, NDIM - tp%xh(j, tpind(1:nenc(i))) = tpenci%xh(j,1:nenc(i)) + pl%inner(NTPHENC)%x(j, i) - tp%vh(j, tpind(1:nenc(i))) = tpenci%vh(j,1:nenc(i)) + pl%inner(NTPHENC)%v(j, i) - end do - deallocate(pl%planetocentric(i)%tp) - deallocate(cbenci%inner) - do inner_index = 0, NTPHENC - deallocate(plenci%inner(inner_index)%x) - deallocate(plenci%inner(inner_index)%v) - deallocate(plenci%inner(inner_index)%aobl) - end do - end associate - end do - end associate + select type(cb => system%cb) + class is (rmvs_cb) + select type(pl => system%pl) + class is (rmvs_pl) + select type (tp => system%tp) + class is (rmvs_tp) + associate (npl => pl%nbody, ntp => tp%nbody) + do i = 1, npl + if (pl%nenc(i) == 0) cycle + select type(cbenci => pl%planetocentric(i)%cb) + class is (rmvs_cb) + select type(plenci => pl%planetocentric(i)%pl) + class is (rmvs_pl) + select type(tpenci => pl%planetocentric(i)%tp) + class is (rmvs_tp) + if (allocated(tpind)) deallocate(tpind) + allocate(tpind(pl%nenc(i))) + ! Index array of encountering test particles + if (allocated(encmask)) deallocate(encmask) + allocate(encmask(ntp)) + encmask(:) = tp%plencP(:) == i + tpind(:) = pack([(j,j=1,ntp)], encmask(:)) + + ! Copy the results of the integration back over and shift back to heliocentric reference + tp%status(tpind(1:pl%nenc(i))) = tpenci%status(1:pl%nenc(i)) + do j = 1, NDIM + tp%xh(j, tpind(1:pl%nenc(i))) = tpenci%xh(j,1:pl%nenc(i)) + pl%inner(NTPHENC)%x(j, i) + tp%vh(j, tpind(1:pl%nenc(i))) = tpenci%vh(j,1:pl%nenc(i)) + pl%inner(NTPHENC)%v(j, i) + end do + deallocate(pl%planetocentric(i)%tp) + deallocate(cbenci%inner) + do inner_index = 0, NTPHENC + deallocate(plenci%inner(inner_index)%x) + deallocate(plenci%inner(inner_index)%v) + deallocate(plenci%inner(inner_index)%aobl) + end do + end select + end select + end select + end do + end associate + end select + end select + end select return end subroutine rmvs_end_planetocentric diff --git a/src/setup/setup.f90 b/src/setup/setup.f90 index 55140ac88..03fd6f0f0 100644 --- a/src/setup/setup.f90 +++ b/src/setup/setup.f90 @@ -9,7 +9,7 @@ module subroutine setup_construct_system(system, param) implicit none ! Arguments class(swiftest_nbody_system), allocatable, intent(inout) :: system !! Swiftest system object - type(swiftest_parameters), intent(in) :: param !! Swiftest parameters + type(swiftest_parameters), intent(in) :: param !! Swiftest parameters select case(param%integrator) case (BS) diff --git a/src/user/user_getacch.f90 b/src/user/user_getacch.f90 index e3733303a..da3d00578 100644 --- a/src/user/user_getacch.f90 +++ b/src/user/user_getacch.f90 @@ -1,7 +1,7 @@ submodule(swiftest_classes) s_user_getacch use swiftest contains - module subroutine user_getacch_body(self, cb, param, t) + module subroutine user_getacch_body(self, system, param, t) !! author: David A. Minton !! !! Add user-supplied heliocentric accelerations to planets @@ -9,10 +9,10 @@ module subroutine user_getacch_body(self, cb, param, t) !! Adapted from David E. Kaufmann's Swifter routine whm_user_getacch.f90 implicit none ! Arguments - class(swiftest_body), intent(inout) :: self !! Swiftest massive body particle data structure - class(swiftest_cb), intent(inout) :: cb !! Swiftest central body particle data structuree - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters of user parameters - real(DP), intent(in) :: t !! Current time + class(swiftest_body), intent(inout) :: self !! Swiftest massive body particle data structure + class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody_system_object + class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters of user parameters + real(DP), intent(in) :: t !! Current time return end subroutine user_getacch_body diff --git a/src/whm/whm_drift.f90 b/src/whm/whm_drift.f90 index c18648364..0ac170991 100644 --- a/src/whm/whm_drift.f90 +++ b/src/whm/whm_drift.f90 @@ -10,21 +10,17 @@ module subroutine whm_drift_pl(self, system, param, dt) !! Adapted from David E. Kaufmann's Swifter routine whm_drift.f90 implicit none ! Arguments - class(whm_pl), intent(inout) :: self !! WHM 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) :: dt !! Stepsize + class(whm_pl), intent(inout) :: self !! WHM massive body particle data structure + 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) :: dt !! Stepsize ! Internals integer(I4B) :: i real(DP) :: energy, vmag2, rmag !! Variables used in GR calculation integer(I4B), dimension(:), allocatable :: iflag real(DP), dimension(:), allocatable :: dtp - associate(npl => self%nbody, & - xj => self%xj, & - vj => self%vj, & - status => self%status,& - mu => self%muj) + associate(pl => self, npl => self%nbody) if (npl == 0) return @@ -34,24 +30,24 @@ module subroutine whm_drift_pl(self, system, param, dt) if (param%lgr) then do i = 1,npl - rmag = norm2(xj(:, i)) - vmag2 = dot_product(vj(:, i), vj(:, i)) - energy = 0.5_DP * vmag2 - mu(i) / rmag + rmag = norm2(pl%xj(:, i)) + vmag2 = dot_product(pl%vj(:, i), pl%vj(:, i)) + energy = 0.5_DP * vmag2 - pl%muj(i) / rmag dtp(i) = dt * (1.0_DP + 3 * param%inv_c2 * energy) end do else dtp(:) = dt end if - call drift_one(mu(1:npl), xj(1,1:npl), xj(2,1:npl), xj(3,1:npl), & - vj(1,1:npl), vj(2,1:npl), vj(3,1:npl), & + call drift_one(pl%muj(1:npl), pl%xj(1,1:npl), pl%xj(2,1:npl), pl%xj(3,1:npl), & + pl%vj(1,1:npl), pl%vj(2,1:npl), pl%vj(3,1:npl), & dtp(1:npl), iflag(1:npl)) if (any(iflag(1:npl) /= 0)) then do i = 1, npl if (iflag(i) /= 0) then write(*, *) " Planet ", self%name(i), " is lost!!!!!!!!!!" - write(*, *) xj(:,i) - write(*, *) vj(:,i) + write(*, *) pl%xj(:,i) + write(*, *) pl%vj(:,i) write(*, *) " stopping " call util_exit(FAILURE) end if @@ -73,44 +69,38 @@ module subroutine whm_drift_tp(self, system, param, dt) !! Adapted from David E. Kaufmann's Swifter routine whm_drift_tp.f90 implicit none ! Arguments - class(whm_tp), intent(inout) :: self !! WHM 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 of - real(DP), intent(in) :: dt !! Stepsize + class(whm_tp), intent(inout) :: self !! WHM test particle data structure + 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) :: dt !! Stepsize ! Internals integer(I4B) :: i real(DP) :: energy, vmag2, rmag !! Variables used in GR calculation integer(I4B), dimension(:), allocatable :: iflag real(DP), dimension(:), allocatable :: dtp - - associate(ntp => self%nbody, & - xh => self%xh, & - vh => self%vh, & - status => self%status,& - mu => self%mu, GMcb => system%cb%Gmass) - + associate(tp => self, ntp => self%nbody) if (ntp == 0) return allocate(iflag(ntp)) iflag(:) = 0 allocate(dtp(ntp)) if (param%lgr) then - do i = 1,ntp - rmag = norm2(xh(:, i)) - vmag2 = dot_product(vh(:, i), vh(:, i)) - energy = 0.5_DP * vmag2 - GMcb / rmag + do i = 1, ntp + rmag = norm2(tp%xh(:, i)) + vmag2 = dot_product(tp%vh(:, i), tp%vh(:, i)) + energy = 0.5_DP * vmag2 - tp%mu(i) / rmag dtp(i) = dt * (1.0_DP + 3 * param%inv_c2 * energy) end do else dtp(:) = dt end if - do concurrent(i = 1:ntp, status(i) == ACTIVE) - call drift_one(mu(i), xh(1,i), xh(2,i), xh(3,i), & - vh(1,i), vh(2,i), vh(3,i), & - dtp(i), iflag(i)) + do concurrent(i = 1:ntp, tp%status(i) == ACTIVE) + call drift_one(tp%mu(i), tp%xh(1,i), tp%xh(2,i), tp%xh(3,i), & + tp%vh(1,i), tp%vh(2,i), tp%vh(3,i), & + dtp(i), iflag(i)) end do if (any(iflag(1:ntp) /= 0)) then - where(iflag(1:ntp) /= 0) status(1:ntp) = DISCARDED_DRIFTERR + where(iflag(1:ntp) /= 0) tp%status(1:ntp) = DISCARDED_DRIFTERR do i = 1, ntp if (iflag(i) /= 0) write(*, *) "Particle ", self%name(i), " lost due to error in Danby drift" end do @@ -118,6 +108,5 @@ module subroutine whm_drift_tp(self, system, param, dt) end associate return - end subroutine whm_drift_tp end submodule whm_drift diff --git a/src/whm/whm_getacch.f90 b/src/whm/whm_getacch.f90 index 53cd9b97a..d2e1bdc3e 100644 --- a/src/whm/whm_getacch.f90 +++ b/src/whm/whm_getacch.f90 @@ -10,15 +10,15 @@ module subroutine whm_getacch_pl(self, system, param, t) !! Adapted from David E. Kaufmann's Swifter routine whm_getacch.f90 implicit none ! Arguments - class(whm_pl), intent(inout) :: self !! WHM massive body particle data structure - class(whm_nbody_system), intent(inout) :: system !! Swiftest central body particle data structure - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters of - real(DP), intent(in) :: t !! Current time + class(whm_pl), intent(inout) :: self !! WHM massive body particle data structure + class(swiftest_nbody_system), intent(inout) :: system !! Swiftest central body particle data structure + class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters of + real(DP), intent(in) :: t !! Current time ! Internals - integer(I4B) :: i - real(DP), dimension(NDIM) :: ah0 + integer(I4B) :: i + real(DP), dimension(NDIM) :: ah0 - associate(pl => self, cb => system%cb, npl => self%nbody) + associate(pl => self, npl => self%nbody) if (npl == 0) return call pl%set_ir3() @@ -26,19 +26,23 @@ module subroutine whm_getacch_pl(self, system, param, t) do i = 1, npl pl%ah(:, i) = ah0(:) end do - call whm_getacch_ah1(cb, pl) - call whm_getacch_ah2(cb, pl) - call whm_getacch_ah3(pl) + select type(cb => system%cb) + class is (whm_cb) + + call whm_getacch_ah1(cb, pl) + call whm_getacch_ah2(cb, pl) + call whm_getacch_ah3(pl) - if (param%loblatecb) call pl%obl_acc(cb) - if (param%lextra_force) call pl%user_getacch(cb, param, t) - if (param%lgr) call pl%gr_getacch(param) + 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) + end select end associate return end subroutine whm_getacch_pl - module subroutine whm_getacch_tp(self, system, param, t, xh) + module subroutine whm_getacch_tp(self, system, param, t, xhp) !! author: David A. Minton !! !! Compute heliocentric accelerations of test particles @@ -47,39 +51,38 @@ module subroutine whm_getacch_tp(self, system, param, t, xh) !! Adapted from David E. Kaufmann's Swifter routine whm_getacch_tp.f90 implicit none ! Arguments - class(whm_tp), intent(inout) :: self !! WHM test particle data structure - class(whm_nbody_system), intent(inout) :: system !! Swiftest central body particle data structure - class(whm_pl), intent(inout) :: pl !! Generic Swiftest massive body particle data structure. - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters of - real(DP), intent(in) :: t !! Current time - real(DP), dimension(:,:), intent(in) :: xh !! Heliocentric positions of planets + class(whm_tp), intent(inout) :: self !! WHM test particle data structure + class(swiftest_nbody_system), intent(inout) :: system !! Swiftest 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), dimension(:,:), intent(in) :: xhp !! Heliocentric positions of planets at the current substep ! Internals - integer(I4B) :: i - real(DP), dimension(NDIM) :: ah0 + integer(I4B) :: i + real(DP), dimension(NDIM) :: ah0 - associate(tp => self, ntp => self%nbody, pl => system%pl, npl => system%pl%nbody) + associate(tp => self, ntp => self%nbody, pl => system%pl, cb => system%cb, npl => system%pl%nbody) if (ntp == 0 .or. npl == 0) return - ah0 = whm_getacch_ah0(pl%Gmass(:), xh(:,:), npl) + ah0 = whm_getacch_ah0(pl%Gmass(:), xhp(:,:), npl) do i = 1, ntp tp%ah(:, i) = ah0(:) end do - call whm_getacch_ah3_tp(system, xh) + call whm_getacch_ah3_tp(system) if (param%loblatecb) call tp%obl_acc(cb) - if (param%lextra_force) call tp%user_getacch(cb, param, t) + if (param%lextra_force) call tp%user_getacch(system, param, t) if (param%lgr) call tp%gr_getacch(param) end associate return end subroutine whm_getacch_tp - function whm_getacch_ah0(mu, xh, n) result(ah0) + function whm_getacch_ah0(mu, xhp, n) result(ah0) !! author: David A. Minton !! !! Compute zeroth term heliocentric accelerations of planets implicit none ! Arguments real(DP), dimension(:), intent(in) :: mu - real(DP), dimension(:,:), intent(in) :: xh + real(DP), dimension(:,:), intent(in) :: xhp integer(I4B), intent(in) :: n ! Result real(DP), dimension(NDIM) :: ah0 @@ -89,10 +92,10 @@ function whm_getacch_ah0(mu, xh, n) result(ah0) ah0(:) = 0.0_DP do i = 1, n - r2 = dot_product(xh(:, i), xh(:, i)) + r2 = dot_product(xhp(:, i), xhp(:, i)) ir3h = 1.0_DP / (r2 * sqrt(r2)) fac = mu(i) * ir3h - ah0(:) = ah0(:) - fac * xh(:, i) + ah0(:) = ah0(:) - fac * xhp(:, i) end do return @@ -107,17 +110,17 @@ pure subroutine whm_getacch_ah1(cb, pl) !! Adapted from David E. Kaufmann's Swifter routine whm_getacch_ah1.f90 implicit none ! Arguments - class(swiftest_cb), intent(in) :: cb !! Swiftest central body object - class(whm_pl), intent(inout) :: pl !! WHM massive body object + class(whm_cb), intent(in) :: cb !! WHM central body object + class(whm_pl), intent(inout) :: pl !! WHM massive body object ! Internals - integer(I4B) :: i - real(DP), dimension(NDIM) :: ah1h, ah1j + integer(I4B) :: i + real(DP), dimension(NDIM) :: ah1h, ah1j - associate(npl => pl%nbody, msun => cb%Gmass, xh => pl%xh, xj => pl%xj, ir3j => pl%ir3j, ir3h => pl%ir3h ) + associate(npl => pl%nbody) do i = 2, npl - ah1j(:) = xj(:, i) * ir3j(i) - ah1h(:) = xh(:, i) * ir3h(i) - pl%ah(:, i) = pl%ah(:, i) + msun * (ah1j(:) - ah1h(:)) + ah1j(:) = pl%xj(:, i) * pl%ir3j(i) + ah1h(:) = pl%xh(:, i) * pl%ir3h(i) + pl%ah(:, i) = pl%ah(:, i) + cb%Gmass * (ah1j(:) - ah1h(:)) end do end associate @@ -134,21 +137,21 @@ pure subroutine whm_getacch_ah2(cb, pl) !! Adapted from David E. Kaufmann's Swifter routine whm_getacch_ah2.f90 implicit none ! Arguments - class(swiftest_cb), intent(in) :: cb !! Swiftest central body object - class(whm_pl), intent(inout) :: pl !! WHM massive body object + class(whm_cb), intent(in) :: cb !! Swiftest central body object + class(whm_pl), intent(inout) :: pl !! WHM massive body object ! Internals - integer(I4B) :: i - real(DP) :: etaj, fac - real(DP), dimension(NDIM) :: ah2, ah2o + integer(I4B) :: i + real(DP) :: etaj, fac + real(DP), dimension(NDIM) :: ah2, ah2o - associate(npl => pl%nbody, Gmsun => cb%Gmass, xh => pl%xh, xj => pl%xj, Gmpl => pl%Gmass, ir3j => pl%ir3j) + associate(npl => pl%nbody, GMcb => cb%Gmass) ah2(:) = 0.0_DP ah2o(:) = 0.0_DP - etaj = Gmsun + etaj = GMcb do i = 2, npl - etaj = etaj + Gmpl(i - 1) - fac = Gmpl(i) * Gmsun * ir3j(i) / etaj - ah2(:) = ah2o + fac * xj(:, i) + etaj = etaj + pl%Gmass(i - 1) + fac = pl%Gmass(i) * GMcb * pl%ir3j(i) / etaj + ah2(:) = ah2o + fac * pl%xj(:, i) pl%ah(:,i) = pl%ah(:, i) + ah2(:) ah2o(:) = ah2(:) end do @@ -172,17 +175,17 @@ pure subroutine whm_getacch_ah3(pl) real(DP), dimension(NDIM) :: dx real(DP), dimension(:,:), allocatable :: ah3 - associate(npl => pl%nbody, xh => pl%xh, Gmpl => pl%Gmass) + associate(npl => pl%nbody) allocate(ah3, mold=pl%ah) ah3(:, 1:npl) = 0.0_DP do i = 1, npl - 1 do j = i + 1, npl - dx(:) = xh(:, j) - xh(:, i) + dx(:) = pl%xh(:, j) - pl%xh(:, i) rji2 = dot_product(dx(:), dx(:)) irij3 = 1.0_DP / (rji2 * sqrt(rji2)) - faci = Gmpl(i) * irij3 - facj = Gmpl(j) * irij3 + faci = pl%Gmass(i) * irij3 + facj = pl%Gmass(j) * irij3 ah3(:, i) = ah3(:, i) + facj * dx(:) ah3(:, j) = ah3(:, j) - faci * dx(:) end do @@ -196,7 +199,7 @@ pure subroutine whm_getacch_ah3(pl) return end subroutine whm_getacch_ah3 - pure subroutine whm_getacch_ah3_tp(system, xh) + pure subroutine whm_getacch_ah3_tp(system) !! author: David A. Minton !! !! Compute direct cross (third) term heliocentric accelerations of test particles @@ -205,25 +208,24 @@ pure subroutine whm_getacch_ah3_tp(system, xh) !! Adapted from David E. Kaufmann's Swifter routine whm_getacch_ah3.f90 implicit none ! Arguments - class(whm_nbody_system) :: system !! WHM nbody system object - real(DP), dimension(:,:), intent(in) :: xh !! Position vector of massive bodies at required point in step + class(swiftest_nbody_system), intent(inout) :: system !! WHM nbody system object ! Internals - integer(I4B) :: i, j - real(DP) :: rji2, irij3, fac - real(DP), dimension(NDIM) :: dx, acc + integer(I4B) :: i, j + real(DP) :: rji2, irij3, fac + real(DP), dimension(NDIM) :: dx, acc - associate(ntp => system%tp%nbody, npl => system%pl%nbody, msun => system%cb%Gmass, GMpl => system%pl%Gmass, xht => system%tp%xh, aht => system%tp%ah) + associate(ntp => system%tp%nbody, npl => system%pl%nbody, tp => system%tp, pl => system%pl) if (ntp == 0) return do i = 1, ntp acc(:) = 0.0_DP do j = 1, npl - dx(:) = xht(:, i) - xh(:, j) + dx(:) = tp%xh(:, i) - pl%xh(:, j) rji2 = dot_product(dx(:), dx(:)) irij3 = 1.0_DP / (rji2 * sqrt(rji2)) - fac = GMpl(j) * irij3 + fac = pl%Gmass(j) * irij3 acc(:) = acc(:) - fac * dx(:) end do - aht(:, i) = aht(:, i) + acc(:) + tp%ah(:, i) = tp%ah(:, i) + acc(:) end do end associate return diff --git a/src/whm/whm_setup.f90 b/src/whm/whm_setup.f90 index cdec4e1fb..e0812d00d 100644 --- a/src/whm/whm_setup.f90 +++ b/src/whm/whm_setup.f90 @@ -133,7 +133,7 @@ module subroutine whm_setup_set_beg_end(self, xbeg, xend, vbeg) !! Sets one or more of the values of xbeg and xend implicit none ! Arguments - class(whm_tp), intent(inout) :: self !! Swiftest test particle object + class(whm_nbody_system), intent(inout) :: self !! WHM nbody system object real(DP), dimension(:,:), intent(in), optional :: xbeg, xend real(DP), dimension(:,:), intent(in), optional :: vbeg ! vbeg is an unused variable to keep this method forward compatible with RMVS diff --git a/src/whm/whm_step.f90 b/src/whm/whm_step.f90 index c3c7a098f..b37d06bbe 100644 --- a/src/whm/whm_step.f90 +++ b/src/whm/whm_step.f90 @@ -2,7 +2,7 @@ use swiftest contains - module subroutine whm_step_system(self, param, dt) + module subroutine whm_step_system(self, param, t, dt) !! author: David A. Minton !! !! Step massive bodies and and active test particles ahead in heliocentric coordinates @@ -11,35 +11,24 @@ module subroutine whm_step_system(self, param, dt) !! Adapted from David E. Kaufmann's Swifter routine whm_step.f90 implicit none ! Arguments - class(whm_nbody_system), intent(inout) :: self !! WHM nbody system object - class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters of on parameters - integer(I4B), intent(in) :: dt !! Current stepsize + class(whm_nbody_system), intent(inout) :: self !! WHM nbody system object + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters of on parameters + real(DP), intent(in) :: t !! Current simulation time + real(DP), intent(in) :: dt !! Current stepsize - select type (system => self) - class is (whm_nbody_system) - select type(cb => self%cb) - class is (whm_cb) - select type(pl => self%pl) - class is (whm_pl) - select type(tp => self%tp) - class is (whm_tp) - associate(ntp => tp%nbody, npl => pl%nbody, t => param%t, dt => param%dt) + associate(system => self, cb => self%cb, pl => self%pl, tp => self%tp, ntp => self%tp%nbody, npl => self%pl%nbody) call pl%set_rhill(cb) - call tp%set_beg_end(xbeg = pl%xh) - call pl%step(system, param, dt) + call self%set_beg_end(xbeg = pl%xh) + call pl%step(system, param, t, dt) if (ntp > 0) then - call tp%set_beg_end(xend = pl%xh) - call tp%step(system, param, dt) + call self%set_beg_end(xend = pl%xh) + call tp%step(system, param, t, dt) end if end associate - end select - end select - end select - end select return end subroutine whm_step_system - module subroutine whm_step_pl(self, system, param, dt) + module subroutine whm_step_pl(self, system, param, t, dt) !! author: David A. Minton !! !! Step planets ahead using kick-drift-kick algorithm @@ -52,23 +41,23 @@ module subroutine whm_step_pl(self, system, param, dt) class(whm_pl), intent(inout) :: self !! WHM massive body particle data structure class(swiftest_nbody_system), intent(inout) :: system !! Swiftest system object class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters - integer(I4B), intent(in) :: dt !! Current stepsize + real(DP), intent(in) :: t !! Current simulation time + real(DP), intent(in) :: dt !! Current stepsize ! Internals - real(DP) :: dth + real(DP) :: dth - associate(cb => system%cb, t => param%t) + associate(pl => self, cb => system%cb) dth = 0.5_DP * dt if (pl%lfirst) then call pl%h2j(cb) call pl%getacch(system, param, t) pl%lfirst = .false. end if - call pl%kickvh(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(cb, param, dt) + 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) @@ -77,7 +66,7 @@ module subroutine whm_step_pl(self, system, param, dt) return end subroutine whm_step_pl - module subroutine whm_step_tp(self, system, param, dt) + module subroutine whm_step_tp(self, system, param, t, dt) !! author: David A. Minton !! !! Step active test particles ahead using kick-drift-kick algorithm @@ -89,24 +78,28 @@ module subroutine whm_step_tp(self, system, param, dt) class(whm_tp), intent(inout) :: self !! WHM test particle data structure class(swiftest_nbody_system), intent(inout) :: system !! Swiftest system object class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters - integer(I4B), intent(in) :: dt !! Current stepsize + real(DP), intent(in) :: t !! Current simulation time + real(DP), intent(in) :: dt !! Current stepsize ! Internals - real(DP) :: dth + real(DP) :: dth - associate(cb => system%cb, pl => system%pl, t => param%t, xbeg => self%xbeg, xend => self%xend) - dth = 0.5_DP * dt - if (tp%lfirst) then - call tp%getacch(system, param, t, xbeg) - tp%lfirst = .false. - end if - call tp%kickvh(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(cb, param, dt) - if (param%lgr) call tp%gr_p4(param, dth) - call tp%getacch(system, param, t + dt, xend) - call tp%kickvh(dth) - end associate + select type(system) + 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%getacch(system, param, t, system%xbeg) + tp%lfirst = .false. + end if + call tp%kickvh(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, system%xend) + call tp%kickvh(dth) + end associate + end select return end subroutine whm_step_tp