diff --git a/src/discard/discard.f90 b/src/discard/discard.f90 index 1d9922207..cf38c2bbd 100644 --- a/src/discard/discard.f90 +++ b/src/discard/discard.f90 @@ -236,13 +236,12 @@ subroutine discard_pl_tp(tp, system, param) tp%status(i) = DISCARDED_PLR tp%lmask(i) = .false. pl%ldiscard(j) = .true. - write(*, *) "Particle ", tp%id(i), " too close to massive body ", pl%id(j), " at t = ", t write(idstri, *) tp%id(i) write(idstrj, *) pl%id(j) write(timestr, *) param%t - write(*, *) "Particle " // trim(adjustl(tp%info(i)%name)) // " (" // trim(adjustl(idstri)) // ")" & - // " too close to massive body " // trim(adjustl(pl%info(i)%name)) // " (" // trim(adjustl(idstrj)) & - // " at t = " // trim(adjustl(timestr)) + write(*, *) "Test particle " // trim(adjustl(tp%info(i)%name)) // " (" // trim(adjustl(idstri)) // ")" & + // " too close to massive body " // trim(adjustl(pl%info(i)%name)) // " (" // trim(adjustl(idstrj)) & + // " at t = " // trim(adjustl(timestr)) tp%ldiscard(i) = .true. exit end if diff --git a/src/fragmentation/fragmentation.f90 b/src/fragmentation/fragmentation.f90 index 43ed02599..c62d337ca 100644 --- a/src/fragmentation/fragmentation.f90 +++ b/src/fragmentation/fragmentation.f90 @@ -152,12 +152,12 @@ module subroutine fragmentation_initialize(system, param, family, x, v, L_spin, ! write(*, "(' -------------------------------------------------------------------------------------')") ! call calculate_system_energy(linclude_fragments=.true.) if (lfailure) then - ! write(*,*) "symba_frag_pos failed after: ",try," tries" + write(*,*) "symba_frag_pos failed after: ",try," tries" do ii = 1, nfrag vb_frag(:, ii) = vcom(:) end do - ! else - ! write(*,*) "symba_frag_pos succeeded after: ",try," tries" + else + write(*,*) "symba_frag_pos succeeded after: ",try," tries" ! write(*, "(' dL_tot should be very small' )") ! write(*,fmtlabel) ' dL_tot |', dLmag / Lmag_before ! write(*, "(' dE_tot should be negative and equal to Qloss' )") diff --git a/src/rmvs/rmvs_discard.f90 b/src/rmvs/rmvs_discard.f90 index bcdb9f902..724c107be 100644 --- a/src/rmvs/rmvs_discard.f90 +++ b/src/rmvs/rmvs_discard.f90 @@ -16,6 +16,7 @@ module subroutine rmvs_discard_tp(self, system, param) class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters ! Internals integer(I4B) :: i + character(len=STRMAX) :: timestr, idstri, idstrj if (self%nbody == 0) return @@ -25,7 +26,12 @@ module subroutine rmvs_discard_tp(self, system, param) if ((tp%status(i) == ACTIVE) .and. (tp%lperi(i))) then if ((tp%peri(i) < pl%radius(iplperP))) then tp%status(i) = DISCARDED_PLQ - write(*, *) "Particle ",tp%id(i)," q with respect to Planet ",pl%id(iplperP)," is too small at t = ",t + write(idstri, *) tp%id(i) + write(idstrj, *) pl%id(iplperP) + write(timestr, *) t + write(*, *) "Particle " // trim(adjustl(tp%info(i)%name)) // " (" // trim(adjustl(idstri)) & + // ") q with respect to massive body " // trim(adjustl(pl%info(iplperP)%name)) // " (" // trim(adjustl(idstrj)) & + // ") is too small at t = " // trim(adjustl(timestr)) tp%ldiscard(i) = .true. tp%lmask(i) = .false. end if diff --git a/src/rmvs/rmvs_util.f90 b/src/rmvs/rmvs_util.f90 index 140c43426..3b823c190 100644 --- a/src/rmvs/rmvs_util.f90 +++ b/src/rmvs/rmvs_util.f90 @@ -189,9 +189,11 @@ module subroutine rmvs_util_sort_pl(self, sortby, ascending) character(*), intent(in) :: sortby !! Sorting attribute logical, intent(in) :: ascending !! Logical flag indicating whether or not the sorting should be in ascending or descending order ! Internals - integer(I4B), dimension(self%nbody) :: ind + integer(I4B), dimension(:), allocatable :: ind integer(I4B) :: direction + if (self%nbody == 0) return + if (ascending) then direction = 1 else @@ -199,6 +201,7 @@ module subroutine rmvs_util_sort_pl(self, sortby, ascending) end if associate(pl => self, npl => self%nbody) + allocate(ind(npl)) select case(sortby) case("nenc") call util_sort(direction * pl%nenc(1:npl), ind(1:npl)) @@ -231,8 +234,10 @@ module subroutine rmvs_util_sort_tp(self, sortby, ascending) character(*), intent(in) :: sortby !! Sorting attribute logical, intent(in) :: ascending !! Logical flag indicating whether or not the sorting should be in ascending or descending order ! Internals - integer(I4B), dimension(self%nbody) :: ind - integer(I4B) :: direction + integer(I4B), dimension(:), allocatable :: ind + integer(I4B) :: direction + + if (self%nbody == 0) return if (ascending) then direction = 1 @@ -241,6 +246,7 @@ module subroutine rmvs_util_sort_tp(self, sortby, ascending) end if associate(tp => self, ntp => self%nbody) + allocate(ind(ntp)) select case(sortby) case("plperP") call util_sort(direction * tp%plperP(1:ntp), ind(1:ntp)) diff --git a/src/symba/symba_collision.f90 b/src/symba/symba_collision.f90 index 7a68a56cd..df215665f 100644 --- a/src/symba/symba_collision.f90 +++ b/src/symba/symba_collision.f90 @@ -740,7 +740,7 @@ module subroutine symba_collision_encounter_extract_collisions(self, system, par class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters ! Internals - logical, dimension(self%nenc) :: lplpl_collision + logical, dimension(:), allocatable :: lplpl_collision logical, dimension(:), allocatable :: lplpl_unique_parent integer(I4B), dimension(:), pointer :: plparent integer(I4B), dimension(:), allocatable :: collision_idx, unique_parent_idx @@ -749,6 +749,7 @@ module subroutine symba_collision_encounter_extract_collisions(self, system, par select type (pl => system%pl) class is (symba_pl) associate(plplenc_list => self, nplplenc => self%nenc, idx1 => self%index1, idx2 => self%index2, plparent => pl%kin%parent) + allocate(lplpl_collision(nplplenc)) lplpl_collision(:) = plplenc_list%status(1:nplplenc) == COLLISION if (.not.any(lplpl_collision)) return ! Collisions have been detected in this step. So we need to determine which of them are between unique bodies. @@ -1139,6 +1140,7 @@ module subroutine symba_collision_resolve_plplenc(self, system, param, t, dt, ir ! Internals real(DP) :: Eorbit_before, Eorbit_after logical :: lplpl_collision + character(len=STRMAX) :: timestr associate(plplenc_list => self, plplcollision_list => system%plplcollision_list) select type(pl => system%pl) @@ -1157,7 +1159,8 @@ module subroutine symba_collision_resolve_plplenc(self, system, param, t, dt, ir end if do - write(*, *) "Collision between massive bodies detected at time t = ", t + write(timestr,*) t + write(*, *) "Collision between massive bodies detected at time t = " // trim(adjustl(timestr)) if (param%lfragmentation) then call plplcollision_list%resolve_fragmentations(system, param) else diff --git a/src/symba/symba_discard.f90 b/src/symba/symba_discard.f90 index b1a433fca..f812dc664 100644 --- a/src/symba/symba_discard.f90 +++ b/src/symba/symba_discard.f90 @@ -36,14 +36,14 @@ subroutine symba_discard_cb_pl(pl, system, param) pl%status(i) = DISCARDED_RMAX write(idstr, *) pl%id(i) write(timestr, *) param%t - write(*, *) "Massive body " // trim(adjustl(pl%info(i)%name)) // " (" // trim(adjustl(idstr)) // ")" // " too far from the central body at t = " // trim(adjustl(timestr)) + write(*, *) trim(adjustl(pl%info(i)%name)) // " (" // trim(adjustl(idstr)) // ")" // " too far from the central body at t = " // trim(adjustl(timestr)) else if ((param%rmin >= 0.0_DP) .and. (rh2 < rmin2)) then pl%ldiscard(i) = .true. pl%lcollision(i) = .false. pl%status(i) = DISCARDED_RMIN write(idstr, *) pl%id(i) write(timestr, *) param%t - write(*, *) "Massive body " // trim(adjustl(pl%info(i)%name)) // " (" // trim(adjustl(idstr)) // ")" // " too close to the central body at t = " // trim(adjustl(timestr)) + write(*, *) trim(adjustl(pl%info(i)%name)) // " (" // trim(adjustl(idstr)) // ")" // " too close to the central body at t = " // trim(adjustl(timestr)) else if (param%rmaxu >= 0.0_DP) then rb2 = dot_product(pl%xb(:,i), pl%xb(:,i)) vb2 = dot_product(pl%vb(:,i), pl%vb(:,i)) @@ -54,7 +54,7 @@ subroutine symba_discard_cb_pl(pl, system, param) pl%status(i) = DISCARDED_RMAXU write(idstr, *) pl%id(i) write(timestr, *) param%t - write(*, *) "Massive body " // trim(adjustl(pl%info(i)%name)) // " (" // trim(adjustl(idstr)) // ")" // " is unbound and too far from barycenter at t = " // trim(adjustl(timestr)) + write(*, *) trim(adjustl(pl%info(i)%name)) // " (" // trim(adjustl(idstr)) // ")" // " is unbound and too far from barycenter at t = " // trim(adjustl(timestr)) end if end if end if @@ -255,6 +255,7 @@ subroutine symba_discard_peri_pl(pl, system, param) logical, save :: lfirst = .true. logical :: lfirst_orig integer(I4B) :: i + character(len=STRMAX) :: timestr, idstr lfirst_orig = pl%lfirst @@ -271,7 +272,9 @@ subroutine symba_discard_peri_pl(pl, system, param) pl%ldiscard(i) = .true. pl%lcollision(i) = .false. pl%status(i) = DISCARDED_PERI - write(*, *) "Particle ", pl%id(i), " perihelion distance too small at t = ", param%t + write(timestr, *) param%t + write(idstr, *) pl%id(i) + write(*, *) trim(adjustl(pl%info(i)%name)) // " (" // trim(adjustl(idstr)) // ") perihelion distance too small at t = " // trim(adjustl(timestr)) end if end if end if diff --git a/src/symba/symba_encounter_check.f90 b/src/symba/symba_encounter_check.f90 index 5aa0af850..78403c348 100644 --- a/src/symba/symba_encounter_check.f90 +++ b/src/symba/symba_encounter_check.f90 @@ -99,7 +99,7 @@ module function symba_encounter_check(self, system, dt, irec) result(lany_encoun integer(I4B), intent(in) :: irec !! Current recursion level logical :: lany_encounter !! Returns true if there is at least one close encounter ! Internals - integer(I4B) :: k + integer(I4B) :: i, j,k real(DP), dimension(NDIM) :: xr, vr logical :: lencounter, isplpl real(DP) :: rlim2, rji2 @@ -122,40 +122,40 @@ module function symba_encounter_check(self, system, dt, irec) result(lany_encoun allocate(lencmask(self%nenc)) lencmask(:) = (self%status(1:self%nenc) == ACTIVE) .and. (self%level(1:self%nenc) == irec - 1) if (.not.any(lencmask(:))) return - associate(ind1 => self%index1, ind2 => self%index2) - do concurrent(k = 1:self%nenc, lencmask(k)) + do concurrent(k = 1:self%nenc, lencmask(k)) + i = self%index1(k) + j = self%index2(k) + if (isplpl) then + xr(:) = pl%xh(:,j) - pl%xh(:,i) + vr(:) = pl%vb(:,j) - pl%vb(:,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, self%lvdotr(k)) + else + xr(:) = tp%xh(:,j) - pl%xh(:,i) + vr(:) = tp%vb(:,j) - pl%vb(:,i) + call symba_encounter_check_one(xr(1), xr(2), xr(3), vr(1), vr(2), vr(3), pl%rhill(i), 0.0_DP, dt, irec, lencounter, self%lvdotr(k)) + end if + if (lencounter) then if (isplpl) then - xr(:) = pl%xh(:,ind2(k)) - pl%xh(:,ind1(k)) - vr(:) = pl%vb(:,ind2(k)) - pl%vb(:,ind1(k)) - call symba_encounter_check_one(xr(1), xr(2), xr(3), vr(1), vr(2), vr(3), pl%rhill(ind1(k)), pl%rhill(ind2(k)), dt, irec, lencounter, self%lvdotr(k)) + rlim2 = (pl%radius(i) + pl%radius(j))**2 else - xr(:) = tp%xh(:,ind2(k)) - pl%xh(:,ind1(k)) - vr(:) = tp%vb(:,ind2(k)) - pl%vb(:,ind1(k)) - call symba_encounter_check_one(xr(1), xr(2), xr(3), vr(1), vr(2), vr(3), pl%rhill(ind1(k)), 0.0_DP, dt, irec, lencounter, self%lvdotr(k)) + rlim2 = (pl%radius(i))**2 end if - if (lencounter) then + rji2 = dot_product(xr(:), xr(:))! Check to see if these are physically overlapping bodies first, which we should ignore + if (rji2 > rlim2) then + lany_encounter = .true. + pl%levelg(i) = irec + pl%levelm(i) = MAX(irec, pl%levelm(i)) if (isplpl) then - rlim2 = (pl%radius(ind1(k)) + pl%radius(ind2(k)))**2 + pl%levelg(j) = irec + pl%levelm(j) = MAX(irec, pl%levelm(j)) else - rlim2 = (pl%radius(ind1(k)))**2 + tp%levelg(j) = irec + tp%levelm(j) = MAX(irec, tp%levelm(j)) end if - rji2 = dot_product(xr(:), xr(:))! Check to see if these are physically overlapping bodies first, which we should ignore - if (rji2 > rlim2) then - lany_encounter = .true. - pl%levelg(ind1(k)) = irec - pl%levelm(ind1(k)) = MAX(irec, pl%levelm(ind1(k))) - if (isplpl) then - pl%levelg(ind2(k)) = irec - pl%levelm(ind2(k)) = MAX(irec, pl%levelm(ind2(k))) - else - tp%levelg(ind2(k)) = irec - tp%levelm(ind2(k)) = MAX(irec, tp%levelm(ind2(k))) - end if - self%level(k) = irec - end if - end if - end do - end associate + self%level(k) = irec + end if + end if + end do end select end select diff --git a/src/symba/symba_kick.f90 b/src/symba/symba_kick.f90 index 9a5b11229..c80f53536 100644 --- a/src/symba/symba_kick.f90 +++ b/src/symba/symba_kick.f90 @@ -36,7 +36,7 @@ module subroutine symba_kick_getacch_pl(self, system, param, t, lbeg) ! Internals integer(I4B) :: i, j integer(I8B) :: k, nplplenc - real(DP) :: rji2, rlim2, xr, yr, zr + real(DP) :: rjj, rlim2, xr, yr, zr real(DP), dimension(NDIM,self%nbody) :: ah_enc integer(I4B), dimension(:,:), allocatable :: k_plpl_enc @@ -78,7 +78,7 @@ module subroutine symba_kick_getacch_tp(self, system, param, t, lbeg) 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) :: rji2, fac, rlim2 + real(DP) :: rjj, fac, rlim2 real(DP), dimension(NDIM) :: dx if (self%nbody == 0) return @@ -95,8 +95,8 @@ module subroutine symba_kick_getacch_tp(self, system, param, t, lbeg) else dx(:) = tp%xh(:,j) - pl%xend(:,i) end if - rji2 = dot_product(dx(:), dx(:)) - fac = pl%Gmass(i) / (rji2 * sqrt(rji2)) + rjj = dot_product(dx(:), dx(:)) + fac = pl%Gmass(i) / (rjj * sqrt(rjj)) tp%ah(:,j) = tp%ah(:,j) + fac * dx(:) end IF end associate @@ -123,15 +123,17 @@ module subroutine symba_kick_encounter(self, system, dt, irec, sgn) integer(I4B), intent(in) :: irec !! Current recursion level integer(I4B), intent(in) :: sgn !! sign to be applied to acceleration ! Internals - integer(I4B) :: i, irm1, irecl, ngood + integer(I4B) :: i, j, irm1, irecl, ngood + integer(I8B) :: k real(DP) :: r, rr, ri, ris, rim1, r2, ir3, fac, faci, facj real(DP), dimension(NDIM) :: dx logical :: isplpl - logical, dimension(self%nenc) :: lgoodlevel + logical, dimension(:), allocatable :: lgoodlevel integer(I4B), dimension(:), allocatable :: good_idx if (self%nenc == 0) return + select type(self) class is (symba_plplenc) isplpl = .true. @@ -142,9 +144,10 @@ module subroutine symba_kick_encounter(self, system, dt, irec, sgn) class is (symba_pl) select type(tp => system%tp) class is (symba_tp) - associate(ind1 => self%index1, ind2 => self%index2, npl => pl%nbody, ntp => tp%nbody, nenc => self%nenc) + associate(npl => pl%nbody, ntp => tp%nbody, nenc => self%nenc) if (npl > 0) pl%lmask(1:npl) = pl%status(1:npl) /= INACTIVE if (ntp > 0) tp%lmask(1:ntp) = tp%status(1:ntp) /= INACTIVE + allocate(lgoodlevel(nenc)) irm1 = irec - 1 @@ -154,89 +157,89 @@ module subroutine symba_kick_encounter(self, system, dt, irec, sgn) irecl = irec end if - do i = 1, nenc + do k = 1, nenc + i = self%index1(k) + j = self%index2(k) if (isplpl) then - lgoodlevel(i) = (pl%levelg(ind1(i)) >= irm1) .and. (pl%levelg(ind2(i)) >= irm1) + lgoodlevel(k) = (pl%levelg(i) >= irm1) .and. (pl%levelg(j) >= irm1) else - lgoodlevel(i) = (pl%levelg(ind1(i)) >= irm1) .and. (tp%levelg(ind2(i)) >= irm1) + lgoodlevel(k) = (pl%levelg(i) >= irm1) .and. (tp%levelg(j) >= irm1) end if - lgoodlevel(i) = (self%status(i) == ACTIVE) .and. lgoodlevel(i) + lgoodlevel(k) = (self%status(k) == ACTIVE) .and. lgoodlevel(k) end do ngood = count(lgoodlevel(:)) - if (ngood > 0) then + if (ngood > 0_I8B) then allocate(good_idx(ngood)) good_idx(:) = pack([(i, i = 1, nenc)], lgoodlevel(:)) if (isplpl) then - do i = 1, ngood - associate(i1 => ind1(good_idx(i)), i2 => ind2(good_idx(i))) - pl%ah(:,i1) = 0.0_DP - pl%ah(:,i2) = 0.0_DP - end associate + do concurrent (k = 1:ngood) + i = self%index1(good_idx(k)) + j = self%index2(good_idx(k)) + pl%ah(:,i) = 0.0_DP + pl%ah(:,j) = 0.0_DP end do else - do i = 1, ngood - associate(i2 => ind2(good_idx(i))) - tp%ah(:,i2) = 0.0_DP - end associate + do concurrent (k = 1_I8B:ngood) + j = self%index2(good_idx(k)) + tp%ah(:,j) = 0.0_DP end do end if - do i = 1, ngood - associate(k => good_idx(i), i1 => ind1(good_idx(i)), i2 => ind2(good_idx(i))) - if (isplpl) then - ri = ((pl%rhill(i1) + pl%rhill(i2))**2) * (RHSCALE**2) * (RSHELL**(2*irecl)) - rim1 = ri * (RSHELL**2) - dx(:) = pl%xh(:,i2) - pl%xh(:,i1) - else - ri = ((pl%rhill(i1))**2) * (RHSCALE**2) * (RSHELL**(2*irecl)) - rim1 = ri * (RSHELL**2) - dx(:) = tp%xh(:,i2) - pl%xh(:,i1) - end if - r2 = dot_product(dx(:), dx(:)) - if (r2 < rim1) then - fac = 0.0_DP - lgoodlevel(k) = .false. - cycle - end if - if (r2 < ri) then - ris = sqrt(ri) - r = sqrt(r2) - rr = (ris - r) / (ris * (1.0_DP - RSHELL)) - fac = (r2**(-1.5_DP)) * (1.0_DP - 3 * (rr**2) + 2 * (rr**3)) - else - ir3 = 1.0_DP / (r2 * sqrt(r2)) - fac = ir3 - end if - faci = fac * pl%Gmass(i1) - if (isplpl) then - facj = fac * pl%Gmass(i2) - pl%ah(:, i1) = pl%ah(:, i1) + facj * dx(:) - pl%ah(:, i2) = pl%ah(:, i2) - faci * dx(:) - else - tp%ah(:, i2) = tp%ah(:, i2) - faci * dx(:) - end if - end associate + do k = 1, ngood + i = self%index1(good_idx(k)) + j = self%index2(good_idx(k)) + if (isplpl) then + ri = ((pl%rhill(i) + pl%rhill(j))**2) * (RHSCALE**2) * (RSHELL**(2*irecl)) + rim1 = ri * (RSHELL**2) + dx(:) = pl%xh(:,j) - pl%xh(:,i) + else + ri = ((pl%rhill(i))**2) * (RHSCALE**2) * (RSHELL**(2*irecl)) + rim1 = ri * (RSHELL**2) + dx(:) = tp%xh(:,j) - pl%xh(:,i) + end if + r2 = dot_product(dx(:), dx(:)) + if (r2 < rim1) then + fac = 0.0_DP + lgoodlevel(good_idx(k)) = .false. + cycle + end if + if (r2 < ri) then + ris = sqrt(ri) + r = sqrt(r2) + rr = (ris - r) / (ris * (1.0_DP - RSHELL)) + fac = (r2**(-1.5_DP)) * (1.0_DP - 3 * (rr**2) + 2 * (rr**3)) + else + ir3 = 1.0_DP / (r2 * sqrt(r2)) + fac = ir3 + end if + faci = fac * pl%Gmass(i) + if (isplpl) then + facj = fac * pl%Gmass(j) + pl%ah(:, i) = pl%ah(:, i) + facj * dx(:) + pl%ah(:, j) = pl%ah(:, j) - faci * dx(:) + else + tp%ah(:, j) = tp%ah(:, j) - faci * dx(:) + end if end do ngood = count(lgoodlevel(:)) - if (ngood == 0) return + if (ngood == 0_I8B) return good_idx(1:ngood) = pack([(i, i = 1, nenc)], lgoodlevel(:)) if (isplpl) then - do i = 1, ngood - associate(i1 => ind1(good_idx(i)), i2 => ind2(good_idx(i))) - pl%vb(:,i1) = pl%vb(:,i1) + sgn * dt * pl%ah(:,i1) - pl%vb(:,i2) = pl%vb(:,i2) + sgn * dt * pl%ah(:,i2) - pl%ah(:,i1) = 0.0_DP - pl%ah(:,i2) = 0.0_DP - end associate + do k = 1, ngood + i = self%index1(good_idx(k)) + j = self%index2(good_idx(k)) + pl%vb(:,i) = pl%vb(:,i) + sgn * dt * pl%ah(:,i) + pl%vb(:,j) = pl%vb(:,j) + sgn * dt * pl%ah(:,j) + pl%ah(:,i) = 0.0_DP + pl%ah(:,j) = 0.0_DP end do else - do i = 1, ngood - associate(i2 => ind2(good_idx(i))) - tp%vb(:,i2) = tp%vb(:,i2) + sgn * dt * tp%ah(:,i2) - tp%ah(:,i2) = 0.0_DP - end associate + do k = 1, ngood + j = self%index2(good_idx(k)) + tp%vb(:,j) = tp%vb(:,j) + sgn * dt * tp%ah(:,j) + tp%ah(:,j) = 0.0_DP end do end if end if diff --git a/src/symba/symba_util.f90 b/src/symba/symba_util.f90 index 150bf7825..91a81091a 100644 --- a/src/symba/symba_util.f90 +++ b/src/symba/symba_util.f90 @@ -646,7 +646,7 @@ module subroutine symba_util_sort_pl(self, sortby, ascending) character(*), intent(in) :: sortby !! Sorting attribute logical, intent(in) :: ascending !! Logical flag indicating whether or not the sorting should be in ascending or descending order ! Internals - integer(I4B), dimension(self%nbody) :: ind + integer(I4B), dimension(:), allocatable :: ind integer(I4B) :: direction if (self%nbody == 0) return @@ -658,6 +658,7 @@ module subroutine symba_util_sort_pl(self, sortby, ascending) end if associate(pl => self, npl => self%nbody) + allocate(ind(npl)) select case(sortby) case("nplenc") call util_sort(direction * pl%nplenc(1:npl), ind(1:npl)) @@ -696,7 +697,7 @@ module subroutine symba_util_sort_tp(self, sortby, ascending) character(*), intent(in) :: sortby !! Sorting attribute logical, intent(in) :: ascending !! Logical flag indicating whether or not the sorting should be in ascending or descending order ! Internals - integer(I4B), dimension(self%nbody) :: ind + integer(I4B), dimension(:), allocatable :: ind integer(I4B) :: direction if (self%nbody == 0) return @@ -708,6 +709,7 @@ module subroutine symba_util_sort_tp(self, sortby, ascending) end if associate(tp => self, ntp => self%nbody) + allocate(ind(ntp)) select case(sortby) case("nplenc") call util_sort(direction * tp%nplenc(1:ntp), ind(1:ntp)) diff --git a/src/util/util_append.f90 b/src/util/util_append.f90 index 6abbc0e9a..ecefb1572 100644 --- a/src/util/util_append.f90 +++ b/src/util/util_append.f90 @@ -258,7 +258,6 @@ module subroutine util_append_pl(self, source, lsource_mask) class(swiftest_body), intent(in) :: source !! Source object to append logical, dimension(:), intent(in) :: lsource_mask !! Logical mask indicating which elements to append to - select type(source) class is (swiftest_pl) associate(nold => self%nbody, nsrc => source%nbody) @@ -276,6 +275,8 @@ module subroutine util_append_pl(self, source, lsource_mask) call util_append(self%Q, source%Q, nold, nsrc, lsource_mask) call util_append(self%tlag, source%tlag, nold, nsrc, lsource_mask) + if (allocated(self%k_plpl)) deallocate(self%k_plpl) + call util_append_body(self, source, lsource_mask) end associate class default diff --git a/src/util/util_fill.f90 b/src/util/util_fill.f90 index 541ee40cd..26b21310e 100644 --- a/src/util/util_fill.f90 +++ b/src/util/util_fill.f90 @@ -206,6 +206,8 @@ module subroutine util_fill_pl(self, inserts, lfill_list) call util_fill(keeps%vbeg, inserts%vbeg, lfill_list) call util_fill(keeps%Ip, inserts%Ip, lfill_list) call util_fill(keeps%rot, inserts%rot, lfill_list) + + if (allocated(keeps%k_plpl)) deallocate(keeps%k_plpl) call util_fill_body(keeps, inserts, lfill_list) class default diff --git a/src/util/util_get_energy_momentum.f90 b/src/util/util_get_energy_momentum.f90 index 4243070ac..70fbb16df 100644 --- a/src/util/util_get_energy_momentum.f90 +++ b/src/util/util_get_energy_momentum.f90 @@ -14,18 +14,19 @@ module subroutine util_get_energy_momentum_system(self, param) class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters ! Internals integer(I4B) :: i, j - integer(I8B) :: k + integer(I8B) :: k, nplpl real(DP) :: oblpot, kecb, kespincb real(DP), dimension(self%pl%nbody) :: irh, kepl, kespinpl, pecb real(DP), dimension(self%pl%nbody) :: Lplorbitx, Lplorbity, Lplorbitz real(DP), dimension(self%pl%nbody) :: Lplspinx, Lplspiny, Lplspinz - real(DP), dimension(self%pl%nplpl) :: pepl + logical, dimension(self%pl%nbody) :: lstatus real(DP), dimension(NDIM) :: Lcborbit, Lcbspin - logical, dimension(self%pl%nplpl) :: lstatpl - logical, dimension(self%pl%nbody) :: lstatus + real(DP), dimension(:), allocatable :: pepl + logical, dimension(:), allocatable :: lstatpl real(DP) :: hx, hy, hz associate(system => self, pl => self%pl, npl => self%pl%nbody, cb => self%cb) + nplpl = pl%nplpl system%Lorbit(:) = 0.0_DP system%Lspin(:) = 0.0_DP system%Ltot(:) = 0.0_DP @@ -82,24 +83,35 @@ module subroutine util_get_energy_momentum_system(self, param) end if ! Do the central body potential energy component first - associate(px => pl%xb(1,:), py => pl%xb(2,:), pz => pl%xb(3,:)) - do concurrent(i = 1:npl, lstatus(i)) - pecb(i) = -cb%Gmass * pl%mass(i) / sqrt(px(i)**2 + py(i)**2 + pz(i)**2) - end do - end associate + where(.not. lstatus(1:npl)) + pecb(1:npl) = 0.0_DP + end where + + do concurrent(i = 1:npl, lstatus(i)) + pecb(i) = -cb%Gmass * pl%mass(i) / norm2(pl%xb(:,i)) + end do ! Do the potential energy between pairs of massive bodies - associate(indi => pl%k_plpl(1, :), indj => pl%k_plpl(2, :)) - do concurrent (k = 1:pl%nplpl) - lstatpl(k) = (lstatus(indi(k)) .and. lstatus(indj(k))) - end do + allocate(lstatpl(nplpl)) + allocate(pepl(nplpl)) + do concurrent (k = 1:nplpl) + i = pl%k_plpl(1,k) + j = pl%k_plpl(2,k) + lstatpl(k) = (lstatus(i) .and. lstatus(j)) + end do - do concurrent (k = 1:pl%nplpl, lstatpl(k)) - pepl(k) = -(pl%Gmass(indi(k)) * pl%mass(indj(k))) / norm2(pl%xb(:, indi(k)) - pl%xb(:, indj(k))) - end do - end associate + where(.not.lstatpl(1:nplpl)) + pepl(1:nplpl) = 0.0_DP + end where + + do concurrent (k = 1:nplpl, lstatpl(k)) + i = pl%k_plpl(1,k) + j = pl%k_plpl(2,k) + pepl(k) = -(pl%Gmass(i) * pl%mass(j)) / norm2(pl%xb(:, i) - pl%xb(:, j)) + end do system%pe = sum(pepl(:), lstatpl(:)) + sum(pecb(1:npl), lstatus(1:npl)) + deallocate(lstatpl, pepl) system%ke_orbit = 0.5_DP * (kecb + sum(kepl(1:npl), lstatus(:))) if (param%lrotation) system%ke_spin = 0.5_DP * (kespincb + sum(kespinpl(1:npl), lstatus(:))) diff --git a/src/util/util_resize.f90 b/src/util/util_resize.f90 index ed28d1007..889703ac8 100644 --- a/src/util/util_resize.f90 +++ b/src/util/util_resize.f90 @@ -307,6 +307,8 @@ module subroutine util_resize_pl(self, nnew) call util_resize(self%Q, nnew) call util_resize(self%tlag, nnew) + if (allocated(self%k_plpl)) deallocate(self%k_plpl) + return end subroutine util_resize_pl diff --git a/src/util/util_sort.f90 b/src/util/util_sort.f90 index 9717f1e3e..b69627b0b 100644 --- a/src/util/util_sort.f90 +++ b/src/util/util_sort.f90 @@ -13,9 +13,11 @@ module subroutine util_sort_body(self, sortby, ascending) character(*), intent(in) :: sortby !! Sorting attribute logical, intent(in) :: ascending !! Logical flag indicating whether or not the sorting should be in ascending or descending order ! Internals - integer(I4B), dimension(self%nbody) :: ind + integer(I4B), dimension(:), allocatable :: ind integer(I4B) :: direction + if (self%nbody == 0) return + if (ascending) then direction = 1 else @@ -23,6 +25,7 @@ module subroutine util_sort_body(self, sortby, ascending) end if associate(body => self, n => self%nbody) + allocate(ind(n)) select case(sortby) case("id") call util_sort(direction * body%id(1:n), ind(1:n)) @@ -55,7 +58,6 @@ module subroutine util_sort_body(self, sortby, ascending) end subroutine util_sort_body - module subroutine util_sort_dp(arr) !! author: David A. Minton !! @@ -235,9 +237,11 @@ module subroutine util_sort_pl(self, sortby, ascending) character(*), intent(in) :: sortby !! Sorting attribute logical, intent(in) :: ascending !! Logical flag indicating whether or not the sorting should be in ascending or descending order ! Internals - integer(I4B), dimension(self%nbody) :: ind + integer(I4B), dimension(:), allocatable :: ind integer(I4B) :: direction + if (self%nbody == 0) return + if (ascending) then direction = 1 else @@ -245,6 +249,7 @@ module subroutine util_sort_pl(self, sortby, ascending) end if associate(pl => self, npl => self%nbody) + allocate(ind(npl)) select case(sortby) case("Gmass","mass") call util_sort(direction * pl%Gmass(1:npl), ind(1:npl)) @@ -286,9 +291,11 @@ module subroutine util_sort_tp(self, sortby, ascending) character(*), intent(in) :: sortby !! Sorting attribute logical, intent(in) :: ascending !! Logical flag indicating whether or not the sorting should be in ascending or descending order ! Internals - integer(I4B), dimension(self%nbody) :: ind + integer(I4B), dimension(:), allocatable :: ind integer(I4B) :: direction + if (self%nbody == 0) return + if (ascending) then direction = 1 else @@ -296,6 +303,7 @@ module subroutine util_sort_tp(self, sortby, ascending) end if associate(tp => self, ntp => self%nbody) + allocate(ind(ntp)) select case(sortby) case("peri") call util_sort(direction * tp%peri(1:ntp), ind(1:ntp)) @@ -506,6 +514,8 @@ module subroutine util_sort_rearrange_pl(self, ind) call util_sort_rearrange(pl%Q, ind, npl) call util_sort_rearrange(pl%tlag, ind, npl) + if (allocated(pl%k_plpl)) deallocate(pl%k_plpl) + call util_sort_rearrange_body(pl, ind) end associate diff --git a/src/util/util_spill.f90 b/src/util/util_spill.f90 index b4d43e019..8f9567ef9 100644 --- a/src/util/util_spill.f90 +++ b/src/util/util_spill.f90 @@ -417,6 +417,8 @@ module subroutine util_spill_pl(self, discards, lspill_list, ldestructive) call util_spill(keeps%Ip, discards%Ip, lspill_list, ldestructive) call util_spill(keeps%rot, discards%rot, lspill_list, ldestructive) + if (ldestructive .and. allocated(keeps%k_plpl)) deallocate(keeps%k_plpl) + call util_spill_body(keeps, discards, lspill_list, ldestructive) class default write(*,*) 'Error! spill method called for incompatible return type on swiftest_pl' diff --git a/src/whm/whm_util.f90 b/src/whm/whm_util.f90 index 6689e78cb..bedf359ae 100644 --- a/src/whm/whm_util.f90 +++ b/src/whm/whm_util.f90 @@ -141,9 +141,11 @@ module subroutine whm_util_sort_pl(self, sortby, ascending) character(*), intent(in) :: sortby !! Sorting attribute logical, intent(in) :: ascending !! Logical flag indicating whether or not the sorting should be in ascending or descending order ! Internals - integer(I4B), dimension(self%nbody) :: ind + integer(I4B), dimension(:), allocatable :: ind integer(I4B) :: direction + if (self%nbody == 0) return + if (ascending) then direction = 1 else @@ -151,6 +153,7 @@ module subroutine whm_util_sort_pl(self, sortby, ascending) end if associate(pl => self, npl => self%nbody) + allocate(ind(npl)) select case(sortby) case("eta") call util_sort(direction * pl%eta(1:npl), ind(1:npl))