From 22be9947586cfed521d230103f9fcc6eeafb1f06 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Mon, 20 Sep 2021 12:25:43 -0400 Subject: [PATCH] Refactored util_index to util_flatten as it is more descriptive and consistent with the rest of the terminology. Put in safety check in case allocating k_plpl fails --- src/fraggle/fraggle_generate.f90 | 4 +- src/fraggle/fraggle_util.f90 | 4 +- src/helio/helio_setup.f90 | 2 +- src/modules/fraggle_classes.f90 | 4 +- src/modules/swiftest_classes.f90 | 22 +++++----- src/modules/symba_classes.f90 | 8 ++-- src/modules/walltime_classes.f90 | 2 + src/modules/whm_classes.f90 | 7 ---- src/setup/setup.f90 | 2 +- src/symba/symba_util.f90 | 10 ++--- src/util/{util_index.f90 => util_flatten.f90} | 40 ++++++++++--------- src/whm/whm_setup.f90 | 2 +- 12 files changed, 53 insertions(+), 54 deletions(-) rename src/util/{util_index.f90 => util_flatten.f90} (78%) diff --git a/src/fraggle/fraggle_generate.f90 b/src/fraggle/fraggle_generate.f90 index 1c189c5ca..486aed2ea 100644 --- a/src/fraggle/fraggle_generate.f90 +++ b/src/fraggle/fraggle_generate.f90 @@ -18,7 +18,7 @@ module subroutine fraggle_generate_fragments(self, colliders, system, param, lfa class(fraggle_fragments), intent(inout) :: self !! Fraggle system object the outputs will be the fragmentation class(fraggle_colliders), intent(inout) :: colliders !! Fraggle colliders object containing the two-body equivalent values of the colliding bodies class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters logical, intent(out) :: lfailure !! Answers the question: Should this have been a merger instead? ! Internals integer(I4B) :: i @@ -140,7 +140,7 @@ module subroutine fraggle_generate_fragments(self, colliders, system, param, lfa call frag%set_original_scale(colliders) ! Restore the big array - if (lk_plpl) call pl%index(param) + if (lk_plpl) call pl%flatten(param) end associate call ieee_set_halting_mode(IEEE_ALL,fpe_halting_modes) ! Save the current halting modes so we can turn them off temporarily diff --git a/src/fraggle/fraggle_util.f90 b/src/fraggle/fraggle_util.f90 index df8b89c41..f65bd81c5 100644 --- a/src/fraggle/fraggle_util.f90 +++ b/src/fraggle/fraggle_util.f90 @@ -135,7 +135,7 @@ module subroutine fraggle_util_get_energy_momentum(self, colliders, system, para class(fraggle_fragments), intent(inout) :: self !! Fraggle fragment system object class(fraggle_colliders), intent(inout) :: colliders !! Fraggle collider system object class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object - class(swiftest_parameters), intent(in) :: param !! Current swiftest run configuration parameters + class(swiftest_parameters), intent(inout) :: param !! Current swiftest run configuration parameters logical, intent(in) :: lbefore !! Flag indicating that this the "before" state of the system, with colliders included and fragments excluded or vice versa ! Internals integer(I4B) :: i, nplm @@ -172,7 +172,7 @@ module subroutine fraggle_util_get_energy_momentum(self, colliders, system, para call fraggle_util_add_fragments_to_system(frag, colliders, tmpsys, tmpparam) end if - call tmpsys%pl%index(param) + call tmpsys%pl%flatten(param) call tmpsys%get_energy_and_momentum(param) diff --git a/src/helio/helio_setup.f90 b/src/helio/helio_setup.f90 index cf5f57d8a..5839628a0 100644 --- a/src/helio/helio_setup.f90 +++ b/src/helio/helio_setup.f90 @@ -16,7 +16,7 @@ module subroutine helio_setup_initialize_system(self, param) call self%pl%h2b(self%cb) call self%tp%h2b(self%cb) call self%pl%sort("mass", ascending=.false.) - call self%pl%index(param) + call self%pl%flatten(param) return end subroutine helio_setup_initialize_system diff --git a/src/modules/fraggle_classes.f90 b/src/modules/fraggle_classes.f90 index 5a8a4a392..2f6d75a84 100644 --- a/src/modules/fraggle_classes.f90 +++ b/src/modules/fraggle_classes.f90 @@ -104,7 +104,7 @@ module subroutine fraggle_generate_fragments(self, colliders, system, param, lfa class(fraggle_fragments), intent(inout) :: self !! Fraggle fragment system object class(fraggle_colliders), intent(inout) :: colliders !! Fraggle colliders object containing the two-body equivalent values of the colliding bodies class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters logical, intent(out) :: lfailure !! Answers the question: Should this have been a merger instead? end subroutine fraggle_generate_fragments @@ -249,7 +249,7 @@ module subroutine fraggle_util_get_energy_momentum(self, colliders, system, para class(fraggle_fragments), intent(inout) :: self !! Fraggle fragment system object class(fraggle_colliders), intent(inout) :: colliders !! Fraggle collider system object class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object - class(swiftest_parameters), intent(in) :: param !! Current swiftest run configuration parameters + class(swiftest_parameters), intent(inout) :: param !! Current swiftest run configuration parameters logical, intent(in) :: lbefore !! Flag indicating that this the "before" state of the system, with colliders included and fragments excluded or vice versa end subroutine fraggle_util_get_energy_momentum diff --git a/src/modules/swiftest_classes.f90 b/src/modules/swiftest_classes.f90 index af367c5ed..dc5eacc3c 100644 --- a/src/modules/swiftest_classes.f90 +++ b/src/modules/swiftest_classes.f90 @@ -133,7 +133,7 @@ module swiftest_classes logical :: loblatecb = .false. !! Calculate acceleration from oblate central body (automatically turns true if nonzero J2 is input) logical :: lrotation = .false. !! Include rotation states of big bodies logical :: ltides = .false. !! Include tidal dissipation - logical :: lflatten_interactions = .true. !! Use the flattened upper triangular matrix for pl-pl interactions (turning this on improves the speed but uses more memory) + logical :: lflatten_interactions = .false. !! Use the flattened upper triangular matrix for pl-pl interactions (turning this on improves the speed but uses more memory) ! Initial values to pass to the energy report subroutine (usually only used in the case of a restart, otherwise these will be updated with initial conditions values) real(DP) :: Eorbit_orig = 0.0_DP !! Initial orbital energy @@ -324,7 +324,7 @@ module swiftest_classes ! Massive body-specific concrete methods ! These are concrete because they are the same implemenation for all integrators procedure :: discard => discard_pl !! Placeholder method for discarding massive bodies - procedure :: index => util_index_eucl_plpl !! Sets up the (i, j) -> k indexing used for the single-loop blocking Euclidean distance matrix + procedure :: flatten => util_flatten_eucl_plpl !! Sets up the (i, j) -> k indexing used for the single-loop blocking Euclidean distance matrix procedure :: accel_int => kick_getacch_int_pl !! Compute direct cross (third) term heliocentric accelerations of massive bodies procedure :: accel_obl => obl_acc_pl !! Compute the barycentric accelerations of bodies due to the oblateness of the central body procedure :: setup => setup_pl !! A base constructor that sets the number of bodies and allocates and initializes all arrays @@ -561,34 +561,34 @@ module pure elemental subroutine drift_one(mu, px, py, pz, vx, vy, vz, dt, iflag integer(I4B), intent(out) :: iflag !! iflag : error status flag for Danby drift (0 = OK, nonzero = ERROR) end subroutine drift_one - module pure subroutine util_index_eucl_ij_to_k(n, i, j, k) - !$omp declare simd(util_index_eucl_ij_to_k) + module pure subroutine util_flatten_eucl_ij_to_k(n, i, j, k) + !$omp declare simd(util_flatten_eucl_ij_to_k) implicit none integer(I4B), intent(in) :: n !! Number of bodies integer(I4B), intent(in) :: i !! Index of the ith body integer(I4B), intent(in) :: j !! Index of the jth body integer(I8B), intent(out) :: k !! Index of the flattened matrix - end subroutine util_index_eucl_ij_to_k + end subroutine util_flatten_eucl_ij_to_k - module pure subroutine util_index_eucl_k_to_ij(n, k, i, j) + module pure subroutine util_flatten_eucl_k_to_ij(n, k, i, j) implicit none integer(I4B), intent(in) :: n !! Number of bodies integer(I8B), intent(in) :: k !! Index of the flattened matrix integer(I4B), intent(out) :: i !! Index of the ith body integer(I4B), intent(out) :: j !! Index of the jth body - end subroutine util_index_eucl_k_to_ij + end subroutine util_flatten_eucl_k_to_ij - module subroutine util_index_eucl_plpl(self, param) + module subroutine util_flatten_eucl_plpl(self, param) implicit none class(swiftest_pl), intent(inout) :: self !! Swiftest massive body object - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters end subroutine - module subroutine util_index_eucl_pltp(self, pl, param) + module subroutine util_flatten_eucl_pltp(self, pl, param) implicit none class(swiftest_tp), intent(inout) :: self !! Swiftest test particle object class(swiftest_pl), intent(in) :: pl !! Swiftest massive body object - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters end subroutine module pure subroutine gr_kick_getaccb_ns_body(self, system, param) diff --git a/src/modules/symba_classes.f90 b/src/modules/symba_classes.f90 index 0faefc3d1..52c66d293 100644 --- a/src/modules/symba_classes.f90 +++ b/src/modules/symba_classes.f90 @@ -67,7 +67,7 @@ module symba_classes type(symba_kinship), dimension(:), allocatable :: kin !! Array of merger relationship structures that can account for multiple pairwise mergers in a single step contains procedure :: make_colliders => symba_collision_make_colliders_pl !! When a single body is involved in more than one collision in a single step, it becomes part of a family - procedure :: index => symba_util_index_eucl_plpl !! Sets up the (i, j) -> k indexing used for the single-loop blocking Euclidean distance matrix + procedure :: flatten => symba_util_flatten_eucl_plpl !! Sets up the (i, j) -> k indexing used for the single-loop blocking Euclidean distance matrix procedure :: discard => symba_discard_pl !! Process massive body discards procedure :: drift => symba_drift_pl !! Method for Danby drift in Democratic Heliocentric coordinates. Sets the mask to the current recursion level procedure :: encounter_check => symba_encounter_check_pl !! Checks if massive bodies are going through close encounters with each other @@ -316,12 +316,12 @@ module function symba_collision_casemerge(system, param, colliders, frag) result integer(I4B) :: status !! Status flag assigned to this outcome end function symba_collision_casemerge - module subroutine symba_util_index_eucl_plpl(self, param) + module subroutine symba_util_flatten_eucl_plpl(self, param) use swiftest_classes, only : swiftest_parameters implicit none class(symba_pl), intent(inout) :: self !! SyMBA massive body object - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters - end subroutine symba_util_index_eucl_plpl + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters + end subroutine symba_util_flatten_eucl_plpl module subroutine symba_io_param_reader(self, unit, iotype, v_list, iostat, iomsg) implicit none diff --git a/src/modules/walltime_classes.f90 b/src/modules/walltime_classes.f90 index 82fda6886..c27a5b06a 100644 --- a/src/modules/walltime_classes.f90 +++ b/src/modules/walltime_classes.f90 @@ -6,6 +6,8 @@ module walltime_classes implicit none public + integer(I4B) :: INTERACTION_TIMER_CADENCE = 1000 !! Minimum number of steps to wait before timing an interaction loop in ADAPTIVE mode + type :: walltimer integer(I8B) :: count_rate !! Rate at wich the clock ticks integer(I8B) :: count_max !! Maximum value of the clock ticker diff --git a/src/modules/whm_classes.f90 b/src/modules/whm_classes.f90 index bd17f98e0..0bae1f654 100644 --- a/src/modules/whm_classes.f90 +++ b/src/modules/whm_classes.f90 @@ -108,13 +108,6 @@ module subroutine whm_drift_pl(self, system, param, dt) real(DP), intent(in) :: dt !! Stepsize end subroutine whm_drift_pl - module subroutine whm_util_index_eucl_plpl(self, param) - use swiftest_classes, only : swiftest_parameters - implicit none - class(whm_pl), intent(inout) :: self !! WHM massive body object - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters - end subroutine whm_util_index_eucl_plpl - !> Get heliocentric accelration of massive bodies module subroutine whm_kick_getacch_pl(self, system, param, t, lbeg) use swiftest_classes, only : swiftest_cb, swiftest_parameters diff --git a/src/setup/setup.f90 b/src/setup/setup.f90 index 18d8f2189..2dc57d952 100644 --- a/src/setup/setup.f90 +++ b/src/setup/setup.f90 @@ -175,7 +175,7 @@ module subroutine setup_initialize_system(self, param) call self%pl%el2xv(self%cb) call self%tp%el2xv(self%cb) end if - call self%pl%index(param) + call self%pl%flatten(param) if (.not.param%lrhill_present) call self%pl%set_rhill(self%cb) self%pl%lfirst = param%lfirstkick self%tp%lfirst = param%lfirstkick diff --git a/src/symba/symba_util.f90 b/src/symba/symba_util.f90 index 0d062894f..4b20cf0bf 100644 --- a/src/symba/symba_util.f90 +++ b/src/symba/symba_util.f90 @@ -269,7 +269,7 @@ module subroutine symba_util_fill_tp(self, inserts, lfill_list) end subroutine symba_util_fill_tp - module subroutine symba_util_index_eucl_plpl(self, param) + module subroutine symba_util_flatten_eucl_plpl(self, param) !! author: Jacob R. Elliott and David A. Minton !! !! Turns i,j indices into k index for use in the Euclidean distance matrix. This also sets the lmtiny flag and computes the @@ -283,7 +283,7 @@ module subroutine symba_util_index_eucl_plpl(self, param) implicit none ! Arguments class(symba_pl), intent(inout) :: self !! SyMBA massive body object - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters ! Internals integer(I8B) :: k, nplpl, nplplm integer(I4B) :: i, j, npl, nplm, ip, jp @@ -299,7 +299,7 @@ module subroutine symba_util_index_eucl_plpl(self, param) allocate(self%k_plpl(2, pl%nplpl)) do concurrent (i = 1:npl) do concurrent (j = i+1:npl) - call util_index_eucl_ij_to_k(npl, i, j, k) + call util_flatten_eucl_ij_to_k(npl, i, j, k) self%k_plpl(1, k) = i self%k_plpl(2, k) = j end do @@ -308,7 +308,7 @@ module subroutine symba_util_index_eucl_plpl(self, param) end associate return - end subroutine symba_util_index_eucl_plpl + end subroutine symba_util_flatten_eucl_plpl module subroutine symba_util_peri_pl(self, system, param) @@ -478,7 +478,7 @@ module subroutine symba_util_rearray_pl(self, system, param) ! Reindex the new list of bodies call pl%sort("mass", ascending=.false.) - call pl%index(param) + call pl%flatten(param) ! Reset the kinship trackers call pl%reset_kinship([(i, i=1, npl)]) diff --git a/src/util/util_index.f90 b/src/util/util_flatten.f90 similarity index 78% rename from src/util/util_index.f90 rename to src/util/util_flatten.f90 index d0dd83896..e3e97c1e7 100644 --- a/src/util/util_index.f90 +++ b/src/util/util_flatten.f90 @@ -2,8 +2,8 @@ use swiftest contains - module pure subroutine util_index_eucl_ij_to_k(n, i, j, k) - !$omp declare simd(util_index_eucl_ij_to_k) + module pure subroutine util_flatten_eucl_ij_to_k(n, i, j, k) + !$omp declare simd(util_flatten_eucl_ij_to_k) !! author: Jacob R. Elliott and David A. Minton !! !! Turns i,j indices into k index for use in the Euclidean distance matrix for pl-pl interactions. @@ -27,10 +27,10 @@ module pure subroutine util_index_eucl_ij_to_k(n, i, j, k) k = (i8 - 1_I8B) * n8 - i8 * (i8 - 1_I8B) / 2_I8B + (j8 - i8) return - end subroutine util_index_eucl_ij_to_k + end subroutine util_flatten_eucl_ij_to_k - module pure subroutine util_index_eucl_k_to_ij(n, k, i, j) + module pure subroutine util_flatten_eucl_k_to_ij(n, k, i, j) !! author: Jacob R. Elliott and David A. Minton !! !! Turns k index into i,j indices for use in the Euclidean distance matrix for pl-pl interactions. @@ -59,10 +59,10 @@ module pure subroutine util_index_eucl_k_to_ij(n, k, i, j) j = int(j8, kind=I4B) return - end subroutine util_index_eucl_k_to_ij + end subroutine util_flatten_eucl_k_to_ij - module subroutine util_index_eucl_plpl(self, param) + module subroutine util_flatten_eucl_plpl(self, param) !! author: Jacob R. Elliott and David A. Minton !! !! Turns i,j indices into k index for use in the Euclidean distance matrix for pl-pl interactions for a Swiftest massive body object @@ -74,9 +74,9 @@ module subroutine util_index_eucl_plpl(self, param) implicit none ! Arguments class(swiftest_pl), intent(inout) :: self !! Swiftest massive body object - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters ! Internals - integer(I4B) :: i, j, npl + integer(I4B) :: i, j, npl, err integer(I8B) :: k npl = int(self%nbody, kind=I8B) @@ -84,20 +84,24 @@ module subroutine util_index_eucl_plpl(self, param) nplpl = (npl * (npl - 1) / 2) ! number of entries in a strict lower triangle, npl x npl 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)) - do concurrent (i=1:npl, j=1:npl, j>i) - call util_index_eucl_ij_to_k(npl, i, j, k) - self%k_plpl(1, k) = i - self%k_plpl(2, k) = j - end do + 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(npl, i, j, k) + self%k_plpl(1, k) = i + self%k_plpl(2, k) = j + end do + end if end if end associate return - end subroutine util_index_eucl_plpl + end subroutine util_flatten_eucl_plpl - module subroutine util_index_eucl_pltp(self, pl, param) + module subroutine util_flatten_eucl_pltp(self, pl, param) !! author: Jacob R. Elliott and David A. Minton !! !! Turns i,j indices into k index for use in the Euclidean distance matrix for pl-tp interactions @@ -110,7 +114,7 @@ module subroutine util_index_eucl_pltp(self, pl, param) ! Arguments class(swiftest_tp), intent(inout) :: self !! Swiftest test particle object class(swiftest_pl), intent(in) :: pl !! Swiftest massive body object - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters + class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters ! Internals integer(I8B) :: i, j, counter, npl, ntp @@ -131,6 +135,6 @@ module subroutine util_index_eucl_pltp(self, pl, param) end associate return - end subroutine util_index_eucl_pltp + end subroutine util_flatten_eucl_pltp end submodule s_util_index diff --git a/src/whm/whm_setup.f90 b/src/whm/whm_setup.f90 index 3c5cdeec3..60be140fd 100644 --- a/src/whm/whm_setup.f90 +++ b/src/whm/whm_setup.f90 @@ -82,7 +82,7 @@ module subroutine whm_setup_initialize_system(self, param) ! First we need to make sure that the massive bodies are sorted by heliocentric distance before computing jacobies call util_set_ir3h(self%pl) call self%pl%sort("ir3h", ascending=.false.) - call self%pl%index(param) + call self%pl%flatten(param) ! Make sure that the discard list gets allocated initially call self%tp_discards%setup(0, param)