diff --git a/python/swiftest/swiftest/simulation_class.py b/python/swiftest/swiftest/simulation_class.py index c908d322d..991758441 100644 --- a/python/swiftest/swiftest/simulation_class.py +++ b/python/swiftest/swiftest/simulation_class.py @@ -283,29 +283,19 @@ def __init__(self,read_param: bool = False, execute. If false, will start a new run. If the file given by `output_file_name` exists, it will be replaced when the run is executed. Parameter input file equivalent: `OUT_STAT` - interaction_loops : {"TRIANGULAR","FLAT","ADAPTIVE"}, default "TRIANGULAR" + interaction_loops : {"TRIANGULAR","FLAT"}, default "TRIANGULAR" > *Swiftest Experimental feature* Specifies which algorithm to use for the computation of body-body gravitational forces. * "TRIANGULAR" : Upper-triangular double-loops . * "FLAT" : Body-body interation pairs are flattened into a 1-D array. - * "ADAPTIVE" : Periodically times the TRIANGULAR and FLAT methods and determines which one to use based on - the wall time to complete the loop. *Note:* Using ADAPTIVE means that bit-identical repeatability cannot - be assured, as the choice of algorithm depends on possible external factors that can influence the wall - time calculation. The exact floating-point results of the interaction will be different between the two - algorithm types. Parameter input file equivalent: `INTERACTION_LOOPS` - encounter_check_loops : {"TRIANGULAR","SORTSWEEP","ADAPTIVE"}, default "TRIANGULAR" + encounter_check_loops : {"TRIANGULAR","SORTSWEEP"}, default "TRIANGULAR" > *Swiftest Experimental feature* Specifies which algorithm to use for checking whether bodies are in a close encounter state or not. * "TRIANGULAR" : Upper-triangular double-loops. * "SORTSWEEP" : A Sort-Sweep algorithm is used to reduce the population of potential close encounter bodies. This algorithm is still in development, and does not necessarily speed up the encounter checking. Use with caution. - * "ADAPTIVE" : Periodically times the TRIANGULAR and SORTSWEEP methods and determines which one to use based - on the wall time to complete the loop. *Note:* Using ADAPTIVE means that bit-identical repeatability cannot - be assured, as the choice of algorithm depends on possible external factors that can influence the wall - time calculation. The exact floating-point results of the interaction will be different between the two - algorithm types. Parameter input file equivalent: `ENCOUNTER_CHECK` dask : bool, default False Use Dask to lazily load data (useful for very large datasets) @@ -1072,8 +1062,8 @@ def set_feature(self, rhill_present: bool | None = None, restart: bool | None = None, tides: bool | None = None, - interaction_loops: Literal["TRIANGULAR", "FLAT", "ADAPTIVE"] | None = None, - encounter_check_loops: Literal["TRIANGULAR", "SORTSWEEP", "ADAPTIVE"] | None = None, + interaction_loops: Literal["TRIANGULAR", "FLAT"] | None = None, + encounter_check_loops: Literal["TRIANGULAR", "SORTSWEEP"] | None = None, encounter_save: Literal["NONE", "TRAJECTORY", "CLOSEST", "BOTH"] | None = None, verbose: bool | None = None, simdir: str | os.PathLike = None, @@ -1128,28 +1118,18 @@ def set_feature(self, Includes big bodies when performing a discard (Swifter only) rhill_present: bool, optional Include the Hill's radius with the input files. - interaction_loops : {"TRIANGULAR","FLAT","ADAPTIVE"}, default "TRIANGULAR" + interaction_loops : {"TRIANGULAR","FLAT"}, default "TRIANGULAR" *Swiftest Experimental feature* Specifies which algorithm to use for the computation of body-body gravitational forces. * "TRIANGULAR" : Upper-triangular double-loops . * "FLAT" : Body-body interation pairs are flattened into a 1-D array. - * "ADAPTIVE" : Periodically times the TRIANGULAR and FLAT methods and determines which one to use based on - the wall time to complete the loop. *Note:* Using ADAPTIVE means that bit-identical repeatability cannot - be assured, as the choice of algorithm depends on possible external factors that can influence the wall - time calculation. The exact floating-point results of the interaction will be different between the two - algorithm types. - encounter_check_loops : {"TRIANGULAR","SORTSWEEP","ADAPTIVE"}, default "TRIANGULAR" + encounter_check_loops : {"TRIANGULAR","SORTSWEEP"}, default "TRIANGULAR" *Swiftest Experimental feature* Specifies which algorithm to use for checking whether bodies are in a close encounter state or not. * "TRIANGULAR" : Upper-triangular double-loops. * "SORTSWEEP" : A Sort-Sweep algorithm is used to reduce the population of potential close encounter bodies. This algorithm is still in development, and does not necessarily speed up the encounter checking. Use with caution. - * "ADAPTIVE" : Periodically times the TRIANGULAR and SORTSWEEP methods and determines which one to use based - on the wall time to complete the loop. *Note:* Using ADAPTIVE means that bit-identical repeatability cannot - be assured, as the choice of algorithm depends on possible external factors that can influence the wall - time calculation. The exact floating-point results of the interaction will be different between the two - algorithm types. tides: bool, optional Turns on tidal model (IN DEVELOPMENT - IGNORED) Yarkovsky: bool, optional @@ -1250,7 +1230,8 @@ def set_feature(self, update_list.append("restart") if interaction_loops is not None: - valid_vals = ["TRIANGULAR", "FLAT", "ADAPTIVE"] + valid_vals = ["TRIANGULAR", "FLAT"] + interaction_loops = interaction_loops.upper() if interaction_loops not in valid_vals: msg = f"{interaction_loops} is not a valid option for interaction loops." msg += f"\nMust be one of {valid_vals}" @@ -1262,7 +1243,8 @@ def set_feature(self, update_list.append("interaction_loops") if encounter_check_loops is not None: - valid_vals = ["TRIANGULAR", "SORTSWEEP", "ADAPTIVE"] + valid_vals = ["TRIANGULAR", "SORTSWEEP"] + encounter_check_loops = encounter_check_loops.upper() if encounter_check_loops not in valid_vals: msg = f"{encounter_check_loops} is not a valid option for interaction loops." msg += f"\nMust be one of {valid_vals}" diff --git a/src/base/base_module.f90 b/src/base/base_module.f90 index 480587b04..0a9429dbc 100644 --- a/src/base/base_module.f90 +++ b/src/base/base_module.f90 @@ -70,11 +70,8 @@ module base ! The following are used internally, and are not set by the user, but instead are determined by the input value of INTERACTION_LOOPS logical :: lflatten_interactions = .false. !! Use the flattened upper triangular matrix for pl-pl interaction loops - logical :: ladaptive_interactions = .false. !! Adaptive interaction loop is turned on (choose between TRIANGULAR and FLAT based on periodic timing tests) logical :: lencounter_sas_plpl = .false. !! Use the Sort and Sweep algorithm to prune the encounter list before checking for close encounters logical :: lencounter_sas_pltp = .false. !! Use the Sort and Sweep algorithm to prune the encounter list before checking for close encounters - logical :: ladaptive_encounters_plpl = .false. !! Adaptive encounter checking is turned on (choose between TRIANGULAR or SORTSWEEP based on periodic timing tests) - logical :: ladaptive_encounters_pltp = .false. !! Adaptive encounter checking is turned on (choose between TRIANGULAR or SORTSWEEP based on periodic timing tests) ! Logical flags to turn on or off various features of the code logical :: lrhill_present = .false. !! Hill radii are given as an input rather than calculated by the code (can be used to inflate close encounter regions manually) diff --git a/src/encounter/encounter_check.f90 b/src/encounter/encounter_check.f90 index 917fcae45..105aaf7fe 100644 --- a/src/encounter/encounter_check.f90 +++ b/src/encounter/encounter_check.f90 @@ -11,31 +11,29 @@ use swiftest contains - module subroutine encounter_check_all_plpl(param, npl, x, v, renc, dt, nenc, index1, index2, lvdotr) + module subroutine encounter_check_all_plpl(param, npl, r, 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 !! implicit none ! Arguments - class(base_parameters), intent(inout) :: param !! Current Swiftest run configuration parameter5s + class(base_parameters), intent(inout) :: param !! Current Swiftest run configuration parameter5s integer(I4B), intent(in) :: npl !! Total number of massive bodies - real(DP), dimension(:,:), intent(in) :: x !! Position vectors of massive bodies - 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), dimension(:,:), intent(in) :: r !! Position vectors of massive bodies + real(DP), dimension(:,:), intent(in) :: v !! Velocity vectors of massive bodies real(DP), intent(in) :: dt !! Step size 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 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 - - call encounter_check_all_triangular_plpl(npl, x, v, renc, dt, nenc, index1, index2, lvdotr) + if (param%lencounter_sas_plpl) then + call encounter_check_all_sort_and_sweep_plpl(npl, r, v, renc, dt, nenc, index1, index2, lvdotr) + else + call encounter_check_all_triangular_plpl(npl, r, v, renc, dt, nenc, index1, index2, lvdotr) + end if return end subroutine encounter_check_all_plpl @@ -48,7 +46,7 @@ module subroutine encounter_check_all_plplm(param, nplm, nplt, rplm, vplm, rplt, !! implicit none ! Arguments - class(base_parameters), intent(inout) :: param !! Current Swiftest run configuration parameter5s + class(base_parameters), intent(inout) :: param !! Current Swiftest run configuration parameter5s integer(I4B), intent(in) :: nplm !! Total number of fully interacting massive bodies integer(I4B), intent(in) :: nplt !! Total number of partially interacting masive bodies (GM < GMTINY) real(DP), dimension(:,:), intent(in) :: rplm !! Position vectors of fully interacting massive bodies @@ -79,13 +77,16 @@ module subroutine encounter_check_all_plplm(param, nplm, nplt, rplm, vplm, rplt, allocate(tmp_param, source=param) - ! Turn off adaptive encounter checks for the pl-pl group - tmp_param%ladaptive_encounters_plpl = .false. - ! Start with the pl-pl group call encounter_check_all_plpl(tmp_param, nplm, rplm, vplm, rencm, dt, nenc, index1, index2, lvdotr) - call encounter_check_all_triangular_plplm(nplm, nplt, rplm, vplm, rplt, vplt, rencm, renct, dt, plmplt_nenc, plmplt_index1, plmplt_index2, plmplt_lvdotr) + if (param%lencounter_sas_plpl) then + call encounter_check_all_sort_and_sweep_plplm(nplm, nplt, rplm, vplm, rplt, vplt, rencm, renct, dt, & + plmplt_nenc, plmplt_index1, plmplt_index2, plmplt_lvdotr) + else + call encounter_check_all_triangular_plplm(nplm, nplt, rplm, vplm, rplt, vplt, rencm, renct, dt, & + plmplt_nenc, plmplt_index1, plmplt_index2, plmplt_lvdotr) + end if if (plmplt_nenc > 0) then ! Consolidate the two lists allocate(itmp(nenc+plmplt_nenc)) @@ -113,19 +114,19 @@ module subroutine encounter_check_all_plplm(param, nplm, nplt, rplm, vplm, rplt, end subroutine encounter_check_all_plplm - module subroutine encounter_check_all_pltp(param, npl, ntp, rpl, vpl, xtp, vtp, renc, dt, nenc, index1, index2, lvdotr) + module subroutine encounter_check_all_pltp(param, npl, ntp, rpl, vpl, rtp, 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 !! implicit none ! Arguments - class(base_parameters), intent(inout) :: param !! Current Swiftest run configuration parameter5s + class(base_parameters), intent(inout) :: param !! Current Swiftest run configuration parameter5s integer(I4B), intent(in) :: npl !! Total number of massive bodies integer(I4B), intent(in) :: ntp !! Total number of test particles real(DP), dimension(:,:), intent(in) :: rpl !! Position vectors of massive bodies real(DP), dimension(:,:), intent(in) :: vpl !! Velocity vectors of massive bodies - real(DP), dimension(:,:), intent(in) :: xtp !! Position vectors of test particlse + real(DP), dimension(:,:), intent(in) :: rtp !! Position vectors of test particlse 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 @@ -133,19 +134,18 @@ module subroutine encounter_check_all_pltp(param, npl, ntp, rpl, vpl, xtp, vtp, 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 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 :: lsecond = .false. - integer(I8B) :: npltp = 0_I8B - call encounter_check_all_triangular_pltp(npl, ntp, rpl, vpl, xtp, vtp, renc, dt, nenc, index1, index2, lvdotr) + if (param%lencounter_sas_pltp) then + call encounter_check_all_sort_and_sweep_pltp(npl, ntp, rpl, vpl, rtp, vtp, renc, dt, nenc, index1, index2, lvdotr) + else + call encounter_check_all_triangular_pltp(npl, ntp, rpl, vpl, rtp, vtp, renc, dt, nenc, index1, index2, lvdotr) + end if return end subroutine encounter_check_all_pltp - subroutine encounter_check_all_sort_and_sweep_plpl(npl, x, v, renc, dt, nenc, index1, index2, lvdotr) + subroutine encounter_check_all_sort_and_sweep_plpl(npl, r, v, renc, dt, nenc, index1, index2, lvdotr) !! author: David A. Minton !! !! Check for encounters between massive bodies. @@ -155,7 +155,7 @@ subroutine encounter_check_all_sort_and_sweep_plpl(npl, x, v, renc, dt, nenc, in implicit none ! Arguments integer(I4B), intent(in) :: npl !! Total number of massive bodies - real(DP), dimension(:,:), intent(in) :: x !! Position vectors of massive bodies + real(DP), dimension(:,:), intent(in) :: r !! Position vectors of massive bodies 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 @@ -164,10 +164,11 @@ subroutine encounter_check_all_sort_and_sweep_plpl(npl, x, v, renc, dt, nenc, in integer(I4B), dimension(:), allocatable, intent(out) :: index2 !! List of indices for body 2 in each encounter logical, dimension(:), allocatable, intent(out) :: lvdotr !! Logical flag indicating the sign of v .dot. x ! Internals - integer(I4B) :: dim, n + integer(I4B) :: i, n integer(I4B), save :: npl_last = 0 type(encounter_bounding_box), save :: boundingbox - integer(I2B), dimension(npl) :: vshift_min, vshift_max + real(DP), dimension(npl) :: rmin,rmax + real(DP) :: rmag if (npl == 0) return @@ -178,23 +179,15 @@ subroutine encounter_check_all_sort_and_sweep_plpl(npl, x, v, renc, dt, nenc, in npl_last = npl end if - !$omp parallel do default(private) schedule(static) & - !$omp shared(x, v, renc, boundingbox) & - !$omp firstprivate(dt, npl, n) - do dim = 1, SWEEPDIM - where(v(dim,1:npl) < 0.0_DP) - vshift_min(1:npl) = 1 - vshift_max(1:npl) = 0 - elsewhere - vshift_min(1:npl) = 0 - vshift_max(1:npl) = 1 - end where - call boundingbox%aabb(dim)%sort(npl, [x(dim,1:npl) - renc(1:npl) + vshift_min(1:npl) * v(dim,1:npl) * dt, & - x(dim,1:npl) + renc(1:npl) + vshift_max(1:npl) * v(dim,1:npl) * dt]) + do concurrent (i = 1:npl) + rmag = .mag.r(:,i) + rmax(i) = rmag + RSWEEP_FACTOR * renc(i) + rmin(i) = rmag - RSWEEP_FACTOR * renc(i) end do - !$omp end parallel do - call boundingbox%sweep(npl, x, v, renc, dt, nenc, index1, index2, lvdotr) + call boundingbox%aabb%sort(npl, [rmin,rmax]) + + call boundingbox%sweep(npl, r, v, renc, dt, nenc, index1, index2, lvdotr) return end subroutine encounter_check_all_sort_and_sweep_plpl @@ -224,10 +217,10 @@ subroutine encounter_check_all_sort_and_sweep_plplm(nplm, nplt, rplm, vplm, rplt logical, dimension(:), allocatable, intent(out) :: lvdotr !! Logical flag indicating the sign of v .dot. x ! Internals type(encounter_bounding_box), save :: boundingbox - integer(I4B) :: dim, n, ntot + integer(I4B) :: i, n, ntot integer(I4B), save :: ntot_last = 0 - integer(I2B), dimension(nplm) :: vplmshift_min, vplmshift_max - integer(I2B), dimension(nplt) :: vpltshift_min, vpltshift_max + real(DP), dimension(nplm+nplt) :: rmin,rmax + real(DP) :: rmag ! If this is the first time through, build the index lists if ((nplm == 0) .or. (nplt == 0)) return @@ -239,32 +232,18 @@ subroutine encounter_check_all_sort_and_sweep_plplm(nplm, nplt, rplm, vplm, rplt ntot_last = ntot end if - !$omp parallel do default(private) schedule(static) & - !$omp shared(rplm, rplt, vplm, vplt, rencm, renct, boundingbox) & - !$omp firstprivate(dt, nplm, nplt, ntot) - do dim = 1, SWEEPDIM - where(vplm(dim,1:nplm) < 0.0_DP) - vplmshift_min(1:nplm) = 1 - vplmshift_max(1:nplm) = 0 - elsewhere - vplmshift_min(1:nplm) = 0 - vplmshift_max(1:nplm) = 1 - end where - - where(vplt(dim,1:nplt) < 0.0_DP) - vpltshift_min(1:nplt) = 1 - vpltshift_max(1:nplt) = 0 - elsewhere - vpltshift_min(1:nplt) = 0 - vpltshift_max(1:nplt) = 1 - end where - - call boundingbox%aabb(dim)%sort(ntot, [rplm(dim,1:nplm) - rencm(1:nplm) + vplmshift_min(1:nplm) * vplm(dim,1:nplm) * dt, & - rplt(dim,1:nplt) - renct(1:nplt) + vpltshift_min(1:nplt) * vplt(dim,1:nplt) * dt, & - rplm(dim,1:nplm) + rencm(1:nplm) + vplmshift_max(1:nplm) * vplm(dim,1:nplm) * dt, & - rplt(dim,1:nplt) + renct(1:nplt) + vpltshift_max(1:nplt) * vplt(dim,1:nplt) * dt]) + do concurrent (i = 1:nplm) + rmag = .mag.rplm(:,i) + rmax(i) = rmag + RSWEEP_FACTOR * rencm(i) + rmin(i) = rmag - RSWEEP_FACTOR * rencm(i) end do - !$omp end parallel do + do concurrent (i = 1:nplt) + rmag = .mag.rplt(:,i) + rmax(nplm+i) = rmag + RSWEEP_FACTOR * renct(i) + rmin(nplm+i) = rmag - RSWEEP_FACTOR * renct(i) + end do + + call boundingbox%aabb%sort(ntot, [rmin, rmax]) call boundingbox%sweep(nplm, nplt, rplm, vplm, rplt, vplt, rencm, renct, dt, nenc, index1, index2, lvdotr) @@ -272,7 +251,7 @@ subroutine encounter_check_all_sort_and_sweep_plplm(nplm, nplt, rplm, vplm, rplt end subroutine encounter_check_all_sort_and_sweep_plplm - subroutine encounter_check_all_sort_and_sweep_pltp(npl, ntp, rpl, vpl, xtp, vtp, rencpl, dt, nenc, index1, index2, lvdotr) + subroutine encounter_check_all_sort_and_sweep_pltp(npl, ntp, rpl, vpl, rtp, vtp, rencpl, dt, nenc, index1, index2, lvdotr) !! author: David A. Minton !! !! Check for encounters between massive bodies and test particles. @@ -285,7 +264,7 @@ subroutine encounter_check_all_sort_and_sweep_pltp(npl, ntp, rpl, vpl, xtp, vtp, integer(I4B), intent(in) :: ntp !! Total number of test particles real(DP), dimension(:,:), intent(in) :: rpl !! Position vectors of massive bodies 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) :: rtp !! Position vectors of massive bodies real(DP), dimension(:,:), intent(in) :: vtp !! Velocity vectors of massive bodies real(DP), dimension(:), intent(in) :: rencpl !! Critical radii of massive bodies that defines an encounter real(DP), intent(in) :: dt !! Step size @@ -295,11 +274,11 @@ subroutine encounter_check_all_sort_and_sweep_pltp(npl, ntp, rpl, vpl, xtp, vtp, logical, dimension(:), allocatable, intent(out) :: lvdotr !! Logical flag indicating the sign of v .dot. x ! Internals type(encounter_bounding_box), save :: boundingbox - integer(I4B) :: dim, n, ntot + integer(I4B) :: i, n, ntot integer(I4B), save :: ntot_last = 0 - integer(I2B), dimension(npl) :: vplshift_min, vplshift_max - integer(I2B), dimension(ntp) :: vtpshift_min, vtpshift_max - real(DP), dimension(ntp) :: renctp + real(DP), dimension(npl+ntp) :: rmin,rmax + real(DP), dimension(ntp) :: renctp + real(DP) :: rmag ! If this is the first time through, build the index lists if ((ntp == 0) .or. (npl == 0)) return @@ -312,35 +291,21 @@ subroutine encounter_check_all_sort_and_sweep_pltp(npl, ntp, rpl, vpl, xtp, vtp, end if renctp(:) = 0.0_DP - - !$omp parallel do default(private) schedule(static) & - !$omp shared(rpl, xtp, vpl, vtp, rencpl, renctp, boundingbox) & - !$omp firstprivate(dt, npl, ntp, ntot) - do dim = 1, SWEEPDIM - where(vpl(dim,1:npl) < 0.0_DP) - vplshift_min(1:npl) = 1 - vplshift_max(1:npl) = 0 - elsewhere - vplshift_min(1:npl) = 0 - vplshift_max(1:npl) = 1 - end where - - where(vtp(dim,1:ntp) < 0.0_DP) - vtpshift_min(1:ntp) = 1 - vtpshift_max(1:ntp) = 0 - elsewhere - vtpshift_min(1:ntp) = 0 - vtpshift_max(1:ntp) = 1 - end where - - call boundingbox%aabb(dim)%sort(ntot, [rpl(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, & - rpl(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]) + + do concurrent (i = 1:npl) + rmag = .mag.rpl(:,i) + rmax(i) = rmag + RSWEEP_FACTOR * rencpl(i) + rmin(i) = rmag - RSWEEP_FACTOR * rencpl(i) end do - !$omp end parallel do + do concurrent (i = 1:ntp) + rmag = .mag.rtp(:,i) + rmax(npl+i) = rmag + RSWEEP_FACTOR * renctp(i) + rmin(npl+i) = rmag - RSWEEP_FACTOR * renctp(i) + end do + + call boundingbox%aabb%sort(ntot, [rmin, rmax]) - call boundingbox%sweep(npl, ntp, rpl, vpl, xtp, vtp, rencpl, renctp, dt, nenc, index1, index2, lvdotr) + call boundingbox%sweep(npl, ntp, rpl, vpl, rtp, vtp, rencpl, renctp, dt, nenc, index1, index2, lvdotr) return end subroutine encounter_check_all_sort_and_sweep_pltp @@ -445,7 +410,7 @@ subroutine encounter_check_all_triangular_one(i, n, xi, yi, zi, vxi, vyi, vzi, x end subroutine encounter_check_all_triangular_one - subroutine encounter_check_all_triangular_plpl(npl, x, v, renc, dt, nenc, index1, index2, lvdotr) + subroutine encounter_check_all_triangular_plpl(npl, r, 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 @@ -453,9 +418,9 @@ subroutine encounter_check_all_triangular_plpl(npl, x, v, renc, dt, nenc, index1 implicit none ! Arguments integer(I4B), intent(in) :: npl !! Total number of massive bodies - real(DP), dimension(:,:), intent(in) :: x !! Position vectors of massive bodies + real(DP), dimension(:,:), intent(in) :: r !! Position vectors of massive bodies 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), dimension(:), intent(in) :: renc !! Critical radii of massive bodies that defines an encounter real(DP), intent(in) :: dt !! Step size 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 @@ -469,12 +434,12 @@ subroutine encounter_check_all_triangular_plpl(npl, x, v, renc, dt, nenc, index1 call swiftest_util_index_array(ind_arr, npl) !$omp parallel do default(private) schedule(static)& - !$omp shared(x, v, renc, lenc, ind_arr) & + !$omp shared(r, v, renc, lenc, ind_arr) & !$omp firstprivate(npl, dt) do i = 1,npl - call encounter_check_all_triangular_one(i, npl, x(1,i), x(2,i), x(3,i), & + call encounter_check_all_triangular_one(i, npl, r(1,i), r(2,i), r(3,i), & v(1,i), v(2,i), v(3,i), & - x(1,:), x(2,:), x(3,:), & + r(1,:), r(2,:), r(3,:), & v(1,:), v(2,:), v(3,:), & renc(i), renc(:), dt, ind_arr(:), lenc(i)) if (lenc(i)%nenc > 0) lenc(i)%index1(:) = i @@ -483,7 +448,6 @@ subroutine encounter_check_all_triangular_plpl(npl, x, v, renc, dt, nenc, index1 call encounter_check_collapse_ragged_list(lenc, npl, nenc, index1, index2, lvdotr) - nenc = 0 return end subroutine encounter_check_all_triangular_plpl @@ -820,9 +784,8 @@ module subroutine encounter_check_sweep_aabb_double_list(self, n1, n2, r1, v1, r ! Internals integer(I4B) :: ii, i, ntot, nbox, dim logical, dimension(n1+n2) :: loverlap - logical, dimension(SWEEPDIM,n1+n2) :: loverlap_by_dimension - logical, dimension(SWEEPDIM,2*(n1+n2)) :: llist1 - integer(I4B), dimension(SWEEPDIM,2*(n1+n2)) :: ext_ind + logical, dimension(2*(n1+n2)) :: llist1 + integer(I4B), dimension(2*(n1+n2)) :: ext_ind type(collision_list_pltp), dimension(n1+n2) :: lenc !! Array of encounter lists (one encounter list per body) integer(I4B), dimension(:), allocatable, save :: ind_arr integer(I8B) :: ibeg, iend @@ -831,39 +794,33 @@ module subroutine encounter_check_sweep_aabb_double_list(self, n1, n2, r1, v1, r ntot = n1 + n2 call swiftest_util_index_array(ind_arr, ntot) - 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 + loverlap(:) = (self%aabb%ibeg(:) + 1_I8B) < (self%aabb%iend(:) - 1_I8B) + where(self%aabb%ind(:) > ntot) + ext_ind(:) = self%aabb%ind(:) - ntot + elsewhere + ext_ind(:) = self%aabb%ind(:) + endwhere + llist1(:) = ext_ind(:) <= n1 + where(.not.llist1(:)) ext_ind(:) = ext_ind(:) - n1 - loverlap(:) = loverlap_by_dimension(1,:) - do dim = 2, SWEEPDIM - loverlap(:) = loverlap(:) .and. loverlap_by_dimension(dim,:) - end do dim = 1 - where(llist1(dim,:)) - xind(:) = r1(1,ext_ind(dim,:)) - yind(:) = r1(2,ext_ind(dim,:)) - zind(:) = r1(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,:)) + where(llist1(:)) + xind(:) = r1(1,ext_ind(:)) + yind(:) = r1(2,ext_ind(:)) + zind(:) = r1(3,ext_ind(:)) + vxind(:) = v1(1,ext_ind(:)) + vyind(:) = v1(2,ext_ind(:)) + vzind(:) = v1(3,ext_ind(:)) + rencind(:) = renc1(ext_ind(:)) elsewhere - xind(:) = r2(1,ext_ind(dim,:)) - yind(:) = r2(2,ext_ind(dim,:)) - zind(:) = r2(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,:)) + xind(:) = r2(1,ext_ind(:)) + yind(:) = r2(2,ext_ind(:)) + zind(:) = r2(3,ext_ind(:)) + vxind(:) = v2(1,ext_ind(:)) + vyind(:) = v2(2,ext_ind(:)) + vzind(:) = v2(3,ext_ind(:)) + rencind(:) = renc2(ext_ind(:)) endwhere where(.not.loverlap(:)) lenc(:)%nenc = 0 @@ -875,14 +832,14 @@ module subroutine encounter_check_sweep_aabb_double_list(self, n1, n2, r1, v1, r !$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 + ibeg = self%aabb%ibeg(i) + 1_I8B + iend = self%aabb%iend(i) - 1_I8B nbox = iend - ibeg + 1 call encounter_check_all_sweep_one(i, nbox, r1(1,i), r1(2,i), r1(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) + renc1(i), rencind(ibeg:iend), dt, ext_ind(ibeg:iend), & + .not.llist1(ibeg:iend), lenc(i)%nenc, lenc(i)%index1, lenc(i)%index2, lenc(i)%lvdotr) end if end do !$omp end do nowait @@ -891,15 +848,15 @@ module subroutine encounter_check_sweep_aabb_double_list(self, n1, n2, r1, v1, r !$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 + ibeg = self%aabb%ibeg(i) + 1_I8B + iend = self%aabb%iend(i) - 1_I8B nbox = iend - ibeg + 1 ii = i - n1 - call encounter_check_all_sweep_one(ii, nbox, r1(1,ii), r1(2,ii), r1(3,ii), v2(1,ii), v2(2,ii), v2(3,ii), & + call encounter_check_all_sweep_one(ii, nbox, r2(1,ii), r2(2,ii), r2(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) + renc2(ii), rencind(ibeg:iend), dt, ext_ind(ibeg:iend), & + llist1(ibeg:iend), lenc(i)%nenc, lenc(i)%index2, lenc(i)%index1, lenc(i)%lvdotr) end if end do !$omp end do nowait @@ -914,7 +871,7 @@ module subroutine encounter_check_sweep_aabb_double_list(self, n1, n2, r1, v1, r end subroutine encounter_check_sweep_aabb_double_list - module subroutine encounter_check_sweep_aabb_single_list(self, n, x, v, renc, dt, nenc, index1, index2, lvdotr) + module subroutine encounter_check_sweep_aabb_single_list(self, n, r, v, renc, dt, nenc, index1, index2, lvdotr) !! author: David A. Minton !! !! Sweeps the sorted bounding box extents and returns the true encounters (combines broad and narrow phases) @@ -923,7 +880,7 @@ module subroutine encounter_check_sweep_aabb_single_list(self, n, x, v, renc, dt ! Arguments 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) :: r, 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 @@ -936,7 +893,7 @@ module subroutine encounter_check_sweep_aabb_single_list(self, n, x, v, renc, dt 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 + integer(I4B), dimension(2*n) :: ext_ind type(collision_list_plpl), dimension(n) :: lenc !! Array of encounter lists (one encounter list per body) integer(I4B), dimension(:), allocatable, save :: ind_arr integer(I8B) :: ibeg, iend @@ -946,36 +903,36 @@ module subroutine encounter_check_sweep_aabb_single_list(self, n, x, v, renc, dt ! 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 + where(self%aabb%ind(:) > n) + ext_ind(:) = self%aabb%ind(:) - n elsewhere - ext_ind(dim,:) = self%aabb(1)%ind(:) + ext_ind(:) = self%aabb%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,:)) + xind(:) = r(1,ext_ind(:)) + yind(:) = r(2,ext_ind(:)) + zind(:) = r(3,ext_ind(:)) + vxind(:) = v(1,ext_ind(:)) + vyind(:) = v(2,ext_ind(:)) + vzind(:) = v(3,ext_ind(:)) + rencind(:) = renc(ext_ind(:)) - loverlap(:) = (self%aabb(dim)%ibeg(:) + 1_I8B) < (self%aabb(dim)%iend(:) - 1_I8B) + loverlap(:) = (self%aabb%ibeg(:) + 1_I8B) < (self%aabb%iend(:) - 1_I8B) where(.not.loverlap(:)) lenc(:)%nenc = 0 !$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 shared(self, ext_ind, lenc, loverlap, r, v, renc, xind, yind, zind, vxind, vyind, vzind, rencind) & !$omp firstprivate(n, dt, dim) do i = 1, n if (loverlap(i)) then - ibeg = self%aabb(dim)%ibeg(i) + 1_I8B - iend = self%aabb(dim)%iend(i) - 1_I8B + ibeg = self%aabb%ibeg(i) + 1_I8B + iend = self%aabb%iend(i) - 1_I8B nbox = int(iend - ibeg + 1, kind=I4B) 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), & + call encounter_check_all_sweep_one(i, nbox, r(1,i), r(2,i), r(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), & + renc(i), rencind(ibeg:iend), dt, ext_ind(ibeg:iend), & lencounteri(ibeg:iend), lenc(i)%nenc, lenc(i)%index1, lenc(i)%index2, lenc(i)%lvdotr) end if end do diff --git a/src/encounter/encounter_module.f90 b/src/encounter/encounter_module.f90 index 720627485..f3c24c763 100644 --- a/src/encounter/encounter_module.f90 +++ b/src/encounter/encounter_module.f90 @@ -17,8 +17,8 @@ module encounter implicit none public - integer(I4B), parameter :: SWEEPDIM = 3 character(len=*), parameter :: ENCOUNTER_OUTFILE = 'encounters.nc' !! Name of NetCDF output file for encounter information + real(DP), parameter :: RSWEEP_FACTOR = 1.1_DP type, abstract :: encounter_list integer(I8B) :: nenc = 0 !! Total number of encounters @@ -95,7 +95,7 @@ module encounter type encounter_bounding_box - type(encounter_bounding_box_1D), dimension(SWEEPDIM) :: aabb + type(encounter_bounding_box_1D) :: aabb contains procedure :: dealloc => encounter_util_dealloc_bounding_box !! Deallocates all allocatables procedure :: setup => encounter_util_setup_aabb !! Setup a new axis-aligned bounding box structure @@ -107,12 +107,12 @@ module encounter interface - module subroutine encounter_check_all_plpl(param, npl, x, v, renc, dt, nenc, index1, index2, lvdotr) + module subroutine encounter_check_all_plpl(param, npl, r, v, renc, dt, nenc, index1, index2, lvdotr) use base, only: base_parameters implicit none class(base_parameters), intent(inout) :: param !! Current Swiftest run configuration parameter5s integer(I4B), intent(in) :: npl !! Total number of massive bodies - real(DP), dimension(:,:), intent(in) :: x !! Position vectors of massive bodies + real(DP), dimension(:,:), intent(in) :: r !! Position vectors of massive bodies 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 @@ -142,7 +142,7 @@ module subroutine encounter_check_all_plplm(param, nplm, nplt, rplm, vplm, rplt, 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, rpl, vpl, xtp, vtp, renc, dt, nenc, index1, index2, lvdotr) + module subroutine encounter_check_all_pltp(param, npl, ntp, rpl, vpl, rtp, vtp, renc, dt, nenc, index1, index2, lvdotr) use base, only: base_parameters implicit none class(base_parameters), intent(inout) :: param !! Current Swiftest run configuration parameter5s @@ -150,7 +150,7 @@ module subroutine encounter_check_all_pltp(param, npl, ntp, rpl, vpl, xtp, vtp, integer(I4B), intent(in) :: ntp !! Total number of test particles real(DP), dimension(:,:), intent(in) :: rpl !! Position vectors of massive bodies 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) :: rtp !! 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), intent(in) :: dt !! Step size @@ -205,11 +205,11 @@ module subroutine encounter_check_sweep_aabb_double_list(self, n1, n2, r1, v1, r 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, x, v, renc, dt, nenc, index1, index2, lvdotr) + module subroutine encounter_check_sweep_aabb_single_list(self, n, r, 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) :: r, 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 diff --git a/src/encounter/encounter_util.f90 b/src/encounter/encounter_util.f90 index 941ed6e16..01848d571 100644 --- a/src/encounter/encounter_util.f90 +++ b/src/encounter/encounter_util.f90 @@ -104,9 +104,7 @@ module subroutine encounter_util_dealloc_bounding_box(self) ! Internals integer(I4B) :: i - do i = 1,NDIM - call self%aabb(i)%dealloc() - end do + call self%aabb%dealloc() return end subroutine encounter_util_dealloc_bounding_box @@ -351,26 +349,20 @@ module subroutine encounter_util_setup_aabb(self, n, n_last) next_last = 2 * n_last if (n > n_last) then ! The number of bodies has grown. Resize and append the new bodies - do dim = 1, SWEEPDIM - allocate(itmp(next)) - if (n_last > 0) itmp(1:next_last) = self%aabb(dim)%ind(1:next_last) - call move_alloc(itmp, self%aabb(dim)%ind) - self%aabb(dim)%ind(next_last+1:next) = [(k, k = next_last+1, next)] - end do + allocate(itmp(next)) + if (n_last > 0) itmp(1:next_last) = self%aabb%ind(1:next_last) + call move_alloc(itmp, self%aabb%ind) + self%aabb%ind(next_last+1:next) = [(k, k = next_last+1, next)] else ! The number of bodies has gone down. Resize and chop of the old indices - do dim = 1, SWEEPDIM - allocate(itmp(next)) - itmp(1:next) = pack(self%aabb(dim)%ind(1:next_last), self%aabb(dim)%ind(1:next_last) <= next) - call move_alloc(itmp, self%aabb(dim)%ind) - end do + allocate(itmp(next)) + itmp(1:next) = pack(self%aabb%ind(1:next_last), self%aabb%ind(1:next_last) <= next) + call move_alloc(itmp, self%aabb%ind) end if - do dim = 1, SWEEPDIM - if (allocated(self%aabb(dim)%ibeg)) deallocate(self%aabb(dim)%ibeg) - allocate(self%aabb(dim)%ibeg(n)) - if (allocated(self%aabb(dim)%iend)) deallocate(self%aabb(dim)%iend) - allocate(self%aabb(dim)%iend(n)) - end do + if (allocated(self%aabb%ibeg)) deallocate(self%aabb%ibeg) + allocate(self%aabb%ibeg(n)) + if (allocated(self%aabb%iend)) deallocate(self%aabb%iend) + allocate(self%aabb%iend(n)) return end subroutine encounter_util_setup_aabb diff --git a/src/rmvs/rmvs_encounter_check.f90 b/src/rmvs/rmvs_encounter_check.f90 index 00aafd1fb..0880cb136 100644 --- a/src/rmvs/rmvs_encounter_check.f90 +++ b/src/rmvs/rmvs_encounter_check.f90 @@ -40,7 +40,7 @@ module function rmvs_encounter_check_tp(self, param, nbody_system, dt) result(le select type(pl => nbody_system%pl) class is (rmvs_pl) - associate(tp => self, ntp => self%nbody, npl => pl%nbody) + associate(tp => self, ntp => self%nbody, npl => pl%nbody, cb => nbody_system%cb) tp%plencP(1:ntp) = 0 call encounter_check_all_pltp(param, npl, ntp, pl%rbeg, pl%vbeg, tp%rh, tp%vh, pl%renc, dt, & nenc, index1, index2, lvdotr) diff --git a/src/swiftest/swiftest_io.f90 b/src/swiftest/swiftest_io.f90 index 3e5efa606..9406e44c2 100644 --- a/src/swiftest/swiftest_io.f90 +++ b/src/swiftest/swiftest_io.f90 @@ -2270,72 +2270,42 @@ module subroutine swiftest_io_param_reader(self, unit, iotype, v_list, iostat, i end if select case(trim(adjustl(param%interaction_loops))) - case("ADAPTIVE") - param%ladaptive_interactions = .true. - param%lflatten_interactions = .true. - call swiftest_io_log_start(param, INTERACTION_TIMER_LOG_OUT, "Interaction loop timer logfile") - call swiftest_io_log_one_message(INTERACTION_TIMER_LOG_OUT, "Diagnostic values: loop style, time count, nplpl, metric") case("TRIANGULAR") - param%ladaptive_interactions = .false. param%lflatten_interactions = .false. case("FLAT") - param%ladaptive_interactions = .false. param%lflatten_interactions = .true. case default write(*,*) "Unknown value for parameter INTERACTION_LOOPS: -> ",trim(adjustl(param%interaction_loops)) - write(*,*) "Must be one of the following: TRIANGULAR, FLAT, or ADAPTIVE" - write(*,*) "Using default value of ADAPTIVE" - param%interaction_loops = "ADAPTIVE" - param%ladaptive_interactions = .true. - param%lflatten_interactions = .true. - call swiftest_io_log_start(param, INTERACTION_TIMER_LOG_OUT, "Interaction loop timer logfile") - call swiftest_io_log_one_message(INTERACTION_TIMER_LOG_OUT, "Diagnostic values: loop style, time count, nplpl, metric") + write(*,*) "Must be one of the following: TRIANGULAR or FLAT" + write(*,*) "Using default value of TRIANGULAR" + param%interaction_loops = "TRIANGULAR" + param%lflatten_interactions = .false. end select select case(trim(adjustl(param%encounter_check_plpl))) - case("ADAPTIVE") - param%ladaptive_encounters_plpl = .true. - param%lencounter_sas_plpl = .true. - call swiftest_io_log_start(param, ENCOUNTER_PLPL_TIMER_LOG_OUT, "Encounter check loop timer logfile") - call swiftest_io_log_one_message(ENCOUNTER_PLPL_TIMER_LOG_OUT, "Diagnostic values: loop style, time count, nplpl, metric") case("TRIANGULAR") - param%ladaptive_encounters_plpl = .false. param%lencounter_sas_plpl = .false. case("SORTSWEEP") - param%ladaptive_encounters_plpl = .false. param%lencounter_sas_plpl = .true. case default write(*,*) "Unknown value for parameter ENCOUNTER_CHECK_PLPL: -> ",trim(adjustl(param%encounter_check_plpl)) - write(*,*) "Must be one of the following: TRIANGULAR, SORTSWEEP, or ADAPTIVE" - write(*,*) "Using default value of ADAPTIVE" - param%encounter_check_plpl = "ADAPTIVE" - param%ladaptive_encounters_plpl = .true. - param%lencounter_sas_plpl = .true. - call swiftest_io_log_start(param, ENCOUNTER_PLPL_TIMER_LOG_OUT, "Encounter check loop timer logfile") - call swiftest_io_log_one_message(ENCOUNTER_PLPL_TIMER_LOG_OUT, "Diagnostic values: loop style, time count, nplpl, metric") + write(*,*) "Must be one of the following: TRIANGULAR or SORTSWEEP" + write(*,*) "Using default value of TRIANGULAR" + param%encounter_check_plpl = "TRIANGULAR" + param%lencounter_sas_plpl = .false. end select select case(trim(adjustl(param%encounter_check_pltp))) - case("ADAPTIVE") - param%ladaptive_encounters_pltp = .true. - param%lencounter_sas_pltp = .true. - call swiftest_io_log_start(param, ENCOUNTER_PLTP_TIMER_LOG_OUT, "Encounter check loop timer logfile") - call swiftest_io_log_one_message(ENCOUNTER_PLTP_TIMER_LOG_OUT, "Diagnostic values: loop style, time count, npltp, metric") case("TRIANGULAR") - param%ladaptive_encounters_pltp = .false. param%lencounter_sas_pltp = .false. case("SORTSWEEP") - param%ladaptive_encounters_pltp = .false. param%lencounter_sas_pltp = .true. case default write(*,*) "Unknown value for parameter ENCOUNTER_CHECK_PLTP: -> ",trim(adjustl(param%encounter_check_pltp)) - write(*,*) "Must be one of the following: TRIANGULAR, SORTSWEEP, or ADAPTIVE" - write(*,*) "Using default value of ADAPTIVE" - param%encounter_check_pltp = "ADAPTIVE" - param%ladaptive_encounters_pltp = .true. - param%lencounter_sas_pltp = .true. - call swiftest_io_log_start(param, ENCOUNTER_PLTP_TIMER_LOG_OUT, "Encounter check loop timer logfile") - call swiftest_io_log_one_message(ENCOUNTER_PLTP_TIMER_LOG_OUT, "Diagnostic values: loop style, time count, npltp, metric") + write(*,*) "Must be one of the following: TRIANGULAR or SORTSWEEP" + write(*,*) "Using default value of TRIANGULAR" + param%encounter_check_pltp = "TRIANGULAR" + param%lencounter_sas_pltp = .false. end select iostat = 0 diff --git a/src/swiftest/swiftest_module.f90 b/src/swiftest/swiftest_module.f90 index d360db56b..7913fd1ae 100644 --- a/src/swiftest/swiftest_module.f90 +++ b/src/swiftest/swiftest_module.f90 @@ -140,8 +140,8 @@ module swiftest procedure :: write_info => swiftest_io_netcdf_write_info_body !! Dump contents of particle information metadata to file procedure :: accel_obl => swiftest_obl_acc_body !! Compute the barycentric accelerations of bodies due to the oblateness of the central body procedure :: el2xv => swiftest_orbel_el2xv_vec !! Convert orbital elements to position and velocity vectors - procedure :: xv2el => swiftest_orbel_xv2el_vec !! Convert position and velocity vectors to orbital elements - procedure :: setup => swiftest_util_setup_body !! A constructor that sets the number of bodies and allocates all allocatable arrays + procedure :: xv2el => swiftest_orbel_xv2el_vec !! Convert position and velocity vectors to orbital elements + procedure :: setup => swiftest_util_setup_body !! A constructor that sets the number of bodies and allocates all allocatable arrays procedure :: accel_user => swiftest_user_kick_getacch_body !! Add user-supplied heliocentric accelerations to planets procedure :: append => swiftest_util_append_body !! Appends elements from one structure to another procedure :: dealloc => swiftest_util_dealloc_body !! Deallocates all allocatable arrays diff --git a/src/swiftest/swiftest_orbel.f90 b/src/swiftest/swiftest_orbel.f90 index bf72c91b2..8e351e911 100644 --- a/src/swiftest/swiftest_orbel.f90 +++ b/src/swiftest/swiftest_orbel.f90 @@ -680,7 +680,7 @@ real(DP) pure function swiftest_orbel_fhybrid(e,n) abn = n if(n < 0._DP) abn = -abn - if(abn < 0.636_DP * e -0.6_DP) then + if(abn < 0.636_DP * e - 0.6_DP) then swiftest_orbel_fhybrid = swiftest_orbel_flon(e,n) else swiftest_orbel_fhybrid = swiftest_orbel_fget(e,n) @@ -715,7 +715,7 @@ pure elemental module subroutine swiftest_orbel_xv2aeq(mu, px, py, pz, vx, vy, v q = 0.0_DP x = [px, py, pz] v = [vx, vy, vz] - r = sqrt(dot_product(x(:), x(:))) + r = .mag.x(:) v2 = dot_product(v(:), v(:)) hvec(:) = x(:) .cross. v(:) h2 = dot_product(hvec(:), hvec(:)) diff --git a/src/symba/symba_encounter_check.f90 b/src/symba/symba_encounter_check.f90 index ca7761ff0..895a9f34d 100644 --- a/src/symba/symba_encounter_check.f90 +++ b/src/symba/symba_encounter_check.f90 @@ -249,7 +249,7 @@ module function symba_encounter_check_tp(self, param, nbody_system, dt, irec) re lany_encounter = .false. if (self%nbody == 0) return - associate(tp => self, ntp => self%nbody, pl => nbody_system%pl, npl => nbody_system%pl%nbody) + associate(tp => self, ntp => self%nbody, pl => nbody_system%pl, npl => nbody_system%pl%nbody, cb => nbody_system%cb) call pl%set_renc(irec) call encounter_check_all_pltp(param, npl, ntp, pl%rh, pl%vb, tp%rh, tp%vb, pl%renc, dt, nenc, index1, index2, lvdotr) diff --git a/src/symba/symba_step.f90 b/src/symba/symba_step.f90 index c811b1968..6e72c3e95 100644 --- a/src/symba/symba_step.f90 +++ b/src/symba/symba_step.f90 @@ -27,6 +27,7 @@ module subroutine symba_step_system(self, param, t, dt) real(DP), intent(in) :: dt !! Current stepsize ! Internals logical :: lencounter + type(walltimer) :: timer1,timer2,timer3 !! Object used for computing elapsed wall time select type(pl => self%pl) class is (symba_pl) @@ -48,7 +49,6 @@ module subroutine symba_step_system(self, param, t, dt) end select end select end select - return end subroutine symba_step_system