From e80355fb5058ce7913e3fb3402d79c37a24a6042 Mon Sep 17 00:00:00 2001 From: David Minton Date: Mon, 19 Apr 2021 11:59:40 -0400 Subject: [PATCH] Numerous small changes and refactoring for clarity and hunting for source of divergence with original Swifter RMVS --- src/driftkick/drift.f90 | 5 +- src/helio/helio_drift.f90 | 72 +++++++------ src/modules/rmvs_classes.f90 | 17 +-- src/modules/swiftest_classes.f90 | 3 +- src/rmvs/rmvs_getacch.f90 | 17 +-- src/rmvs/rmvs_setup.f90 | 4 +- src/rmvs/rmvs_step.f90 | 180 +++++++++++++++++-------------- src/whm/whm_drift.f90 | 33 +++--- 8 files changed, 184 insertions(+), 147 deletions(-) diff --git a/src/driftkick/drift.f90 b/src/driftkick/drift.f90 index 68d9d2822..bea9ef2d8 100644 --- a/src/driftkick/drift.f90 +++ b/src/driftkick/drift.f90 @@ -9,8 +9,7 @@ integer(I2B), parameter :: NLAG2 = 40 contains - module pure subroutine drift_one(mu, px, py, pz, vx, vy, vz, dt, iflag) - !$omp declare simd + module pure elemental subroutine drift_one(mu, px, py, pz, vx, vy, vz, dt, iflag) !! author: The Purdue Swiftest Team - David A. Minton, Carlisle A. Wishard, Jennifer L.L. Pouplin, and Jacob R. Elliott !! !! Perform Danby drift for one body, redoing drift with smaller substeps if original accuracy is insufficient @@ -36,7 +35,7 @@ module pure subroutine drift_one(mu, px, py, pz, vx, vy, vz, dt, iflag) dttmp = 0.1_DP * dt do i = 1, 10 call drift_dan(mu, x(:), v(:), dttmp, iflag) - if (iflag /= 0) return + if (iflag /= 0) exit end do end if px = x(1); py = x(2); pz = x(3) diff --git a/src/helio/helio_drift.f90 b/src/helio/helio_drift.f90 index 522b3403b..e03ff98bf 100644 --- a/src/helio/helio_drift.f90 +++ b/src/helio/helio_drift.f90 @@ -16,9 +16,10 @@ module subroutine helio_drift_pl(self, cb, config, dt) class(swiftest_configuration), intent(in) :: config !! Input collection of real(DP), intent(in) :: dt !! Stepsize ! Internals - integer(I4B), dimension(:),allocatable :: iflag !! Vectorized error code flag integer(I4B) :: i !! Loop counter - real(DP) :: rmag, vmag2, energy, dtp + 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, & @@ -30,23 +31,28 @@ module subroutine helio_drift_pl(self, cb, config, dt) allocate(iflag(npl)) iflag(:) = 0 - - do i = 1, npl - if (config%lgr) then + allocate(dtp(npl)) + + if (config%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 - dtp = dt * (1.0_DP + 3 * config%inv_c2 * energy) - else - dtp = dt - end if - call drift_one(mu(i), xh(1, i), xh(2, i), xh(3, i), vb(1, i), vb(2, i), vb(3, i), dtp, iflag(i)) - end do + dtp(i) = dt * (1.0_DP + 3 * config%inv_c2 * energy) + end do + else + dtp(:) = dt + end if + + !!$omp simd + 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)) if (any(iflag(1:npl) /= 0)) then do i = 1, npl write(*, *) " Planet ", self%name(i), " is lost!!!!!!!!!!" - write(*, *) xh - write(*, *) vb + write(*, *) xh(:,i) + write(*, *) vb(:,i) write(*, *) " stopping " call util_exit(FAILURE) end do @@ -103,33 +109,39 @@ module subroutine helio_drift_tp(self, cb, config, dt) class(swiftest_configuration), intent(in) :: config !! Input collection of real(DP), intent(in) :: dt !! Stepsize ! Internals - integer(I4B), dimension(:),allocatable :: iflag !! Vectorized error code flag integer(I4B) :: i !! Loop counter - real(DP) :: rmag, vmag2, energy, dtp + 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, & - vb => self%vb, & + vh => self%vh, & status => self%status,& mu => self%mu) if (ntp == 0) return allocate(iflag(ntp)) + allocate(dtp(ntp)) iflag(:) = 0 - do i = 1,ntp - if (status(i) == ACTIVE) then - if (config%lgr) then - rmag = norm2(xh(:, i)) - vmag2 = dot_product(vb(:, i), vb(:, i)) - energy = 0.5_DP * vmag2 - cb%Gmass / rmag - dtp = dt * (1.0_DP + 3 * config%inv_c2 * energy) - else - dtp = dt - end if - call drift_one(mu(i), xh(1, i), xh(2, i), xh(3, i), vb(1, i), vb(2, i), vb(3, i), dtp, iflag(i)) - if (iflag(i) /= 0) status = DISCARDED_DRIFTERR - end if - end do + + iflag(:) = 0 + allocate(dtp(ntp)) + + if (config%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 + dtp(i) = dt * (1.0_DP + 3 * config%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)) if (any(iflag(1:ntp) /= 0)) then + status = DISCARDED_DRIFTERR do i = 1, ntp if (iflag(i) /= 0) write(*, *) "Particle ", self%name(i), " lost due to error in Danby drift" end do diff --git a/src/modules/rmvs_classes.f90 b/src/modules/rmvs_classes.f90 index f32122d2a..b5445c156 100644 --- a/src/modules/rmvs_classes.f90 +++ b/src/modules/rmvs_classes.f90 @@ -29,14 +29,20 @@ module rmvs_classes procedure, public :: step => rmvs_step_system end type rmvs_nbody_system + type, private :: rmvs_interp + real(DP), dimension(:, :), allocatable :: x !! interpolated heliocentric planet position for outer encounter + real(DP), dimension(:, :), allocatable :: v !! interpolated heliocentric planet velocity for outer encounter + real(DP), dimension(:, :), allocatable :: aobl !! Encountering planet's oblateness acceleration value + end type rmvs_interp + !******************************************************************************************************************************** ! rmvs_cb class definitions and method interfaces !******************************************************************************************************************************* !> 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 - real(DP), dimension(NDIM) :: xin = [0.0_DP, 0.0_DP, 0.0_DP] - real(DP), dimension(NDIM) :: vin = [0.0_DP, 0.0_DP, 0.0_DP] + 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 end type rmvs_cb !******************************************************************************************************************************** @@ -73,12 +79,7 @@ module rmvs_classes !******************************************************************************************************************************** ! rmvs_pl class definitions and method interfaces !******************************************************************************************************************************* - - type, private :: rmvs_interp - real(DP), dimension(:, :), allocatable :: x !! interpolated heliocentric planet position for outer encounter - real(DP), dimension(:, :), allocatable :: v !! interpolated heliocentric planet velocity for outer encounter - real(DP), dimension(:, :), allocatable :: aobl !! Encountering planet's oblateness acceleration value - end type rmvs_interp + !> RMVS massive body particle class type, private, extends(whm_pl) :: rmvs_base_pl diff --git a/src/modules/swiftest_classes.f90 b/src/modules/swiftest_classes.f90 index b0e7c4b67..3a8ae9bfd 100644 --- a/src/modules/swiftest_classes.f90 +++ b/src/modules/swiftest_classes.f90 @@ -356,8 +356,7 @@ module subroutine discard_system(self, config) class(swiftest_configuration), intent(in) :: config !! Input collection of configuration parameters end subroutine discard_system - module pure subroutine drift_one(mu, px, py, pz, vx, vy, vz, dt, iflag) - !$omp declare simd(drift_one) + module pure elemental subroutine drift_one(mu, px, py, pz, vx, vy, vz, dt, iflag) implicit none real(DP), intent(in) :: mu !! G * (Mcb + m), G = gravitational constant, Mcb = mass of central body, m = mass of body to drift real(DP), intent(inout) :: px, py, pz, vx, vy, vz !! Position and velocity of body to drift diff --git a/src/rmvs/rmvs_getacch.f90 b/src/rmvs/rmvs_getacch.f90 index 1628fc8d2..21e80a716 100644 --- a/src/rmvs/rmvs_getacch.f90 +++ b/src/rmvs/rmvs_getacch.f90 @@ -28,7 +28,10 @@ module subroutine rmvs_getacch_tp(self, cb, pl, config, t, xh) 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) + class is (rmvs_pl) + select type (cb) + class is (rmvs_cb) + associate(xpc => pl%xh, xpct => self%xh) allocate(xh_original, source=tp%xh) config_planetocen = config ! Temporarily turn off the heliocentric-dependent acceleration terms during an inner encounter @@ -36,16 +39,16 @@ module subroutine rmvs_getacch_tp(self, cb, pl, config, t, xh) config_planetocen%lextra_force = .false. config_planetocen%lgr = .false. ! Now compute the planetocentric values of acceleration - call whm_getacch_tp(tp, cb, pl%planetocentric(ipleP)%pl, config_planetocen, t, pl%xh) + call whm_getacch_tp(tp, cb, pl, config_planetocen, t, pl%xh) ! Now compute any heliocentric values of acceleration if (tp%lfirst) then do i = 1, ntp - tp%xheliocentric(:,i) = tp%xh(:,i) + pl%inner(enc_index - 1)%x(:,ipleP) + tp%xheliocentric(:,i) = tp%xh(:,i) + cb%inner(enc_index - 1)%x(:,1) end do else do i = 1, ntp - tp%xheliocentric(:,i) = tp%xh(:,i) + pl%inner(enc_index )%x(:,ipleP) + tp%xheliocentric(:,i) = tp%xh(:,i) + cb%inner(enc_index )%x(:,1) end do end if ! Swap the planetocentric and heliocentric position vectors @@ -53,9 +56,9 @@ module subroutine rmvs_getacch_tp(self, cb, pl, config, t, xh) if (config%loblatecb) then ! Put in the current encountering planet's oblateness acceleration as the central body's if (tp%lfirst) then - cb_heliocentric%aobl(:) = pl%inner(enc_index - 1)%aobl(:, ipleP) + cb_heliocentric%aobl(:) = cb%inner(enc_index - 1)%aobl(:,1) else - cb_heliocentric%aobl(:) = pl%inner(enc_index )%aobl(:, ipleP) + cb_heliocentric%aobl(:) = cb%inner(enc_index )%aobl(:,1) end if call tp%obl_acc(cb_heliocentric) end if @@ -64,6 +67,8 @@ module subroutine rmvs_getacch_tp(self, cb, pl, config, t, xh) if (config%lgr) call tp%gr_getacch(cb_heliocentric, config) call move_alloc(xh_original, tp%xh) + 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, config, t, pl%xh) diff --git a/src/rmvs/rmvs_setup.f90 b/src/rmvs/rmvs_setup.f90 index 23a00cce9..c135f8a1a 100644 --- a/src/rmvs/rmvs_setup.f90 +++ b/src/rmvs/rmvs_setup.f90 @@ -20,17 +20,17 @@ module subroutine rmvs_setup_pl(self,n) call whm_setup_pl(pl, n) if (n <= 0) return + allocate(pl%outer(0:NTENC)) + allocate(pl%inner(0:NTPHENC)) if (.not.pl%lplanetocentric) then allocate(pl%nenc(n)) pl%nenc(:) = 0 ! Set up inner and outer planet interpolation vector storage containers - allocate(pl%outer(0:NTENC)) do i = 0, NTENC allocate(pl%outer(i)%x(NDIM, n)) allocate(pl%outer(i)%v(NDIM, n)) end do - allocate(pl%inner(0:NTPHENC)) do i = 0, NTPHENC allocate(pl%inner(i)%x(NDIM, n)) allocate(pl%inner(i)%v(NDIM, n)) diff --git a/src/rmvs/rmvs_step.f90 b/src/rmvs/rmvs_step.f90 index 01dff9e9c..ae4712efb 100644 --- a/src/rmvs/rmvs_step.f90 +++ b/src/rmvs/rmvs_step.f90 @@ -24,7 +24,9 @@ module subroutine rmvs_step_system(self, config) class is (rmvs_pl) select type(tp => self%tp) class is (rmvs_tp) - associate(ntp => tp%nbody, npl => pl%nbody, t => config%t, dt => config%dt) + associate(ntp => tp%nbody, npl => pl%nbody, t => config%t, dt => config%dt, & + xhpl => pl%xh, vhpl => pl%vh, xjpl => pl%xj, vjpl => pl%vj, & + xhtp => tp%xh, vhtp => tp%vh) allocate(xbeg, source=pl%xh) allocate(vbeg, source=pl%vh) call pl%set_rhill(cb) @@ -61,7 +63,7 @@ module subroutine rmvs_step_system(self, config) end subroutine rmvs_step_system - subroutine rmvs_interp_in(pl, cb, dt, enc_index, config) + subroutine rmvs_interp_in(pl, cb, dt, outer_index, config) !! author: David A. Minton !! !! Interpolate planet positions between two Keplerian orbits in inner encounter regio @@ -74,10 +76,10 @@ subroutine rmvs_interp_in(pl, cb, dt, enc_index, config) 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) :: enc_index !! Outer substep number within current se + integer(I4B), intent(in) :: outer_index !! Outer substep number within current se class(swiftest_configuration), intent(in) :: config !! Swiftest configuration file ! Internals - integer(I4B) :: i, j + integer(I4B) :: i, inner_index real(DP) :: frac, dntphenc real(DP), dimension(:,:), allocatable :: xtmp, vtmp, xh_original real(DP), dimension(:), allocatable :: msun, dti @@ -87,10 +89,10 @@ subroutine rmvs_interp_in(pl, cb, dt, enc_index, config) dntphenc = real(NTPHENC, kind=DP) ! Set the endpoints of the inner region from the outer region values in the current outer step index - pl%inner(0)%x(:,:) = pl%outer(enc_index - 1)%x(:, :) - pl%inner(0)%v(:,:) = pl%outer(enc_index - 1)%v(:, :) - pl%inner(NTPHENC)%x(:,:) = pl%outer(enc_index)%x(:, :) - pl%inner(NTPHENC)%v(:,:) = pl%outer(enc_index)%v(:, :) + 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) @@ -108,12 +110,11 @@ subroutine rmvs_interp_in(pl, cb, dt, enc_index, config) pl%inner(0)%aobl(:, :) = pl%aobl(:, :) ! Save the oblateness acceleration on the planet for this substep end if - do j = 1, NTPHENC - 1 - !!$omp simd - do i = 1, npl - call drift_one(msun(i), xtmp(1,i), xtmp(2,i), xtmp(3,i), vtmp(1,i), vtmp(2,i), vtmp(3,i), dti(i), iflag(i)) - end do - if (any(iflag(:) /= 0)) then + 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!!!!!!!!!!" @@ -125,20 +126,19 @@ subroutine rmvs_interp_in(pl, cb, dt, enc_index, config) end if end do end if - frac = 1.0_DP - j / dntphenc - pl%inner(j)%x(:, :) = frac * xtmp(:,:) - pl%inner(j)%v(:, :) = frac * vtmp(:,:) + 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 j = NTPHENC - 1, 1, -1 - !!$omp simd - do i = 1, npl - call drift_one(msun(i), xtmp(1,i), xtmp(2,i), xtmp(3,i), vtmp(1,i), vtmp(2,i), vtmp(3,i), -dti(i), iflag(i)) - end do - if (any(iflag(:) /= 0)) then + 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!!!!!!!!!!" @@ -150,14 +150,14 @@ subroutine rmvs_interp_in(pl, cb, dt, enc_index, config) end if end do end if - frac = j / dntphenc - pl%inner(j)%x(:, :) = pl%inner(j)%x(:, :) + frac * xtmp(:, :) - pl%inner(j)%v(:, :) = pl%inner(j)%v(:, :) + frac * vtmp(:, :) + 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 (config%loblatecb) then - pl%xh(:,:) = pl%inner(j)%x(:, :) + pl%xh(:,:) = pl%inner(inner_index)%x(:, :) call pl%obl_acc(cb) - pl%inner(j)%aobl(:, :) = pl%aobl(:, :) + pl%inner(inner_index)%aobl(:, :) = pl%aobl(:, :) end if end do ! Put the planet positions back into place @@ -182,7 +182,7 @@ subroutine rmvs_interp_out(pl, cb, dt, config) real(DP), intent(in) :: dt !! Step size class(swiftest_configuration), intent(in) :: config !! Swiftest configuration file ! Internals - integer(I4B) :: i, j + integer(I4B) :: i, outer_index real(DP) :: frac, dntenc real(DP), dimension(:,:), allocatable :: xtmp, vtmp real(DP), dimension(:), allocatable :: msun, dto @@ -200,12 +200,12 @@ subroutine rmvs_interp_out(pl, cb, dt, config) msun(:) = cb%Gmass xtmp(:,:) = pl%outer(0)%x(:, :) vtmp(:,:) = pl%outer(0)%v(:, :) - do j = 1, NTENC - 1 - !!$omp simd - do i = 1, npl - call drift_one(msun(i), xtmp(1,i), xtmp(2,i), xtmp(3,i), vtmp(1,i), vtmp(2,i), vtmp(3,i), dto(i), iflag(i)) - end do - if (any(iflag(:) /= 0)) then + 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)) + + if (any(iflag(1:npl) /= 0)) then do i = 1, npl if (iflag(i) /= 0) then write(*, *) " Planet ", pl%name(i), " is lost!!!!!!!!!!" @@ -217,18 +217,17 @@ subroutine rmvs_interp_out(pl, cb, dt, config) end if end do end if - frac = 1.0_DP - j / dntenc - pl%outer(j)%x(:, :) = frac * xtmp(:,:) - pl%outer(j)%v(:, :) = frac * vtmp(:,:) + 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 j = NTENC - 1, 1, -1 - !!$omp simd - do i = 1, npl - call drift_one(msun(i), xtmp(1,i), xtmp(2,i), xtmp(3,i), vtmp(1,i), vtmp(2,i), vtmp(3,i), -dto(i), iflag(i)) - end do - if (any(iflag(:) /= 0)) then + 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!!!!!!!!!!" @@ -240,9 +239,9 @@ subroutine rmvs_interp_out(pl, cb, dt, config) end if end do end if - frac = j / dntenc - pl%outer(j)%x(:, :) = pl%outer(j)%x(:, :) + frac * xtmp(:,:) - pl%outer(j)%v(:, :) = pl%outer(j)%v(:, :) + frac * vtmp(:,:) + 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 associate @@ -250,7 +249,7 @@ subroutine rmvs_interp_out(pl, cb, dt, config) end subroutine rmvs_interp_out - subroutine rmvs_peri_tp(tp, pl, t, dt, lfirst, enc_index, ipleP, config) + subroutine rmvs_peri_tp(tp, pl, t, dt, lfirst, inner_index, ipleP, config) !! author: David A. Minton !! !! Determine planetocentric pericenter passages for test particles in close encounters with a planet @@ -264,7 +263,7 @@ subroutine rmvs_peri_tp(tp, pl, t, dt, lfirst, enc_index, ipleP, config) real(DP), intent(in) :: t !! current time real(DP), intent(in) :: dt !! step size logical, intent(in) :: lfirst !! Logical flag indicating whether current invocation is the first - integer(I4B), intent(in) :: enc_index !! Outer substep number within current set + integer(I4B), intent(in) :: inner_index !! Outer substep number within current set integer(I4B), intent(in) :: ipleP !! index of RMVS planet being closely encountered class(swiftest_configuration), intent(in) :: config !! Input collection of configuration parameters ! Internals @@ -299,8 +298,8 @@ subroutine rmvs_peri_tp(tp, pl, t, dt, lfirst, enc_index, ipleP, config) if (config%encounter_file /= "") then id1 = pl%name(ipleP) rpl = pl%radius(ipleP) - xh1(:) = pl%inner(enc_index)%x(:, ipleP) - vh1(:) = pl%inner(enc_index)%v(:, ipleP) + xh1(:) = pl%inner(inner_index)%x(:, ipleP) + vh1(:) = pl%inner(inner_index)%v(:, ipleP) id2 = tp%name(i) xh2(:) = xpc(:, i) + xh1(:) vh2(:) = xpc(:, i) + vh1(:) @@ -349,7 +348,7 @@ subroutine rmvs_step_out(pl, cb, tp, dt, config) real(DP), intent(in) :: dt !! Step size class(swiftest_configuration), intent(in) :: config !! Input collection of configuration parameters ! Internals - integer(I4B) :: enc_index, j, k + integer(I4B) :: outer_index, j, k real(DP) :: dto, outer_time, rts logical :: lencounter, lfirsttp @@ -360,16 +359,16 @@ subroutine rmvs_step_out(pl, cb, tp, dt, config) elsewhere tp%lperi(:) = .false. end where - do enc_index = 1, NTENC - outer_time = t + (enc_index - 1) * dto - call tp%set_beg_end(xbeg = pl%outer(enc_index - 1)%x(:, :), & - vbeg = pl%outer(enc_index - 1)%v(:, :), & - xend = pl%outer(enc_index )%x(:, :)) + 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(pl, cb, dto, enc_index, config) + call rmvs_interp_in(pl, cb, dto, outer_index, config) ! Step through the inner region call rmvs_step_in(pl, cb, tp, config, outer_time, dto) lfirsttp = tp%lfirst @@ -422,32 +421,23 @@ subroutine rmvs_step_in(pl, cb, tp, config, outer_time, dto) if (nenc(i) == 0) cycle associate(cbenci => pl%planetocentric(i)%cb, plenci => pl%planetocentric(i)%pl, & tpenci => pl%planetocentric(i)%tp) - associate(enc_index => tpenci%index, plind => plenci%plind) - ! There are inner encounters with this planet...switch to planetocentric coordinates to proceed + 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, config) ! now step the encountering test particles fully through the inner encounter lfirsttp = .true. - do enc_index = 1, NTPHENC ! Integrate over the encounter region, using the "substitute" planetocentric systems at each level - xbeg(:,1) = cb%xin(:) - pl%inner(enc_index - 1)%x(:, i) - xend(:,1) = cb%xin(:) - pl%inner(enc_index )%x(:, i) - vbeg(:,1) = cb%vin(:) - pl%inner(enc_index - 1)%v(:, i) - do j = 2, npl - ipleP = plind(j) - xbeg(:,j) = pl%inner(enc_index - 1)%x(:,ipleP) - pl%inner(enc_index - 1)%x(:,i) - xend(:,j) = pl%inner(enc_index )%x(:,ipleP) - pl%inner(enc_index )%x(:,i) - vbeg(:,j) = pl%inner(enc_index - 1)%v(:,ipleP) - pl%inner(enc_index - 1)%v(:,i) - end do - plenci%xh(:,:) = xbeg(:,:) - plenci%vh(:,:) = vbeg(:,:) - call tpenci%set_beg_end(xbeg = xbeg, xend = xend) - call tpenci%step(cbenci, pl, config, inner_time, dti) + 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, config, inner_time, dti) do j = 1, nenc(i) - tpenci%xheliocentric(:, j) = tpenci%xh(:, j) + pl%inner(enc_index )%x(:,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., enc_index, i, config) + call rmvs_peri_tp(tpenci, pl, inner_time, dti, .false., inner_index, i, config) end do where(tpenci%status(:) == ACTIVE) tpenci%status(:) = INACTIVE end associate @@ -472,7 +462,7 @@ subroutine rmvs_make_planetocentric(pl, cb, tp, config) class(rmvs_tp), intent(inout) :: tp !! RMVS test particle object class(swiftest_configuration), intent(in) :: config !! Input collection of configuration parameters ! Internals - integer(I4B) :: i, j + integer(I4B) :: i, j, inner_index, ipc2hc logical, dimension(:), allocatable :: encmask associate(npl => pl%nbody, ntp => tp%nbody, GMpl => pl%Gmass, nenc => pl%nenc) @@ -485,7 +475,9 @@ subroutine rmvs_make_planetocentric(pl, cb, tp, config) ! Create encountering test particle structure allocate(rmvs_tp :: pl%planetocentric(i)%tp) - associate(tpenci => pl%planetocentric(i)%tp, cbenci => pl%planetocentric(i)%cb, plenci => pl%planetocentric(i)%pl) + 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%ipleP = i @@ -498,6 +490,26 @@ subroutine rmvs_make_planetocentric(pl, cb, tp, config) 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 associate end do @@ -516,14 +528,16 @@ subroutine rmvs_end_planetocentric(pl, cb, tp) class(rmvs_cb), intent(inout) :: cb !! RMVS central body object class(rmvs_tp), intent(inout) :: tp !! RMVS test particle object ! Internals - integer(I4B) :: i, j + 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(tpenci => pl%planetocentric(i)%tp) + associate(cbenci => pl%planetocentric(i)%cb, & + tpenci => pl%planetocentric(i)%tp, & + plenci => pl%planetocentric(i)%pl) allocate(tpind(nenc(i))) ! Index array of encountering test particles if (allocated(encmask)) deallocate(encmask) @@ -538,6 +552,12 @@ subroutine rmvs_end_planetocentric(pl, cb, tp) 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 diff --git a/src/whm/whm_drift.f90 b/src/whm/whm_drift.f90 index 082953963..dd8cbe403 100644 --- a/src/whm/whm_drift.f90 +++ b/src/whm/whm_drift.f90 @@ -44,16 +44,18 @@ module subroutine whm_drift_pl(self, cb, config, dt) end if !!$omp simd - do i = 1, npl - call drift_one(mu(i), xj(1, i), xj(2, i), xj(3, i), vj(1, i), vj(2, i), vj(3, i), dtp(i), iflag(i)) - end do - if (any(iflag(:) /= 0)) then + 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), & + 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(*, *) xj(:,i) - write(*, *) vj(:,i) - write(*, *) " stopping " - call util_exit(FAILURE) + if (iflag(i) /= 0) then + write(*, *) " Planet ", self%name(i), " is lost!!!!!!!!!!" + write(*, *) xj(:,i) + write(*, *) vj(:,i) + write(*, *) " stopping " + call util_exit(FAILURE) + end if end do end if end associate @@ -103,14 +105,13 @@ module subroutine whm_drift_tp(self, cb, config, dt) else dtp(:) = dt end if - !!$omp simd - do i = 1,ntp - if (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, 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)) end do - if (any(iflag(:) /= 0)) then - where(iflag(:) /= 0) status(:) = DISCARDED_DRIFTERR + if (any(iflag(1:ntp) /= 0)) then + where(iflag(1:ntp) /= 0) 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