diff --git a/Makefile.Defines b/Makefile.Defines index 90a2ceb58..9264c39c6 100644 --- a/Makefile.Defines +++ b/Makefile.Defines @@ -48,12 +48,14 @@ COLLRESOLVE_HOME = $(ROOT_DIR)/collresolve/ # DO NOT include in FFLAGS the "-c" option to compile object only # this is done explicitly as needed in the Makefile -ADVIXE_FLAGS = -g -O2 -qopt-report=5 -vecabi=cmdtarget -simd -shared-intel -debug inline-debug-info -DTBB_DEBUG -DTBB_USE_THREADING_TOOLS -xhost -traceback -parallel-source-info=2 +ADVIXE_INCLUDE = -I$(ADVISOR_2019_DIR)/include/intel64 +ADVIXE_LINK = -L$(ADVISOR_2019_DIR)/lib64 -ladvisor +ADVIXE_FLAGS = -g -O2 -qopt-report=5 -vecabi=cmdtarget -simd -shared-intel -debug inline-debug-info -DTBB_DEBUG -DTBB_USE_THREADING_TOOLS -xhost -traceback #Be sure to set the environment variable KMP_FORKJOIN_FRAMES=1 for OpenMP debuging in vtune -IDEBUG = -O0 -init=snan,arrays -nogen-interfaces -no-pie -no-ftz -fpe-all=0 -g -traceback -mp1 -qopt-matmul -fp-model strict -fpe0 -debug all -align all -pad -ip -prec-div -prec-sqrt -assume protect-parens -CB -no-wrap-margin +IDEBUG = -O0 -init=snan,arrays -nogen-interfaces -no-pie -no-ftz -fpe-all=0 -g -traceback -mp1 -qopt-matmul -fp-model strict -fpe0 -debug all -align all -pad -ip -prec-div -prec-sqrt -assume protect-parens -CB -no-wrap-margin STRICTREAL = -fp-model=precise -prec-div -prec-sqrt -assume protect-parens -SIMDVEC = -simd -xhost -align all -assume contiguous_assumed_shape -vecabi=cmdtarget -fp-model no-except -fma +SIMDVEC = -simd -xhost -align all -pad -assume contiguous_assumed_shape -vecabi=cmdtarget -fp-model no-except -fma PAR = -qopenmp -parallel -parallel-source-info=2 HEAPARR = -heap-arrays 4194304 OPTREPORT = -qopt-report=5 @@ -66,14 +68,14 @@ GMEM = -fsanitize-address-use-after-scope -fstack-check -fsanitize=bounds-stri GWARNINGS = -Wall -Warray-bounds -Wimplicit-interface -Wextra -Warray-temporaries GPRODUCTION = -O3 -ffree-line-length-none $(GPAR) -INCLUDES = -I$(SWIFTEST_HOME)/include -I$(NETCDF_FORTRAN_HOME)/include -I$(MKLROOT)/include -LINKS = -L$(MKLROOT)/lib/intel64 -L$(NETCDF_FORTRAN_HOME)/lib -lswiftest -lnetcdf -lnetcdff -qopt-matmul #$(PAR) +INCLUDES = -I$(SWIFTEST_HOME)/include -I$(NETCDF_FORTRAN_HOME)/include -I$(MKLROOT)/include +LINKS = -L$(MKLROOT)/lib/intel64 -L$(NETCDF_FORTRAN_HOME)/lib -lswiftest -lnetcdf -lnetcdff -qopt-matmul $(PAR) -FSTRICTFLAGS = $(IDEBUG) #$(SIMDVEC) $(PAR) $(HEAPARR) -FFLAGS = $(IDEBUG) #$(SIMDVEC) $(PAR) $(HEAPARR) -#FSTRICTFLAGS = $(IPRODUCTION) $(STRICTREAL) -g -traceback -#FFLAGS = $(IPRODUCTION) -fp-model=fast -g -traceback +#FSTRICTFLAGS = $(IDEBUG) $(SIMDVEC) $(PAR) +#FFLAGS = $(IDEBUG) $(SIMDVEC) $(PAR) +FSTRICTFLAGS = $(ADVIXE_FLAGS) $(PAR) $(SIMDVEC) $(STRICTREAL) +FFLAGS = $(ADVIXE_FLAGS) $(PAR) $(SIMDVEC) -fp-model=fast FORTRAN = ifort AR = xiar CC = icc diff --git a/src/discard/discard.f90 b/src/discard/discard.f90 index d791274ef..0fe86e8fa 100644 --- a/src/discard/discard.f90 +++ b/src/discard/discard.f90 @@ -128,18 +128,22 @@ subroutine discard_cb_tp(tp, system, param) tp%status(i) = DISCARDED_RMAX write(idstr, *) tp%id(i) write(timestr, *) param%t - write(*, *) "Particle " // trim(adjustl(tp%info(i)%name)) // " (" // trim(adjustl(idstr)) // ")" // " too far from the central body at t = " // trim(adjustl(timestr)) + write(*, *) "Particle " // trim(adjustl(tp%info(i)%name)) // " (" // trim(adjustl(idstr)) // ")" // & + " too far from the central body at t = " // trim(adjustl(timestr)) tp%ldiscard(i) = .true. tp%lmask(i) = .false. - call tp%info(i)%set_value(status="DISCARDED_RMAX", discard_time=param%t, discard_xh=tp%xh(:,i), discard_vh=tp%vh(:,i)) + call tp%info(i)%set_value(status="DISCARDED_RMAX", discard_time=param%t, discard_xh=tp%xh(:,i), & + discard_vh=tp%vh(:,i)) else if ((param%rmin >= 0.0_DP) .and. (rh2 < rmin2)) then tp%status(i) = DISCARDED_RMIN write(idstr, *) tp%id(i) write(timestr, *) param%t - write(*, *) "Particle " // trim(adjustl(tp%info(i)%name)) // " (" // trim(adjustl(idstr)) // ")" // " too close to the central body at t = " // trim(adjustl(timestr)) + write(*, *) "Particle " // trim(adjustl(tp%info(i)%name)) // " (" // trim(adjustl(idstr)) // ")" // & + " too close to the central body at t = " // trim(adjustl(timestr)) tp%ldiscard(i) = .true. tp%lmask(i) = .false. - call tp%info(i)%set_value(status="DISCARDED_RMIN", discard_time=param%t, discard_xh=tp%xh(:,i), discard_vh=tp%vh(:,i), discard_body_id=cb%id) + call tp%info(i)%set_value(status="DISCARDED_RMIN", discard_time=param%t, discard_xh=tp%xh(:,i), & + discard_vh=tp%vh(:,i), discard_body_id=cb%id) else if (param%rmaxu >= 0.0_DP) then rb2 = dot_product(tp%xb(:, i), tp%xb(:, i)) vb2 = dot_product(tp%vb(:, i), tp%vb(:, i)) @@ -148,10 +152,12 @@ subroutine discard_cb_tp(tp, system, param) tp%status(i) = DISCARDED_RMAXU write(idstr, *) tp%id(i) write(timestr, *) param%t - write(*, *) "Particle " // trim(adjustl(tp%info(i)%name)) // " (" // trim(adjustl(idstr)) // ")" // " is unbound and too far from barycenter at t = " // trim(adjustl(timestr)) + write(*, *) "Particle " // trim(adjustl(tp%info(i)%name)) // " (" // trim(adjustl(idstr)) // ")" // & + " is unbound and too far from barycenter at t = " // trim(adjustl(timestr)) tp%ldiscard(i) = .true. tp%lmask(i) = .false. - call tp%info(i)%set_value(status="DISCARDED_RMAXU", discard_time=param%t, discard_xh=tp%xh(:,i), discard_vh=tp%vh(:,i)) + call tp%info(i)%set_value(status="DISCARDED_RMAXU", discard_time=param%t, discard_xh=tp%xh(:,i), & + discard_vh=tp%vh(:,i)) end if end if end if @@ -199,9 +205,11 @@ subroutine discard_peri_tp(tp, system, param) tp%status(i) = DISCARDED_PERI write(idstr, *) tp%id(i) write(timestr, *) param%t - write(*, *) "Particle " // trim(adjustl(tp%info(i)%name)) // " (" // trim(adjustl(idstr)) // ")" // " perihelion distance too small at t = " // trim(adjustl(timestr)) + write(*, *) "Particle " // trim(adjustl(tp%info(i)%name)) // " (" // trim(adjustl(idstr)) // ")" // & + " perihelion distance too small at t = " // trim(adjustl(timestr)) tp%ldiscard(i) = .true. - call tp%info(i)%set_value(status="DISCARDED_PERI", discard_time=param%t, discard_xh=tp%xh(:,i), discard_vh=tp%vh(:,i), discard_body_id=pl%id(j)) + call tp%info(i)%set_value(status="DISCARDED_PERI", discard_time=param%t, discard_xh=tp%xh(:,i), & + discard_vh=tp%vh(:,i), discard_body_id=pl%id(j)) end if end if end if @@ -250,7 +258,8 @@ subroutine discard_pl_tp(tp, system, param) // " too close to massive body " // trim(adjustl(pl%info(j)%name)) // " (" // trim(adjustl(idstrj)) // ")" & // " at t = " // trim(adjustl(timestr)) tp%ldiscard(i) = .true. - call tp%info(i)%set_value(status="DISCARDED_PLR", discard_time=param%t, discard_xh=tp%xh(:,i), discard_vh=tp%vh(:,i), discard_body_id=pl%id(j)) + call tp%info(i)%set_value(status="DISCARDED_PLR", discard_time=param%t, discard_xh=tp%xh(:,i), & + discard_vh=tp%vh(:,i), discard_body_id=pl%id(j)) exit end if end do diff --git a/src/encounter/encounter_check.f90 b/src/encounter/encounter_check.f90 index b125d0ff9..20fcb3bd5 100644 --- a/src/encounter/encounter_check.f90 +++ b/src/encounter/encounter_check.f90 @@ -2,50 +2,8 @@ use swiftest contains - module subroutine encounter_check_all(nenc, index1, index2, x1, v1, x2, v2, renc1, renc2, dt, lencounter, lvdotr) - !! author: David A. Minton - !! - !! Check for encounters between n pairs of bodies. - !! This implementation is general for any type of body. For instance, for massive bodies, you would pass x1 = x2, for test particles renc2 is an array of zeros, etc. - !! - implicit none - ! Arguments - integer(I4B), intent(in) :: nenc !! Number of encounters in the encounter lists - integer(I4B), dimension(:), intent(in) :: index1 !! List of indices for body 1 in each encounter - integer(I4B), dimension(:), intent(in) :: index2 !! List of indices for body 2 in each encounter1 - real(DP), dimension(:,:), intent(in) :: x1, v1 !! Array of indices of bodies 1 - real(DP), dimension(:,:), intent(in) :: x2, v2 !! Array of indices of bodies 2 - real(DP), dimension(:), intent(in) :: renc1 !! Radius of encounter regions of bodies 1 - real(DP), dimension(:), intent(in) :: renc2 !! Radius of encounter regions of bodies 2 - real(DP), intent(in) :: dt !! Step size - logical, dimension(:), intent(out) :: lencounter !! Logical array indicating which pairs are in a close encounter state - logical, dimension(:), intent(out) :: lvdotr !! Logical array indicating which pairs are approaching - ! Internals - integer(I4B) :: i, j, k - real(DP) :: xr, yr, zr, vxr, vyr, vzr, renc12 - - !$omp parallel do simd default(firstprivate) schedule(static)& - !$omp shared(lencounter, lvdotr, index1, index2, x1, v1, x2, v2) & - !$omp lastprivate(i, j, xr, yr, zr, vxr, vyr, vzr, renc12) - do k = 1, nenc - i = index1(k) - j = index2(k) - xr = x2(1, j) - x1(1, i) - yr = x2(2, j) - x1(2, i) - zr = x2(3, j) - x1(3, i) - vxr = v2(1, j) - v1(1, i) - vyr = v2(2, j) - v1(2, i) - vzr = v2(3, j) - v1(3, i) - renc12 = renc1(i) + renc2(j) - call encounter_check_one(xr, yr, zr, vxr, vyr, vzr, renc12, dt, lencounter(k), lvdotr(k)) - end do - !$omp end parallel do simd - - return - end subroutine encounter_check_all - - - module subroutine encounter_check_all_plpl(param, npl, x, v, renc, dt, lvdotr, index1, index2, nenc) + module subroutine encounter_check_all_plpl(param, npl, x, v, renc, dt, & + nenc, index1, index2, lvdotr) !! author: David A. Minton !! !! Check for encounters between massive bodies. Choose between the standard triangular or the Sort & Sweep method based on user inputs @@ -58,15 +16,16 @@ module subroutine encounter_check_all_plpl(param, npl, x, v, renc, dt, lvdotr, i real(DP), dimension(:,:), intent(in) :: v !! Velocity vectors of massive bodies real(DP), dimension(:), intent(in) :: renc !! Critical radii of massive bodies that defines an encounter real(DP), intent(in) :: dt !! Step size - logical, dimension(:), allocatable, intent(out) :: lvdotr !! Logical flag indicating the sign of v .dot. x + integer(I8B), intent(out) :: nenc !! Total number of encounters integer(I4B), dimension(:), allocatable, intent(out) :: index1 !! List of indices for body 1 in each encounter integer(I4B), dimension(:), allocatable, intent(out) :: index2 !! List of indices for body 2 in each encounter - integer(I4B), intent(out) :: nenc !! Total number of encounters + logical, dimension(:), allocatable, intent(out) :: lvdotr !! Logical flag indicating the sign of v .dot. x ! Internals type(interaction_timer), save :: itimer logical, save :: lfirst = .true. logical, save :: skipit = .false. ! This will be used to ensure that the sort & sweep subroutine gets called at least once before timing it so that the extent array is nearly sorted when it is timed integer(I8B) :: nplpl = 0_I8B + integer(I8B) :: k if (param%ladaptive_encounters_plpl .and. (.not. skipit)) then nplpl = (npl * (npl - 1) / 2) @@ -83,11 +42,11 @@ module subroutine encounter_check_all_plpl(param, npl, x, v, renc, dt, lvdotr, i param%lencounter_sas_plpl = .false. end if end if - + if (param%lencounter_sas_plpl) then - call encounter_check_all_sort_and_sweep_plpl(npl, x, v, renc, dt, lvdotr, index1, index2, nenc) + call encounter_check_all_sort_and_sweep_plpl(npl, x, v, renc, dt, nenc, index1, index2, lvdotr) else - call encounter_check_all_triangular_plpl(npl, x, v, renc, dt, lvdotr, index1, index2, nenc) + call encounter_check_all_triangular_plpl(npl, x, v, renc, dt, nenc, index1, index2, lvdotr) end if if (skipit) then @@ -105,7 +64,8 @@ module subroutine encounter_check_all_plpl(param, npl, x, v, renc, dt, lvdotr, i end subroutine encounter_check_all_plpl - module subroutine encounter_check_all_plplm(param, nplm, nplt, xplm, vplm, xplt, vplt, rencm, renct, dt, lvdotr, index1, index2, nenc) + module subroutine encounter_check_all_plplm(param, nplm, nplt, xplm, vplm, xplt, vplt, rencm, renct, dt, & + nenc, index1, index2, lvdotr) !! author: David A. Minton !! !! Check for encounters between fully interacting massive bodies partially interacting massive bodies. Choose between the standard triangular or the Sort & Sweep method based on user inputs @@ -122,22 +82,24 @@ module subroutine encounter_check_all_plplm(param, nplm, nplt, xplm, vplm, xplt, real(DP), dimension(:), intent(in) :: rencm !! Critical radii of fully interacting massive bodies that defines an encounter real(DP), dimension(:), intent(in) :: renct !! Critical radii of partially interacting massive bodies that defines an encounter real(DP), intent(in) :: dt !! Step size - logical, dimension(:), allocatable, intent(out) :: lvdotr !! Logical flag indicating the sign of v .dot. x + integer(I8B), intent(out) :: nenc !! Total number of encounters integer(I4B), dimension(:), allocatable, intent(out) :: index1 !! List of indices for body 1 in each encounter integer(I4B), dimension(:), allocatable, intent(out) :: index2 !! List of indices for body 2 in each encounter - integer(I4B), intent(out) :: nenc !! Total number of encounters + logical, dimension(:), allocatable, intent(out) :: lvdotr !! Logical flag indicating the sign of v .dot. x ! Internals type(interaction_timer), save :: itimer logical, save :: lfirst = .true. logical, save :: skipit = .false. integer(I8B) :: nplplm = 0_I8B integer(I4B) :: npl, i + integer(I8B) :: k logical, dimension(:), allocatable :: plmplt_lvdotr !! Logical flag indicating the sign of v .dot. x in the plm-plt group integer(I4B), dimension(:), allocatable :: plmplt_index1 !! List of indices for body 1 in each encounter in the plm-plt group integer(I4B), dimension(:), allocatable :: plmplt_index2 !! List of indices for body 2 in each encounter in the plm-lt group - integer(I4B) :: plmplt_nenc !! Number of encounters of the plm-plt group + integer(I8B) :: plmplt_nenc !! Number of encounters of the plm-plt group class(swiftest_parameters), allocatable :: tmp_param !! Temporary parameter structure to turn off adaptive timer for the pl-pl phase if necessary - integer(I4B), dimension(:), allocatable :: itmp, ind + integer(I8B), dimension(:), allocatable :: ind + integer(I4B), dimension(:), allocatable :: itmp logical, dimension(:), allocatable :: ltmp if (param%ladaptive_encounters_plpl .and. (.not. skipit)) then @@ -163,12 +125,14 @@ module subroutine encounter_check_all_plplm(param, nplm, nplt, xplm, vplm, xplt, tmp_param%ladaptive_encounters_plpl = .false. ! Start with the pl-pl group - call encounter_check_all_plpl(tmp_param, nplm, xplm, vplm, rencm, dt, lvdotr, index1, index2, nenc) + call encounter_check_all_plpl(tmp_param, nplm, xplm, vplm, rencm, dt, nenc, index1, index2, lvdotr) if (param%lencounter_sas_plpl) then - call encounter_check_all_sort_and_sweep_plplm(nplm, nplt, xplm, vplm, xplt, vplt, rencm, renct, dt, plmplt_lvdotr, plmplt_index1, plmplt_index2, plmplt_nenc) + call encounter_check_all_sort_and_sweep_plplm(nplm, nplt, xplm, vplm, xplt, vplt, rencm, renct, dt, & + plmplt_nenc, plmplt_index1, plmplt_index2, plmplt_lvdotr) else - call encounter_check_all_triangular_plplm(nplm, nplt, xplm, vplm, xplt, vplt, rencm, renct, dt, plmplt_lvdotr, plmplt_index1, plmplt_index2, plmplt_nenc) + call encounter_check_all_triangular_plplm(nplm, nplt, xplm, vplm, xplt, vplt, rencm, renct, dt, & + plmplt_nenc, plmplt_index1, plmplt_index2, plmplt_lvdotr) end if if (skipit) then @@ -197,17 +161,19 @@ module subroutine encounter_check_all_plplm(param, nplm, nplt, xplm, vplm, xplt, call move_alloc(ltmp, lvdotr) nenc = nenc + plmplt_nenc - call util_sort(index1, itmp) - call util_sort_rearrange(index1, itmp, nenc) - call util_sort_rearrange(index2, itmp, nenc) - call util_sort_rearrange(lvdotr, itmp, nenc) + call util_sort(index1, ind) + call util_sort_rearrange(index1, ind, nenc) + call util_sort_rearrange(index2, ind, nenc) + call util_sort_rearrange(lvdotr, ind, nenc) + end if return end subroutine encounter_check_all_plplm - module subroutine encounter_check_all_pltp(param, npl, ntp, xpl, vpl, xtp, vtp, renc, dt, lvdotr, index1, index2, nenc) + module subroutine encounter_check_all_pltp(param, npl, ntp, xpl, vpl, xtp, vtp, renc, dt, & + nenc, index1, index2, lvdotr) !! author: David A. Minton !! !! Check for encounters between massive bodies and test particles. Choose between the standard triangular or the Sort & Sweep method based on user inputs @@ -223,10 +189,10 @@ module subroutine encounter_check_all_pltp(param, npl, ntp, xpl, vpl, xtp, vtp, real(DP), dimension(:,:), intent(in) :: vtp !! Velocity vectors of test particles real(DP), dimension(:), intent(in) :: renc !! Critical radii of massive bodies that defines an encounter real(DP), intent(in) :: dt !! Step size - logical, dimension(:), allocatable, intent(out) :: lvdotr !! Logical flag indicating the sign of v .dot. x + integer(I8B), intent(out) :: nenc !! Total number of encounters integer(I4B), dimension(:), allocatable, intent(out) :: index1 !! List of indices for body 1 in each encounter integer(I4B), dimension(:), allocatable, intent(out) :: index2 !! List of indices for body 2 in each encounter - integer(I4B), intent(out) :: nenc !! Total number of encounters + logical, dimension(:), allocatable, intent(out) :: lvdotr !! Logical flag indicating the sign of v .dot. x ! Internals type(interaction_timer), save :: itimer logical, save :: lfirst = .true. @@ -256,9 +222,9 @@ module subroutine encounter_check_all_pltp(param, npl, ntp, xpl, vpl, xtp, vtp, end if if (param%lencounter_sas_pltp) then - call encounter_check_all_sort_and_sweep_pltp(npl, ntp, xpl, vpl, xtp, vtp, renc, dt, lvdotr, index1, index2, nenc) + call encounter_check_all_sort_and_sweep_pltp(npl, ntp, xpl, vpl, xtp, vtp, renc, dt, nenc, index1, index2, lvdotr) else - call encounter_check_all_triangular_pltp(npl, ntp, xpl, vpl, xtp, vtp, renc, dt, lvdotr, index1, index2, nenc) + call encounter_check_all_triangular_pltp(npl, ntp, xpl, vpl, xtp, vtp, renc, dt, nenc, index1, index2, lvdotr) end if if (.not.lfirst .and. param%ladaptive_encounters_pltp .and. npltp > 0) then @@ -269,105 +235,7 @@ module subroutine encounter_check_all_pltp(param, npl, ntp, xpl, vpl, xtp, vtp, end subroutine encounter_check_all_pltp - subroutine encounter_check_reduce_broadphase(n, nenc, index1, index2, lencounter, lvdotr) - !! author: David A. Minton - !! - !! Takes the candidate encounter lists that came out of the broad phase and narrow it down to the true encounters - implicit none - ! Arguments - integer(I4B), intent(in) :: n !! Number of bodies - integer(I4B), intent(inout) :: nenc !! Number of encountering bodies (input is the broad phase value, output is the final narrow phase value) - integer(I4B), dimension(:), allocatable, intent(inout) :: index1 !! List of indices for body 1 in each encounter - integer(I4B), dimension(:), allocatable, intent(inout) :: index2 !! List of indices for body 2 in each encounter - logical, dimension(:), allocatable, intent(inout) :: lencounter !! Logical flag indicating which of the pairs identified in the broad phase was selected in the narrow phase - logical, dimension(:), allocatable, intent(inout) :: lvdotr !! Logical flag indicating the sign of v .dot. x - ! Internals - integer(I4B) :: i, i0, j, k, nenci, klo, khi - integer(I4B), dimension(:), allocatable :: itmp - integer(I4B), dimension(:), allocatable :: ind - integer(I4B), dimension(n) :: ibeg, iend - logical, dimension(:), allocatable :: ltmp - - nenc = count(lencounter(:)) ! Count the true number of encounters - if (nenc == 0) then - if (allocated(index1)) deallocate(index1) - if (allocated(index2)) deallocate(index2) - if (allocated(lvdotr)) deallocate(lvdotr) - return - end if - - allocate(itmp(nenc)) - itmp(:) = pack(index1(:), lencounter(:)) - call move_alloc(itmp, index1) - - allocate(itmp(nenc)) - itmp(:) = pack(index2(:), lencounter(:)) - call move_alloc(itmp, index2) - - allocate(ltmp(nenc)) - ltmp(:) = pack(lvdotr(:), lencounter(:)) - call move_alloc(ltmp, lvdotr) - - if (allocated(lencounter)) deallocate(lencounter) - allocate(lencounter(nenc)) - lencounter(:) = .true. - - call util_sort(index1, ind) - call util_sort_rearrange(index1, ind, nenc) - call util_sort_rearrange(index2, ind, nenc) - call util_sort_rearrange(lvdotr, ind, nenc) - - ! Get the bounds on the bodies in the first index - ibeg(:) = n - iend(:) = 1 - i0 = index1(1) - ibeg(i0) = 1 - do k = 2, nenc - i = index1(k) - if (i /= i0) then - iend(i0) = k - 1 - ibeg(i) = k - i0 = i - end if - if (k == nenc) iend(i) = k - end do - - ! Sort on the second index and remove duplicates - if (allocated(itmp)) deallocate(itmp) - allocate(itmp, source=index2) - do concurrent(i = 1:n, iend(i) - ibeg(i) > 0) - klo = ibeg(i) - khi = iend(i) - nenci = khi - klo + 1 - if (allocated(ind)) deallocate(ind) - call util_sort(index2(klo:khi), ind) - index2(klo:khi) = itmp(klo - 1 + ind(:)) - do j = klo + 1, khi - if (index2(j) == index2(j - 1)) lencounter(j) = .false. - end do - end do - - if (count(lencounter(:)) == nenc) return - nenc = count(lencounter(:)) ! Count the true number of encounters - - if (allocated(itmp)) deallocate(itmp) - allocate(itmp(nenc)) - itmp(:) = pack(index1(:), lencounter(:)) - call move_alloc(itmp, index1) - - allocate(itmp(nenc)) - itmp(:) = pack(index2(:), lencounter(:)) - call move_alloc(itmp, index2) - - allocate(ltmp(nenc)) - ltmp(:) = pack(lvdotr(:), lencounter(:)) - call move_alloc(ltmp, lvdotr) - - return - end subroutine encounter_check_reduce_broadphase - - - subroutine encounter_check_all_sort_and_sweep_plpl(npl, x, v, renc, dt, lvdotr, index1, index2, nenc) + subroutine encounter_check_all_sort_and_sweep_plpl(npl, x, v, renc, dt, nenc, index1, index2, lvdotr) !! author: David A. Minton !! !! Check for encounters between massive bodies. @@ -381,17 +249,16 @@ subroutine encounter_check_all_sort_and_sweep_plpl(npl, x, v, renc, dt, lvdotr, real(DP), dimension(:,:), intent(in) :: v !! Velocity vectors of massive bodies real(DP), dimension(:), intent(in) :: renc !! Critical radii of massive bodies that defines an encounter real(DP), intent(in) :: dt !! Step size - logical, dimension(:), allocatable, intent(out) :: lvdotr !! Logical flag indicating the sign of v .dot. x + integer(I8B), intent(out) :: nenc !! Total number of encounters integer(I4B), dimension(:), allocatable, intent(out) :: index1 !! List of indices for body 1 in each encounter integer(I4B), dimension(:), allocatable, intent(out) :: index2 !! List of indices for body 2 in each encounter - integer(I4B), intent(out) :: nenc !! Total number of encounters + logical, dimension(:), allocatable, intent(out) :: lvdotr !! Logical flag indicating the sign of v .dot. x ! Internals integer(I4B) :: i, dim, n integer(I4B), save :: npl_last = 0 type(encounter_bounding_box), save :: boundingbox logical, dimension(:), allocatable :: lencounter integer(I2B), dimension(npl) :: vshift_min, vshift_max - type(walltimer) :: timer if (npl == 0) return @@ -418,24 +285,14 @@ subroutine encounter_check_all_sort_and_sweep_plpl(npl, x, v, renc, dt, lvdotr, end do !$omp end parallel do - call boundingbox%sweep(npl, nenc, index1, index2) - - if (nenc > 0) then - ! Now that we have identified potential pairs, use the narrow-phase process to get the final values - allocate(lencounter(nenc)) - allocate(lvdotr(nenc)) - - call encounter_check_all(nenc, index1, index2, x, v, x, v, renc, renc, dt, lencounter, lvdotr) - - call encounter_check_reduce_broadphase(npl, nenc, index1, index2, lencounter, lvdotr) - deallocate(lencounter) - end if + call boundingbox%sweep(npl, x, v, renc, dt, nenc, index1, index2, lvdotr) return end subroutine encounter_check_all_sort_and_sweep_plpl - subroutine encounter_check_all_sort_and_sweep_plplm(nplm, nplt, xplm, vplm, xplt, vplt, rencm, renct, dt, lvdotr, index1, index2, nenc) + subroutine encounter_check_all_sort_and_sweep_plplm(nplm, nplt, xplm, vplm, xplt, vplt, rencm, renct, dt, & + nenc, index1, index2, lvdotr) !! author: David A. Minton !! !! Check for encounters between massive bodies and test particles. @@ -453,10 +310,10 @@ subroutine encounter_check_all_sort_and_sweep_plplm(nplm, nplt, xplm, vplm, xplt real(DP), dimension(:), intent(in) :: rencm !! Critical radii of fully interacting massive bodies that defines an encounter real(DP), dimension(:), intent(in) :: renct !! Critical radii of partially interacting massive bodies that defines an encounter real(DP), intent(in) :: dt !! Step size - logical, dimension(:), allocatable, intent(out) :: lvdotr !! Logical flag indicating the sign of v .dot. x + integer(I8B), intent(out) :: nenc !! Total number of encounter integer(I4B), dimension(:), allocatable, intent(out) :: index1 !! List of indices for body 1 in each encounter integer(I4B), dimension(:), allocatable, intent(out) :: index2 !! List of indices for body 2 in each encounter - integer(I4B), intent(out) :: nenc !! Total number of encounter + logical, dimension(:), allocatable, intent(out) :: lvdotr !! Logical flag indicating the sign of v .dot. x ! Internals type(encounter_bounding_box), save :: boundingbox integer(I4B) :: i, dim, n, ntot @@ -464,7 +321,7 @@ subroutine encounter_check_all_sort_and_sweep_plplm(nplm, nplt, xplm, vplm, xplt logical, dimension(:), allocatable :: lencounter integer(I2B), dimension(nplm) :: vplmshift_min, vplmshift_max integer(I2B), dimension(nplt) :: vpltshift_min, vpltshift_max - type(walltimer) :: timer + logical, save :: lfirst=.true. ! If this is the first time through, build the index lists if ((nplm == 0) .or. (nplt == 0)) return @@ -475,9 +332,7 @@ subroutine encounter_check_all_sort_and_sweep_plplm(nplm, nplt, xplm, vplm, xplt call boundingbox%setup(ntot, ntot_last) ntot_last = ntot end if - - call timer%reset() - call timer%start() + !$omp parallel do default(private) schedule(static) & !$omp shared(xplm, xplt, vplm, vplt, rencm, renct, boundingbox) & !$omp firstprivate(dt, nplm, nplt, ntot) @@ -504,26 +359,15 @@ subroutine encounter_check_all_sort_and_sweep_plplm(nplm, nplt, xplm, vplm, xplt xplt(dim,1:nplt) + renct(1:nplt) + vpltshift_max(1:nplt) * vplt(dim,1:nplt) * dt]) end do !$omp end parallel do - - call boundingbox%sweep(nplm, nplt, nenc, index1, index2) - - if (nenc > 0) then - ! Shift tiny body indices back into the range of the input position and velocity arrays - index2(:) = index2(:) - nplm - ! Now that we have identified potential pairs, use the narrow-phase process to get the final values - allocate(lencounter(nenc)) - allocate(lvdotr(nenc)) + call boundingbox%sweep(nplm, nplt, xplm, vplm, xplt, vplt, rencm, renct, dt, nenc, index1, index2, lvdotr) - call encounter_check_all(nenc, index1, index2, xplm, vplm, xplt, vplt, rencm, renct, dt, lencounter, lvdotr) - - call encounter_check_reduce_broadphase(ntot, nenc, index1, index2, lencounter, lvdotr) - end if return end subroutine encounter_check_all_sort_and_sweep_plplm - subroutine encounter_check_all_sort_and_sweep_pltp(npl, ntp, xpl, vpl, xtp, vtp, renc, dt, lvdotr, index1, index2, nenc) + subroutine encounter_check_all_sort_and_sweep_pltp(npl, ntp, xpl, vpl, xtp, vtp, rencpl, dt, & + nenc, index1, index2, lvdotr) !! author: David A. Minton !! !! Check for encounters between massive bodies and test particles. @@ -538,12 +382,12 @@ subroutine encounter_check_all_sort_and_sweep_pltp(npl, ntp, xpl, vpl, xtp, vtp, real(DP), dimension(:,:), intent(in) :: vpl !! Velocity vectors of massive bodies real(DP), dimension(:,:), intent(in) :: xtp !! Position vectors of massive bodies real(DP), dimension(:,:), intent(in) :: vtp !! Velocity vectors of massive bodies - real(DP), dimension(:), intent(in) :: renc !! Critical radii of massive bodies that defines an encounter + real(DP), dimension(:), intent(in) :: rencpl !! Critical radii of massive bodies that defines an encounter real(DP), intent(in) :: dt !! Step size - logical, dimension(:), allocatable, intent(out) :: lvdotr !! Logical flag indicating the sign of v .dot. x + integer(I8B), intent(out) :: nenc !! Total number of encounter integer(I4B), dimension(:), allocatable, intent(out) :: index1 !! List of indices for body 1 in each encounter integer(I4B), dimension(:), allocatable, intent(out) :: index2 !! List of indices for body 2 in each encounter - integer(I4B), intent(out) :: nenc !! Total number of encounter + logical, dimension(:), allocatable, intent(out) :: lvdotr !! Logical flag indicating the sign of v .dot. x ! Internals type(encounter_bounding_box), save :: boundingbox integer(I4B) :: i, dim, n, ntot @@ -551,7 +395,7 @@ subroutine encounter_check_all_sort_and_sweep_pltp(npl, ntp, xpl, vpl, xtp, vtp, logical, dimension(:), allocatable :: lencounter integer(I2B), dimension(npl) :: vplshift_min, vplshift_max integer(I2B), dimension(ntp) :: vtpshift_min, vtpshift_max - real(DP), dimension(ntp) :: renctp + real(DP), dimension(ntp) :: renctp ! If this is the first time through, build the index lists if ((ntp == 0) .or. (npl == 0)) return @@ -562,9 +406,11 @@ subroutine encounter_check_all_sort_and_sweep_pltp(npl, ntp, xpl, vpl, xtp, vtp, call boundingbox%setup(ntot, ntot_last) ntot_last = ntot end if + + renctp(:) = 0.0_DP !$omp parallel do default(private) schedule(static) & - !$omp shared(xpl, xtp, vpl, vtp, renc, boundingbox) & + !$omp shared(xpl, xtp, vpl, vtp, rencpl, renctp, boundingbox) & !$omp firstprivate(dt, npl, ntp, ntot) do dim = 1, SWEEPDIM where(vpl(dim,1:npl) < 0.0_DP) @@ -583,50 +429,93 @@ subroutine encounter_check_all_sort_and_sweep_pltp(npl, ntp, xpl, vpl, xtp, vtp, vtpshift_max(1:ntp) = 1 end where - call boundingbox%aabb(dim)%sort(ntot, [xpl(dim,1:npl) - renc(1:npl) + vplshift_min(1:npl) * vpl(dim,1:npl) * dt, & - xtp(dim,1:ntp) + vtpshift_min(1:ntp) * vtp(dim,1:ntp) * dt, & - xpl(dim,1:npl) + renc(1:npl) + vplshift_max(1:npl) * vpl(dim,1:npl) * dt, & - xtp(dim,1:ntp) + vtpshift_max(1:ntp) * vtp(dim,1:ntp) * dt]) + call boundingbox%aabb(dim)%sort(ntot, [xpl(dim,1:npl) - rencpl(1:npl) + vplshift_min(1:npl) * vpl(dim,1:npl) * dt, & + xtp(dim,1:ntp) - renctp(1:ntp) + vtpshift_min(1:ntp) * vtp(dim,1:ntp) * dt, & + xpl(dim,1:npl) + rencpl(1:npl) + vplshift_max(1:npl) * vpl(dim,1:npl) * dt, & + xtp(dim,1:ntp) + renctp(1:ntp) + vtpshift_max(1:ntp) * vtp(dim,1:ntp) * dt]) end do !$omp end parallel do - call boundingbox%sweep(npl, ntp, nenc, index1, index2) - - if (nenc > 0) then - ! Shift test particle indices back into the proper range - index2(:) = index2(:) - npl + call boundingbox%sweep(npl, ntp, xpl, vpl, xtp, vtp, rencpl, renctp, dt, nenc, index1, index2, lvdotr) - ! Now that we have identified potential pairs, use the narrow-phase process to get the final values - allocate(lencounter(nenc)) - allocate(lvdotr(nenc)) - renctp(:) = 0.0_DP + return + end subroutine encounter_check_all_sort_and_sweep_pltp - call encounter_check_all(nenc, index1, index2, xpl, vpl, xtp, vtp, renc, renctp, dt, lencounter, lvdotr) - call encounter_check_reduce_broadphase(ntot, nenc, index1, index2, lencounter, lvdotr) + pure subroutine encounter_check_all_sweep_one(i, n, xi, yi, zi, vxi, vyi, vzi, x, y, z, vx, vy, vz, renci, renc, dt, & + ind_arr, lgood, nenci, index1, index2, lvdotr) + !! author: David A. Minton + !! + !! Check for encounters between the ith body and all other bodies. + !! This is used in the narrow phase of the sort & sweep algorithm + implicit none + ! Arguments + integer(I4B), intent(in) :: i !! Index of the ith body that is being checked + integer(I4B), intent(in) :: n !! Total number of bodies being checked + real(DP), intent(in) :: xi, yi, zi !! Position vector components of the ith body + real(DP), intent(in) :: vxi, vyi, vzi !! Velocity vector components of the ith body + real(DP), dimension(:), intent(in) :: x, y, z !! Arrays of position vector components of all bodies + real(DP), dimension(:), intent(in) :: vx, vy, vz !! Arrays of velocity vector components of all bodies + real(DP), intent(in) :: renci !! Encounter radius of the ith body + real(DP), dimension(:), intent(in) :: renc !! Array of encounter radii of all bodies + real(DP), intent(in) :: dt !! Step size + integer(I4B), dimension(:), intent(in) :: ind_arr !! Index array [1, 2, ..., n] + logical, dimension(:), intent(in) :: lgood !! Logical array mask where true values correspond to bodies selected in the broad phase + integer(I8B), intent(out) :: nenci !! Total number of encountering bodies + integer(I4B), dimension(:), allocatable, intent(inout) :: index1 !! Array of indices of the ith body of size nenci [i, i, ..., i] + integer(I4B), dimension(:), allocatable, intent(inout) :: index2 !! Array of indices of the encountering bodies of size nenci + logical, dimension(:), allocatable, intent(inout) :: lvdotr !! v.dot.r direction array + ! Internals + integer(I4B) :: j + integer(I8B) :: k + real(DP) :: xr, yr, zr, vxr, vyr, vzr, renc12 + logical, dimension(n) :: lencounteri, lvdotri + lencounteri(:) = .false. + do concurrent(j = 1:n, lgood(j)) + xr = x(j) - xi + yr = y(j) - yi + zr = z(j) - zi + vxr = vx(j) - vxi + vyr = vy(j) - vyi + vzr = vz(j) - vzi + renc12 = renci + renc(j) + call encounter_check_one(xr, yr, zr, vxr, vyr, vzr, renc12, dt, lencounteri(j), lvdotri(j)) + end do + if (any(lencounteri(:))) then + nenci = count(lencounteri(:)) + allocate(lvdotr(nenci), index1(nenci), index2(nenci)) + index1(:) = i + index2(:) = pack(ind_arr(1:n), lencounteri(1:n)) + lvdotr(:) = pack(lvdotri(1:n), lencounteri(1:n)) end if return - end subroutine encounter_check_all_sort_and_sweep_pltp + end subroutine encounter_check_all_sweep_one - pure subroutine encounter_check_all_triangular_one(i, n, xi, yi, zi, vxi, vyi, vzi, x, y, z, vx, vy, vz, renci, renc, dt, ind_arr, lenci) + pure subroutine encounter_check_all_triangular_one(i, n, xi, yi, zi, vxi, vyi, vzi, x, y, z, vx, vy, vz, renci, renc, & + dt, ind_arr, lenci) + !! author: David A. Minton + !! + !! Check for encounters between the ith body and all other bodies. + !! This is the upper triangular (double loop) version. implicit none ! Arguments - integer(I4B), intent(in) :: i - integer(I4B), intent(in) :: n - real(DP), intent(in) :: xi, yi, zi - real(DP), intent(in) :: vxi, vyi, vzi - real(DP), dimension(:), intent(in) :: x, y, z - real(DP), dimension(:), intent(in) :: vx, vy, vz - real(DP), intent(in) :: renci - real(DP), dimension(:), intent(in) :: renc - real(DP), intent(in) :: dt - integer(I4B), dimension(:), intent(in) :: ind_arr - type(encounter_list), intent(out) :: lenci + integer(I4B), intent(in) :: i !! Index of the ith body that is being checked + integer(I4B), intent(in) :: n !! Total number of bodies being checked + real(DP), intent(in) :: xi, yi, zi !! Position vector components of the ith body + real(DP), intent(in) :: vxi, vyi, vzi !! Velocity vector components of the ith body + real(DP), dimension(:), intent(in) :: x, y, z !! Arrays of position vector components of all bodies + real(DP), dimension(:), intent(in) :: vx, vy, vz !! Arrays of velocity vector components of all bodies + real(DP), intent(in) :: renci !! Encounter radius of the ith body + real(DP), dimension(:), intent(in) :: renc !! Array of encounter radii of all bodies + real(DP), intent(in) :: dt !! Step size + integer(I4B), dimension(:), intent(in) :: ind_arr !! Index array [1, 2, ..., n] + type(encounter_list), intent(out) :: lenci !! Output encounter lists containing number of encounters, the v.dot.r direction array, and the index list of encountering bodies ! Internals - integer(I4B) :: j, nenci + integer(I4B) :: j + integer(I8B) :: k, nenci real(DP) :: xr, yr, zr, vxr, vyr, vzr, renc12 logical, dimension(n) :: lencounteri, lvdotri @@ -641,18 +530,19 @@ pure subroutine encounter_check_all_triangular_one(i, n, xi, yi, zi, vxi, vyi, v call encounter_check_one(xr, yr, zr, vxr, vyr, vzr, renc12, dt, lencounteri(j), lvdotri(j)) end do nenci = count(lencounteri(i+1:n)) - if (nenci > 0) then - allocate(lenci%lvdotr(nenci), lenci%index2(nenci)) + if (nenci > 0_I8B) then + allocate(lenci%lvdotr(nenci), lenci%index1(nenci), lenci%index2(nenci)) lenci%nenc = nenci - lenci%lvdotr(:) = pack(lvdotri(i+1:n), lencounteri(i+1:n)) + lenci%index1(:) = i lenci%index2(:) = pack(ind_arr(i+1:n), lencounteri(i+1:n)) + lenci%lvdotr(:) = pack(lvdotri(i+1:n), lencounteri(i+1:n)) end if return end subroutine encounter_check_all_triangular_one - subroutine encounter_check_all_triangular_plpl(npl, x, v, renc, dt, lvdotr, index1, index2, nenc) + subroutine encounter_check_all_triangular_plpl(npl, x, v, renc, dt, nenc, index1, index2, lvdotr) !! author: David A. Minton !! !! Check for encounters between massive bodies. Split off from the main subroutine for performance @@ -664,17 +554,16 @@ subroutine encounter_check_all_triangular_plpl(npl, x, v, renc, dt, lvdotr, inde real(DP), dimension(:,:), intent(in) :: v !! Velocity vectors of massive bodies real(DP), dimension(:), intent(in) :: renc !! Critical radii of massive bodies that defines an encounter real(DP), intent(in) :: dt !! Step size - logical, dimension(:), allocatable, intent(out) :: lvdotr !! Logical flag indicating the sign of v .dot. x + integer(I8B), intent(out) :: nenc !! Total number of encounters integer(I4B), dimension(:), allocatable, intent(out) :: index1 !! List of indices for body 1 in each encounter integer(I4B), dimension(:), allocatable, intent(out) :: index2 !! List of indices for body 2 in each encounter - integer(I4B), intent(out) :: nenc !! Total number of encounters + logical, dimension(:), allocatable, intent(out) :: lvdotr !! Logical flag indicating the sign of v .dot. x ! Internals integer(I4B) :: i, j, k, nenci, j0, j1 real(DP) :: xr, yr, zr, vxr, vyr, vzr, renc12 logical, dimension(npl) :: lencounteri, lvdotri integer(I4B), dimension(:), allocatable, save :: ind_arr type(encounter_list), dimension(npl) :: lenc - type(walltimer) :: timer call util_index_array(ind_arr, npl) @@ -687,6 +576,7 @@ subroutine encounter_check_all_triangular_plpl(npl, x, v, renc, dt, lvdotr, inde x(1,:), x(2,:), x(3,:), & v(1,:), v(2,:), v(3,:), & renc(i), renc(:), dt, ind_arr(:), lenc(i)) + if (lenc(i)%nenc > 0) lenc(i)%index1(:) = i end do !$omp end parallel do @@ -696,7 +586,8 @@ subroutine encounter_check_all_triangular_plpl(npl, x, v, renc, dt, lvdotr, inde end subroutine encounter_check_all_triangular_plpl - subroutine encounter_check_all_triangular_plplm(nplm, nplt, xplm, vplm, xplt, vplt, rencm, renct, dt, lvdotr, index1, index2, nenc) + subroutine encounter_check_all_triangular_plplm(nplm, nplt, xplm, vplm, xplt, vplt, rencm, renct, dt, & + nenc, index1, index2, lvdotr) !! author: David A. Minton !! !! Check for encounters between massive bodies and test particles. Split off from the main subroutine for performance @@ -712,21 +603,19 @@ subroutine encounter_check_all_triangular_plplm(nplm, nplt, xplm, vplm, xplt, vp real(DP), dimension(:), intent(in) :: rencm !! Critical radii of fully interacting massive bodies that defines an encounter real(DP), dimension(:), intent(in) :: renct !! Critical radii of partially interacting massive bodies that defines an encounter real(DP), intent(in) :: dt !! Step size - logical, dimension(:), allocatable, intent(out) :: lvdotr !! Logical flag indicating the sign of v .dot. x + integer(I8B), intent(out) :: nenc !! Total number of encounters integer(I4B), dimension(:), allocatable, intent(out) :: index1 !! List of indices for body 1 in each encounter integer(I4B), dimension(:), allocatable, intent(out) :: index2 !! List of indices for body 2 in each encounter - integer(I4B), intent(out) :: nenc !! Total number of encounters + logical, dimension(:), allocatable, intent(out) :: lvdotr !! Logical flag indicating the sign of v .dot. x ! Internals - integer(I4B) :: i, j, nenci, j0, j1 - real(DP) :: xr, yr, zr, vxr, vyr, vzr, renc12 + integer(I4B) :: i logical, dimension(nplt) :: lencounteri, lvdotri integer(I4B), dimension(:), allocatable, save :: ind_arr type(encounter_list), dimension(nplm) :: lenc - type(walltimer) :: timer call util_index_array(ind_arr, nplt) - !$omp parallel do default(private) schedule(static)& + !$omp parallel do default(private) schedule(dynamic)& !$omp shared(xplm, vplm, xplt, vplt, rencm, renct, lenc, ind_arr) & !$omp firstprivate(nplm, nplt, dt) do i = 1, nplm @@ -735,6 +624,7 @@ subroutine encounter_check_all_triangular_plplm(nplm, nplt, xplm, vplm, xplt, vp xplt(1,:), xplt(2,:), xplt(3,:), & vplt(1,:), vplt(2,:), vplt(3,:), & rencm(i), renct(:), dt, ind_arr(:), lenc(i)) + if (lenc(i)%nenc > 0) lenc(i)%index1(:) = i end do !$omp end parallel do @@ -744,7 +634,8 @@ subroutine encounter_check_all_triangular_plplm(nplm, nplt, xplm, vplm, xplt, vp end subroutine encounter_check_all_triangular_plplm - subroutine encounter_check_all_triangular_pltp(npl, ntp, xpl, vpl, xtp, vtp, renc, dt, lvdotr, index1, index2, nenc) + subroutine encounter_check_all_triangular_pltp(npl, ntp, xpl, vpl, xtp, vtp, renc, dt, & + nenc, index1, index2, lvdotr) !! author: David A. Minton !! !! Check for encounters between massive bodies and test particles. Split off from the main subroutine for performance @@ -759,13 +650,12 @@ subroutine encounter_check_all_triangular_pltp(npl, ntp, xpl, vpl, xtp, vtp, ren real(DP), dimension(:,:), intent(in) :: vtp !! Velocity vectors of massive bodies real(DP), dimension(:), intent(in) :: renc !! Critical radii of massive bodies that defines an encounter real(DP), intent(in) :: dt !! Step size - logical, dimension(:), allocatable, intent(out) :: lvdotr !! Logical flag indicating the sign of v .dot. x + integer(I8B), intent(out) :: nenc !! Total number of encounters integer(I4B), dimension(:), allocatable, intent(out) :: index1 !! List of indices for body 1 in each encounter integer(I4B), dimension(:), allocatable, intent(out) :: index2 !! List of indices for body 2 in each encounter - integer(I4B), intent(out) :: nenc !! Total number of encounters + logical, dimension(:), allocatable, intent(out) :: lvdotr !! Logical flag indicating the sign of v .dot. x ! Internals - integer(I4B) :: i, j, nenci, j0, j1 - real(DP) :: xr, yr, zr, vxr, vyr, vzr, renc1, renc2 + integer(I4B) :: i logical, dimension(ntp) :: lencounteri, lvdotri integer(I4B), dimension(:), allocatable, save :: ind_arr type(encounter_list), dimension(npl) :: lenc @@ -774,7 +664,7 @@ subroutine encounter_check_all_triangular_pltp(npl, ntp, xpl, vpl, xtp, vtp, ren call util_index_array(ind_arr, ntp) renct(:) = 0.0_DP - !$omp parallel do default(private) schedule(static)& + !$omp parallel do default(private) schedule(dynamic)& !$omp shared(xpl, vpl, xtp, vtp, renc, renct, lenc, ind_arr) & !$omp firstprivate(npl, ntp, dt) do i = 1, npl @@ -783,6 +673,7 @@ subroutine encounter_check_all_triangular_pltp(npl, ntp, xpl, vpl, xtp, vtp, ren xtp(1,:), xtp(2,:), xtp(3,:), & vtp(1,:), vtp(2,:), vtp(3,:), & renc(i), renct(:), dt, ind_arr(:), lenc(i)) + if (lenc(i)%nenc > 0) lenc(i)%index1(:) = i end do !$omp end parallel do @@ -792,7 +683,7 @@ subroutine encounter_check_all_triangular_pltp(npl, ntp, xpl, vpl, xtp, vtp, ren end subroutine encounter_check_all_triangular_pltp - module pure subroutine encounter_check_one(xr, yr, zr, vxr, vyr, vzr, renc, dt, lencounter, lvdotr) + module elemental subroutine encounter_check_one(xr, yr, zr, vxr, vyr, vzr, renc, dt, lencounter, lvdotr) !! author: David A. Minton !! !! Determine whether a test particle and planet are having or will have an encounter within the next time step @@ -847,12 +738,13 @@ module subroutine encounter_check_collapse_ragged_list(ragged_list, n1, nenc, in ! Arguments type(encounter_list), dimension(:), intent(in) :: ragged_list !! The ragged encounter list integer(I4B), intent(in) :: n1 !! Number of bodies 1 - integer(I4B), intent(out) :: nenc !! Total number of encountersj + integer(I8B), intent(out) :: nenc !! Total number of encountersj integer(I4B), dimension(:), allocatable, intent(out) :: index1 !! Array of indices for body 1 integer(I4B), dimension(:), allocatable, intent(out) :: index2 !! Array of indices for body 1 logical, dimension(:), allocatable, intent(out), optional :: lvdotr !! Array indicating which bodies are approaching ! Internals - integer(I4B) :: i, j0, j1, nenci + integer(I4B) :: i + integer(I8B) :: j1, j0, nenci integer(I4B), dimension(n1) :: ibeg associate(nenc_arr => ragged_list(:)%nenc) @@ -872,43 +764,107 @@ module subroutine encounter_check_collapse_ragged_list(ragged_list, n1, nenc, in j0 = j0 + nenci end do - !$omp parallel do default(private) & + !$omp parallel do simd default(private) schedule(simd: static)& !$omp shared(ragged_list, index1, index2, ibeg, lvdotr) & !$omp firstprivate(n1) do i = 1,n1 - if (ragged_list(i)%nenc == 0) cycle + if (ragged_list(i)%nenc == 0_I8B) cycle nenci = ragged_list(i)%nenc j0 = ibeg(i) - j1 = j0 + nenci - 1 - index1(j0:j1) = i + j1 = j0 + nenci - 1_I8B + index1(j0:j1) = ragged_list(i)%index1(:) index2(j0:j1) = ragged_list(i)%index2(:) if (present(lvdotr)) lvdotr(j0:j1) = ragged_list(i)%lvdotr(:) end do - !$omp end parallel do + !$omp end parallel do simd return end subroutine encounter_check_collapse_ragged_list - pure subroutine encounter_check_make_ragged_list(lencounteri, ind_arr, nenc,index2) + subroutine encounter_check_remove_duplicates(n, nenc, index1, index2, lvdotr) + !! author: David A. Minton + !! + !! Takes the candidate encounter lists that came out of the sort & sweep method and remove any duplicates. implicit none ! Arguments - logical, dimension(:), intent(in) :: lencounteri - integer(I4B), dimension(:), intent(in) :: ind_arr - integer(I4B), intent(out) :: nenc - integer(I4B), dimension(:), allocatable, intent(out) :: index2 + integer(I4B), intent(in) :: n !! Number of bodies + integer(I8B), intent(inout) :: nenc !! Number of encountering bodies (input is the broad phase value, output is the final narrow phase value) + integer(I4B), dimension(:), allocatable, intent(inout) :: index1 !! List of indices for body 1 in each encounter + integer(I4B), dimension(:), allocatable, intent(inout) :: index2 !! List of indices for body 2 in each encounter + logical, dimension(:), allocatable, intent(inout) :: lvdotr !! Logical flag indicating the sign of v .dot. x ! Internals + integer(I4B) :: i, i0 + integer(I8B) :: j, k, klo, khi, nenci + integer(I4B), dimension(:), allocatable :: itmp + integer(I8B), dimension(:), allocatable :: ind + integer(I8B), dimension(n) :: ibeg, iend + logical, dimension(:), allocatable :: ltmp + logical, dimension(nenc) :: lencounter - nenc = count(lencounteri(:)) - if (nenc > 0) then - allocate(index2(nenc)) - index2(:) = pack(ind_arr(:), lencounteri(:)) + if (nenc == 0) then + if (allocated(index1)) deallocate(index1) + if (allocated(index2)) deallocate(index2) + if (allocated(lvdotr)) deallocate(lvdotr) + return end if + call util_sort(index1, ind) + call util_sort_rearrange(index1, ind, nenc) + call util_sort_rearrange(index2, ind, nenc) + call util_sort_rearrange(lvdotr, ind, nenc) + + ! Get the bounds on the bodies in the first index + ibeg(:) = n + iend(:) = 1_I8B + i0 = index1(1) + ibeg(i0) = 1_I8B + do k = 2_I8B, nenc + i = index1(k) + if (i /= i0) then + iend(i0) = k - 1_I8B + ibeg(i) = k + i0 = i + end if + if (k == nenc) iend(i) = k + end do + + lencounter(:) = .true. + ! Sort on the second index and remove duplicates + if (allocated(itmp)) deallocate(itmp) + allocate(itmp, source=index2) + do concurrent(i = 1:n, iend(i) - ibeg(i) > 0_I8B) + klo = ibeg(i) + khi = iend(i) + nenci = khi - klo + 1_I8B + if (allocated(ind)) deallocate(ind) + call util_sort(index2(klo:khi), ind) + index2(klo:khi) = itmp(klo - 1_I8B + ind(:)) + do j = klo + 1_I8B, khi + if (index2(j) == index2(j - 1_I8B)) lencounter(j) = .false. + end do + end do + + if (count(lencounter(:)) == nenc) return + nenc = count(lencounter(:)) ! Count the true number of encounters + + if (allocated(itmp)) deallocate(itmp) + allocate(itmp(nenc)) + itmp(:) = pack(index1(:), lencounter(:)) + call move_alloc(itmp, index1) + + allocate(itmp(nenc)) + itmp(:) = pack(index2(:), lencounter(:)) + call move_alloc(itmp, index2) + + allocate(ltmp(nenc)) + ltmp(:) = pack(lvdotr(:), lencounter(:)) + call move_alloc(ltmp, lvdotr) + return - end subroutine encounter_check_make_ragged_list + end subroutine encounter_check_remove_duplicates + - module pure subroutine encounter_check_sort_aabb_1D(self, n, extent_arr) !! author: David A. Minton !! @@ -920,12 +876,11 @@ module pure subroutine encounter_check_sort_aabb_1D(self, n, extent_arr) integer(I4B), intent(in) :: n !! Number of bodies with extents real(DP), dimension(:), intent(in) :: extent_arr !! Array of extents of size 2*n ! Internals - integer(I4B) :: i, j, k, ibox, jbox - logical, dimension(:), allocatable :: lfresh + integer(I8B) :: i, k call util_sort(extent_arr, self%ind) - do concurrent(k = 1:2*n) + do concurrent(k = 1_I8B:2_I8B * n) i = self%ind(k) if (i <= n) then self%ibeg(i) = k @@ -937,246 +892,207 @@ module pure subroutine encounter_check_sort_aabb_1D(self, n, extent_arr) return end subroutine encounter_check_sort_aabb_1D - - module subroutine encounter_check_sweep_aabb_double_list(self, n1, n2, nenc, index1, index2) + + module subroutine encounter_check_sweep_aabb_double_list(self, n1, n2, x1, v1, x2, v2, renc1, renc2, dt, & + nenc, index1, index2, lvdotr) !! author: David A. Minton !! - !! Sweeps the sorted bounding box extents and returns the encounter candidates + !! Sweeps the sorted bounding box extents and returns the true encounters (combines broad and narrow phases) + !! Double list version (e.g. pl-tp or plm-plt) implicit none ! Arguments - class(encounter_bounding_box), intent(inout) :: self !! Multi-dimensional bounding box structure - integer(I4B), intent(in) :: n1 !! Number of bodies 1 - integer(I4B), intent(in) :: n2 !! Number of bodies 2 - integer(I4B), intent(out) :: nenc !! Total number of encounter candidates - integer(I4B), dimension(:), allocatable, intent(out) :: index1 !! List of indices for body 1 in each encounter candidate pair - integer(I4B), dimension(:), allocatable, intent(out) :: index2 !! List of indices for body 2 in each encounter candidate pair - !Internals - integer(I4B) :: i, k, ntot, dim + class(encounter_bounding_box), intent(inout) :: self !! Multi-dimensional bounding box structure + integer(I4B), intent(in) :: n1 !! Number of bodies 1 + integer(I4B), intent(in) :: n2 !! Number of bodies 2 + real(DP), dimension(:,:), intent(in) :: x1, v1 !! Array of position and velocity vectorrs for bodies 1 + real(DP), dimension(:,:), intent(in) :: x2, v2 !! Array of position and velocity vectorrs for bodies 2 + real(DP), dimension(:), intent(in) :: renc1 !! Radius of encounter regions of bodies 1 + real(DP), dimension(:), intent(in) :: renc2 !! Radius of encounter regions of bodies 2 + real(DP), intent(in) :: dt !! Step size + integer(I8B), intent(out) :: nenc !! Total number of encounter candidates + integer(I4B), dimension(:), allocatable, intent(out) :: index1 !! List of indices for body 1 in each encounter candidate pair + integer(I4B), dimension(:), allocatable, intent(out) :: index2 !! List of indices for body 2 in each encounter candidate pair + logical, dimension(:), allocatable, intent(out) :: lvdotr !! Logical array indicating which pairs are approaching + ! Internals + integer(I4B) :: ii, i, j, ntot, nbox, dim + integer(I8B) :: k + logical, dimension(n1+n2) :: loverlap + logical, dimension(SWEEPDIM,n1+n2) :: loverlap_by_dimension + integer(I4B), dimension(SWEEPDIM) :: noverlap + integer(I4B), dimension(SWEEPDIM,n1+n2) :: nbox_arr + logical, dimension(SWEEPDIM,2*(n1+n2)) :: llist1 + integer(I4B), dimension(SWEEPDIM,2*(n1+n2)) :: ext_ind + integer(I4B), dimension(:), allocatable :: x_ind type(encounter_list), dimension(n1+n2) :: lenc !! Array of encounter lists (one encounter list per body) integer(I4B), dimension(:), allocatable, save :: ind_arr - integer(I4B), dimension(:), allocatable :: ibeg, iend + integer(I8B) :: ibeg, iend + real(DP), dimension(2*(n1+n2)) :: xind, yind, zind, vxind, vyind, vzind, rencind ntot = n1 + n2 call util_index_array(ind_arr, ntot) - allocate(ibeg(SWEEPDIM * ntot)) - allocate(iend(SWEEPDIM * ntot)) - do i = 1, ntot - do dim = 1, SWEEPDIM - ibeg((i - 1) * SWEEPDIM + dim) = self%aabb(dim)%ibeg(i) - iend((i - 1) * SWEEPDIM + dim) = self%aabb(dim)%iend(i) - end do - end do - ! Sweep the intervals for each of the massive bodies along one dimension - ! This will build a ragged pair of index lists inside of the lenc data structure - call encounter_check_sweep_aabb_all_double_list(n1, n2, self%aabb(1)%ind(:), reshape(ibeg(:), [SWEEPDIM, ntot]), reshape(iend(:), [SWEEPDIM, ntot]), ind_arr(:), lenc(:)) - - call encounter_check_collapse_ragged_list(lenc, ntot, nenc, index1, index2) - ! Reorder the pairs and sort the first index in order to remove any duplicates - do concurrent(k = 1:nenc, index2(k) < index1(k)) - i = index2(k) - index2(k) = index1(k) - index1(k) = i + do concurrent(dim = 1:SWEEPDIM) + loverlap_by_dimension(dim,:) = (self%aabb(dim)%ibeg(:) + 1_I8B) < (self%aabb(dim)%iend(:) - 1_I8B) + where(self%aabb(dim)%ind(:) > ntot) + ext_ind(dim,:) = self%aabb(dim)%ind(:) - ntot + elsewhere + ext_ind(dim,:) = self%aabb(dim)%ind(:) + endwhere end do + llist1(:,:) = ext_ind(:,:) <= n1 + where(.not.llist1(:,:)) ext_ind(:,:) = ext_ind(:,:) - n1 - return - end subroutine encounter_check_sweep_aabb_double_list - + loverlap(:) = loverlap_by_dimension(1,:) + do dim = 2, SWEEPDIM + loverlap(:) = loverlap(:) .and. loverlap_by_dimension(dim,:) + end do - module subroutine encounter_check_sweep_aabb_single_list(self, n, nenc, index1, index2) - !! author: David A. Minton - !! - !! Sweeps the sorted bounding box extents and returns the encounter candidates. Mutual encounters - !! allowed. That is, all bodies are from the same list - implicit none - ! Arguments - class(encounter_bounding_box), intent(inout) :: self !! Multi-dimensional bounding box structure - integer(I4B), intent(in) :: n !! Number of bodies 1 - integer(I4B), intent(out) :: nenc !! Total number of encounter candidates - integer(I4B), dimension(:), allocatable, intent(out) :: index1 !! List of indices for body 1 in each encounter candidate pair - integer(I4B), dimension(:), allocatable, intent(out) :: index2 !! List of indices for body 2 in each encounter candidate pair - !Internals - Integer(I4B) :: i, k, dim - type(encounter_list), dimension(n) :: lenc !! Array of encounter lists (one encounter list per body) - integer(I4B), dimension(:), allocatable, save :: ind_arr - type(walltimer) :: timer - integer(I4B), dimension(:), allocatable :: ibeg, iend + dim = 1 + where(llist1(dim,:)) + xind(:) = x1(1,ext_ind(dim,:)) + yind(:) = x1(2,ext_ind(dim,:)) + zind(:) = x1(3,ext_ind(dim,:)) + vxind(:) = v1(1,ext_ind(dim,:)) + vyind(:) = v1(2,ext_ind(dim,:)) + vzind(:) = v1(3,ext_ind(dim,:)) + rencind(:) = renc1(ext_ind(dim,:)) + elsewhere + xind(:) = x2(1,ext_ind(dim,:)) + yind(:) = x2(2,ext_ind(dim,:)) + zind(:) = x2(3,ext_ind(dim,:)) + vxind(:) = v2(1,ext_ind(dim,:)) + vyind(:) = v2(2,ext_ind(dim,:)) + vzind(:) = v2(3,ext_ind(dim,:)) + rencind(:) = renc2(ext_ind(dim,:)) + endwhere - call util_index_array(ind_arr, n) - allocate(ibeg(SWEEPDIM * n)) - allocate(iend(SWEEPDIM * n)) - do i = 1, n - do dim = 1, SWEEPDIM - ibeg((i - 1) * SWEEPDIM + dim) = self%aabb(dim)%ibeg(i) - iend((i - 1) * SWEEPDIM + dim) = self%aabb(dim)%iend(i) - end do + where(.not.loverlap(:)) lenc(:)%nenc = 0 + !$omp parallel default(private) & + !$omp shared(self, ext_ind, lenc, loverlap, x1, v1, x2, v2, renc1, renc2, xind, yind, zind, vxind, vyind, vzind, rencind, llist1) & + !$omp firstprivate(ntot, n1, n2, dt, dim) + + ! Do the first group of bodies (i is in list 1, all the others are from list 2) + !$omp do schedule(static) + do i = 1, n1 + if (loverlap(i)) then + ibeg = self%aabb(dim)%ibeg(i) + 1_I8B + iend = self%aabb(dim)%iend(i) - 1_I8B + nbox = iend - ibeg + 1 + call encounter_check_all_sweep_one(i, nbox, x1(1,i), x1(2,i), x1(3,i), v1(1,i), v1(2,i), v1(3,i), & + xind(ibeg:iend), yind(ibeg:iend), zind(ibeg:iend),& + vxind(ibeg:iend), vyind(ibeg:iend), vzind(ibeg:iend), & + renc1(i), rencind(ibeg:iend), dt, ext_ind(dim,ibeg:iend), & + .not.llist1(dim,ibeg:iend), lenc(i)%nenc, lenc(i)%index1, lenc(i)%index2, lenc(i)%lvdotr) + end if + end do + !$omp end do nowait + + ! Do the second group of bodies (i is in list 2, all the others are in list 1) + !$omp do schedule(static) + do i = n1+1, ntot + if (loverlap(i)) then + ibeg = self%aabb(dim)%ibeg(i) + 1_I8B + iend = self%aabb(dim)%iend(i) - 1_I8B + nbox = iend - ibeg + 1 + ii = i - n1 + call encounter_check_all_sweep_one(ii, nbox, x2(1,ii), x2(2,ii), x2(3,ii), v2(1,ii), v2(2,ii), v2(3,ii), & + xind(ibeg:iend), yind(ibeg:iend), zind(ibeg:iend),& + vxind(ibeg:iend), vyind(ibeg:iend), vzind(ibeg:iend), & + renc2(ii), rencind(ibeg:iend), dt, ext_ind(dim,ibeg:iend), & + llist1(dim,ibeg:iend), lenc(i)%nenc, lenc(i)%index2, lenc(i)%index1, lenc(i)%lvdotr) + end if end do + !$omp end do nowait - ! Sweep the intervals for each of the massive bodies along one dimension - ! This will build a ragged pair of index lists inside of the lenc data structure - call encounter_check_sweep_aabb_all_single_list(n, self%aabb(1)%ind(:), reshape(ibeg(:), [SWEEPDIM, n]), reshape(iend(:), [SWEEPDIM, n]), ind_arr(:), lenc(:)) + !$omp end parallel - call encounter_check_collapse_ragged_list(lenc, n, nenc, index1, index2) + call encounter_check_collapse_ragged_list(lenc, ntot, nenc, index1, index2, lvdotr) - ! Reorder the pairs and sort the first index in order to remove any duplicates - do concurrent(k = 1:nenc, index2(k) < index1(k)) - i = index2(k) - index2(k) = index1(k) - index1(k) = i - end do + call encounter_check_remove_duplicates(ntot, nenc, index1, index2, lvdotr) return - end subroutine encounter_check_sweep_aabb_single_list + end subroutine encounter_check_sweep_aabb_double_list - subroutine encounter_check_sweep_aabb_all_double_list(n1, n2, ext_ind, ibeg, iend, ind_arr, lenc) + module subroutine encounter_check_sweep_aabb_single_list(self, n, x, v, renc, dt, nenc, index1, index2, lvdotr) !! author: David A. Minton !! - !! Performs the loop part of the sweep operation. Double list version (e.g. pl-tp or plm-plt) + !! Sweeps the sorted bounding box extents and returns the true encounters (combines broad and narrow phases) + !! Double list version (e.g. pl-tp or plm-plt) implicit none ! Arguments - integer(I4B), intent(in) :: n1, n2 !! Number of bodies - integer(I4B), dimension(:), intent(in) :: ext_ind !! Sorted index array of extents - integer(I4B), dimension(:,:), intent(in) :: ibeg, iend !! Beginning and ending index lists in the n-dimensions - integer(I4B), dimension(:), intent(in) :: ind_arr !! index array for mapping body 2 indexes - type(encounter_list), dimension(:), intent(inout) :: lenc !! Encounter list for the ith body + class(encounter_bounding_box), intent(inout) :: self !! Multi-dimensional bounding box structure + integer(I4B), intent(in) :: n !! Number of bodies + real(DP), dimension(:,:), intent(in) :: x, v !! Array of position and velocity vectors + real(DP), dimension(:), intent(in) :: renc !! Radius of encounter regions of bodies 1 + real(DP), intent(in) :: dt !! Step size + integer(I8B), intent(out) :: nenc !! Total number of encounter candidates + integer(I4B), dimension(:), allocatable, intent(out) :: index1 !! List of indices for one body in each encounter candidate pair + integer(I4B), dimension(:), allocatable, intent(out) :: index2 !! List of indices for the other body in each encounter candidate pair + logical, dimension(:), allocatable, intent(out) :: lvdotr !! Logical array indicating which pairs are approaching ! Internals - integer(I4B) :: i, ntot - logical, dimension(n1+n2) :: lencounteri - integer(I4B), dimension(SWEEPDIM) :: ibegi, iendi + integer(I4B) :: i, nbox, dim + integer(I8B) :: k, itmp + logical, dimension(n) :: loverlap + logical, dimension(2*n) :: lencounteri + real(DP), dimension(2*n) :: xind, yind, zind, vxind, vyind, vzind, rencind + integer(I4B), dimension(SWEEPDIM,2*n) :: ext_ind + type(encounter_list), dimension(n) :: lenc !! Array of encounter lists (one encounter list per body) + integer(I4B), dimension(:), allocatable, save :: ind_arr + integer(I8B) :: ibeg, iend - ntot = n1 + n2 - !$omp parallel do default(private) schedule(dynamic)& - !$omp shared(ext_ind, ibeg, iend, ind_arr, lenc) & - !$omp firstprivate(ntot, n1, n2) - do i = 1, ntot - ibegi(1) = ibeg(1,i) + 1 - iendi(1) = iend(1,i) - 1 - if (iendi(1) >= ibegi(1)) then - ibegi(2:SWEEPDIM) = ibeg(2:SWEEPDIM,i) - iendi(2:SWEEPDIM) = iend(2:SWEEPDIM,i) - call encounter_check_sweep_aabb_one_double_list(i, n1, n2, ntot, ext_ind(:), ibegi(:), iendi(:), ibeg(:,:), iend(:,:), lencounteri(:)) - call encounter_check_make_ragged_list(lencounteri(:), ind_arr(:), lenc(i)%nenc, lenc(i)%index2) - else - lenc(i)%nenc = 0 - end if - end do - !$omp end parallel do + call util_index_array(ind_arr, n) + dim = 1 - return - end subroutine encounter_check_sweep_aabb_all_double_list + ! Sweep the intervals for each of the massive bodies along one dimension + ! This will build a ragged pair of index lists inside of the lenc data structure + where(self%aabb(dim)%ind(:) > n) + ext_ind(dim,:) = self%aabb(1)%ind(:) - n + elsewhere + ext_ind(dim,:) = self%aabb(1)%ind(:) + endwhere + xind(:) = x(1,ext_ind(dim,:)) + yind(:) = x(2,ext_ind(dim,:)) + zind(:) = x(3,ext_ind(dim,:)) + vxind(:) = v(1,ext_ind(dim,:)) + vyind(:) = v(2,ext_ind(dim,:)) + vzind(:) = v(3,ext_ind(dim,:)) + rencind(:) = renc(ext_ind(dim,:)) - subroutine encounter_check_sweep_aabb_all_single_list(n, ext_ind, ibeg, iend, ind_arr, lenc) - !! author: David A. Minton - !! - !! Performs the loop part of the sweep operation. Single list version (e.g. pl-pl) - implicit none - ! Arguments - integer(I4B), intent(in) :: n !! Number of bodies - integer(I4B), dimension(:), intent(in) :: ext_ind !! Sorted index array of extents - integer(I4B), dimension(:,:), intent(in) :: ibeg, iend !! Beginning and ending index lists in the n-dimension - integer(I4B), dimension(:), intent(in) :: ind_arr !! index array for mapping body 2 indexes - type(encounter_list), dimension(:), intent(inout) :: lenc !! Encounter list for the ith body - ! Internals - integer(I4B) :: i - logical, dimension(n) :: lencounteri - integer(I4B), dimension(SWEEPDIM) :: ibegi, iendi + loverlap(:) = (self%aabb(dim)%ibeg(:) + 1_I8B) < (self%aabb(dim)%iend(:) - 1_I8B) + where(.not.loverlap(:)) lenc(:)%nenc = 0 - !$omp parallel do default(private) schedule(dynamic)& - !$omp shared(ext_ind, ibeg, iend, ind_arr, lenc) & - !$omp firstprivate(n) & - !$omp lastprivate(ibegi, iendi, lencounteri) + !$omp parallel do default(private) schedule(static)& + !$omp shared(self, ext_ind, lenc, loverlap, x, v, renc, xind, yind, zind, vxind, vyind, vzind, rencind) & + !$omp firstprivate(n, dt, dim) do i = 1, n - ibegi(1) = ibeg(1,i) + 1 - iendi(1) = iend(1,i) - 1 - if (iendi(1) >= ibegi(1)) then - ibegi(2:SWEEPDIM) = ibeg(2:SWEEPDIM,i) - iendi(2:SWEEPDIM) = iend(2:SWEEPDIM,i) - call encounter_check_sweep_aabb_one_single_list(n, ext_ind(:), ibegi(:), iendi(:), ibeg(:,:), iend(:,:), lencounteri(:)) - call encounter_check_make_ragged_list(lencounteri(:), ind_arr(:), lenc(i)%nenc, lenc(i)%index2) - else - lenc(i)%nenc = 0 - end if + if (loverlap(i)) then + ibeg = self%aabb(dim)%ibeg(i) + 1_I8B + iend = self%aabb(dim)%iend(i) - 1_I8B + nbox = iend - ibeg + 1 + lencounteri(ibeg:iend) = .true. + call encounter_check_all_sweep_one(i, nbox, x(1,i), x(2,i), x(3,i), v(1,i), v(2,i), v(3,i), & + xind(ibeg:iend), yind(ibeg:iend), zind(ibeg:iend),& + vxind(ibeg:iend), vyind(ibeg:iend), vzind(ibeg:iend), & + renc(i), rencind(ibeg:iend), dt, ext_ind(dim,ibeg:iend), & + lencounteri(ibeg:iend), lenc(i)%nenc, lenc(i)%index1, lenc(i)%index2, lenc(i)%lvdotr) + end if end do !$omp end parallel do - return - end subroutine encounter_check_sweep_aabb_all_single_list - + call encounter_check_collapse_ragged_list(lenc, n, nenc, index1, index2, lvdotr) - pure subroutine encounter_check_sweep_aabb_one_double_list(i, n1, n2, ntot, ext_ind, ibegi, iendi, ibeg, iend, lencounteri) - !! author: David A. Minton - !! - !! Performs a sweep operation on a single body. Encounters from the same lists not allowed (e.g. pl-tp encounters only) - implicit none - ! Arguments - integer(I4B), intent(in) :: i !! The current index of the ith body - integer(I4B), intent(in) :: n1 !! Number of bodies 1 - integer(I4B), intent(in) :: n2 !! Number of bodies 2 - integer(I4B), intent(in) :: ntot !! n1 + n2 - integer(I4B), dimension(:), intent(in) :: ext_ind !! Sorted index array of extents - integer(I4B), dimension(:), intent(in) :: ibegi, iendi !! The beginning and ending indices of the ith bounding box in the n-dimensions - integer(I4B), dimension(:,:), intent(in) :: ibeg, iend !! Beginning and ending index lists in the n-dimensions - logical, dimension(:), intent(out) :: lencounteri !! Encounter list for the ith body - ! Internals - integer(I4B) :: j, jbox, dim, jlo, jhi - integer(I4B), dimension(:), allocatable :: box - - lencounteri(:) = .false. - jlo = ibegi(1) - jhi = iendi(1) - allocate(box(jlo:jhi)) - box(:) = ext_ind(jlo:jhi) - where(box(:) > ntot) - box(:) = box(:) - ntot - endwhere - - do concurrent (jbox = jlo:jhi) - j = box(jbox) - if (((i <= n1) .and. (j <= n1)) .or. ((i > n1) .and. (j > n1))) cycle ! only pairs from the two different lists allowed - ! Check the other dimensions - lencounteri(j) = all((iend(2:SWEEPDIM,j) > ibegi(2:SWEEPDIM)) .and. (ibeg(2:SWEEPDIM,j) < iendi(2:SWEEPDIM))) + ! By convention, we always assume that index1 < index2, and so we must swap any that are out of order + do concurrent(k = 1_I8B:nenc, index1(k) > index2(k)) + itmp = index1(k) + index1(k) = index2(k) + index2(k) = itmp end do - return - end subroutine encounter_check_sweep_aabb_one_double_list - - - pure subroutine encounter_check_sweep_aabb_one_single_list(n, ext_ind, ibegi, iendi, ibeg, iend, lencounteri) - !! author: David A. Minton - !! - !! Performs a sweep operation on a single body. Mutual encounters allowed (e.g. pl-pl) - implicit none - ! Arguments - integer(I4B), intent(in) :: n !! Number of bodies - integer(I4B), dimension(:), intent(in) :: ext_ind !! Sorted index array of extents - integer(I4B), dimension(:), intent(in) :: ibegi, iendi !! The beginning and ending indices of the ith bounding box in the n-dimensions - integer(I4B), dimension(:,:), intent(in) :: ibeg, iend !! Beginning and ending index lists in the x-dimension - logical, dimension(:), intent(out) :: lencounteri !! Encounter list for the ith body - ! Internals - integer(I4B) :: j, jbox, dim, jlo, jhi - integer(I4B), dimension(:), allocatable :: box - - lencounteri(:) = .false. - - jlo = ibegi(1) - jhi = iendi(1) - - allocate(box(jlo:jhi)) - box(:) = ext_ind(jlo:jhi) - - where(box(:) > n) - box(:) = box(:) - n - endwhere - - do concurrent(jbox = jlo:jhi) ! Sweep forward until the end of the interval - j = box(jbox) - ! Check the other dimensions - lencounteri(j) = all((iend(2:SWEEPDIM,j) > ibegi(2:SWEEPDIM)) .and. (ibeg(2:SWEEPDIM,j) < iendi(2:SWEEPDIM))) - end do + call encounter_check_remove_duplicates(n, nenc, index1, index2, lvdotr) return - end subroutine encounter_check_sweep_aabb_one_single_list + end subroutine encounter_check_sweep_aabb_single_list end submodule s_encounter_check diff --git a/src/encounter/encounter_io.f90 b/src/encounter/encounter_io.f90 index 60fdf5304..d7785a71a 100644 --- a/src/encounter/encounter_io.f90 +++ b/src/encounter/encounter_io.f90 @@ -22,9 +22,9 @@ module subroutine encounter_io_write_frame(iu, t, id1, id2, Gmass1, Gmass2, radi ! Internals character(len=STRMAX) :: errmsg - write(iu, err = 667, iomsg = errmsg) t - write(iu, err = 667, iomsg = errmsg) id1, xh1(1), xh1(2), xh1(3), vh1(1), vh1(2), Gmass1, radius1 - write(iu, err = 667, iomsg = errmsg) id2, xh2(1), xh2(2), xh2(3), vh2(1), vh2(2), Gmass2, radius2 + write(iu, err=667, iomsg=errmsg) t + write(iu, err=667, iomsg=errmsg) id1, xh1(1), xh1(2), xh1(3), vh1(1), vh1(2), Gmass1, radius1 + write(iu, err=667, iomsg=errmsg) id2, xh2(1), xh2(2), xh2(3), vh2(1), vh2(2), Gmass2, radius2 return 667 continue @@ -46,10 +46,10 @@ module subroutine encounter_io_write_list(self, pl, encbody, param) if (param%enc_out == "" .or. self%nenc == 0) return - open(unit = LUN, file = param%enc_out, status = 'OLD', position = 'APPEND', form = 'UNFORMATTED', iostat = ierr, iomsg = errmsg) + open(unit=LUN, file=param%enc_out, status='OLD', position='APPEND', form='UNFORMATTED', iostat=ierr, iomsg=errmsg) if (ierr /= 0) then if (lfirst) then - open(unit = LUN, file = param%enc_out, status = 'NEW', form = 'UNFORMATTED', err = 667, iomsg = errmsg) + open(unit=LUN, file=param%enc_out, status='NEW', form='UNFORMATTED', err=667, iomsg=errmsg) else goto 667 end if diff --git a/src/encounter/encounter_setup.f90 b/src/encounter/encounter_setup.f90 index cb2e2f3d8..85439bd30 100644 --- a/src/encounter/encounter_setup.f90 +++ b/src/encounter/encounter_setup.f90 @@ -52,7 +52,7 @@ module subroutine encounter_setup_list(self, n) implicit none ! Arguments class(encounter_list), intent(inout) :: self !! Swiftest encounter structure - integer(I4B), intent(in) :: n !! Number of encounters to allocate space for + integer(I8B), intent(in) :: n !! Number of encounters to allocate space for if (n < 0) return @@ -69,7 +69,7 @@ module subroutine encounter_setup_list(self, n) if (allocated(self%t)) deallocate(self%t) self%nenc = n - if (n == 0) return + if (n == 0_I8B) return allocate(self%lvdotr(n)) allocate(self%status(n)) diff --git a/src/encounter/encounter_util.f90 b/src/encounter/encounter_util.f90 index c6dc03d7b..fb971f403 100644 --- a/src/encounter/encounter_util.f90 +++ b/src/encounter/encounter_util.f90 @@ -140,27 +140,27 @@ module subroutine encounter_util_resize_list(self, nnew) implicit none ! Arguments class(encounter_list), intent(inout) :: self !! Swiftest encounter list - integer(I4B), intent(in) :: nnew !! New size of list needed + integer(I8B), intent(in) :: nnew !! New size of list needed ! Internals class(encounter_list), allocatable :: enc_temp - integer(I4B) :: nold + integer(I8B) :: nold logical :: lmalloc lmalloc = allocated(self%status) if (lmalloc) then nold = size(self%status) else - nold = 0 + nold = 0_I8B end if if (nnew > nold) then if (lmalloc) allocate(enc_temp, source=self) - call self%setup(2 * nnew) + call self%setup(2_I8B * nnew) if (lmalloc) then call self%copy(enc_temp) deallocate(enc_temp) end if else - self%status(nnew+1:nold) = INACTIVE + self%status(nnew+1_I8B:nold) = INACTIVE end if self%nenc = nnew @@ -179,7 +179,7 @@ module subroutine encounter_util_spill_list(self, discards, lspill_list, ldestru logical, dimension(:), intent(in) :: lspill_list !! Logical array of bodies to spill into the discards logical, intent(in) :: ldestructive !! Logical flag indicating whether or not this operation should alter body by removing the discard list ! Internals - integer(I4B) :: nenc_old + integer(I8B) :: nenc_old associate(keeps => self) call util_spill(keeps%lvdotr, discards%lvdotr, lspill_list, ldestructive) diff --git a/src/fraggle/fraggle_generate.f90 b/src/fraggle/fraggle_generate.f90 index 0b0126702..2b443bebe 100644 --- a/src/fraggle/fraggle_generate.f90 +++ b/src/fraggle/fraggle_generate.f90 @@ -110,14 +110,16 @@ module subroutine fraggle_generate_fragments(self, colliders, system, param, lfa lfailure = ((abs(dEtot + frag%Qloss) > FRAGGLE_ETOL) .or. (dEtot > 0.0_DP)) if (lfailure) then write(message, *) dEtot, abs(dEtot + frag%Qloss) / FRAGGLE_ETOL - call io_log_one_message(FRAGGLE_LOG_OUT, "Fraggle failed due to high energy error: " // trim(adjustl(message))) + call io_log_one_message(FRAGGLE_LOG_OUT, "Fraggle failed due to high energy error: " // & + trim(adjustl(message))) cycle end if lfailure = ((abs(dLmag) / (.mag.frag%Ltot_before)) > FRAGGLE_LTOL) if (lfailure) then write(message,*) dLmag / (.mag.frag%Ltot_before(:)) - call io_log_one_message(FRAGGLE_LOG_OUT, "Fraggle failed due to high angular momentum error: " // trim(adjustl(message))) + call io_log_one_message(FRAGGLE_LOG_OUT, "Fraggle failed due to high angular momentum error: " // & + trim(adjustl(message))) cycle end if @@ -131,9 +133,11 @@ module subroutine fraggle_generate_fragments(self, colliders, system, param, lfa write(message,*) try if (lfailure) then - call io_log_one_message(FRAGGLE_LOG_OUT, "Fraggle fragment generation failed after " // trim(adjustl(message)) // " tries") + call io_log_one_message(FRAGGLE_LOG_OUT, "Fraggle fragment generation failed after " // & + trim(adjustl(message)) // " tries") else - call io_log_one_message(FRAGGLE_LOG_OUT, "Fraggle fragment generation succeeded after " // trim(adjustl(message)) // " tries") + call io_log_one_message(FRAGGLE_LOG_OUT, "Fraggle fragment generation succeeded after " // & + trim(adjustl(message)) // " tries") call fraggle_io_log_generate(frag) end if @@ -240,14 +244,16 @@ subroutine fraggle_generate_spins(frag, colliders, f_spin, lfailure) frag%ke_spin = 0.0_DP do i = 1, nfrag ! Convert a fraction (f_spin) of either the remaining angular momentum or kinetic energy budget into spin, whichever gives the smaller rotation so as not to blow any budgets - rot_ke(:) = sqrt(2 * f_spin * frag%ke_budget / (nfrag * frag%mass(i) * frag%radius(i)**2 * frag%Ip(3, i))) * L_remainder(:) / norm2(L_remainder(:)) + rot_ke(:) = sqrt(2 * f_spin * frag%ke_budget / (nfrag * frag%mass(i) * frag%radius(i)**2 * frag%Ip(3, i))) & + * L_remainder(:) / norm2(L_remainder(:)) rot_L(:) = f_spin * L_remainder(:) / (nfrag * frag%mass(i) * frag%radius(i)**2 * frag%Ip(3, i)) if (norm2(rot_ke) < norm2(rot_L)) then frag%rot(:,i) = rot_ke(:) else frag%rot(:, i) = rot_L(:) end if - frag%ke_spin = frag%ke_spin + frag%mass(i) * frag%Ip(3, i) * frag%radius(i)**2 * dot_product(frag%rot(:, i), frag%rot(:, i)) + frag%ke_spin = frag%ke_spin + frag%mass(i) * frag%Ip(3, i) * frag%radius(i)**2 & + * dot_product(frag%rot(:, i), frag%rot(:, i)) end do frag%ke_spin = 0.5_DP * frag%ke_spin @@ -340,7 +346,8 @@ subroutine fraggle_generate_tan_vel(frag, colliders, lfailure) frag%v_t_mag(1:nfrag) = solve_fragment_tan_vel(v_t_mag_input=v_t_initial(7:nfrag), lfailure=lfailure) ! Perform one final shift of the radial velocity vectors to align with the center of mass of the collisional system (the origin) - frag%vb(:,1:nfrag) = fraggle_util_vmag_to_vb(frag%v_r_mag(1:nfrag), frag%v_r_unit(:,1:nfrag), frag%v_t_mag(1:nfrag), frag%v_t_unit(:,1:nfrag), frag%mass(1:nfrag), frag%vbcom(:)) + frag%vb(:,1:nfrag) = fraggle_util_vmag_to_vb(frag%v_r_mag(1:nfrag), frag%v_r_unit(:,1:nfrag), frag%v_t_mag(1:nfrag), & + frag%v_t_unit(:,1:nfrag), frag%mass(1:nfrag), frag%vbcom(:)) do concurrent (i = 1:nfrag) frag%v_coll(:,i) = frag%vb(:,i) - frag%vbcom(:) end do @@ -506,7 +513,8 @@ subroutine fraggle_generate_rad_vel(frag, colliders, lfailure) ! Shift the radial velocity vectors to align with the center of mass of the collisional system (the origin) frag%ke_orbit = 0.0_DP - frag%vb(:,1:nfrag) = fraggle_util_vmag_to_vb(frag%v_r_mag(1:nfrag), frag%v_r_unit(:,1:nfrag), frag%v_t_mag(1:nfrag), frag%v_t_unit(:,1:nfrag), frag%mass(1:nfrag), frag%vbcom(:)) + frag%vb(:,1:nfrag) = fraggle_util_vmag_to_vb(frag%v_r_mag(1:nfrag), frag%v_r_unit(:,1:nfrag), & + frag%v_t_mag(1:nfrag), frag%v_t_unit(:,1:nfrag), frag%mass(1:nfrag), frag%vbcom(:)) do i = 1, nfrag frag%v_coll(:, i) = frag%vb(:, i) - frag%vbcom(:) frag%ke_orbit = frag%ke_orbit + frag%mass(i) * dot_product(frag%vb(:, i), frag%vb(:, i)) diff --git a/src/fraggle/fraggle_regime.f90 b/src/fraggle/fraggle_regime.f90 index f53888df5..abe5766f0 100644 --- a/src/fraggle/fraggle_regime.f90 +++ b/src/fraggle/fraggle_regime.f90 @@ -50,8 +50,9 @@ module subroutine fraggle_regime_colliders(self, frag, system, param) dentot = sum(mass_si(:) * density_si(:)) / mtot !! Use the positions and velocities of the parents from indside the step (at collision) to calculate the collisional regime - call fraggle_regime_collresolve(Mcb_si, mass_si(jtarg), mass_si(jproj), radius_si(jtarg), radius_si(jproj), x1_si(:), x2_si(:),& - v1_si(:), v2_si(:), density_si(jtarg), density_si(jproj), min_mfrag_si, frag%regime, mlr, mslr, frag%Qloss) + call fraggle_regime_collresolve(Mcb_si, mass_si(jtarg), mass_si(jproj), radius_si(jtarg), radius_si(jproj), & + x1_si(:), x2_si(:), v1_si(:), v2_si(:), density_si(jtarg), density_si(jproj), & + min_mfrag_si, frag%regime, mlr, mslr, frag%Qloss) frag%mass_dist(1) = min(max(mlr, 0.0_DP), mtot) frag%mass_dist(2) = min(max(mslr, 0.0_DP), mtot) @@ -73,7 +74,8 @@ module subroutine fraggle_regime_colliders(self, frag, system, param) end subroutine fraggle_regime_colliders - subroutine fraggle_regime_collresolve(Mcb, m1, m2, rad1, rad2, xh1, xh2, vb1, vb2, den1, den2, min_mfrag, regime, Mlr, Mslr, Qloss) + subroutine fraggle_regime_collresolve(Mcb, m1, m2, rad1, rad2, xh1, xh2, vb1, vb2, den1, den2, min_mfrag, & + regime, Mlr, Mslr, Qloss) !! Author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton !! !! Determine the collisional regime of two colliding bodies. @@ -182,7 +184,8 @@ subroutine fraggle_regime_collresolve(Mcb, m1, m2, rad1, rad2, xh1, xh2, vb1, vb Mlr = Mtot Mslr = 0.0_DP Qloss = 0.0_DP - call io_log_one_message(FRAGGLE_LOG_OUT, "Fragments would have mass below the minimum. Converting this collision into a merger.") + call io_log_one_message(FRAGGLE_LOG_OUT, & + "Fragments would have mass below the minimum. Converting this collision into a merger.") else if( Vimp < Vescp) then regime = COLLRESOLVE_REGIME_MERGE !perfect merging regime diff --git a/src/fraggle/fraggle_set.f90 b/src/fraggle/fraggle_set.f90 index 4012eda34..6902bfb48 100644 --- a/src/fraggle/fraggle_set.f90 +++ b/src/fraggle/fraggle_set.f90 @@ -212,7 +212,8 @@ module subroutine fraggle_set_natural_scale_factors(self, colliders) associate(frag => self) ! Set scale factors - frag%Escale = 0.5_DP * (colliders%mass(1) * dot_product(colliders%vb(:,1), colliders%vb(:,1)) + colliders%mass(2) * dot_product(colliders%vb(:,2), colliders%vb(:,2))) + frag%Escale = 0.5_DP * (colliders%mass(1) * dot_product(colliders%vb(:,1), colliders%vb(:,1)) & + + colliders%mass(2) * dot_product(colliders%vb(:,2), colliders%vb(:,2))) frag%dscale = sum(colliders%radius(:)) frag%mscale = frag%mtot frag%vscale = sqrt(frag%Escale / frag%mscale) diff --git a/src/fraggle/fraggle_util.f90 b/src/fraggle/fraggle_util.f90 index f65bd81c5..f2c082242 100644 --- a/src/fraggle/fraggle_util.f90 +++ b/src/fraggle/fraggle_util.f90 @@ -166,7 +166,8 @@ module subroutine fraggle_util_get_energy_momentum(self, colliders, system, para lexclude(1:npl_after) = .false. lexclude(colliders%idx(1:colliders%ncoll)) = .true. if (.not.allocated(tmpsys)) then - write(*,*) "Error in fraggle_util_get_energy_momentum. This must be called with lbefore=.true. at least once before calling it with lbefore=.false." + write(*,*) "Error in fraggle_util_get_energy_momentum. " // & + " This must be called with lbefore=.true. at least once before calling it with lbefore=.false." call util_exit(FAILURE) end if call fraggle_util_add_fragments_to_system(frag, colliders, tmpsys, tmpparam) @@ -232,7 +233,8 @@ module subroutine fraggle_util_restructure(self, colliders, try, f_spin, r_max_s ! Linearly interpolate the last two failed solution ke deficits to find a new distance value to try ke_tot_deficit = ke_tot_deficit - (frag%ke_budget - frag%ke_orbit - frag%ke_spin) ke_avg_deficit = ke_tot_deficit / try - delta_r = (r_max_start - r_max_start_old) * (ke_avg_deficit_target - ke_avg_deficit_old) / (ke_avg_deficit - ke_avg_deficit_old) + delta_r = (r_max_start - r_max_start_old) * (ke_avg_deficit_target - ke_avg_deficit_old) & + / (ke_avg_deficit - ke_avg_deficit_old) if (abs(delta_r) > delta_r_max) delta_r = sign(delta_r_max, delta_r) end if r_max_start_old = r_max_start diff --git a/src/gr/gr.f90 b/src/gr/gr.f90 index 8b6bb1120..b05465e98 100644 --- a/src/gr/gr.f90 +++ b/src/gr/gr.f90 @@ -28,7 +28,8 @@ module pure subroutine gr_kick_getaccb_ns_body(self, system, param) rmag = norm2(self%xh(:,i)) vmag2 = dot_product(self%vh(:,i), self%vh(:,i)) rdotv = dot_product(self%xh(:,i), self%vh(:,i)) - self%agr(:, i) = self%mu * inv_c2 / rmag**3 * ((4 * self%mu(i) / rmag - vmag2) * self%xh(:,i) + 4 * rdotv * self%vh(:,i)) + self%agr(:, i) = self%mu * inv_c2 / rmag**3 * ((4 * self%mu(i) / rmag - vmag2) & + * self%xh(:,i) + 4 * rdotv * self%vh(:,i)) end do select type(self) diff --git a/src/io/io.f90 b/src/io/io.f90 index 7ec71d9ae..18a22ccb9 100644 --- a/src/io/io.f90 +++ b/src/io/io.f90 @@ -29,10 +29,11 @@ module subroutine io_conservation_report(self, param, lterminal) associate(system => self, pl => self%pl, cb => self%cb, npl => self%pl%nbody) if (((param%out_type == REAL4_TYPE) .or. (param%out_type == REAL8_TYPE)) .and. (param%energy_out /= "")) then if (param%lfirstenergy .and. (param%out_stat /= "OLD")) then - open(unit = EGYIU, file = param%energy_out, form = "formatted", status = "replace", action = "write", err = 667, iomsg = errmsg) - write(EGYIU,EGYHEADER, err = 667, iomsg = errmsg) + open(unit=EGYIU, file=param%energy_out, form="formatted", status="replace", action="write", err=667, iomsg=errmsg) + write(EGYIU,EGYHEADER, err=667, iomsg=errmsg) else - open(unit = EGYIU, file = param%energy_out, form = "formatted", status = "old", action = "write", position = "append", err = 667, iomsg = errmsg) + open(unit=EGYIU, file=param%energy_out, form="formatted", status="old", action="write", & + position="append", err=667, iomsg=errmsg) end if end if @@ -177,9 +178,9 @@ module subroutine io_dump_particle_info_base(self, param, idx) if (lfirst) then select case(param%out_stat) case('APPEND') - open(unit = LUN, file = param%particle_out, status = 'OLD', position = 'APPEND', form = 'UNFORMATTED', err = 667, iomsg = errmsg) + open(unit=LUN, file=param%particle_out, status='OLD', position='APPEND', form='UNFORMATTED', err=667, iomsg=errmsg) case('NEW', 'UNKNOWN', 'REPLACE') - open(unit = LUN, file = param%particle_out, status = param%out_stat, form = 'UNFORMATTED', err = 667, iomsg = errmsg) + open(unit=LUN, file=param%particle_out, status=param%out_stat, form='UNFORMATTED', err=667, iomsg=errmsg) case default write(*,*) 'Invalid status code',trim(adjustl(param%out_stat)) call util_exit(FAILURE) @@ -187,7 +188,7 @@ module subroutine io_dump_particle_info_base(self, param, idx) lfirst = .false. else - open(unit = LUN, file = param%particle_out, status = 'OLD', position = 'APPEND', form = 'UNFORMATTED', err = 667, iomsg = errmsg) + open(unit=LUN, file=param%particle_out, status='OLD', position= 'APPEND', form='UNFORMATTED', err=667, iomsg=errmsg) end if select type(self) @@ -760,7 +761,8 @@ module subroutine io_param_reader(self, unit, iotype, v_list, iostat, iomsg) iostat = -1 return end if - if ((param%in_type /= REAL8_TYPE) .and. (param%in_type /= "ASCII") .and. (param%in_type /= NETCDF_FLOAT_TYPE) .and. (param%in_type /= NETCDF_DOUBLE_TYPE)) then + if ((param%in_type /= REAL8_TYPE) .and. (param%in_type /= "ASCII") & + .and. (param%in_type /= NETCDF_FLOAT_TYPE) .and. (param%in_type /= NETCDF_DOUBLE_TYPE)) then write(iomsg,*) 'Invalid input file type:',trim(adjustl(param%in_type)) iostat = -1 return @@ -788,7 +790,8 @@ module subroutine io_param_reader(self, unit, iotype, v_list, iostat, iomsg) iostat = -1 return end if - if ((param%out_stat /= "NEW") .and. (param%out_stat /= "REPLACE") .and. (param%out_stat /= "APPEND") .and. (param%out_stat /= "UNKNOWN")) then + if ((param%out_stat /= "NEW") .and. (param%out_stat /= "REPLACE") .and. (param%out_stat /= "APPEND") & + .and. (param%out_stat /= "UNKNOWN")) then write(iomsg,*) 'Invalid out_stat: ',trim(adjustl(param%out_stat)) iostat = -1 return @@ -1956,9 +1959,9 @@ module subroutine io_write_discard(self, param) end if select case(out_stat) case('APPEND') - open(unit = LUN, file = param%discard_out, status = 'OLD', position = 'APPEND', form = 'FORMATTED', err = 667, iomsg = errmsg) + open(unit=LUN, file=param%discard_out, status='OLD', position='APPEND', form='FORMATTED', err=667, iomsg=errmsg) case('NEW', 'REPLACE', 'UNKNOWN') - open(unit = LUN, file = param%discard_out, status = param%out_stat, form = 'FORMATTED', err = 667, iomsg = errmsg) + open(unit=LUN, file=param%discard_out, status=param%out_stat, form='FORMATTED', err=667, iomsg=errmsg) case default write(*,*) 'Invalid status code for OUT_STAT: ',trim(adjustl(param%out_stat)) call util_exit(FAILURE) @@ -2137,9 +2140,9 @@ module subroutine io_write_frame_system(self, param) if (lfirst) then select case(param%out_stat) case('APPEND') - open(unit = iu, file = param%outfile, status = 'OLD', position = 'APPEND', form = 'UNFORMATTED', err = 667, iomsg = errmsg) + open(unit=iu, file=param%outfile, status='OLD', position='APPEND', form='UNFORMATTED', err=667, iomsg=errmsg) case('NEW', 'REPLACE', 'UNKNOWN') - open(unit = iu, file = param%outfile, status = param%out_stat, form = 'UNFORMATTED', err = 667, iomsg = errmsg) + open(unit=iu, file=param%outfile, status=param%out_stat, form='UNFORMATTED', err=667, iomsg=errmsg) case default write(*,*) 'Invalid status code for OUT_STAT: ',trim(adjustl(param%out_stat)) call util_exit(FAILURE) @@ -2147,7 +2150,7 @@ module subroutine io_write_frame_system(self, param) lfirst = .false. else - open(unit = iu, file = param%outfile, status = 'OLD', position = 'APPEND', form = 'UNFORMATTED', err = 667, iomsg = errmsg) + open(unit=iu, file=param%outfile, status='OLD', position= 'APPEND', form='UNFORMATTED', err=667, iomsg=errmsg) end if else if ((param%out_type == NETCDF_FLOAT_TYPE) .or. (param%out_type == NETCDF_DOUBLE_TYPE)) then diff --git a/src/kick/kick.f90 b/src/kick/kick.f90 index c3e37d927..e680fde98 100644 --- a/src/kick/kick.f90 +++ b/src/kick/kick.f90 @@ -190,7 +190,7 @@ module subroutine kick_getacch_int_all_triangular_pl(npl, nplm, x, Gmass, radius !$omp end parallel do else !$omp parallel do default(private) schedule(static)& - !$omp shared(npl, nplm, x, Gmass, radius) & + !$omp shared(npl, nplm, x, Gmass) & !$omp lastprivate(rji2, xr, yr, zr) & !$omp reduction(+:ahi) & !$omp reduction(-:ahj) diff --git a/src/main/swiftest_driver.f90 b/src/main/swiftest_driver.f90 index d4e1f2dc8..afc62fddf 100644 --- a/src/main/swiftest_driver.f90 +++ b/src/main/swiftest_driver.f90 @@ -22,8 +22,10 @@ program swiftest_driver real(DP) :: old_t_final = 0.0_DP !! Output time at which writing should start, in order to prevent duplicate lines being written for restarts type(walltimer) :: integration_timer !! Object used for computing elapsed wall time real(DP) :: tfrac - character(*), parameter :: statusfmt = '("Time = ", ES12.5, "; fraction done = ", F6.3, "; Number of active pl, tp = ", I5, ", ", I5)' - character(*), parameter :: symbastatfmt = '("Time = ", ES12.5, "; fraction done = ", F6.3, "; Number of active plm, pl, tp = ", I5, ", ", I5, ", ", I5)' + character(*), parameter :: statusfmt = '("Time = ", ES12.5, "; fraction done = ", F6.3, ' // & + '"; Number of active pl, tp = ", I5, ", ", I5)' + character(*), parameter :: symbastatfmt = '("Time = ", ES12.5, "; fraction done = ", F6.3, ' // & + '"; Number of active plm, pl, tp = ", I5, ", ", I5, ", ", I5)' ierr = io_get_args(integrator, param_file_name) if (ierr /= 0) then @@ -77,7 +79,6 @@ program swiftest_driver !$ write(*,'(a,i3,/)') ' Number of threads = ', nthreads write(*, *) " *************** Main Loop *************** " if (param%lrestart .and. param%lenergy) call nbody_system%conservation_report(param, lterminal=.true.) - call integration_timer%reset() do iloop = 1, nloops !> Step the system forward in time call integration_timer%start() @@ -105,7 +106,7 @@ program swiftest_driver write(*, statusfmt) param%t, tfrac, pl%nbody, nbody_system%tp%nbody end select if (param%lenergy) call nbody_system%conservation_report(param, lterminal=.true.) - call integration_timer%report(nsubsteps=istep_dump, message="Integration steps:") + call integration_timer%report(message="Integration steps:", nsubsteps=istep_dump) call integration_timer%reset() iout = istep_out diff --git a/src/modules/encounter_classes.f90 b/src/modules/encounter_classes.f90 index 1c182619b..f7758ef54 100644 --- a/src/modules/encounter_classes.f90 +++ b/src/modules/encounter_classes.f90 @@ -10,7 +10,7 @@ module encounter_classes integer(I4B), parameter :: SWEEPDIM = 3 type :: encounter_list - integer(I4B) :: nenc = 0 !! Total number of encounters + integer(I8B) :: nenc = 0 !! Total number of encounters logical, dimension(:), allocatable :: lvdotr !! relative vdotr flag integer(I4B), dimension(:), allocatable :: status !! status of the interaction integer(I4B), dimension(:), allocatable :: index1 !! position of the first body in the encounter @@ -35,9 +35,9 @@ module encounter_classes type encounter_bounding_box_1D integer(I4B) :: n !! Number of bodies with extents - integer(I4B), dimension(:), allocatable :: ind !! Sorted minimum/maximum extent indices - integer(I4B), dimension(:), allocatable :: ibeg !! Beginning index for box - integer(I4B), dimension(:), allocatable :: iend !! Ending index for box + integer(I4B), dimension(:), allocatable :: ind !! Sorted minimum/maximum extent indices (value > n indicates an ending index) + integer(I8B), dimension(:), allocatable :: ibeg !! Beginning index for box + integer(I8B), dimension(:), allocatable :: iend !! Ending index for box contains procedure :: sort => encounter_check_sort_aabb_1D !! Sorts the bounding box extents along a single dimension prior to the sweep phase procedure :: dealloc => encounter_util_dealloc_aabb !! Deallocates all allocatables @@ -54,21 +54,7 @@ module encounter_classes end type interface - module subroutine encounter_check_all(nenc, index1, index2, x1, v1, x2, v2, renc1, renc2, dt, lencounter, lvdotr) - implicit none - integer(I4B), intent(in) :: nenc !! Number of encounters in the encounter lists - integer(I4B), dimension(:), intent(in) :: index1 !! List of indices for body 1 in each encounter - integer(I4B), dimension(:), intent(in) :: index2 !! List of indices for body 2 in each encounter1 - real(DP), dimension(:,:), intent(in) :: x1, v1 !! Array of indices of bodies 1 - real(DP), dimension(:,:), intent(in) :: x2, v2 !! Array of indices of bodies 2 - real(DP), dimension(:), intent(in) :: renc1 !! Radius of encounter regions of bodies 1 - real(DP), dimension(:), intent(in) :: renc2 !! Radius of encounter regions of bodies 2 - real(DP), intent(in) :: dt !! Step size - logical, dimension(:), intent(out) :: lencounter !! Logical array indicating which pairs are in a close encounter state - logical, dimension(:), intent(out) :: lvdotr !! Logical array indicating which pairs are approaching - end subroutine encounter_check_all - - module subroutine encounter_check_all_plpl(param, npl, x, v, renc, dt, lvdotr, index1, index2, nenc) + module subroutine encounter_check_all_plpl(param, npl, x, v, renc, dt, nenc, index1, index2, lvdotr) use swiftest_classes, only: swiftest_parameters implicit none class(swiftest_parameters), intent(inout) :: param !! Current Swiftest run configuration parameter5s @@ -80,10 +66,11 @@ module subroutine encounter_check_all_plpl(param, npl, x, v, renc, dt, lvdotr, i logical, dimension(:), allocatable, intent(out) :: lvdotr !! Logical flag indicating the sign of v .dot. x integer(I4B), dimension(:), allocatable, intent(out) :: index1 !! List of indices for body 1 in each encounter integer(I4B), dimension(:), allocatable, intent(out) :: index2 !! List of indices for body 2 in each encounter - integer(I4B), intent(out) :: nenc !! Total number of encounters + integer(I8B), intent(out) :: nenc !! Total number of encounters end subroutine encounter_check_all_plpl - module subroutine encounter_check_all_plplm(param, nplm, nplt, xplm, vplm, xplt, vplt, rencm, renct, dt, lvdotr, index1, index2, nenc) + module subroutine encounter_check_all_plplm(param, nplm, nplt, xplm, vplm, xplt, vplt, rencm, renct, dt, & + nenc, index1, index2, lvdotr) use swiftest_classes, only: swiftest_parameters implicit none class(swiftest_parameters), intent(inout) :: param !! Current Swiftest run configuration parameter5s @@ -96,13 +83,13 @@ module subroutine encounter_check_all_plplm(param, nplm, nplt, xplm, vplm, xplt, real(DP), dimension(:), intent(in) :: rencm !! Critical radii of fully interacting massive bodies that defines an encounter real(DP), dimension(:), intent(in) :: renct !! Critical radii of partially interacting massive bodies that defines an encounter real(DP), intent(in) :: dt !! Step size - logical, dimension(:), allocatable, intent(out) :: lvdotr !! Logical flag indicating the sign of v .dot. x + integer(I8B), intent(out) :: nenc !! Total number of encounters integer(I4B), dimension(:), allocatable, intent(out) :: index1 !! List of indices for body 1 in each encounter integer(I4B), dimension(:), allocatable, intent(out) :: index2 !! List of indices for body 2 in each encounter - integer(I4B), intent(out) :: nenc !! Total number of encounters + logical, dimension(:), allocatable, intent(out) :: lvdotr !! Logical flag indicating the sign of v .dot. x end subroutine encounter_check_all_plplm - module subroutine encounter_check_all_pltp(param, npl, ntp, xpl, vpl, xtp, vtp, renc, dt, lvdotr, index1, index2, nenc) + module subroutine encounter_check_all_pltp(param, npl, ntp, xpl, vpl, xtp, vtp, renc, dt, nenc, index1, index2, lvdotr) use swiftest_classes, only: swiftest_parameters implicit none class(swiftest_parameters), intent(inout) :: param !! Current Swiftest run configuration parameter5s @@ -114,13 +101,13 @@ module subroutine encounter_check_all_pltp(param, npl, ntp, xpl, vpl, xtp, vtp, real(DP), dimension(:,:), intent(in) :: vtp !! Velocity vectors of massive bodies real(DP), dimension(:), intent(in) :: renc !! Critical radii of massive bodies that defines an encounter real(DP), intent(in) :: dt !! Step size - logical, dimension(:), allocatable, intent(out) :: lvdotr !! Logical flag indicating the sign of v .dot. x + integer(I8B), intent(out) :: nenc !! Total number of encounters integer(I4B), dimension(:), allocatable, intent(out) :: index1 !! List of indices for body 1 in each encounter integer(I4B), dimension(:), allocatable, intent(out) :: index2 !! List of indices for body 2 in each encounter - integer(I4B), intent(out) :: nenc !! Total number of encounters + logical, dimension(:), allocatable, intent(out) :: lvdotr !! Logical flag indicating the sign of v .dot. x end subroutine encounter_check_all_pltp - module pure subroutine encounter_check_one(xr, yr, zr, vxr, vyr, vzr, renc, dt, lencounter, lvdotr) + module elemental subroutine encounter_check_one(xr, yr, zr, vxr, vyr, vzr, renc, dt, lencounter, lvdotr) !$omp declare simd(encounter_check_one) implicit none real(DP), intent(in) :: xr, yr, zr !! Relative distance vector components @@ -135,7 +122,7 @@ module subroutine encounter_check_collapse_ragged_list(ragged_list, n1, nenc, in implicit none type(encounter_list), dimension(:), intent(in) :: ragged_list !! The ragged encounter list integer(I4B), intent(in) :: n1 !! Number of bodies 1 - integer(I4B), intent(out) :: nenc !! Total number of encountersj + integer(I8B), intent(out) :: nenc !! Total number of encountersj integer(I4B), dimension(:), allocatable, intent(out) :: index1 !! Array of indices for body 1 integer(I4B), dimension(:), allocatable, intent(out) :: index2 !! Array of indices for body 1 logical, dimension(:), allocatable, intent(out), optional :: lvdotr !! Array indicating which bodies are approaching @@ -148,23 +135,34 @@ module pure subroutine encounter_check_sort_aabb_1D(self, n, extent_arr) real(DP), dimension(:), intent(in) :: extent_arr !! Array of extents of size 2*n end subroutine encounter_check_sort_aabb_1D - module subroutine encounter_check_sweep_aabb_double_list(self, n1, n2, nenc, index1, index2) - implicit none - class(encounter_bounding_box), intent(inout) :: self !! Multi-dimensional bounding box structure - integer(I4B), intent(in) :: n1 !! Number of bodies 1 - integer(I4B), intent(in) :: n2 !! Number of bodies 2 - integer(I4B), intent(out) :: nenc !! Total number of encounter candidates - integer(I4B), dimension(:), allocatable, intent(out) :: index1 !! List of indices for body 1 in each encounter candidate pair - integer(I4B), dimension(:), allocatable, intent(out) :: index2 !! List of indices for body 2 in each encounter candidate pair + module subroutine encounter_check_sweep_aabb_double_list(self, n1, n2, x1, v1, x2, v2, renc1, renc2, dt, & + nenc, index1, index2, lvdotr) + implicit none + class(encounter_bounding_box), intent(inout) :: self !! Multi-dimensional bounding box structure + integer(I4B), intent(in) :: n1 !! Number of bodies 1 + integer(I4B), intent(in) :: n2 !! Number of bodies 2 + real(DP), dimension(:,:), intent(in) :: x1, v1 !! Array of indices of bodies 1 + real(DP), dimension(:,:), intent(in) :: x2, v2 !! Array of indices of bodies 2 + real(DP), dimension(:), intent(in) :: renc1 !! Radius of encounter regions of bodies 1 + real(DP), dimension(:), intent(in) :: renc2 !! Radius of encounter regions of bodies 2 + real(DP), intent(in) :: dt !! Step size + integer(I8B), intent(out) :: nenc !! Total number of encounter candidates + integer(I4B), dimension(:), allocatable, intent(out) :: index1 !! List of indices for body 1 in each encounter candidate pair + integer(I4B), dimension(:), allocatable, intent(out) :: index2 !! List of indices for body 2 in each encounter candidate pair + logical, dimension(:), allocatable, intent(out) :: lvdotr !! Logical array indicating which pairs are approaching end subroutine encounter_check_sweep_aabb_double_list - module subroutine encounter_check_sweep_aabb_single_list(self, n, nenc, index1, index2) - implicit none - class(encounter_bounding_box), intent(inout) :: self !! Multi-dimensional bounding box structure - integer(I4B), intent(in) :: n !! Number of bodies 1 - integer(I4B), intent(out) :: nenc !! Total number of encounter candidates - integer(I4B), dimension(:), allocatable, intent(out) :: index1 !! List of indices for body 1 in each encounter candidate pair - integer(I4B), dimension(:), allocatable, intent(out) :: index2 !! List of indices for body 2 in each encounter candidate pair + module subroutine encounter_check_sweep_aabb_single_list(self, n, x, v, renc, dt, nenc, index1, index2, lvdotr) + implicit none + class(encounter_bounding_box), intent(inout) :: self !! Multi-dimensional bounding box structure + integer(I4B), intent(in) :: n !! Number of bodies + real(DP), dimension(:,:), intent(in) :: x, v !! Array of position and velocity vectors + real(DP), dimension(:), intent(in) :: renc !! Radius of encounter regions of bodies 1 + real(DP), intent(in) :: dt !! Step size + integer(I8B), intent(out) :: nenc !! Total number of encounter candidates + integer(I4B), dimension(:), allocatable, intent(out) :: index1 !! List of indices for one body in each encounter candidate pair + integer(I4B), dimension(:), allocatable, intent(out) :: index2 !! List of indices for the other body in each encounter candidate pair + logical, dimension(:), allocatable, intent(out) :: lvdotr !! Logical array indicating which pairs are approaching end subroutine encounter_check_sweep_aabb_single_list module subroutine encounter_io_write_frame(iu, t, id1, id2, Gmass1, Gmass2, radius1, radius2, xh1, xh2, vh1, vh2) @@ -181,7 +179,7 @@ end subroutine encounter_io_write_frame module subroutine encounter_io_write_list(self, pl, encbody, param) use swiftest_classes, only : swiftest_pl, swiftest_body, swiftest_parameters implicit none - class(encounter_list), intent(in) :: self !! Swiftest encounter list object + class(encounter_list), intent(in) :: self !! Swiftest encounter list object class(swiftest_pl), intent(in) :: pl !! Swiftest massive body object class(swiftest_body), intent(in) :: encbody !! Encountering body - Swiftest generic body object (pl or tp) class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters @@ -197,14 +195,14 @@ end subroutine encounter_setup_aabb module subroutine encounter_setup_list(self, n) implicit none class(encounter_list), intent(inout) :: self !! Swiftest encounter structure - integer(I4B), intent(in) :: n !! Number of encounters to allocate space for + integer(I8B), intent(in) :: n !! Number of encounters to allocate space for end subroutine encounter_setup_list module subroutine encounter_util_append_list(self, source, lsource_mask) implicit none class(encounter_list), intent(inout) :: self !! Swiftest encounter list object class(encounter_list), intent(in) :: source !! Source object to append - logical, dimension(:), 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 encounter_util_append_list module subroutine encounter_util_copy_list(self, source) @@ -236,15 +234,15 @@ end subroutine encounter_util_final_list module subroutine encounter_util_resize_list(self, nnew) implicit none class(encounter_list), intent(inout) :: self !! Swiftest encounter list - integer(I4B), intent(in) :: nnew !! New size of list needed + integer(I8B), intent(in) :: nnew !! New size of list needed end subroutine encounter_util_resize_list module subroutine encounter_util_spill_list(self, discards, lspill_list, ldestructive) implicit none class(encounter_list), intent(inout) :: self !! Swiftest encounter list class(encounter_list), intent(inout) :: discards !! Discarded object - logical, dimension(:), intent(in) :: lspill_list !! Logical array of bodies to spill into the discards - logical, intent(in) :: ldestructive !! Logical flag indicating whether or not this operation should alter body by removing the discard list + logical, dimension(:), intent(in) :: lspill_list !! Logical array of bodies to spill into the discards + logical, intent(in) :: ldestructive !! Logical flag indicating whether or not this operation should alter body by removing the discard list end subroutine encounter_util_spill_list end interface diff --git a/src/modules/swiftest_classes.f90 b/src/modules/swiftest_classes.f90 index a1489f82c..0cc6eb4f2 100644 --- a/src/modules/swiftest_classes.f90 +++ b/src/modules/swiftest_classes.f90 @@ -1576,7 +1576,8 @@ module subroutine util_set_mu_tp(self, cb) class(swiftest_cb), intent(inout) :: cb !! Swiftest central body object end subroutine util_set_mu_tp - module subroutine util_set_particle_info(self, name, particle_type, status, origin_type, origin_time, collision_id, origin_xh, origin_vh, discard_time, discard_xh, discard_vh, discard_body_id) + module subroutine util_set_particle_info(self, name, particle_type, status, origin_type, origin_time, collision_id, & + origin_xh, origin_vh, discard_time, discard_xh, discard_vh, discard_body_id) implicit none class(swiftest_particle_info), intent(inout) :: self character(len=*), intent(in), optional :: name !! Non-unique name @@ -1663,6 +1664,18 @@ module pure subroutine util_sort_index_i4b(arr,ind) integer(I4B), dimension(:), allocatable, intent(inout) :: ind end subroutine util_sort_index_i4b + module pure subroutine util_sort_index_I4B_I8Bind(arr,ind) + implicit none + integer(I4B), dimension(:), intent(in) :: arr + integer(I8B), dimension(:), allocatable, intent(inout) :: ind + end subroutine util_sort_index_I4b_I8Bind + + module pure subroutine util_sort_index_I8B_I8Bind(arr,ind) + implicit none + integer(I8B), dimension(:), intent(in) :: arr + integer(I8B), dimension(:), allocatable, intent(inout) :: ind + end subroutine util_sort_index_I8B_I8Bind + module pure subroutine util_sort_sp(arr) implicit none real(SP), dimension(:), intent(inout) :: arr @@ -1715,6 +1728,13 @@ module pure subroutine util_sort_rearrange_arr_I4B(arr, ind, n) integer(I4B), intent(in) :: n !! Number of elements in arr and ind to rearrange end subroutine util_sort_rearrange_arr_I4B + module pure subroutine util_sort_rearrange_arr_I4B_I8Bind(arr, ind, n) + implicit none + integer(I4B), dimension(:), allocatable, intent(inout) :: arr !! Destination array + integer(I8B), dimension(:), intent(in) :: ind !! Index to rearrange against + integer(I8B), intent(in) :: n !! Number of elements in arr and ind to rearrange + end subroutine util_sort_rearrange_arr_I4B_I8Bind + module subroutine util_sort_rearrange_arr_info(arr, ind, n) implicit none type(swiftest_particle_info), dimension(:), allocatable, intent(inout) :: arr !! Destination array @@ -1728,6 +1748,13 @@ module pure subroutine util_sort_rearrange_arr_logical(arr, ind, n) integer(I4B), dimension(:), intent(in) :: ind !! Index to rearrange against integer(I4B), intent(in) :: n !! Number of elements in arr and ind to rearrange end subroutine util_sort_rearrange_arr_logical + + module pure subroutine util_sort_rearrange_arr_logical_I8Bind(arr, ind, n) + implicit none + logical, dimension(:), allocatable, intent(inout) :: arr !! Destination array + integer(I8B), dimension(:), intent(in) :: ind !! Index to rearrange against + integer(I8B), intent(in) :: n !! Number of elements in arr and ind to rearrange + end subroutine util_sort_rearrange_arr_logical_I8Bind end interface util_sort_rearrange interface diff --git a/src/modules/symba_classes.f90 b/src/modules/symba_classes.f90 index 99c527143..a01d60e6a 100644 --- a/src/modules/symba_classes.f90 +++ b/src/modules/symba_classes.f90 @@ -450,7 +450,7 @@ end subroutine symba_setup_pl module subroutine symba_setup_encounter_list(self,n) implicit none class(symba_encounter), intent(inout) :: self !! SyMBA pl-tp encounter structure - integer(I4B), intent(in) :: n !! Number of encounters to allocate space for + integer(I8B), intent(in) :: n !! Number of encounters to allocate space for end subroutine symba_setup_encounter_list module subroutine symba_setup_tp(self, n, param) diff --git a/src/modules/walltime_classes.f90 b/src/modules/walltime_classes.f90 index 7e6c971b4..4bd337ef2 100644 --- a/src/modules/walltime_classes.f90 +++ b/src/modules/walltime_classes.f90 @@ -51,11 +51,11 @@ module walltime_classes end type interaction_timer interface - module subroutine walltime_report(self, nsubsteps, message) + module subroutine walltime_report(self, message, nsubsteps) implicit none class(walltimer), intent(inout) :: self !! Walltimer object - integer(I4B), intent(in) :: nsubsteps !! Number of substeps used to compute the time per step character(len=*), intent(in) :: message !! Message to prepend to the wall time terminal output + integer(I4B), optional, intent(in) :: nsubsteps !! Number of substeps used to compute the time per step end subroutine walltime_report module subroutine walltime_reset(self) diff --git a/src/rmvs/rmvs_discard.f90 b/src/rmvs/rmvs_discard.f90 index b179178dc..8e3b1d7f0 100644 --- a/src/rmvs/rmvs_discard.f90 +++ b/src/rmvs/rmvs_discard.f90 @@ -30,11 +30,12 @@ module subroutine rmvs_discard_tp(self, system, param) 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)) + // ") 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. - call tp%info(i)%set_value(status="DISCARDED_PLQ", discard_time=t, discard_xh=tp%xh(:,i), discard_vh=tp%vh(:,i), discard_body_id=pl%id(iplperP)) + call tp%info(i)%set_value(status="DISCARDED_PLQ", discard_time=t, discard_xh=tp%xh(:,i), & + discard_vh=tp%vh(:,i), discard_body_id=pl%id(iplperP)) end if end if end associate diff --git a/src/rmvs/rmvs_encounter_check.f90 b/src/rmvs/rmvs_encounter_check.f90 index 68fe7a2ef..bbb5a4382 100644 --- a/src/rmvs/rmvs_encounter_check.f90 +++ b/src/rmvs/rmvs_encounter_check.f90 @@ -18,7 +18,8 @@ module function rmvs_encounter_check_tp(self, param, system, dt) result(lencount ! Result logical :: lencounter !! Returns true if there is at least one close encounter ! Internals - integer(I4B) :: i, j, nenc + integer(I4B) :: i, j + integer(I8B) :: nenc real(DP) :: xr, yr, zr, vxr, vyr, vzr real(DP), dimension(system%pl%nbody) :: rcrit logical :: lflag @@ -35,11 +36,12 @@ module function rmvs_encounter_check_tp(self, param, system, dt) result(lencount class is (rmvs_pl) associate(tp => self, ntp => self%nbody, npl => pl%nbody) tp%plencP(1:ntp) = 0 - call encounter_check_all_pltp(param, npl, ntp, pl%xbeg, pl%vbeg, tp%xh, tp%vh, pl%renc, dt, lvdotr, index1, index2, nenc) + call encounter_check_all_pltp(param, npl, ntp, pl%xbeg, pl%vbeg, tp%xh, tp%vh, pl%renc, dt, & + nenc, index1, index2, lvdotr) - lencounter = (nenc > 0) + lencounter = (nenc > 0_I8B) if (lencounter) then - tp%plencP(index2(1:nenc)) = index1(1:nenc) + tp%plencP(index2(1_I8B:nenc)) = index1(1_I8B:nenc) do j = 1, npl pl%nenc(j) = count(tp%plencP(1:ntp) == j) end do diff --git a/src/rmvs/rmvs_step.f90 b/src/rmvs/rmvs_step.f90 index 673188f0e..76c5a802c 100644 --- a/src/rmvs/rmvs_step.f90 +++ b/src/rmvs/rmvs_step.f90 @@ -483,8 +483,10 @@ subroutine rmvs_make_planetocentric(param, cb, pl, tp) 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) + 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) @@ -554,7 +556,8 @@ subroutine rmvs_peri_tp(tp, pl, t, dt, lfirst, inner_index, ipleP, param) id2 = tp%id(i) xh2(:) = xpc(:, i) + xh1(:) vh2(:) = xpc(:, i) + vh1(:) - call rmvs_io_write_encounter(t, id1, id2, mu, 0.0_DP, rpl, 0.0_DP, xh1(:), xh2(:), vh1(:), vh2(:), param%enc_out) + call rmvs_io_write_encounter(t, id1, id2, mu, 0.0_DP, rpl, 0.0_DP, & + xh1(:), xh2(:), vh1(:), vh2(:), param%enc_out) end if if (tp%lperi(i)) then if (peri < tp%peri(i)) then diff --git a/src/setup/setup.f90 b/src/setup/setup.f90 index 3f1932d24..e6678effa 100644 --- a/src/setup/setup.f90 +++ b/src/setup/setup.f90 @@ -106,12 +106,15 @@ module subroutine setup_initialize_particle_info_system(self, param) associate(cb => self%cb, pl => self%pl, npl => self%pl%nbody, tp => self%tp, ntp => self%tp%nbody) - call cb%info%set_value(particle_type=CB_TYPE_NAME, status="ACTIVE", origin_type="Initial conditions", origin_time=param%t0, origin_xh=[0.0_DP, 0.0_DP, 0.0_DP], origin_vh=[0.0_DP, 0.0_DP, 0.0_DP]) + call cb%info%set_value(particle_type=CB_TYPE_NAME, status="ACTIVE", origin_type="Initial conditions", & + origin_time=param%t0, origin_xh=[0.0_DP, 0.0_DP, 0.0_DP], origin_vh=[0.0_DP, 0.0_DP, 0.0_DP]) do i = 1, self%pl%nbody - call pl%info(i)%set_value(particle_type=PL_TYPE_NAME, status="ACTIVE", origin_type="Initial conditions", origin_time=param%t0, origin_xh=self%pl%xh(:,i), origin_vh=self%pl%vh(:,i)) + call pl%info(i)%set_value(particle_type=PL_TYPE_NAME, status="ACTIVE", origin_type="Initial conditions", & + origin_time=param%t0, origin_xh=self%pl%xh(:,i), origin_vh=self%pl%vh(:,i)) end do do i = 1, self%tp%nbody - call tp%info(i)%set_value(particle_type=TP_TYPE_NAME, status="ACTIVE", origin_type="Initial conditions", origin_time=param%t0, origin_xh=self%tp%xh(:,i), origin_vh=self%tp%vh(:,i)) + call tp%info(i)%set_value(particle_type=TP_TYPE_NAME, status="ACTIVE", origin_type="Initial conditions", & + origin_time=param%t0, origin_xh=self%tp%xh(:,i), origin_vh=self%tp%vh(:,i)) end do end associate diff --git a/src/symba/symba_collision.f90 b/src/symba/symba_collision.f90 index fd0d569d1..31be12b0b 100644 --- a/src/symba/symba_collision.f90 +++ b/src/symba/symba_collision.f90 @@ -312,7 +312,8 @@ module function symba_collision_check_encounter(self, system, param, t, dt, irec vr(:) = pl%vb(:, i) - pl%vb(:, j) rlim = pl%radius(i) + pl%radius(j) Gmtot = pl%Gmass(i) + pl%Gmass(j) - lcollision(k) = symba_collision_check_one(xr(1), xr(2), xr(3), vr(1), vr(2), vr(3), Gmtot, rlim, dt, self%lvdotr(k)) + lcollision(k) = symba_collision_check_one(xr(1), xr(2), xr(3), vr(1), vr(2), vr(3), & + Gmtot, rlim, dt, self%lvdotr(k)) end do else do concurrent(k = 1:nenc, lmask(k)) @@ -320,7 +321,8 @@ module function symba_collision_check_encounter(self, system, param, t, dt, irec j = self%index2(k) xr(:) = pl%xh(:, i) - tp%xh(:, j) vr(:) = pl%vb(:, i) - tp%vb(:, j) - lcollision(k) = symba_collision_check_one(xr(1), xr(2), xr(3), vr(1), vr(2), vr(3), pl%Gmass(i), pl%radius(i), dt, self%lvdotr(k)) + lcollision(k) = symba_collision_check_one(xr(1), xr(2), xr(3), vr(1), vr(2), vr(3), & + pl%Gmass(i), pl%radius(i), dt, self%lvdotr(k)) end do end if @@ -532,7 +534,9 @@ function symba_collision_consolidate_colliders(pl, cb, param, idx_parent, collid xcrossv(:) = xc(:) .cross. vc(:) colliders%L_spin(:, j) = colliders%L_spin(:, j) + mchild * xcrossv(:) - colliders%L_spin(:, j) = colliders%L_spin(:, j) + mchild * pl%Ip(3, idx_child) * pl%radius(idx_child)**2 * pl%rot(:, idx_child) + colliders%L_spin(:, j) = colliders%L_spin(:, j) + mchild * pl%Ip(3, idx_child) & + * pl%radius(idx_child)**2 & + * pl%rot(:, idx_child) colliders%Ip(:, j) = colliders%Ip(:, j) + mchild * pl%Ip(:, idx_child) end if @@ -752,7 +756,9 @@ subroutine symba_collision_mergeaddsub(system, param, colliders, frag, status) plnew%status(1:nfrag) = NEW_PARTICLE do i = 1, nfrag write(newname, FRAGFMT) frag%id(i) - call plnew%info(i)%set_value(origin_type="Disruption", origin_time=param%t, name=newname, origin_xh=plnew%xh(:,i), origin_vh=plnew%vh(:,i), collision_id=param%maxid_collision) + call plnew%info(i)%set_value(origin_type="Disruption", origin_time=param%t, name=newname, & + origin_xh=plnew%xh(:,i), & + origin_vh=plnew%vh(:,i), collision_id=param%maxid_collision) end do do i = 1, ncolliders if (colliders%idx(i) == ibiggest) then @@ -760,13 +766,16 @@ subroutine symba_collision_mergeaddsub(system, param, colliders, frag, status) else iother = ibiggest end if - call pl%info(colliders%idx(i))%set_value(status="Disruption", discard_time=param%t, discard_xh=pl%xh(:,i), discard_vh=pl%vh(:,i), discard_body_id=iother) + call pl%info(colliders%idx(i))%set_value(status="Disruption", discard_time=param%t, & + discard_xh=pl%xh(:,i), discard_vh=pl%vh(:,i), discard_body_id=iother) end do case(SUPERCATASTROPHIC) plnew%status(1:nfrag) = NEW_PARTICLE do i = 1, nfrag write(newname, FRAGFMT) frag%id(i) - call plnew%info(i)%set_value(origin_type="Supercatastrophic", origin_time=param%t, name=newname, origin_xh=plnew%xh(:,i), origin_vh=plnew%vh(:,i), collision_id=param%maxid_collision) + call plnew%info(i)%set_value(origin_type="Supercatastrophic", origin_time=param%t, name=newname, & + origin_xh=plnew%xh(:,i), origin_vh=plnew%vh(:,i), & + collision_id=param%maxid_collision) end do do i = 1, ncolliders if (colliders%idx(i) == ibiggest) then @@ -774,19 +783,25 @@ subroutine symba_collision_mergeaddsub(system, param, colliders, frag, status) else iother = ibiggest end if - call pl%info(colliders%idx(i))%set_value(status="Supercatastrophic", discard_time=param%t, discard_xh=pl%xh(:,i), discard_vh=pl%vh(:,i), discard_body_id=iother) + call pl%info(colliders%idx(i))%set_value(status="Supercatastrophic", discard_time=param%t, & + discard_xh=pl%xh(:,i), discard_vh=pl%vh(:,i), & + discard_body_id=iother) end do case(HIT_AND_RUN_DISRUPT) call plnew%info(1)%copy(pl%info(ibiggest)) plnew%status(1) = OLD_PARTICLE do i = 2, nfrag write(newname, FRAGFMT) frag%id(i) - call plnew%info(i)%set_value(origin_type="Hit and run fragment", origin_time=param%t, name=newname, origin_xh=plnew%xh(:,i), origin_vh=plnew%vh(:,i), collision_id=param%maxid_collision) + call plnew%info(i)%set_value(origin_type="Hit and run fragment", origin_time=param%t, name=newname, & + origin_xh=plnew%xh(:,i), origin_vh=plnew%vh(:,i), & + collision_id=param%maxid_collision) end do do i = 1, ncolliders if (colliders%idx(i) == ibiggest) cycle iother = ibiggest - call pl%info(colliders%idx(i))%set_value(status="Hit and run fragmention", discard_time=param%t, discard_xh=pl%xh(:,i), discard_vh=pl%vh(:,i), discard_body_id=iother) + call pl%info(colliders%idx(i))%set_value(status="Hit and run fragmention", discard_time=param%t, & + discard_xh=pl%xh(:,i), discard_vh=pl%vh(:,i), & + discard_body_id=iother) end do case(MERGED) call plnew%info(1)%copy(pl%info(ibiggest)) @@ -795,7 +810,8 @@ subroutine symba_collision_mergeaddsub(system, param, colliders, frag, status) if (colliders%idx(i) == ibiggest) cycle iother = ibiggest - call pl%info(colliders%idx(i))%set_value(status="MERGED", discard_time=param%t, discard_xh=pl%xh(:,i), discard_vh=pl%vh(:,i), discard_body_id=iother) + call pl%info(colliders%idx(i))%set_value(status="MERGED", discard_time=param%t, discard_xh=pl%xh(:,i), & + discard_vh=pl%vh(:,i), discard_body_id=iother) end do end select @@ -995,9 +1011,12 @@ module subroutine symba_collision_resolve_plplenc(self, system, param, t, dt, ir do write(timestr,*) t call io_log_one_message(FRAGGLE_LOG_OUT, "") - call io_log_one_message(FRAGGLE_LOG_OUT, "***********************************************************************************************************************") - call io_log_one_message(FRAGGLE_LOG_OUT, "Collision between massive bodies detected at time t = " // trim(adjustl(timestr))) - call io_log_one_message(FRAGGLE_LOG_OUT, "***********************************************************************************************************************") + call io_log_one_message(FRAGGLE_LOG_OUT, "***********************************************************" // & + "***********************************************************") + call io_log_one_message(FRAGGLE_LOG_OUT, "Collision between massive bodies detected at time t = " // & + trim(adjustl(timestr))) + call io_log_one_message(FRAGGLE_LOG_OUT, "***********************************************************" // & + "***********************************************************") allocate(tmp_param, source=param) tmp_param%t = t if (param%lfragmentation) then @@ -1007,7 +1026,7 @@ module subroutine symba_collision_resolve_plplenc(self, system, param, t, dt, ir end if ! Destroy the collision list now that the collisions are resolved - call plplcollision_list%setup(0) + call plplcollision_list%setup(0_I8B) if ((system%pl_adds%nbody == 0) .and. (system%pl_discards%nbody == 0)) exit diff --git a/src/symba/symba_discard.f90 b/src/symba/symba_discard.f90 index 0ae89587f..385f6f25b 100644 --- a/src/symba/symba_discard.f90 +++ b/src/symba/symba_discard.f90 @@ -36,26 +36,34 @@ subroutine symba_discard_cb_pl(pl, system, param) pl%status(i) = DISCARDED_RMAX write(idstr, *) pl%id(i) write(timestr, *) param%t - write(message, *) trim(adjustl(pl%info(i)%name)) // " (" // trim(adjustl(idstr)) // ")" // " too far from the central body at t = " // trim(adjustl(timestr)) + write(message, *) trim(adjustl(pl%info(i)%name)) // " (" // trim(adjustl(idstr)) // ")" // & + " too far from the central body at t = " // trim(adjustl(timestr)) call io_log_one_message(FRAGGLE_LOG_OUT, "") - call io_log_one_message(FRAGGLE_LOG_OUT, "***********************************************************************************************************************") + call io_log_one_message(FRAGGLE_LOG_OUT, "***********************************************************" // & + "***********************************************************") call io_log_one_message(FRAGGLE_LOG_OUT, message) - call io_log_one_message(FRAGGLE_LOG_OUT, "***********************************************************************************************************************") + call io_log_one_message(FRAGGLE_LOG_OUT, "***********************************************************" // & + "***********************************************************") call io_log_one_message(FRAGGLE_LOG_OUT, "") - call pl%info(i)%set_value(status="DISCARDED_RMAX", discard_time=param%t, discard_xh=pl%xh(:,i), discard_vh=pl%vh(:,i)) + call pl%info(i)%set_value(status="DISCARDED_RMAX", discard_time=param%t, discard_xh=pl%xh(:,i), & + discard_vh=pl%vh(:,i)) 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(message, *) trim(adjustl(pl%info(i)%name)) // " (" // trim(adjustl(idstr)) // ")" // " too close to the central body at t = " // trim(adjustl(timestr)) + write(message, *) trim(adjustl(pl%info(i)%name)) // " (" // trim(adjustl(idstr)) // ")" // & + " too close to the central body at t = " // trim(adjustl(timestr)) call io_log_one_message(FRAGGLE_LOG_OUT, "") - call io_log_one_message(FRAGGLE_LOG_OUT, "***********************************************************************************************************************") + call io_log_one_message(FRAGGLE_LOG_OUT, "************************************************************" // & + "************************************************************") call io_log_one_message(FRAGGLE_LOG_OUT, message) - call io_log_one_message(FRAGGLE_LOG_OUT, "***********************************************************************************************************************") + call io_log_one_message(FRAGGLE_LOG_OUT, "************************************************************" // & + "************************************************************") call io_log_one_message(FRAGGLE_LOG_OUT, "") - call pl%info(i)%set_value(status="DISCARDED_RMIN", discard_time=param%t, discard_xh=pl%xh(:,i), discard_vh=pl%vh(:,i), discard_body_id=cb%id) + call pl%info(i)%set_value(status="DISCARDED_RMIN", discard_time=param%t, discard_xh=pl%xh(:,i), & + discard_vh=pl%vh(:,i), discard_body_id=cb%id) 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)) @@ -66,13 +74,17 @@ subroutine symba_discard_cb_pl(pl, system, param) pl%status(i) = DISCARDED_RMAXU write(idstr, *) pl%id(i) write(timestr, *) param%t - write(message, *) trim(adjustl(pl%info(i)%name)) // " (" // trim(adjustl(idstr)) // ")" // " is unbound and too far from barycenter at t = " // trim(adjustl(timestr)) + write(message, *) trim(adjustl(pl%info(i)%name)) // " (" // trim(adjustl(idstr)) // ")" // & + " is unbound and too far from barycenter at t = " // trim(adjustl(timestr)) call io_log_one_message(FRAGGLE_LOG_OUT, "") - call io_log_one_message(FRAGGLE_LOG_OUT, "***********************************************************************************************************************") + call io_log_one_message(FRAGGLE_LOG_OUT, "************************************************************" // & + "************************************************************") call io_log_one_message(FRAGGLE_LOG_OUT, message) - call io_log_one_message(FRAGGLE_LOG_OUT, "***********************************************************************************************************************") + call io_log_one_message(FRAGGLE_LOG_OUT, "************************************************************" // & + "************************************************************") call io_log_one_message(FRAGGLE_LOG_OUT, "") - call pl%info(i)%set_value(status="DISCARDED_RMAXU", discard_time=param%t, discard_xh=pl%xh(:,i), discard_vh=pl%vh(:,i)) + call pl%info(i)%set_value(status="DISCARDED_RMAXU", discard_time=param%t, discard_xh=pl%xh(:,i), & + discard_vh=pl%vh(:,i)) end if end if end if @@ -138,7 +150,8 @@ subroutine symba_discard_conserve_mtm(pl, system, param, ipl, lescape_body) end do Ltot(:) = Ltot(:) - cb%mass * (cb%xb(:) .cross. cb%vb(:)) system%Lescape(:) = system%Lescape(:) + Ltot(:) - if (param%lrotation) system%Lescape(:) = system%Lescape + pl%mass(ipl) * pl%radius(ipl)**2 * pl%Ip(3, ipl) * pl%rot(:, ipl) + if (param%lrotation) system%Lescape(:) = system%Lescape + pl%mass(ipl) * pl%radius(ipl)**2 & + * pl%Ip(3, ipl) * pl%rot(:, ipl) else xcom(:) = (pl%mass(ipl) * pl%xb(:, ipl) + cb%mass * cb%xb(:)) / (cb%mass + pl%mass(ipl)) @@ -305,8 +318,10 @@ subroutine symba_discard_peri_pl(pl, system, param) pl%status(i) = DISCARDED_PERI 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)) - call pl%info(i)%set_value(status="DISCARDED_PERI", discard_time=param%t, discard_xh=pl%xh(:,i), discard_vh=pl%vh(:,i), discard_body_id=system%cb%id) + write(*, *) trim(adjustl(pl%info(i)%name)) // " (" // trim(adjustl(idstr)) // & + ") perihelion distance too small at t = " // trim(adjustl(timestr)) + call pl%info(i)%set_value(status="DISCARDED_PERI", discard_time=param%t, & + discard_xh=pl%xh(:,i), discard_vh=pl%vh(:,i), discard_body_id=system%cb%id) end if end if end if diff --git a/src/symba/symba_encounter_check.f90 b/src/symba/symba_encounter_check.f90 index f75fc49f2..7281916b2 100644 --- a/src/symba/symba_encounter_check.f90 +++ b/src/symba/symba_encounter_check.f90 @@ -17,8 +17,8 @@ module function symba_encounter_check_pl(self, param, system, dt, irec) result(l ! Result logical :: lany_encounter !! Returns true if there is at least one close encounter ! Internals - integer(I8B) :: k, nplplm, kenc - integer(I4B) :: i, j, nenc, npl, nplm, nplt + integer(I8B) :: k, nplplm, kenc, nenc + integer(I4B) :: i, j, npl, nplm, nplt logical, dimension(:), allocatable :: lencounter, loc_lvdotr, lvdotr integer(I4B), dimension(:), allocatable :: index1, index2 integer(I4B), dimension(:,:), allocatable :: k_plpl_enc @@ -38,12 +38,13 @@ module function symba_encounter_check_pl(self, param, system, dt, irec) result(l call pl%set_renc(irec) if (nplt == 0) then - call encounter_check_all_plpl(param, npl, pl%xh, pl%vh, pl%renc, dt, lvdotr, index1, index2, nenc) + call encounter_check_all_plpl(param, npl, pl%xh, pl%vh, pl%renc, dt, nenc, index1, index2, lvdotr) else - call encounter_check_all_plplm(param, nplm, nplt, pl%xh(:,1:nplm), pl%vh(:,1:nplm), pl%xh(:,nplm+1:npl), pl%vh(:,nplm+1:npl), pl%renc(1:nplm), pl%renc(nplm+1:npl), dt, lvdotr, index1, index2, nenc) + call encounter_check_all_plplm(param, nplm, nplt, pl%xh(:,1:nplm), pl%vh(:,1:nplm), pl%xh(:,nplm+1:npl), & + pl%vh(:,nplm+1:npl), pl%renc(1:nplm), pl%renc(nplm+1:npl), dt, nenc, index1, index2, lvdotr) end if - lany_encounter = nenc > 0 + lany_encounter = nenc > 0_I8B if (lany_encounter) then call plplenc_list%resize(nenc) call move_alloc(lvdotr, plplenc_list%lvdotr) @@ -52,7 +53,7 @@ module function symba_encounter_check_pl(self, param, system, dt, irec) result(l end if if (lany_encounter) then - do k = 1, nenc + do k = 1_I8B, nenc i = plplenc_list%index1(k) j = plplenc_list%index2(k) plplenc_list%id1(k) = pl%id(i) @@ -146,7 +147,8 @@ module function symba_encounter_check(self, param, system, dt, irec) result(lany j = self%index2(k) xr(:) = tp%xh(:,j) - pl%xh(:,i) vr(:) = tp%vb(:,j) - pl%vb(:,i) - call encounter_check_one(xr(1), xr(2), xr(3), vr(1), vr(2), vr(3), pl%renc(i), dt, lencounter(lidx), self%lvdotr(k)) + call encounter_check_one(xr(1), xr(2), xr(3), vr(1), vr(2), vr(3), pl%renc(i), dt, & + lencounter(lidx), self%lvdotr(k)) if (lencounter(lidx)) then rlim2 = (pl%radius(i))**2 rji2 = dot_product(xr(:), xr(:))! Check to see if these are physically overlapping bodies first, which we should ignore @@ -197,7 +199,8 @@ module function symba_encounter_check_tp(self, param, system, dt, irec) result(l logical :: lany_encounter !! Returns true if there is at least one close encounter ! Internals real(DP) :: r2crit, vdotr, r2, v2, tmin, r2min, term2 - integer(I4B) :: i, j, k,nenc, plind, tpind + integer(I4B) :: i, j, plind, tpind + integer(I8B) :: k, nenc real(DP), dimension(NDIM) :: xr, vr real(DP) :: rshell_irec logical, dimension(:), allocatable :: lvdotr @@ -208,7 +211,7 @@ module function symba_encounter_check_tp(self, param, system, dt, irec) result(l associate(tp => self, ntp => self%nbody, pl => system%pl, npl => system%pl%nbody) call pl%set_renc(irec) - call encounter_check_all_pltp(param, npl, ntp, pl%xh, pl%vh, tp%xh, tp%vh, pl%renc, dt, lvdotr, index1, index2, nenc) + call encounter_check_all_pltp(param, npl, ntp, pl%xh, pl%vh, tp%xh, tp%vh, pl%renc, dt, nenc, index1, index2, lvdotr) lany_encounter = nenc > 0 if (lany_encounter) then @@ -224,7 +227,7 @@ module function symba_encounter_check_tp(self, param, system, dt, irec) result(l select type(pl) class is (symba_pl) pl%lencounter(1:npl) = .false. - do k = 1, nenc + do k = 1_I8B, nenc plind = pltpenc_list%index1(k) tpind = pltpenc_list%index2(k) pl%lencounter(plind) = .true. diff --git a/src/symba/symba_io.f90 b/src/symba/symba_io.f90 index 3e217bc6a..19e8635db 100644 --- a/src/symba/symba_io.f90 +++ b/src/symba/symba_io.f90 @@ -74,7 +74,8 @@ module subroutine symba_io_param_reader(self, unit, iotype, v_list, iostat, ioms param_value = io_get_token(line, ifirst, ilast, iostat) read(param_value, *) param%seed(i) end do - param%seed(nseeds_from_file+1:nseeds) = [(param%seed(1) - param%seed(nseeds_from_file) + i, i=nseeds_from_file+1, nseeds)] + param%seed(nseeds_from_file+1:nseeds) = [(param%seed(1) - param%seed(nseeds_from_file) + i, & + i=nseeds_from_file+1, nseeds)] end if seed_set = .true. end select @@ -192,9 +193,9 @@ module subroutine symba_io_write_discard(self, param) end if select case(out_stat) case('APPEND') - open(unit = LUN, file = param%discard_out, status = 'OLD', position = 'APPEND', form = 'FORMATTED', err = 667, iomsg = errmsg) + open(unit=LUN, file=param%discard_out, status='OLD', position='APPEND', form='FORMATTED', err=667, iomsg=errmsg) case('NEW', 'REPLACE', 'UNKNOWN') - open(unit = LUN, file = param%discard_out, status = param%out_stat, form = 'FORMATTED', err = 667, iomsg = errmsg) + open(unit=LUN, file=param%discard_out, status=param%out_stat, form='FORMATTED', err=667, iomsg=errmsg) case default write(*,*) 'Invalid status code for OUT_STAT: ',trim(adjustl(param%out_stat)) call util_exit(FAILURE) @@ -205,7 +206,7 @@ module subroutine symba_io_write_discard(self, param) call pl_adds%pv2v(param) end if - write(LUN, HDRFMT, err = 667, iomsg = errmsg) param%t, pl_discards%nbody, param%lbig_discard + write(LUN, HDRFMT, err=667, iomsg=errmsg) param%t, pl_discards%nbody, param%lbig_discard iadd = 1 isub = 1 do while (iadd <= pl_adds%nbody) @@ -213,9 +214,9 @@ module subroutine symba_io_write_discard(self, param) nsub = pl_discards%ncomp(isub) do j = 1, nadd if (iadd <= pl_adds%nbody) then - write(LUN, NAMEFMT, err = 667, iomsg = errmsg) ADD, pl_adds%id(iadd), pl_adds%status(iadd) - write(LUN, VECFMT, err = 667, iomsg = errmsg) pl_adds%xh(1, iadd), pl_adds%xh(2, iadd), pl_adds%xh(3, iadd) - write(LUN, VECFMT, err = 667, iomsg = errmsg) pl_adds%vh(1, iadd), pl_adds%vh(2, iadd), pl_adds%vh(3, iadd) + write(LUN, NAMEFMT, err=667, iomsg=errmsg) ADD, pl_adds%id(iadd), pl_adds%status(iadd) + write(LUN, VECFMT, err=667, iomsg=errmsg) pl_adds%xh(1, iadd), pl_adds%xh(2, iadd), pl_adds%xh(3, iadd) + write(LUN, VECFMT, err=667, iomsg=errmsg) pl_adds%vh(1, iadd), pl_adds%vh(2, iadd), pl_adds%vh(3, iadd) else exit end if @@ -223,9 +224,9 @@ module subroutine symba_io_write_discard(self, param) end do do j = 1, nsub if (isub <= pl_discards%nbody) then - write(LUN, NAMEFMT, err = 667, iomsg = errmsg) SUB, pl_discards%id(isub), pl_discards%status(isub) - write(LUN, VECFMT, err = 667, iomsg = errmsg) pl_discards%xh(1, isub), pl_discards%xh(2, isub), pl_discards%xh(3, isub) - write(LUN, VECFMT, err = 667, iomsg = errmsg) pl_discards%vh(1, isub), pl_discards%vh(2, isub), pl_discards%vh(3, isub) + write(LUN,NAMEFMT,err=667,iomsg=errmsg) SUB, pl_discards%id(isub), pl_discards%status(isub) + write(LUN,VECFMT,err=667,iomsg=errmsg) pl_discards%xh(1,isub), pl_discards%xh(2,isub), pl_discards%xh(3,isub) + write(LUN,VECFMT,err=667,iomsg=errmsg) pl_discards%vh(1,isub), pl_discards%vh(2,isub), pl_discards%vh(3,isub) else exit end if diff --git a/src/symba/symba_setup.f90 b/src/symba/symba_setup.f90 index bb300f63d..cda5214d0 100644 --- a/src/symba/symba_setup.f90 +++ b/src/symba/symba_setup.f90 @@ -18,9 +18,9 @@ module subroutine symba_setup_initialize_system(self, param) ! Call parent method associate(system => self) call helio_setup_initialize_system(system, param) - call system%pltpenc_list%setup(0) - call system%plplenc_list%setup(0) - call system%plplcollision_list%setup(0) + call system%pltpenc_list%setup(0_I8B) + call system%plplenc_list%setup(0_I8B) + call system%plplcollision_list%setup(0_I8B) end associate return @@ -109,14 +109,14 @@ module subroutine symba_setup_encounter_list(self, n) implicit none ! Arguments class(symba_encounter), intent(inout) :: self !! SyMBA pl-tp encounter structure - integer(I4B), intent(in) :: n !! Number of encounters to allocate space for + integer(I8B), intent(in) :: n !! Number of encounters to allocate space for call encounter_setup_list(self, n) - if (n < 0) return + if (n < 0_I8B) return call self%dealloc() - if (n ==0) return + if (n ==0_I8B) return allocate(self%level(n)) diff --git a/src/symba/symba_step.f90 b/src/symba/symba_step.f90 index 0d077d22d..dda813b82 100644 --- a/src/symba/symba_step.f90 +++ b/src/symba/symba_step.f90 @@ -116,7 +116,8 @@ module subroutine symba_step_set_recur_levels_system(self, ireci) ! Internals integer(I4B) :: k, irecp - associate(system => self, plplenc_list => self%plplenc_list, pltpenc_list => self%pltpenc_list, npl => self%pl%nbody, ntp => self%tp%nbody) + associate(system => self, plplenc_list => self%plplenc_list, pltpenc_list => self%pltpenc_list, & + npl => self%pl%nbody, ntp => self%tp%nbody) select type(pl => self%pl) class is (symba_pl) select type(tp => self%tp) @@ -125,8 +126,16 @@ module subroutine symba_step_set_recur_levels_system(self, ireci) if (npl >0) where(pl%levelg(1:npl) == irecp) pl%levelg(1:npl) = ireci if (ntp > 0) where(tp%levelg(1:ntp) == irecp) tp%levelg(1:ntp) = ireci - if (plplenc_list%nenc > 0) where(plplenc_list%level(1:plplenc_list%nenc) == irecp) plplenc_list%level(1:plplenc_list%nenc) = ireci - if (pltpenc_list%nenc > 0) where(pltpenc_list%level(1:pltpenc_list%nenc) == irecp) pltpenc_list%level(1:pltpenc_list%nenc) = ireci + if (plplenc_list%nenc > 0) then + where(plplenc_list%level(1:plplenc_list%nenc) == irecp) + plplenc_list%level(1:plplenc_list%nenc) = ireci + endwhere + end if + if (pltpenc_list%nenc > 0) then + where(pltpenc_list%level(1:pltpenc_list%nenc) == irecp) + pltpenc_list%level(1:pltpenc_list%nenc) = ireci + endwhere + end if system%irec = ireci @@ -179,7 +188,8 @@ module recursive subroutine symba_step_recur_system(self, param, t, ireci) nloops = NTENC end if do j = 1, nloops - lencounter = plplenc_list%encounter_check(param, system, dtl, irecp) .or. pltpenc_list%encounter_check(param, system, dtl, irecp) + lencounter = plplenc_list%encounter_check(param, system, dtl, irecp) & + .or. pltpenc_list%encounter_check(param, system, dtl, irecp) call plplenc_list%kick(system, dth, irecp, 1) call pltpenc_list%kick(system, dth, irecp, 1) @@ -242,7 +252,8 @@ module subroutine symba_step_reset_system(self, param) class(symba_nbody_system), intent(inout) :: self !! SyMBA nbody system object class(symba_parameters), intent(in) :: param !! Current run configuration parameters with SyMBA additions ! Internals - integer(I4B) :: i, nenc_old + integer(I4B) :: i + integer(I8B) :: nenc_old associate(system => self) select type(pl => system%pl) diff --git a/src/symba/symba_util.f90 b/src/symba/symba_util.f90 index 5a985faea..0cbdc9a4a 100644 --- a/src/symba/symba_util.f90 +++ b/src/symba/symba_util.f90 @@ -38,17 +38,19 @@ module subroutine symba_util_append_encounter_list(self, source, lsource_mask) !! This method will automatically resize the destination body if it is too small implicit none ! Arguments - class(symba_encounter), intent(inout) :: self !! SyMBA encounter list object - class(encounter_list), intent(in) :: source !! Source object to append - logical, dimension(:), intent(in) :: lsource_mask !! Logical mask indicating which elements to append to + class(symba_encounter), intent(inout) :: self !! SyMBA encounter list object + class(encounter_list), intent(in) :: source !! Source object to append + logical, dimension(:), intent(in) :: lsource_mask !! Logical mask indicating which elements to append to + ! Internals + integer(I4B) :: nold, nsrc - associate(nold => self%nenc, nsrc => source%nenc) - select type(source) - class is (symba_encounter) - call util_append(self%level, source%level, nold, nsrc, lsource_mask) - end select - call encounter_util_append_list(self, source, lsource_mask) - end associate + nold = self%nenc + nsrc = source%nenc + select type(source) + class is (symba_encounter) + call util_append(self%level, source%level, nold, nsrc, lsource_mask) + end select + call encounter_util_append_list(self, source, lsource_mask) return end subroutine symba_util_append_encounter_list @@ -403,7 +405,7 @@ module subroutine symba_util_flatten_eucl_plpl(self, param) integer(I8B) :: k, npl, nplm integer(I4B) :: i, j, err - associate(pl => self, nplpl => self%nplpl, nplplm => self%nplplm) + associate(pl => self, nplplm => self%nplplm) npl = int(self%nbody, kind=I8B) select type(param) class is (symba_parameters) @@ -411,21 +413,9 @@ module subroutine symba_util_flatten_eucl_plpl(self, param) end select nplm = count(.not. pl%lmtiny(1:npl)) pl%nplm = int(nplm, kind=I4B) - nplpl = (npl * (npl - 1_I8B)) / 2_I8B ! number of entries in a strict lower triangle, npl x npl, minus first column nplplm = nplm * npl - nplm * (nplm + 1_I8B) / 2_I8B ! number of entries in a strict lower triangle, npl x npl, minus first column including only mutually interacting bodies - if (param%lflatten_interactions) then - 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), stat=err) - if (err /= 0) then ! An error occurred trying to allocate this big array. This probably means it's too big to fit in memory, and so we will force the run back into triangular mode - param%lflatten_interactions = .false. - else - do concurrent (i=1:npl, j=1:npl, j>i) - call util_flatten_eucl_ij_to_k(self%nbody, i, j, k) - self%k_plpl(1, k) = i - self%k_plpl(2, k) = j - end do - end if - end if + + call util_flatten_eucl_plpl(pl, param) end associate return @@ -579,7 +569,8 @@ module subroutine symba_util_peri_pl(self, system, param) if (pl%isperi(i) == -1) then if (vdotr >= 0.0_DP) then pl%isperi(i) = 0 - CALL orbel_xv2aeq(system%Gmtot, pl%xb(1,i), pl%xb(2,i), pl%xb(3,i), pl%vb(1,i), pl%vb(2,i), pl%vb(3,i), pl%atp(i), e, pl%peri(i)) + CALL orbel_xv2aeq(system%Gmtot, pl%xb(1,i), pl%xb(2,i), pl%xb(3,i), pl%vb(1,i), pl%vb(2,i), pl%vb(3,i),& + pl%atp(i), e, pl%peri(i)) end if else if (vdotr > 0.0_DP) then @@ -613,7 +604,8 @@ module subroutine symba_util_rearray_pl(self, system, param) logical, dimension(:), allocatable :: lmask, ldump_mask class(symba_plplenc), allocatable :: plplenc_old logical :: lencounter - integer(I4B), dimension(:), allocatable :: levelg_orig_pl, levelm_orig_pl, levelg_orig_tp, levelm_orig_tp, nplenc_orig_pl, nplenc_orig_tp, ntpenc_orig_pl + integer(I4B), dimension(:), allocatable :: levelg_orig_pl, levelm_orig_pl, levelg_orig_tp, levelm_orig_tp + integer(I4B), dimension(:), allocatable :: nplenc_orig_pl, nplenc_orig_tp, ntpenc_orig_pl associate(pl => self, pl_adds => system%pl_adds) diff --git a/src/util/util_get_energy_momentum.f90 b/src/util/util_get_energy_momentum.f90 index a86f5efd4..b935a680c 100644 --- a/src/util/util_get_energy_momentum.f90 +++ b/src/util/util_get_energy_momentum.f90 @@ -144,7 +144,8 @@ subroutine util_get_energy_potential_flat(npl, nplpl, k_plpl, lmask, GMcb, Gmass end do !$omp parallel do default(private) schedule(static)& - !$omp shared(nplpl, k_plpl, xb, mass, Gmass, pepl, lstatpl, lmask) + !$omp shared(k_plpl, xb, mass, Gmass, pepl, lstatpl, lmask) & + !$omp firstprivate(nplpl) do k = 1, nplpl i = k_plpl(1,k) j = k_plpl(2,k) @@ -178,7 +179,7 @@ subroutine util_get_energy_potential_triangular(npl, lmask, GMcb, Gmass, mass, x real(DP), intent(out) :: pe ! Internals integer(I4B) :: i, j - real(DP), dimension(npl) :: pecb + real(DP), dimension(npl) :: pecb, pepl ! Do the central body potential energy component first where(.not. lmask(1:npl)) @@ -191,12 +192,14 @@ subroutine util_get_energy_potential_triangular(npl, lmask, GMcb, Gmass, mass, x pe = 0.0_DP !$omp parallel do default(private) schedule(static)& - !$omp shared(npl, lmask, Gmass, mass, xb) & - !$omp reduction(-:pe) + !$omp shared(lmask, Gmass, mass, xb) & + !$omp firstprivate(npl) & + !$omp reduction(+:pe) do i = 1, npl do concurrent(j = i+1:npl, lmask(i) .and. lmask(j)) - pe = pe - (Gmass(i) * mass(j)) / norm2(xb(:, i) - xb(:, j)) + pepl(j) = - (Gmass(i) * mass(j)) / norm2(xb(:, i) - xb(:, j)) end do + pe = pe + sum(pepl(i+1:npl), lmask(i) .and. lmask(j)) end do !$omp end parallel do pe = pe + sum(pecb(1:npl), lmask(1:npl)) diff --git a/src/util/util_peri.f90 b/src/util/util_peri.f90 index a81748d74..d451e3dbd 100644 --- a/src/util/util_peri.f90 +++ b/src/util/util_peri.f90 @@ -29,7 +29,8 @@ module subroutine util_peri_tp(self, system, param) if (tp%isperi(i) == -1) then if (vdotr(i) >= 0.0_DP) then tp%isperi(i) = 0 - call orbel_xv2aeq(tp%mu(i), tp%xh(1,i), tp%xh(2,i), tp%xh(3,i), tp%vh(1,i), tp%vh(2,i), tp%vh(3,i), tp%atp(i), e, tp%peri(i)) + call orbel_xv2aeq(tp%mu(i), tp%xh(1,i), tp%xh(2,i), tp%xh(3,i), tp%vh(1,i), tp%vh(2,i), tp%vh(3,i), & + tp%atp(i), e, tp%peri(i)) end if else if (vdotr(i) > 0.0_DP) then @@ -45,7 +46,8 @@ module subroutine util_peri_tp(self, system, param) if (tp%isperi(i) == -1) then if (vdotr(i) >= 0.0_DP) then tp%isperi(i) = 0 - call orbel_xv2aeq(system%Gmtot, tp%xb(1,i), tp%xb(2,i), tp%xb(3,i), tp%vb(1,i), tp%vb(2,i), tp%vb(3,i), tp%atp(i), e, tp%peri(i)) + call orbel_xv2aeq(system%Gmtot, tp%xb(1,i), tp%xb(2,i), tp%xb(3,i), tp%vb(1,i), tp%vb(2,i), tp%vb(3,i), & + tp%atp(i), e, tp%peri(i)) end if else if (vdotr(i) > 0.0_DP) then diff --git a/src/util/util_set.f90 b/src/util/util_set.f90 index 98ee4f497..99dff1e0d 100644 --- a/src/util/util_set.f90 +++ b/src/util/util_set.f90 @@ -98,7 +98,8 @@ module subroutine util_set_mu_tp(self, cb) return end subroutine util_set_mu_tp - module subroutine util_set_particle_info(self, name, particle_type, status, origin_type, origin_time, collision_id, origin_xh, origin_vh, discard_time, discard_xh, discard_vh, discard_body_id) + module subroutine util_set_particle_info(self, name, particle_type, status, origin_type, origin_time, collision_id, origin_xh,& + origin_vh, discard_time, discard_xh, discard_vh, discard_body_id) !! author: David A. Minton !! !! Sets one or more values of the particle information metadata object diff --git a/src/util/util_sort.f90 b/src/util/util_sort.f90 index 84371b773..453c3a2d3 100644 --- a/src/util/util_sort.f90 +++ b/src/util/util_sort.f90 @@ -202,7 +202,7 @@ module pure subroutine util_sort_i4b(arr) end subroutine util_sort_i4b - module pure subroutine util_sort_index_i4b(arr, ind) + module pure subroutine util_sort_index_I4B(arr, ind) !! author: David A. Minton !! !! Sort input integer array by index in ascending numerical order using quicksort. @@ -227,20 +227,48 @@ module pure subroutine util_sort_index_i4b(arr, ind) call qsort_I4B(tmparr, ind) return - end subroutine util_sort_index_i4b + end subroutine util_sort_index_I4B + + + module pure subroutine util_sort_index_I4B_I8Bind(arr, ind) + !! author: David A. Minton + !! + !! Sort input integer array by index in ascending numerical order using quicksort. + !! If ind is supplied already allocated, we assume it is an existing index array (e.g. a previously + !! sorted array). If it is not allocated, this subroutine allocates it. + !! + implicit none + ! Arguments + integer(I4B), dimension(:), intent(in) :: arr + integer(I8B), dimension(:), allocatable, intent(inout) :: ind + ! Internals + integer(I8B) :: n, i + integer(I4B), dimension(:), allocatable :: tmparr + + n = size(arr) + if (.not.allocated(ind)) then + allocate(ind(n)) + ind = [(i, i=1_I8B, n)] + end if + allocate(tmparr, mold=arr) + tmparr(:) = arr(ind(:)) + call qsort_I4B_I8Bind(tmparr, ind) + + return + end subroutine util_sort_index_I4B_I8Bind recursive pure subroutine qsort_I4B(arr, ind) !! author: David A. Minton !! - !! Sort input DP precision array by index in ascending numerical order using quicksort. + !! Sort input I4B array by index in ascending numerical order using quicksort. !! implicit none ! Arguments integer(I4B), dimension(:), intent(inout) :: arr integer(I4B), dimension(:), intent(out), optional :: ind - !! Internals - integer :: iq + ! Internals + integer(I4B) :: iq if (size(arr) > 1) then if (present(ind)) then @@ -257,6 +285,61 @@ recursive pure subroutine qsort_I4B(arr, ind) return end subroutine qsort_I4B + recursive pure subroutine qsort_I4B_I8Bind(arr, ind) + !! author: David A. Minton + !! + !! Sort input I4B array by index in ascending numerical order using quicksort. + !! + implicit none + ! Arguments + integer(I4B), dimension(:), intent(inout) :: arr + integer(I8B), dimension(:), intent(out), optional :: ind + ! Internals + integer(I8B) :: iq + + if (size(arr) > 1_I8B) then + if (present(ind)) then + call partition_I4B_I8Bind(arr, iq, ind) + call qsort_I4B_I8Bind(arr(:iq-1_I8B),ind(:iq-1_I8B)) + call qsort_I4B_I8Bind(arr(iq:), ind(iq:)) + else + call partition_I4B_I8Bind(arr, iq) + call qsort_I4B_I8Bind(arr(:iq-1_I8B)) + call qsort_I4B_I8Bind(arr(iq:)) + end if + end if + + return + end subroutine qsort_I4B_I8Bind + + + recursive pure subroutine qsort_I8B_I8Bind(arr, ind) + !! author: David A. Minton + !! + !! Sort input I8B array by index in ascending numerical order using quicksort. + !! + implicit none + ! Arguments + integer(I8B), dimension(:), intent(inout) :: arr + integer(I8B), dimension(:), intent(out), optional :: ind + ! Internals + integer(I8B) :: iq + + if (size(arr) > 1_I8B) then + if (present(ind)) then + call partition_I8B_I8Bind(arr, iq, ind) + call qsort_I8B_I8Bind(arr(:iq-1_I8B),ind(:iq-1_I8B)) + call qsort_I8B_I8Bind(arr(iq:), ind(iq:)) + else + call partition_I8B_I8Bind(arr, iq) + call qsort_I8B_I8Bind(arr(:iq-1_I8B)) + call qsort_I8B_I8Bind(arr(iq:)) + end if + end if + + return + end subroutine qsort_I8B_I8Bind + pure subroutine partition_I4B(arr, marker, ind) !! author: David A. Minton @@ -314,6 +397,118 @@ pure subroutine partition_I4B(arr, marker, ind) return end subroutine partition_I4B + pure subroutine partition_I4B_I8Bind(arr, marker, ind) + !! author: David A. Minton + !! + !! Partition function for quicksort on I4B type + !! + implicit none + ! Arguments + integer(I4B), intent(inout), dimension(:) :: arr + integer(I8B), intent(inout), dimension(:), optional :: ind + integer(I8B), intent(out) :: marker + ! Internals + integer(I8B) :: i, j, itmp, narr, ipiv + integer(I4B) :: temp + integer(I8B) :: x ! pivot point + + narr = size(arr) + + ! Get center as pivot, as this is likely partially sorted + ipiv = narr / 2_I8B + x = arr(ipiv) + i = 0_I8B + j = narr + 1_I8B + + do + j = j - 1_I8B + do + if (arr(j) <= x) exit + j = j - 1_I8B + end do + i = i + 1_I8B + do + if (arr(i) >= x) exit + i = i + 1_I8B + end do + if (i < j) then + ! exchange A(i) and A(j) + temp = arr(i) + arr(i) = arr(j) + arr(j) = temp + if (present(ind)) then + itmp = ind(i) + ind(i) = ind(j) + ind(j) = itmp + end if + else if (i == j) then + marker = i + 1_I8B + return + else + marker = i + return + endif + end do + + return + end subroutine partition_I4B_I8Bind + + pure subroutine partition_I8B_I8Bind(arr, marker, ind) + !! author: David A. Minton + !! + !! Partition function for quicksort on I8B type with I8B index + !! + implicit none + ! Arguments + integer(I8B), intent(inout), dimension(:) :: arr + integer(I8B), intent(inout), dimension(:), optional :: ind + integer(I8B), intent(out) :: marker + ! Internals + integer(I8B) :: i, j, itmp, narr, ipiv + integer(I8B) :: temp + integer(I8B) :: x ! pivot point + + narr = size(arr) + + ! Get center as pivot, as this is likely partially sorted + ipiv = narr / 2_I8B + x = arr(ipiv) + i = 0_I8B + j = narr + 1_I8B + + do + j = j - 1_I8B + do + if (arr(j) <= x) exit + j = j - 1_I8B + end do + i = i + 1_I8B + do + if (arr(i) >= x) exit + i = i + 1_I8B + end do + if (i < j) then + ! exchange A(i) and A(j) + temp = arr(i) + arr(i) = arr(j) + arr(j) = temp + if (present(ind)) then + itmp = ind(i) + ind(i) = ind(j) + ind(j) = itmp + end if + else if (i == j) then + marker = i + 1_I8B + return + else + marker = i + return + endif + end do + + return + end subroutine partition_I8B_I8Bind + module pure subroutine util_sort_sp(arr) !! author: David A. Minton @@ -662,6 +857,26 @@ module pure subroutine util_sort_rearrange_arr_I4B(arr, ind, n) return end subroutine util_sort_rearrange_arr_I4B + module pure subroutine util_sort_rearrange_arr_I4B_I8Bind(arr, ind, n) + !! author: David A. Minton + !! + !! Rearrange a single array of integers in-place from an index list. + implicit none + ! Arguments + integer(I4B), dimension(:), allocatable, intent(inout) :: arr !! Destination array + integer(I8B), dimension(:), intent(in) :: ind !! Index to rearrange against + integer(I8B), intent(in) :: n !! Number of elements in arr and ind to rearrange + ! Internals + integer(I4B), dimension(:), allocatable :: tmp !! Temporary copy of array used during rearrange operation + + if (.not. allocated(arr) .or. n <= 0_I8B) return + allocate(tmp, mold=arr) + tmp(1:n) = arr(ind) + call move_alloc(tmp, arr) + + return + end subroutine util_sort_rearrange_arr_I4B_I8Bind + module pure subroutine util_sort_rearrange_arr_logical(arr, ind, n) !! author: David A. Minton @@ -684,6 +899,27 @@ module pure subroutine util_sort_rearrange_arr_logical(arr, ind, n) end subroutine util_sort_rearrange_arr_logical + module pure subroutine util_sort_rearrange_arr_logical_I8Bind(arr, ind, n) + !! author: David A. Minton + !! + !! Rearrange a single array of logicals in-place from an index list. + implicit none + ! Arguments + logical, dimension(:), allocatable, intent(inout) :: arr !! Destination array + integer(I8B), dimension(:), intent(in) :: ind !! Index to rearrange against + integer(I8B), intent(in) :: n !! Number of elements in arr and ind to rearrange + ! Internals + logical, dimension(:), allocatable :: tmp !! Temporary copy of array used during rearrange operation + + if (.not. allocated(arr) .or. n <= 0) return + allocate(tmp, mold=arr) + tmp(1:n) = arr(ind) + call move_alloc(tmp, arr) + + return + end subroutine util_sort_rearrange_arr_logical_I8Bind + + module subroutine util_sort_rearrange_arr_info(arr, ind, n) !! author: David A. Minton !! diff --git a/src/walltime/walltime.f90 b/src/walltime/walltime.f90 index a8d82c2b5..ba9e8ab57 100644 --- a/src/walltime/walltime.f90 +++ b/src/walltime/walltime.f90 @@ -29,17 +29,19 @@ module subroutine walltime_stop(self) end subroutine walltime_stop - module subroutine walltime_report(self, nsubsteps, message) + module subroutine walltime_report(self, message, nsubsteps) !! author: David A. Minton !! !! Prints the elapsed time information to the terminal implicit none ! Arguments class(walltimer), intent(inout) :: self !! Walltimer object - integer(I4B), intent(in) :: nsubsteps !! Number of substeps used to compute the time per step character(len=*), intent(in) :: message !! Message to prepend to the wall time terminal output + integer(I4B), optional, intent(in) :: nsubsteps !! Number of substeps used to compute the time per step ! Internals - character(len=*), parameter :: walltimefmt = '" Total wall time: ", es12.5, "; Interval wall time: ", es12.5, "; Interval wall time/step: ", es12.5' + character(len=*), parameter :: nosubstepfmt = '" Total wall time: ", es12.5, "; Interval wall time: ", es12.5 ' + character(len=*), parameter :: substepfmt = '" Total wall time: ", es12.5, "; Interval wall time: ", es12.5, ";' //& + 'Interval wall time/step: ", es12.5' character(len=STRMAX) :: fmt integer(I8B) :: count_delta_step, count_delta_main, count_now real(DP) :: wall_main !! Value of total elapsed time at the end of a timed step @@ -56,10 +58,15 @@ module subroutine walltime_report(self, nsubsteps, message) count_delta_step = count_now - self%count_start_step wall_main = count_delta_main / (self%count_rate * 1.0_DP) wall_step = count_delta_step / (self%count_rate * 1.0_DP) - wall_per_substep = wall_step / nsubsteps + if (present(nsubsteps)) then + wall_per_substep = wall_step / nsubsteps + fmt = '("' // adjustl(message) // '",' // substepfmt // ')' + write(*,trim(adjustl(fmt))) wall_main, self%wall_step, wall_per_substep + else + fmt = '("' // adjustl(message) // '",' // nosubstepfmt // ')' + write(*,trim(adjustl(fmt))) wall_main, self%wall_step + end if - fmt = '("' // adjustl(message) // '",' // walltimefmt // ')' - write(*,trim(adjustl(fmt))) wall_main, self%wall_step, wall_per_substep return end subroutine walltime_report @@ -110,7 +117,10 @@ module subroutine walltime_start(self) integer(I8B) :: count_resume, count_delta - if (.not.self%main_is_started) call self%start_main() + if (.not.self%main_is_started) then + call self%reset() + call self%start_main() + end if if (self%is_paused) then ! Resume a paused step timer call system_clock(count_resume) @@ -196,7 +206,8 @@ module subroutine walltime_interaction_adapt(self, param, ninteractions, pl) write(cstr,*) self%count_stop_step - self%count_start_step - call io_log_one_message(logfile, adjustl(lstyle) // " " // trim(adjustl(cstr)) // " " // trim(adjustl(nstr)) // " " // trim(adjustl(mstr))) + call io_log_one_message(logfile, adjustl(lstyle) // " " // trim(adjustl(cstr)) // " " // & + trim(adjustl(nstr)) // " " // trim(adjustl(mstr))) if (self%stage == 2) then if (ladvanced_final) then @@ -204,7 +215,8 @@ module subroutine walltime_interaction_adapt(self, param, ninteractions, pl) else lstyle = standardstyle end if - call io_log_one_message(logfile, trim(adjustl(self%loopname)) // ": the fastest loop method tested is " // trim(adjustl(lstyle))) + call io_log_one_message(logfile, trim(adjustl(self%loopname)) // & + ": the fastest loop method tested is " // trim(adjustl(lstyle))) end if return