diff --git a/Makefile.Defines b/Makefile.Defines index 1663e016c..862e4bb0d 100644 --- a/Makefile.Defines +++ b/Makefile.Defines @@ -67,13 +67,13 @@ GWARNINGS = -Wall -Warray-bounds -Wimplicit-interface -Wextra -Warray-temporari GPRODUCTION = -O3 -ffree-line-length-none $(GPAR) #FFLAGS = $(IDEBUG) $(HEAPARR) $(SIMDVEC) $(PAR) -#FFLAGS = $(IPRODUCTION) -#FORTRAN = ifort +FFLAGS = $(IPRODUCTION) $(OPTREPORT) +FORTRAN = ifort #AR = xiar -FORTRAN = gfortran -#FFLAGS = $(GDEBUG) #$(GMEM) #$(GPAR) -FFLAGS = $(GPRODUCTION) +#FORTRAN = gfortran +#FFLAGS = $(GDEBUG) $(GMEM) $(GPAR) +#FFLAGS = $(GPRODUCTION) AR = ar # DO NOT include in CFLAGS the "-c" option to compile object only diff --git a/examples/symba_mars_disk/param.in b/examples/symba_mars_disk/param.in index f5760d345..13f933aa7 100644 --- a/examples/symba_mars_disk/param.in +++ b/examples/symba_mars_disk/param.in @@ -6,8 +6,8 @@ CB_IN cb.in PL_IN mars.in TP_IN tp.in IN_TYPE ASCII -ISTEP_OUT 100 -ISTEP_DUMP 100 +ISTEP_OUT 10 +ISTEP_DUMP 10 BIN_OUT bin.dat PARTICLE_OUT particle.dat OUT_TYPE REAL8 diff --git a/src/kick/kick.f90 b/src/kick/kick.f90 index af704c13c..c47c9f174 100644 --- a/src/kick/kick.f90 +++ b/src/kick/kick.f90 @@ -2,7 +2,7 @@ use swiftest contains - module pure subroutine kick_getacch_int_pl(self) + module subroutine kick_getacch_int_pl(self) !! author: David A. Minton !! !! Compute direct cross (third) term heliocentric accelerations of massive bodies @@ -13,33 +13,35 @@ module pure subroutine kick_getacch_int_pl(self) ! Arguments class(swiftest_pl), intent(inout) :: self ! Internals - integer(I4B) :: k - real(DP) :: rji2, irij3, faci, facj, rlim2 + integer(I8B) :: k, nplpl + real(DP) :: rji2, rlim2 real(DP) :: dx, dy, dz + integer(I4B) :: i, j + real(DP), dimension(:,:), pointer :: ah, xh + real(DP), dimension(NDIM,self%nbody) :: ahi, ahj + integer(I4B), dimension(:,:), pointer :: k_plpl + logical, dimension(:), pointer :: lmask + real(DP), dimension(:), pointer :: Gmass - associate(pl => self, npl => self%nbody, nplpl => self%nplpl) - do k = 1, nplpl - associate(i => pl%k_plpl(1, k), j => pl%k_plpl(2, k)) - if (pl%lmask(i) .and. pl%lmask(j)) then - dx = pl%xh(1, j) - pl%xh(1, i) - dy = pl%xh(2, j) - pl%xh(2, i) - dz = pl%xh(3, j) - pl%xh(3, i) - rji2 = dx**2 + dy**2 + dz**2 - rlim2 = (pl%radius(i) + pl%radius(j))**2 - if (rji2 > rlim2) then - irij3 = 1.0_DP / (rji2 * sqrt(rji2)) - faci = pl%Gmass(i) * irij3 - facj = pl%Gmass(j) * irij3 - pl%ah(1, i) = pl%ah(1, i) + facj * dx - pl%ah(2, i) = pl%ah(2, i) + facj * dy - pl%ah(3, i) = pl%ah(3, i) + facj * dz - pl%ah(1, j) = pl%ah(1, j) - faci * dx - pl%ah(2, j) = pl%ah(2, j) - faci * dy - pl%ah(3, j) = pl%ah(3, j) - faci * dz - end if - end if - end associate + associate(ah => self%ah, xh => self%xh, k_plpl => self%k_plpl, lmask => self%lmask, Gmass => self%Gmass) + nplpl = self%nplpl + ahi(:,:) = 0.0_DP + ahj(:,:) = 0.0_DP + !$omp parallel do default(shared)& + !$omp private(k, i, j, dx, dy, dz, rji2) & + !$omp reduction(+:ahi) & + !$omp reduction(-:ahj) + do k = 1_I8B, nplpl + i = k_plpl(1,k) + j = k_plpl(2,k) + dx = xh(1, j) - xh(1, i) + dy = xh(2, j) - xh(2, i) + dz = xh(3, j) - xh(3, i) + rji2 = dx**2 + dy**2 + dz**2 + if (lmask(i) .and. lmask(j)) call kick_getacch_int_one_pl(rji2, dx, dy, dz, Gmass(i), Gmass(j), ahi(1,i), ahi(2,i), ahi(3,i), ahj(1,j), ahj(2,j), ahj(3,j)) end do + !$omp end parallel do + ah(:,:) = ah(:,:) + ahi(:,:) + ahj(:,:) end associate return @@ -61,23 +63,71 @@ module pure subroutine kick_getacch_int_tp(self, GMpl, xhp, npl) integer(I4B), intent(in) :: npl !! Number of active massive bodies ! Internals integer(I4B) :: i, j - real(DP) :: rji2, irij3, fac, r2 - real(DP), dimension(NDIM) :: dx + !real(DP), dimension(:,:), allocatable :: aht if ((self%nbody == 0) .or. (npl == 0)) return associate(tp => self, ntp => self%nbody) - do concurrent(i = 1:ntp, tp%lmask(i)) - do j = 1, npl - dx(:) = tp%xh(:,i) - xhp(:, j) - r2 = dot_product(dx(:), dx(:)) - fac = GMpl(j) / (r2 * sqrt(r2)) - tp%ah(:, i) = tp%ah(:, i) - fac * dx(:) - end do + do i = 1, ntp + if (tp%lmask(i)) then + do j = 1, npl + block + real(DP) :: rji2 + real(DP) :: dx, dy, dz + dx = tp%xh(1, i) - xhp(1, j) + dy = tp%xh(2, i) - xhp(1, j) + dz = tp%xh(3, i) - xhp(1, j) + rji2 = dx**2 + dy**2 + dz**2 + call kick_getacch_int_one_tp(rji2, dx, dy, dz, GMpl(i), tp%ah(1,i), tp%ah(2,i), tp%ah(3,i)) + end block + end do + end if end do + !call move_alloc(aht, tp%ah) end associate return end subroutine kick_getacch_int_tp + module pure elemental subroutine kick_getacch_int_one_pl(rji2, dx, dy, dz, Gmi, Gmj, axi, ayi, azi, axj, ayj, azj) + implicit none + real(DP), intent(in) :: rji2 + real(DP), intent(in) :: dx, dy, dz + real(DP), intent(in) :: Gmi + real(DP), intent(in) :: Gmj + real(DP), intent(inout) :: axi, ayi, azi + real(DP), intent(inout) :: axj, ayj, azj + ! Internals + real(DP) :: faci, facj, irij3 + + irij3 = 1.0_DP / (rji2 * sqrt(rji2)) + faci = Gmi * irij3 + facj = Gmj * irij3 + axi = axi + facj * dx + ayi = ayi + facj * dy + azi = azi + facj * dz + axj = axj - faci * dx + ayj = ayj - faci * dy + azj = azj - faci * dz + return + end subroutine kick_getacch_int_one_pl + + + !module pure elemental subroutine kick_getacch_int_one_tp(rji2, dx, dy, dz, GMpl, ax, ay, az) + module pure subroutine kick_getacch_int_one_tp(rji2, dx, dy, dz, GMpl, ax, ay, az) + implicit none + real(DP), intent(in) :: rji2 + real(DP), intent(in) :: dx, dy, dz + real(DP), intent(in) :: GMpl + real(DP), intent(inout) :: ax, ay, az + ! Internals + real(DP) :: fac + + fac = GMpl / (rji2 * sqrt(rji2)) + ax = ax - fac * dx + ay = ay - fac * dy + az = az - fac * dz + return + end subroutine kick_getacch_int_one_tp + end submodule s_kick diff --git a/src/modules/swiftest_classes.f90 b/src/modules/swiftest_classes.f90 index bcacef235..3b6c2d591 100644 --- a/src/modules/swiftest_classes.f90 +++ b/src/modules/swiftest_classes.f90 @@ -246,6 +246,8 @@ module swiftest_classes integer(I4B), dimension(:), allocatable :: isperi !! Perihelion passage flag real(DP), dimension(:), allocatable :: peri !! Perihelion distance real(DP), dimension(:), allocatable :: atp !! Semimajor axis following perihelion passage + integer(I4B), dimension(:,:), allocatable :: k_pltp !! Index array used to convert flattened the body-body comparison upper triangular matrix + integer(I8B) :: npltp !! Number of pl-tp comparisons in the flattened upper triangular matrix !! Note to developers: If you add components to this class, be sure to update methods and subroutines that traverse the !! component list, such as setup_tp and util_spill_tp contains @@ -261,6 +263,7 @@ module swiftest_classes procedure :: vb2vh => util_coord_vb2vh_tp !! Convert test particles from barycentric to heliocentric coordinates (velocity only) procedure :: vh2vb => util_coord_vh2vb_tp !! Convert test particles from heliocentric to barycentric coordinates (velocity only) procedure :: xh2xb => util_coord_xh2xb_tp !! Convert test particles from heliocentric to barycentric coordinates (position only) + !procedure :: index => util_index_eucl_pltp !! Sets up the (i, j) -> k indexing used for the single-loop blocking Euclidean distance matrix procedure :: fill => util_fill_tp !! "Fills" bodies from one object into another depending on the results of a mask (uses the UNPACK intrinsic) procedure :: get_peri => util_peri_tp !! Determine system pericenter passages for test particles procedure :: resize => util_resize_tp !! Checks the current size of a Swiftest body against the requested size and resizes it if it is too small. @@ -458,6 +461,13 @@ module subroutine util_index_eucl_plpl(self, param) class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters end subroutine + module subroutine util_index_eucl_pltp(self, pl, param) + implicit none + class(swiftest_tp), intent(inout) :: self !! Swiftest test particle object + class(swiftest_pl), intent(in) :: pl !! Swiftest massive body object + class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters + end subroutine + module subroutine fragmentation_initialize(system, param, family, x, v, L_spin, Ip, mass, radius, & nfrag, Ip_frag, m_frag, rad_frag, xb_frag, vb_frag, rot_frag, Qloss, lfailure) implicit none @@ -701,7 +711,27 @@ 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 pure subroutine kick_getacch_int_pl(self) + module pure elemental subroutine kick_getacch_int_one_pl(rji2, dx, dy, dz, Gmi, Gmj, axi, ayi, azi, axj, ayj, azj) + implicit none + real(DP), intent(in) :: rji2 + real(DP), intent(in) :: dx, dy, dz + real(DP), intent(in) :: Gmi + real(DP), intent(in) :: Gmj + real(DP), intent(inout) :: axi, ayi, azi + real(DP), intent(inout) :: axj, ayj, azj + end subroutine kick_getacch_int_one_pl + + !module pure elemental subroutine kick_getacch_int_one_tp(rji2, dx, dy, dz, Gmpl, ax, ay, az) + module pure subroutine kick_getacch_int_one_tp(rji2, dx, dy, dz, Gmpl, ax, ay, az) + !$omp declare simd(kick_getacch_int_one_tp) + implicit none + real(DP), intent(in) :: rji2 + real(DP), intent(in) :: dx, dy, dz + real(DP), intent(in) :: Gmpl + real(DP), intent(inout) :: ax, ay, az + end subroutine kick_getacch_int_one_tp + + module subroutine kick_getacch_int_pl(self) implicit none class(swiftest_pl), intent(inout) :: self end subroutine kick_getacch_int_pl diff --git a/src/modules/symba_classes.f90 b/src/modules/symba_classes.f90 index 53997f033..25f9b5940 100644 --- a/src/modules/symba_classes.f90 +++ b/src/modules/symba_classes.f90 @@ -406,7 +406,7 @@ module subroutine symba_io_read_particle(system, param) class(symba_parameters), intent(inout) :: param !! Current run configuration parameters with SyMBA extensions end subroutine symba_io_read_particle - module pure subroutine symba_kick_getacch_int_pl(self) + module subroutine symba_kick_getacch_int_pl(self) implicit none class(symba_pl), intent(inout) :: self end subroutine symba_kick_getacch_int_pl diff --git a/src/symba/symba_encounter_check.f90 b/src/symba/symba_encounter_check.f90 index 96e13ebaf..42b67d079 100644 --- a/src/symba/symba_encounter_check.f90 +++ b/src/symba/symba_encounter_check.f90 @@ -16,24 +16,37 @@ module function symba_encounter_check_pl(self, system, dt, irec) result(lany_enc ! Result logical :: lany_encounter !! Returns true if there is at least one close encounter ! Internals - integer(I8B) :: k - integer(I4B) :: nenc - real(DP), dimension(NDIM) :: xr, vr - logical, dimension(:), allocatable :: lencounter, loc_lvdotr + integer(I8B) :: k, nplplm + integer(I4B) :: i, j, nenc + real(DP) :: xr, yr, zr, vxr, vyr, vzr, rhill1, rhill2 + logical, dimension(self%nplplm) :: lencounter, loc_lvdotr + real(DP), dimension(:,:), pointer :: xh, vh + real(DP), dimension(:), pointer :: rhill + integer(I4B), dimension(:,:), pointer :: k_plpl if (self%nbody == 0) return - associate(pl => self, npl => self%nbody, nplplm => self%nplplm) - allocate(lencounter(nplplm), loc_lvdotr(nplplm)) + associate(pl => self, xh => self%xh, vh => self%vh, rhill => self%rhill, npl => self%nbody, k_plpl => self%k_plpl) + nplplm = self%nplplm lencounter(:) = .false. - - do k = 1, nplplm - associate(i => pl%k_plpl(1, k), j => pl%k_plpl(2, k)) - xr(:) = pl%xh(:, j) - pl%xh(:, i) - vr(:) = pl%vh(:, j) - pl%vh(:, i) - call symba_encounter_check_one(xr(1), xr(2), xr(3), vr(1), vr(2), vr(3), pl%rhill(i), pl%rhill(j), dt, irec, lencounter(k), loc_lvdotr(k)) - end associate + loc_lvdotr(:) = .false. + + !$omp parallel do default(shared)& + !$omp private(k, i, j, xr, yr, zr, vxr, vyr, vzr, rhill1, rhill2) + do k = 1_I8B, nplplm + i = k_plpl(1, k) + j = k_plpl(2, k) + xr = xh(1, j) - xh(1, i) + yr = xh(2, j) - xh(2, i) + zr = xh(3, j) - xh(3, i) + vxr = vh(1, j) - vh(1, i) + vyr = vh(2, j) - vh(2, i) + vzr = vh(3, j) - vh(3, i) + rhill1 = rhill(i) + rhill2 = rhill(j) + call symba_encounter_check_one(xr, yr, zr, vxr, vyr, vzr, rhill1, rhill2, dt, irec, lencounter(k), loc_lvdotr(k)) end do + !$omp end parallel do nenc = count(lencounter(:)) lany_encounter = nenc > 0 @@ -46,18 +59,18 @@ module function symba_encounter_check_pl(self, system, dt, irec) result(lany_enc plplenc_list%id1(1:nenc) = pl%id(plplenc_list%index1(1:nenc)) plplenc_list%id2(1:nenc) = pl%id(plplenc_list%index2(1:nenc)) do k = 1, nenc - associate(i => plplenc_list%index1(k), j => plplenc_list%index2(k)) - plplenc_list%status(k) = ACTIVE - plplenc_list%level(k) = irec - pl%lencounter(i) = .true. - pl%lencounter(j) = .true. - pl%levelg(i) = irec - pl%levelm(i) = irec - pl%levelg(j) = irec - pl%levelm(j) = irec - pl%nplenc(i) = pl%nplenc(i) + 1 - pl%nplenc(j) = pl%nplenc(j) + 1 - end associate + i = plplenc_list%index1(k) + j = plplenc_list%index2(k) + plplenc_list%status(k) = ACTIVE + plplenc_list%level(k) = irec + pl%lencounter(i) = .true. + pl%lencounter(j) = .true. + pl%levelg(i) = irec + pl%levelm(i) = irec + pl%levelg(j) = irec + pl%levelm(j) = irec + pl%nplenc(i) = pl%nplenc(i) + 1 + pl%nplenc(j) = pl%nplenc(j) + 1 end do end associate end if diff --git a/src/symba/symba_kick.f90 b/src/symba/symba_kick.f90 index bf7adc60c..e03c12816 100644 --- a/src/symba/symba_kick.f90 +++ b/src/symba/symba_kick.f90 @@ -2,7 +2,7 @@ use swiftest contains - module pure subroutine symba_kick_getacch_int_pl(self) + module subroutine symba_kick_getacch_int_pl(self) !! author: David A. Minton !! !! Compute direct cross (third) term heliocentric accelerations of massive bodies, with no mutual interactions between bodies below GMTINY @@ -13,33 +13,35 @@ module pure subroutine symba_kick_getacch_int_pl(self) ! Arguments class(symba_pl), intent(inout) :: self ! Internals - integer(I4B) :: k - real(DP) :: rji2, irij3, faci, facj, rlim2 + integer(I8B) :: k, nplplm + real(DP) :: rji2, rlim2 real(DP) :: dx, dy, dz + integer(I4B) :: i, j + real(DP), dimension(:,:), pointer :: ah, xh + real(DP), dimension(NDIM,self%nbody) :: ahi, ahj + integer(I4B), dimension(:,:), pointer :: k_plpl + logical, dimension(:), pointer :: lmask + real(DP), dimension(:), pointer :: Gmass - associate(pl => self, npl => self%nbody, nplplm => self%nplplm) - do k = 1, nplplm - associate(i => pl%k_plpl(1, k), j => pl%k_plpl(2, k)) - if (pl%lmask(i) .and. pl%lmask(j)) then - dx = pl%xh(1, j) - pl%xh(1, i) - dy = pl%xh(2, j) - pl%xh(2, i) - dz = pl%xh(3, j) - pl%xh(3, i) - rji2 = dx**2 + dy**2 + dz**2 - rlim2 = (pl%radius(i) + pl%radius(j))**2 - if (rji2 > rlim2) then - irij3 = 1.0_DP / (rji2 * sqrt(rji2)) - faci = pl%Gmass(i) * irij3 - facj = pl%Gmass(j) * irij3 - pl%ah(1, i) = pl%ah(1, i) + facj * dx - pl%ah(2, i) = pl%ah(2, i) + facj * dy - pl%ah(3, i) = pl%ah(3, i) + facj * dz - pl%ah(1, j) = pl%ah(1, j) - faci * dx - pl%ah(2, j) = pl%ah(2, j) - faci * dy - pl%ah(3, j) = pl%ah(3, j) - faci * dz - end if - end if - end associate + associate(ah => self%ah, xh => self%xh, k_plpl => self%k_plpl, lmask => self%lmask, Gmass => self%Gmass) + nplplm = self%nplplm + ahi(:,:) = 0.0_DP + ahj(:,:) = 0.0_DP + !$omp parallel do default(shared)& + !$omp private(k, i, j, dx, dy, dz, rji2) & + !$omp reduction(+:ahi) & + !$omp reduction(-:ahj) + do k = 1_I8B, nplplm + i = k_plpl(1,k) + j = k_plpl(2,k) + dx = xh(1, j) - xh(1, i) + dy = xh(2, j) - xh(2, i) + dz = xh(3, j) - xh(3, i) + rji2 = dx**2 + dy**2 + dz**2 + if (lmask(i) .and. lmask(j)) call kick_getacch_int_one_pl(rji2, dx, dy, dz, Gmass(i), Gmass(j), ahi(1,i), ahi(2,i), ahi(3,i), ahj(1,j), ahj(2,j), ahj(3,j)) end do + !$omp end parallel do + ah(:,:) = ah(:,:) + ahi(:,:) + ahj(:,:) end associate return @@ -61,30 +63,36 @@ module subroutine symba_kick_getacch_pl(self, system, param, t, lbeg) real(DP), intent(in) :: t !! Current simulation time logical, intent(in) :: lbeg !! Logical flag that determines whether or not this is the beginning or end of the step ! Internals - integer(I4B) :: k - real(DP) :: irij3, rji2, rlim2, faci, facj - real(DP), dimension(NDIM) :: dx + integer(I4B) :: i, j + integer(I8B) :: k, nplplenc + real(DP) :: rji2, dx, dy, dz + real(DP), dimension(NDIM,self%nbody) :: ahi, ahj + class(symba_plplenc), pointer :: plplenc_list if (self%nbody == 0) return select type(system) class is (symba_nbody_system) - associate(pl => self, cb => system%cb, plplenc_list => system%plplenc_list, nplplenc => system%plplenc_list%nenc) + associate(pl => self, xh => self%xh, ah => self%ah, Gmass => self%Gmass, plplenc_list => system%plplenc_list) call helio_kick_getacch_pl(pl, system, param, t, lbeg) ! Remove accelerations from encountering pairs - do k = 1, nplplenc - associate(i => plplenc_list%index1(k), j => plplenc_list%index2(k)) - dx(:) = pl%xh(:, j) - pl%xh(:, i) - rji2 = dot_product(dx(:), dx(:)) - rlim2 = (pl%radius(i) + pl%radius(j))**2 - if (rji2 > rlim2) then - irij3 = 1.0_DP / (rji2 * sqrt(rji2)) - faci = pl%Gmass(i) * irij3 - facj = pl%Gmass(j) * irij3 - pl%ah(:, i) = pl%ah(:, i) - facj * dx(:) - pl%ah(:, j) = pl%ah(:, j) + faci * dx(:) - end if - end associate + nplplenc = int(plplenc_list%nenc, kind=I8B) + ahi(:,:) = 0.0_DP + ahj(:,:) = 0.0_DP + !$omp parallel do default(shared)& + !$omp private(k, i, j, dx, dy, dz, rji2) & + !$omp reduction(+:ahi) & + !$omp reduction(-:ahj) + do k = 1_I8B, nplplenc + i = plplenc_list%index1(k) + j = plplenc_list%index2(k) + dx = xh(1, j) - xh(1, i) + dy = xh(2, j) - xh(2, i) + dz = xh(3, j) - xh(3, i) + rji2 = dx**2 + dy**2 + dz**2 + call kick_getacch_int_one_pl(rji2, dx, dy, dz, Gmass(i), Gmass(j), ahi(1,i), ahi(2,i), ahi(3,i), ahj(1,j), ahj(2,j), ahj(3,j)) end do + !$omp end parallel do + ah(:,:) = ah(:,:) - ahi(:,:) - ahj(:,:) end associate end select diff --git a/src/util/util_index.f90 b/src/util/util_index.f90 index f852a492a..0e42ec7c7 100644 --- a/src/util/util_index.f90 +++ b/src/util/util_index.f90 @@ -5,7 +5,7 @@ module subroutine util_index_eucl_plpl(self, param) !! author: Jacob R. Elliott and David A. Minton !! - !! Turns i,j indices into k index for use in the Euclidean distance matrix + !! Turns i,j indices into k index for use in the Euclidean distance matrix for pl-pl interactions. !! !! Reference: !! @@ -23,7 +23,7 @@ module subroutine util_index_eucl_plpl(self, param) nplpl = (npl * (npl - 1) / 2) ! number of entries in a strict lower triangle, npl x npl, minus first column if (allocated(self%k_plpl)) deallocate(self%k_plpl) ! Reset the index array if it's been set previously allocate(self%k_plpl(2, nplpl)) - do i = 1, npl + do i = 1_I8B, npl counter = (i - 1_I8B) * npl - i * (i - 1_I8B) / 2_I8B + 1_I8B do j = i + 1_I8B, npl self%k_plpl(1, counter) = i @@ -36,4 +36,41 @@ module subroutine util_index_eucl_plpl(self, param) return end subroutine util_index_eucl_plpl + + module subroutine util_index_eucl_pltp(self, pl, param) + !! author: Jacob R. Elliott and David A. Minton + !! + !! Turns i,j indices into k index for use in the Euclidean distance matrix for pl-tp interactions + !! + !! Reference: + !! + !! Mélodie Angeletti, Jean-Marie Bonny, Jonas Koko. Parallel Euclidean distance matrix computation on big datasets *. + !! 2019. hal-0204751 + implicit none + ! Arguments + class(swiftest_tp), intent(inout) :: self !! Swiftest test particle object + class(swiftest_pl), intent(in) :: pl !! Swiftest massive body object + class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters + ! Internals + integer(I8B) :: i, j, counter, npl, ntp + + ntp = int(self%nbody, kind=I8B) + npl = int(pl%nbody, kind=I8B) + associate(npltp => self%npltp) + npltp = npl * ntp + if (allocated(self%k_pltp)) deallocate(self%k_pltp) ! Reset the index array if it's been set previously + allocate(self%k_pltp(2, npltp)) + do i = 1_I8B, npl + counter = (i - 1_I8B) * npl + 1_I8B + do j = 1_I8B, ntp + self%k_pltp(1, counter) = i + self%k_pltp(2, counter) = j + counter = counter + 1_I8B + end do + end do + end associate + + return + end subroutine util_index_eucl_pltp + end submodule s_util_index