diff --git a/src/fragmentation/fragmentation.f90 b/src/fragmentation/fragmentation.f90 index 1977d42f3..f4ff3b4ce 100644 --- a/src/fragmentation/fragmentation.f90 +++ b/src/fragmentation/fragmentation.f90 @@ -49,7 +49,6 @@ module subroutine fragmentation_initialize(system, param, family, x, v, L_spin, class(swiftest_nbody_system), allocatable :: tmpsys class(swiftest_parameters), allocatable :: tmpparam - if (nfrag < NFRAG_MIN) then write(*,*) "symba_frag_pos needs at least ",NFRAG_MIN," fragments, but only ",nfrag," were given." lfailure = .true. @@ -114,10 +113,9 @@ module subroutine fragmentation_initialize(system, param, family, x, v, L_spin, if (lfailure) write(*,*) 'Failed to find radial velocities' if (.not.lfailure) then call calculate_system_energy(linclude_fragments=.true.) - - write(*,*) 'Qloss : ',Qloss - write(*,*) '-dEtot: ',-dEtot - write(*,*) 'delta : ',abs((dEtot + Qloss)) + ! write(*,*) 'Qloss : ',Qloss + ! write(*,*) '-dEtot: ',-dEtot + ! write(*,*) 'delta : ',abs((dEtot + Qloss)) if ((abs(dEtot + Qloss) > Etol) .or. (dEtot > 0.0_DP)) then write(*,*) 'Failed due to high energy error: ',dEtot, abs(dEtot + Qloss) / Etol lfailure = .true. @@ -132,6 +130,10 @@ module subroutine fragmentation_initialize(system, param, family, x, v, L_spin, call restructure_failed_fragments() try = try + 1 end do + call restore_scale_factors() + call calculate_system_energy(linclude_fragments=.true.) + + write(*, "(' -------------------------------------------------------------------------------------')") write(*, "(' Final diagnostic')") write(*, "(' -------------------------------------------------------------------------------------')") @@ -151,7 +153,6 @@ module subroutine fragmentation_initialize(system, param, family, x, v, L_spin, end if write(*, "(' -------------------------------------------------------------------------------------')") - call restore_scale_factors() call ieee_set_halting_mode(IEEE_ALL,fpe_halting_modes) ! Save the current halting modes so we can turn them off temporarily return @@ -248,6 +249,11 @@ subroutine restore_scale_factors() Ltot_after = Ltot_after * Lscale Lmag_after = Lmag_after * Lscale + dLmag = norm2(Ltot_after(:) - Ltot_before(:)) + dEtot = Etot_after - Etot_before + + call tmpsys%rescale(tmpparam, mscale**(-1), dscale**(-1), tscale**(-1)) + mscale = 1.0_DP dscale = 1.0_DP vscale = 1.0_DP diff --git a/src/modules/rmvs_classes.f90 b/src/modules/rmvs_classes.f90 index 4f7255237..6ffb7ba1b 100644 --- a/src/modules/rmvs_classes.f90 +++ b/src/modules/rmvs_classes.f90 @@ -163,7 +163,7 @@ module subroutine rmvs_util_append_pl(self, source, lsource_mask) implicit none class(rmvs_pl), intent(inout) :: self !! RMVS massive body object class(swiftest_body), intent(in) :: source !! Source object to append - logical, dimension(:), optional, intent(in) :: lsource_mask !! Logical mask indicating which elements to append to + logical, dimension(:), intent(in) :: lsource_mask !! Logical mask indicating which elements to append to end subroutine rmvs_util_append_pl module subroutine rmvs_util_append_tp(self, source, lsource_mask) @@ -171,7 +171,7 @@ module subroutine rmvs_util_append_tp(self, source, lsource_mask) implicit none class(rmvs_tp), intent(inout) :: self !! RMVS test particle object class(swiftest_body), intent(in) :: source !! Source object to append - logical, dimension(:), optional, intent(in) :: lsource_mask !! Logical mask indicating which elements to append to + logical, dimension(:), intent(in) :: lsource_mask !! Logical mask indicating which elements to append to end subroutine rmvs_util_append_tp module subroutine rmvs_util_fill_pl(self, inserts, lfill_list) diff --git a/src/modules/swiftest_classes.f90 b/src/modules/swiftest_classes.f90 index 4d0e98704..0d7fac843 100644 --- a/src/modules/swiftest_classes.f90 +++ b/src/modules/swiftest_classes.f90 @@ -836,35 +836,35 @@ module subroutine util_append_arr_char_string(arr, source, lsource_mask) implicit none character(len=STRMAX), dimension(:), allocatable, intent(inout) :: arr !! Destination array character(len=STRMAX), dimension(:), allocatable, intent(in) :: source !! Array to append - logical, dimension(:), optional, intent(in) :: lsource_mask !! Logical mask indicating which elements to append to + logical, dimension(:), intent(in) :: lsource_mask !! Logical mask indicating which elements to append to end subroutine util_append_arr_char_string module subroutine util_append_arr_DP(arr, source, lsource_mask) implicit none real(DP), dimension(:), allocatable, intent(inout) :: arr !! Destination array real(DP), dimension(:), allocatable, intent(in) :: source !! Array to append - logical, dimension(:), optional, intent(in) :: lsource_mask !! Logical mask indicating which elements to append to + logical, dimension(:), intent(in) :: lsource_mask !! Logical mask indicating which elements to append to end subroutine util_append_arr_DP module subroutine util_append_arr_DPvec(arr, source, lsource_mask) implicit none real(DP), dimension(:,:), allocatable, intent(inout) :: arr !! Destination array real(DP), dimension(:,:), allocatable, intent(in) :: source !! Array to append - logical, dimension(:), optional, intent(in) :: lsource_mask !! Logical mask indicating which elements to append to + logical, dimension(:), intent(in) :: lsource_mask !! Logical mask indicating which elements to append to end subroutine util_append_arr_DPvec module subroutine util_append_arr_I4B(arr, source, lsource_mask) implicit none integer(I4B), dimension(:), allocatable, intent(inout) :: arr !! Destination array integer(I4B), dimension(:), allocatable, intent(in) :: source !! Array to append - logical, dimension(:), optional, intent(in) :: lsource_mask !! Logical mask indicating which elements to append to + logical, dimension(:), intent(in) :: lsource_mask !! Logical mask indicating which elements to append to end subroutine util_append_arr_I4B module subroutine util_append_arr_logical(arr, source, lsource_mask) implicit none logical, dimension(:), allocatable, intent(inout) :: arr !! Destination array logical, dimension(:), allocatable, intent(in) :: source !! Array to append - logical, dimension(:), optional, intent(in) :: lsource_mask !! Logical mask indicating which elements to append to + logical, dimension(:), intent(in) :: lsource_mask !! Logical mask indicating which elements to append to end subroutine util_append_arr_logical end interface @@ -873,21 +873,21 @@ module subroutine util_append_body(self, source, lsource_mask) implicit none class(swiftest_body), intent(inout) :: self !! Swiftest body object class(swiftest_body), intent(in) :: source !! Source object to append - logical, dimension(:), optional, intent(in) :: lsource_mask !! Logical mask indicating which elements to append to + logical, dimension(:), intent(in) :: lsource_mask !! Logical mask indicating which elements to append to end subroutine util_append_body module subroutine util_append_pl(self, source, lsource_mask) implicit none class(swiftest_pl), intent(inout) :: self !! Swiftest massive body object class(swiftest_body), intent(in) :: source !! Source object to append - logical, dimension(:), optional, intent(in) :: lsource_mask !! Logical mask indicating which elements to append to + logical, dimension(:), intent(in) :: lsource_mask !! Logical mask indicating which elements to append to end subroutine util_append_pl module subroutine util_append_tp(self, source, lsource_mask) implicit none class(swiftest_tp), intent(inout) :: self !! Swiftest test particle object class(swiftest_body), intent(in) :: source !! Source object to append - logical, dimension(:), optional, intent(in) :: lsource_mask !! Logical mask indicating which elements to append to + logical, dimension(:), intent(in) :: lsource_mask !! Logical mask indicating which elements to append to end subroutine util_append_tp module subroutine util_coord_b2h_pl(self, cb) diff --git a/src/modules/symba_classes.f90 b/src/modules/symba_classes.f90 index dfe1c0326..e8b5918b3 100644 --- a/src/modules/symba_classes.f90 +++ b/src/modules/symba_classes.f90 @@ -492,14 +492,14 @@ module subroutine symba_util_append_arr_info(arr, source, lsource_mask) implicit none type(symba_particle_info), dimension(:), allocatable, intent(inout) :: arr !! Destination array type(symba_particle_info), dimension(:), allocatable, intent(in) :: source !! Array to append - logical, dimension(:), optional, intent(in) :: lsource_mask !! Logical mask indicating which elements to append to + logical, dimension(:), intent(in) :: lsource_mask !! Logical mask indicating which elements to append to end subroutine symba_util_append_arr_info module subroutine symba_util_append_arr_kin(arr, source, lsource_mask) implicit none type(symba_kinship), dimension(:), allocatable, intent(inout) :: arr !! Destination array type(symba_kinship), dimension(:), allocatable, intent(in) :: source !! Array to append - logical, dimension(:), optional, intent(in) :: lsource_mask !! Logical mask indicating which elements to append to + logical, dimension(:), intent(in) :: lsource_mask !! Logical mask indicating which elements to append to end subroutine symba_util_append_arr_kin end interface @@ -509,7 +509,7 @@ module subroutine symba_util_append_merger(self, source, lsource_mask) implicit none class(symba_merger), intent(inout) :: self !! SyMBA massive body object class(swiftest_body), intent(in) :: source !! Source object to append - logical, dimension(:), optional, intent(in) :: lsource_mask !! Logical mask indicating which elements to append to + logical, dimension(:), intent(in) :: lsource_mask !! Logical mask indicating which elements to append to end subroutine symba_util_append_merger module subroutine symba_util_append_pl(self, source, lsource_mask) @@ -517,7 +517,7 @@ module subroutine symba_util_append_pl(self, source, lsource_mask) implicit none class(symba_pl), intent(inout) :: self !! SyMBA massive body object class(swiftest_body), intent(in) :: source !! Source object to append - logical, dimension(:), optional, intent(in) :: lsource_mask !! Logical mask indicating which elements to append to + logical, dimension(:), intent(in) :: lsource_mask !! Logical mask indicating which elements to append to end subroutine symba_util_append_pl module subroutine symba_util_append_tp(self, source, lsource_mask) @@ -525,7 +525,7 @@ module subroutine symba_util_append_tp(self, source, lsource_mask) implicit none class(symba_tp), intent(inout) :: self !! SyMBA test particle object class(swiftest_body), intent(in) :: source !! Source object to append - logical, dimension(:), optional, intent(in) :: lsource_mask !! Logical mask indicating which elements to append to + logical, dimension(:), intent(in) :: lsource_mask !! Logical mask indicating which elements to append to end subroutine symba_util_append_tp end interface diff --git a/src/modules/whm_classes.f90 b/src/modules/whm_classes.f90 index a79f52bca..6d5b26394 100644 --- a/src/modules/whm_classes.f90 +++ b/src/modules/whm_classes.f90 @@ -232,7 +232,7 @@ module subroutine whm_util_append_pl(self, source, lsource_mask) implicit none class(whm_pl), intent(inout) :: self !! WHM massive body object class(swiftest_body), intent(in) :: source !! Source object to append - logical, dimension(:), optional, intent(in) :: lsource_mask !! Logical mask indicating which elements to append to + logical, dimension(:), intent(in) :: lsource_mask !! Logical mask indicating which elements to append to end subroutine whm_util_append_pl module subroutine whm_util_spill_pl(self, discards, lspill_list, ldestructive) diff --git a/src/rmvs/rmvs_util.f90 b/src/rmvs/rmvs_util.f90 index 9f9cf0037..67a76acb3 100644 --- a/src/rmvs/rmvs_util.f90 +++ b/src/rmvs/rmvs_util.f90 @@ -11,7 +11,7 @@ module subroutine rmvs_util_append_pl(self, source, lsource_mask) !! Arguments class(rmvs_pl), intent(inout) :: self !! RMVS massive body object class(swiftest_body), intent(in) :: source !! Source object to append - logical, dimension(:), optional, intent(in) :: lsource_mask !! Logical mask indicating which elements to append to + logical, dimension(:), intent(in) :: lsource_mask !! Logical mask indicating which elements to append to select type(source) class is (rmvs_pl) @@ -44,7 +44,7 @@ module subroutine rmvs_util_append_tp(self, source, lsource_mask) !! Arguments class(rmvs_tp), intent(inout) :: self !! RMVS test particle object class(swiftest_body), intent(in) :: source !! Source object to append - logical, dimension(:), optional, intent(in) :: lsource_mask !! Logical mask indicating which elements to append to + logical, dimension(:), intent(in) :: lsource_mask !! Logical mask indicating which elements to append to select type(source) class is (rmvs_tp) diff --git a/src/symba/symba_fragmentation.f90 b/src/symba/symba_fragmentation.f90 index a27d31e12..9b82e145a 100644 --- a/src/symba/symba_fragmentation.f90 +++ b/src/symba/symba_fragmentation.f90 @@ -133,11 +133,13 @@ module function symba_fragmentation_casedisruption(system, param, family, x, v, plnew%k2 = pl%k2(ibiggest) plnew%tlag = pl%tlag(ibiggest) end if + call plnew%set_mu(cb) + pl%lmtiny(:) = pl%Gmass(:) > param%GMTINY ! Append the new merged body to the list and record how many we made nstart = pl_adds%nbody + 1 nend = pl_adds%nbody + plnew%nbody - call pl_adds%append(plnew) + call pl_adds%append(plnew, lsource_mask=[(.true., i=1, plnew%nbody)]) pl_adds%ncomp(nstart:nend) = plnew%nbody call plnew%setup(0, param) @@ -291,11 +293,13 @@ module function symba_fragmentation_casehitandrun(system, param, family, x, v, m plnew%k2 = pl%k2(ibiggest) plnew%tlag = pl%tlag(ibiggest) end if + call plnew%set_mu(cb) + pl%lmtiny(:) = pl%Gmass(:) > param%GMTINY ! Append the new merged body to the list and record how many we made nstart = pl_adds%nbody + 1 nend = pl_adds%nbody + plnew%nbody - call pl_adds%append(plnew) + call pl_adds%append(plnew, lsource_mask=[(.true., i=1, plnew%nbody)]) pl_adds%ncomp(nstart:nend) = plnew%nbody call plnew%setup(0, param) @@ -427,11 +431,13 @@ module function symba_fragmentation_casemerge(system, param, family, x, v, mass, plnew%k2 = pl%k2(ibiggest) plnew%tlag = pl%tlag(ibiggest) end if + call plnew%set_mu(cb) + pl%lmtiny(:) = pl%Gmass(:) > param%GMTINY ! Append the new merged body to the list and record how many we made nstart = pl_adds%nbody + 1 nend = pl_adds%nbody + plnew%nbody - call pl_adds%append(plnew) + call pl_adds%append(plnew, lsource_mask=[(.true., i=1, plnew%nbody)]) pl_adds%ncomp(nstart:nend) = plnew%nbody call plnew%setup(0, param) @@ -570,11 +576,13 @@ module function symba_fragmentation_casesupercatastrophic(system, param, family, plnew%k2 = pl%k2(ibiggest) plnew%tlag = pl%tlag(ibiggest) end if + call plnew%set_mu(cb) + pl%lmtiny(:) = pl%Gmass(:) > param%GMTINY ! Append the new merged body to the list and record how many we made nstart = pl_adds%nbody + 1 nend = pl_adds%nbody + plnew%nbody - call pl_adds%append(plnew) + call pl_adds%append(plnew, lsource_mask=[(.true., i=1, plnew%nbody)]) pl_adds%ncomp(nstart:nend) = plnew%nbody call plnew%setup(0, param) diff --git a/src/symba/symba_util.f90 b/src/symba/symba_util.f90 index 028b0678c..37fdca873 100644 --- a/src/symba/symba_util.f90 +++ b/src/symba/symba_util.f90 @@ -10,17 +10,13 @@ module subroutine symba_util_append_arr_info(arr, source, lsource_mask) ! Arguments type(symba_particle_info), dimension(:), allocatable, intent(inout) :: arr !! Destination array type(symba_particle_info), dimension(:), allocatable, intent(in) :: source !! Array to append - logical, dimension(:), optional, intent(in) :: lsource_mask !! Logical mask indicating which elements to append to + logical, dimension(:), intent(in) :: lsource_mask !! Logical mask indicating which elements to append to ! Internals integer(I4B) :: narr, nsrc if (.not. allocated(source)) return - if (present(lsource_mask)) then - nsrc = count(lsource_mask) - else - nsrc = size(source) - end if + nsrc = count(lsource_mask) if (allocated(arr)) then narr = size(arr) @@ -31,11 +27,7 @@ module subroutine symba_util_append_arr_info(arr, source, lsource_mask) call util_resize(arr, narr + nsrc) - if (present(lsource_mask)) then - arr(narr + 1:narr + nsrc) = pack(source(:), lsource_mask(:)) - else - arr(narr + 1:narr + nsrc) = source(:) - end if + arr(narr + 1:narr + nsrc) = pack(source(:), lsource_mask(:)) return end subroutine symba_util_append_arr_info @@ -49,17 +41,13 @@ module subroutine symba_util_append_arr_kin(arr, source, lsource_mask) ! Arguments type(symba_kinship), dimension(:), allocatable, intent(inout) :: arr !! Destination array type(symba_kinship), dimension(:), allocatable, intent(in) :: source !! Array to append - logical, dimension(:), optional, intent(in) :: lsource_mask !! Logical mask indicating which elements to append to + logical, dimension(:), intent(in) :: lsource_mask !! Logical mask indicating which elements to append to ! Internals integer(I4B) :: narr, nsrc if (.not. allocated(source)) return - if (present(lsource_mask)) then - nsrc = count(lsource_mask) - else - nsrc = size(source) - end if + nsrc = count(lsource_mask) if (allocated(arr)) then narr = size(arr) @@ -70,11 +58,7 @@ module subroutine symba_util_append_arr_kin(arr, source, lsource_mask) call util_resize(arr, narr + nsrc) - if (present(lsource_mask)) then - arr(narr + 1:narr + nsrc) = pack(source(:), lsource_mask(:)) - else - arr(narr + 1:narr + nsrc) = source(:) - end if + arr(narr + 1:narr + nsrc) = pack(source(:), lsource_mask(:)) return end subroutine symba_util_append_arr_kin @@ -89,7 +73,7 @@ module subroutine symba_util_append_pl(self, source, lsource_mask) !! Arguments class(symba_pl), intent(inout) :: self !! SyMBA massive body object class(swiftest_body), intent(in) :: source !! Source object to append - logical, dimension(:), optional, intent(in) :: lsource_mask !! Logical mask indicating which elements to append to + logical, dimension(:), intent(in) :: lsource_mask !! Logical mask indicating which elements to append to select type(source) class is (symba_pl) @@ -125,7 +109,7 @@ module subroutine symba_util_append_merger(self, source, lsource_mask) ! Arguments class(symba_merger), intent(inout) :: self !! SyMBA massive body object class(swiftest_body), intent(in) :: source !! Source object to append - logical, dimension(:), optional, intent(in) :: lsource_mask !! Logical mask indicating which elements to append to + logical, dimension(:), intent(in) :: lsource_mask !! Logical mask indicating which elements to append to ! Internals integer(I4B), dimension(:), allocatable :: ncomp_tmp !! Temporary placeholder for ncomp incase we are appending a symba_pl object to a symba_merger @@ -156,7 +140,7 @@ module subroutine symba_util_append_tp(self, source, lsource_mask) !! Arguments class(symba_tp), intent(inout) :: self !! SyMBA test particle object class(swiftest_body), intent(in) :: source !! Source object to append - logical, dimension(:), optional, intent(in) :: lsource_mask !! Logical mask indicating which elements to append to + logical, dimension(:), intent(in) :: lsource_mask !! Logical mask indicating which elements to append to select type(source) class is (symba_tp) @@ -380,11 +364,15 @@ module subroutine symba_util_rearray_pl(self, system, param) class(symba_parameters), intent(in) :: param !! Current run configuration parameters ! Internals class(symba_pl), allocatable :: tmp !! The discarded body list. + integer(I4B) :: i + logical, dimension(:), allocatable :: lmask associate(pl => self, pl_adds => system%pl_adds) allocate(tmp, mold=pl) ! Remove the discards and destroy the list, as the system already tracks pl_discards elsewhere - call pl%spill(tmp, lspill_list=(pl%ldiscard(:) .or. pl%status(:) == INACTIVE), ldestructive=.true.) + allocate(lmask, source=pl%ldiscard(:)) + lmask(:) = lmask(:) .or. pl%status(:) == INACTIVE + call pl%spill(tmp, lspill_list=lmask, ldestructive=.true.) call tmp%setup(0,param) deallocate(tmp) @@ -393,7 +381,7 @@ module subroutine symba_util_rearray_pl(self, system, param) if (allocated(pl%xend)) deallocate(pl%xend) ! Add in any new bodies - call pl%append(pl_adds) + call pl%append(pl_adds, lsource_mask=[(.true., i=1, pl_adds%nbody)]) ! If there are still bodies in the system, sort by mass in descending order and re-index if (pl%nbody > 0) then diff --git a/src/util/util_append.f90 b/src/util/util_append.f90 index 0f7ac0bde..cf0bb4117 100644 --- a/src/util/util_append.f90 +++ b/src/util/util_append.f90 @@ -10,17 +10,13 @@ module subroutine util_append_arr_char_string(arr, source, lsource_mask) ! Arguments character(len=STRMAX), dimension(:), allocatable, intent(inout) :: arr !! Destination array character(len=STRMAX), dimension(:), allocatable, intent(in) :: source !! Array to append - logical, dimension(:), optional, intent(in) :: lsource_mask !! Logical mask indicating which elements to append to + logical, dimension(:), intent(in) :: lsource_mask !! Logical mask indicating which elements to append to ! Internals integer(I4B) :: narr, nsrc if (.not. allocated(source)) return - if (present(lsource_mask)) then - nsrc = count(lsource_mask) - else - nsrc = size(source) - end if + nsrc = count(lsource_mask) if (allocated(arr)) then narr = size(arr) @@ -31,11 +27,7 @@ module subroutine util_append_arr_char_string(arr, source, lsource_mask) call util_resize(arr, narr + nsrc) - if (present(lsource_mask)) then - arr(narr + 1:narr + nsrc) = pack(source(:), lsource_mask(:)) - else - arr(narr + 1:narr + nsrc) = source(:) - end if + arr(narr + 1:narr + nsrc) = pack(source(:), lsource_mask(:)) return end subroutine util_append_arr_char_string @@ -49,17 +41,13 @@ module subroutine util_append_arr_DP(arr, source, lsource_mask) ! Arguments real(DP), dimension(:), allocatable, intent(inout) :: arr !! Destination array real(DP), dimension(:), allocatable, intent(in) :: source !! Array to append - logical, dimension(:), optional, intent(in) :: lsource_mask !! Logical mask indicating which elements to append to + logical, dimension(:), intent(in) :: lsource_mask !! Logical mask indicating which elements to append to ! Internals integer(I4B) :: narr, nsrc if (.not. allocated(source)) return - if (present(lsource_mask)) then - nsrc = count(lsource_mask) - else - nsrc = size(source) - end if + nsrc = count(lsource_mask) if (allocated(arr)) then narr = size(arr) @@ -70,11 +58,7 @@ module subroutine util_append_arr_DP(arr, source, lsource_mask) call util_resize(arr, narr + nsrc) - if (present(lsource_mask)) then - arr(narr + 1:narr + nsrc) = pack(source(:), lsource_mask(:)) - else - arr(narr + 1:narr + nsrc) = source(:) - end if + arr(narr + 1:narr + nsrc) = pack(source(:), lsource_mask(:)) return end subroutine util_append_arr_DP @@ -88,17 +72,13 @@ module subroutine util_append_arr_DPvec(arr, source, lsource_mask) ! Arguments real(DP), dimension(:,:), allocatable, intent(inout) :: arr !! Destination array real(DP), dimension(:,:), allocatable, intent(in) :: source !! Array to append - logical, dimension(:), optional, intent(in) :: lsource_mask !! Logical mask indicating which elements to append to + logical, dimension(:), intent(in) :: lsource_mask !! Logical mask indicating which elements to append to ! Internals integer(I4B) :: narr, nsrc if (.not. allocated(source)) return - if (present(lsource_mask)) then - nsrc = count(lsource_mask) - else - nsrc = size(source, dim=2) - end if + nsrc = count(lsource_mask) if (allocated(arr)) then narr = size(arr, dim=2) @@ -109,13 +89,9 @@ module subroutine util_append_arr_DPvec(arr, source, lsource_mask) call util_resize(arr, narr + nsrc) - if (present(lsource_mask)) then - arr(1, narr + 1:narr + nsrc) = pack(source(1,:), lsource_mask(:)) - arr(2, narr + 1:narr + nsrc) = pack(source(2,:), lsource_mask(:)) - arr(3, narr + 1:narr + nsrc) = pack(source(3,:), lsource_mask(:)) - else - arr(:, narr + 1:narr + nsrc) = source(:,:) - end if + arr(1, narr + 1:narr + nsrc) = pack(source(1,:), lsource_mask(:)) + arr(2, narr + 1:narr + nsrc) = pack(source(2,:), lsource_mask(:)) + arr(3, narr + 1:narr + nsrc) = pack(source(3,:), lsource_mask(:)) return end subroutine util_append_arr_DPvec @@ -129,17 +105,13 @@ module subroutine util_append_arr_I4B(arr, source, lsource_mask) ! Arguments integer(I4B), dimension(:), allocatable, intent(inout) :: arr !! Destination array integer(I4B), dimension(:), allocatable, intent(in) :: source !! Array to append - logical, dimension(:), optional, intent(in) :: lsource_mask !! Logical mask indicating which elements to append to + logical, dimension(:), intent(in) :: lsource_mask !! Logical mask indicating which elements to append to ! Internals integer(I4B) :: narr, nsrc if (.not. allocated(source)) return - if (present(lsource_mask)) then - nsrc = count(lsource_mask) - else - nsrc = size(source) - end if + nsrc = count(lsource_mask) if (allocated(arr)) then narr = size(arr) @@ -150,11 +122,7 @@ module subroutine util_append_arr_I4B(arr, source, lsource_mask) call util_resize(arr, narr + nsrc) - if (present(lsource_mask)) then - arr(narr + 1:narr + nsrc) = pack(source(:), lsource_mask(:)) - else - arr(narr + 1:narr + nsrc) = source(:) - end if + arr(narr + 1:narr + nsrc) = pack(source(:), lsource_mask(:)) return end subroutine util_append_arr_I4B @@ -168,7 +136,7 @@ module subroutine util_append_arr_logical(arr, source, lsource_mask) ! Arguments logical, dimension(:), allocatable, intent(inout) :: arr !! Destination array logical, dimension(:), allocatable, intent(in) :: source !! Array to append - logical, dimension(:), optional, intent(in) :: lsource_mask !! Logical mask indicating which elements to append to + logical, dimension(:), intent(in) :: lsource_mask !! Logical mask indicating which elements to append to ! Internals integer(I4B) :: narr, nsrc @@ -181,19 +149,11 @@ module subroutine util_append_arr_logical(arr, source, lsource_mask) narr = 0 end if - if (present(lsource_mask)) then - nsrc = count(lsource_mask) - else - nsrc = size(source) - end if + nsrc = count(lsource_mask) call util_resize(arr, narr + nsrc) - if (present(lsource_mask)) then - arr(narr + 1:narr + nsrc) = pack(source(:), lsource_mask(:)) - else - arr(narr + 1:narr + nsrc) = source(:) - end if + arr(narr + 1:narr + nsrc) = pack(source(:), lsource_mask(:)) return end subroutine util_append_arr_logical @@ -208,7 +168,7 @@ module subroutine util_append_body(self, source, lsource_mask) ! Arguments class(swiftest_body), intent(inout) :: self !! Swiftest body object class(swiftest_body), intent(in) :: source !! Source object to append - logical, dimension(:), optional, intent(in) :: lsource_mask !! Logical mask indicating which elements to append to + logical, dimension(:), intent(in) :: lsource_mask !! Logical mask indicating which elements to append to call util_append(self%name, source%name, lsource_mask) call util_append(self%id, source%id, lsource_mask) @@ -247,7 +207,7 @@ module subroutine util_append_pl(self, source, lsource_mask) ! Arguments class(swiftest_pl), intent(inout) :: self !! Swiftest massive body object class(swiftest_body), intent(in) :: source !! Source object to append - logical, dimension(:), optional, intent(in) :: lsource_mask !! Logical mask indicating which elements to append to + logical, dimension(:), intent(in) :: lsource_mask !! Logical mask indicating which elements to append to select type(source) @@ -287,7 +247,7 @@ module subroutine util_append_tp(self, source, lsource_mask) ! Arguments class(swiftest_tp), intent(inout) :: self !! Swiftest test particle object class(swiftest_body), intent(in) :: source !! Source object to append - logical, dimension(:), optional, intent(in) :: lsource_mask !! Logical mask indicating which elements to append to + logical, dimension(:), intent(in) :: lsource_mask !! Logical mask indicating which elements to append to select type(source) class is (swiftest_tp) diff --git a/src/util/util_get_energy_momentum.f90 b/src/util/util_get_energy_momentum.f90 index f8047b6b6..55d05e823 100644 --- a/src/util/util_get_energy_momentum.f90 +++ b/src/util/util_get_energy_momentum.f90 @@ -94,19 +94,22 @@ module subroutine util_get_energy_momentum_system(self, param) ! Do the central body potential energy component first !$omp simd - do i = 1, npl - associate(px => pl%xb(1,i), py => pl%xb(2,i), pz => pl%xb(3,i)) - pecb(i) = -cb%Gmass * pl%mass(i) / sqrt(px**2 + py**2 + pz**2) - end associate - end do + 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 ! Do the potential energy between pairs of massive bodies - do k = 1, pl%nplpl - associate(ik => pl%k_plpl(1, k), jk => pl%k_plpl(2, k)) - pepl(k) = -pl%Gmass(ik) * pl%mass(jk) / norm2(pl%xb(:, jk) - pl%xb(:, ik)) - lstatpl(k) = (lstatus(ik) .and. lstatus(jk)) - end associate - end do + 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 + + 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 system%pe = sum(pepl(:), lstatpl(:)) + sum(pecb(1:npl), lstatus(1:npl)) @@ -116,7 +119,7 @@ module subroutine util_get_energy_momentum_system(self, param) ! Potential energy from the oblateness term if (param%loblatecb) then !$omp simd - do i = 1, npl + do concurrent(i = 1:npl, lstatus(i)) irh(i) = 1.0_DP / norm2(pl%xh(:,i)) end do call obl_pot(npl, cb%Gmass, pl%mass, cb%j2rp2, cb%j4rp4, pl%xh, irh, oblpot) diff --git a/src/whm/whm_util.f90 b/src/whm/whm_util.f90 index f3dc15d3e..a71e4439c 100644 --- a/src/whm/whm_util.f90 +++ b/src/whm/whm_util.f90 @@ -11,7 +11,7 @@ module subroutine whm_util_append_pl(self, source, lsource_mask) !! Arguments class(whm_pl), intent(inout) :: self !! WHM massive body object class(swiftest_body), intent(in) :: source !! Source object to append - logical, dimension(:), optional, intent(in) :: lsource_mask !! Logical mask indicating which elements to append to + logical, dimension(:), intent(in) :: lsource_mask !! Logical mask indicating which elements to append to select type(source) class is (whm_pl)