From bb22d207437b81b8a1364539c1e55c20387e3675 Mon Sep 17 00:00:00 2001 From: "Brian T.N. Gunney" Date: Fri, 17 May 2024 16:52:09 -0700 Subject: [PATCH 01/17] Directly compute object bounding box and remove an unneeded templated method. The object bounding box from the BVH may have special values. We need the true bounding box. After computing the bounding box, gather it right away. --- .../detail/DistributedClosestPointImpl.hpp | 98 ++++++++++++------- 1 file changed, 60 insertions(+), 38 deletions(-) diff --git a/src/axom/quest/detail/DistributedClosestPointImpl.hpp b/src/axom/quest/detail/DistributedClosestPointImpl.hpp index 13d157db6c..187d40176a 100644 --- a/src/axom/quest/detail/DistributedClosestPointImpl.hpp +++ b/src/axom/quest/detail/DistributedClosestPointImpl.hpp @@ -12,6 +12,7 @@ #include "axom/primal.hpp" #include "axom/spin.hpp" #include "axom/core/execution/runtime_policy.hpp" +#include "axom/core/WhereMacro.hpp" #include "axom/fmt.hpp" @@ -639,7 +640,7 @@ class DistributedClosestPointExec : public DistributedClosestPointImpl } // Copy points to internal memory - PointArray coords(ptCount, ptCount); + m_objectPtCoords = PointArray(ptCount, ptCount, m_allocatorID); axom::Array domIds(ptCount, ptCount); std::size_t copiedCount = 0; conduit::Node tmpValues; @@ -673,7 +674,7 @@ class DistributedClosestPointExec : public DistributedClosestPointImpl const int N = internal::extractSize(copySrc); const std::size_t nBytes = sizeof(double) * DIM * N; - axom::copy(coords.data() + copiedCount, + axom::copy(m_objectPtCoords.data() + copiedCount, copySrc.fetch_existing("x").data_ptr(), nBytes); tmpValues.reset(); @@ -682,16 +683,41 @@ class DistributedClosestPointExec : public DistributedClosestPointImpl copiedCount += N; } + + // Compute bounding box + // Don't rely on BVH, whose bounding box is a bit different. +#if defined(AXOM_USE_RAJA) + // Coordinates may be on device but should be compatible with ExecSpace. + axom::ArrayView coordsView = m_objectPtCoords.view(); + PointType minPt, maxPt; + for (int d=0; d minCoord(std::numeric_limits::max()); + RAJA::ReduceMax maxCoord(std::numeric_limits::min()); + RAJA::forall( + RAJA::RangeSegment(0, ptCount), + AXOM_LAMBDA(RAJA::Index_type n) { + for (int d=0; d{m_objectPtCoords.data(), m_objectPtCoords.size()}; +#endif + gatherBoundingBoxes(m_objectBb, m_objectPartitionBbs); + // copy computed data to ExecSpace - m_objectPtCoords = PointArray(coords, m_allocatorID); m_objectPtDomainIds = axom::Array(domIds, m_allocatorID); } bool generateBVHTree() override { - // Delegates to generateBVHTreeImpl<> which uses - // the execution space templated bvh tree - SLIC_ASSERT_MSG(!m_bvh, "BVH tree already initialized"); // In case user changed the allocator after setObjectMesh, @@ -703,18 +729,20 @@ class DistributedClosestPointExec : public DistributedClosestPointImpl } m_bvh = std::make_unique(); - return generateBVHTreeImpl(m_bvh.get()); - } + const int npts = m_objectPtCoords.size(); + axom::Array boxesArray(npts, npts, m_allocatorID); + auto boxesView = boxesArray.view(); + auto pointsView = m_objectPtCoords.view(); - /// Get local copy of all ranks BVH root bounding boxes. - void gatherBVHRoots() - { - SLIC_ASSERT_MSG( - m_bvh, - "BVH tree must be initialized before calling 'gatherBVHRoots"); + axom::for_all( + npts, + AXOM_LAMBDA(axom::IndexType i) { boxesView[i] = BoxType {pointsView[i]}; }); - BoxType local_bb = m_bvh->getBounds(); - gatherBoundingBoxes(local_bb, m_objectPartitionBbs); + // Build bounding volume hierarchy + m_bvh->setAllocatorID(m_allocatorID); + int result = m_bvh->initialize(boxesView, npts); + + return (result == spin::BVH_BUILD_OK); } /// Allgather one bounding box from each rank. @@ -776,6 +804,19 @@ class DistributedClosestPointExec : public DistributedClosestPointImpl * The worst case could incur nranks^2 sends. To avoid excessive * buffer usage, we occasionally check the sends for completion, * using check_send_requests(). + * + * To exclude fruitless searches and communication, we check the + * distance between object partition and query partition, using + * their bounding boxes. If a query partition is too far from an + * object partition, we don't check that specific pair. + * + * TODO: The bounding box for a partition could be excessively big + * if the domains on that partition are spread out far, leading to + * some fruitless communications and checks. Consider having one + * bounding box per domain instead of one per partition. This goes + * for query mesh too. One way simple to implement this is to have + * a xferNode for each query domain and a BVH for each object + * domain. */ void computeClosestPoints(conduit::Node& queryMesh, const std::string& topologyName) const override @@ -940,28 +981,6 @@ class DistributedClosestPointExec : public DistributedClosestPointImpl // Note: following should be private, but nvcc complains about lambdas in private scope public: - /// Templated implementation of generateBVHTree function - bool generateBVHTreeImpl(BVHTreeType* bvh) - { - SLIC_ASSERT(bvh != nullptr); - - const int npts = m_objectPtCoords.size(); - axom::Array boxesArray(npts, npts, m_allocatorID); - auto boxesView = boxesArray.view(); - auto pointsView = m_objectPtCoords.view(); - - axom::for_all( - npts, - AXOM_LAMBDA(axom::IndexType i) { boxesView[i] = BoxType {pointsView[i]}; }); - - // Build bounding volume hierarchy - bvh->setAllocatorID(m_allocatorID); - int result = bvh->initialize(boxesView, npts); - - gatherBVHRoots(); - - return (result == spin::BVH_BUILD_OK); - } void computeLocalClosestPoints(conduit::Node& xferNode) const { @@ -1170,6 +1189,9 @@ class DistributedClosestPointExec : public DistributedClosestPointImpl axom::Array m_objectPtDomainIds; + //!@brief Bounding box for m_objectPtCoords. + BoxType m_objectBb; + /*! @brief Object partition bounding boxes, one per rank. All are in physical space, not index space. */ From 7e80da79969f3eb610c3d86819abf7168787940e Mon Sep 17 00:00:00 2001 From: "Brian T.N. Gunney" Date: Wed, 5 Jun 2024 14:25:47 -0700 Subject: [PATCH 02/17] Rename m_sqDistanceThreshold to m_sqUserDistanceThreshold. --- .../quest/detail/DistributedClosestPointImpl.hpp | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/src/axom/quest/detail/DistributedClosestPointImpl.hpp b/src/axom/quest/detail/DistributedClosestPointImpl.hpp index 187d40176a..9456f199ba 100644 --- a/src/axom/quest/detail/DistributedClosestPointImpl.hpp +++ b/src/axom/quest/detail/DistributedClosestPointImpl.hpp @@ -255,7 +255,7 @@ class DistributedClosestPointImpl , m_mpiComm(MPI_COMM_NULL) , m_rank(-1) , m_nranks(-1) - , m_sqDistanceThreshold(std::numeric_limits::max()) + , m_sqUserDistanceThreshold(std::numeric_limits::max()) { } virtual ~DistributedClosestPointImpl() { } @@ -304,7 +304,7 @@ class DistributedClosestPointImpl { SLIC_ERROR_IF(sqThreshold < 0.0, "Squared distance-threshold must be non-negative."); - m_sqDistanceThreshold = sqThreshold; + m_sqUserDistanceThreshold = sqThreshold; } /*! @@ -551,7 +551,8 @@ class DistributedClosestPointImpl int m_rank; int m_nranks; - double m_sqDistanceThreshold; + //!@brief Distance threshold specified by user. + double m_sqUserDistanceThreshold; bool m_outputRank = true; bool m_outputIndex = true; @@ -855,7 +856,7 @@ class DistributedClosestPointExec : public DistributedClosestPointImpl const auto& otherQueryBb = allQueryBbs[r]; double sqDistance = axom::primal::squared_distance(otherQueryBb, myObjectBb); - if(sqDistance <= m_sqDistanceThreshold) + if(sqDistance <= m_sqUserDistanceThreshold) { ++remainingRecvs; } @@ -971,7 +972,7 @@ class DistributedClosestPointExec : public DistributedClosestPointImpl } double sqDistance = primal::squared_distance(bb, m_objectPartitionBbs[maybeNextRecip]); - if(sqDistance <= m_sqDistanceThreshold) + if(sqDistance <= m_sqUserDistanceThreshold) { return maybeNextRecip; } @@ -1082,7 +1083,7 @@ class DistributedClosestPointExec : public DistributedClosestPointImpl const int rank = m_rank; double* sqDistThresh = axom::allocate(1, m_allocatorID); - *sqDistThresh = m_sqDistanceThreshold; + *sqDistThresh = m_sqUserDistanceThreshold; auto ptCoordsView = m_objectPtCoords.view(); auto ptDomainIdsView = m_objectPtDomainIds.view(); From 71db1c773b02b4f8d30519cd4a4275b1118f05c6 Mon Sep 17 00:00:00 2001 From: "Brian T.N. Gunney" Date: Wed, 5 Jun 2024 14:27:24 -0700 Subject: [PATCH 03/17] Clean up outdated comments. --- src/axom/quest/detail/DistributedClosestPointImpl.hpp | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/src/axom/quest/detail/DistributedClosestPointImpl.hpp b/src/axom/quest/detail/DistributedClosestPointImpl.hpp index 9456f199ba..ffcc5b2725 100644 --- a/src/axom/quest/detail/DistributedClosestPointImpl.hpp +++ b/src/axom/quest/detail/DistributedClosestPointImpl.hpp @@ -685,8 +685,7 @@ class DistributedClosestPointExec : public DistributedClosestPointImpl copiedCount += N; } - // Compute bounding box - // Don't rely on BVH, whose bounding box is a bit different. + // Compute bounding box. #if defined(AXOM_USE_RAJA) // Coordinates may be on device but should be compatible with ExecSpace. axom::ArrayView coordsView = m_objectPtCoords.view(); @@ -753,7 +752,6 @@ class DistributedClosestPointExec : public DistributedClosestPointImpl aabb.getMin().to_array(&sendbuf[0]); aabb.getMax().to_array(&sendbuf[DIM]); axom::Array recvbuf(m_nranks * sendbuf.size()); - // Note: Using axom::Array may reduce clutter a tad. int errf = MPI_Allgather(sendbuf.data(), 2 * DIM, mpi_traits::type, From cee44f7e1c3b5b79456e7aeef937541deb1c4730 Mon Sep 17 00:00:00 2001 From: "Brian T.N. Gunney" Date: Wed, 5 Jun 2024 14:27:49 -0700 Subject: [PATCH 04/17] Compute and gather min-max distance threshold. Not using the computed data yet. --- .../detail/DistributedClosestPointImpl.hpp | 73 +++++++++++++++++++ 1 file changed, 73 insertions(+) diff --git a/src/axom/quest/detail/DistributedClosestPointImpl.hpp b/src/axom/quest/detail/DistributedClosestPointImpl.hpp index ffcc5b2725..2618f075bc 100644 --- a/src/axom/quest/detail/DistributedClosestPointImpl.hpp +++ b/src/axom/quest/detail/DistributedClosestPointImpl.hpp @@ -772,6 +772,22 @@ class DistributedClosestPointExec : public DistributedClosestPointImpl } } + /// Allgather a primitive value. + template + void gatherPrimitiveValue(const T& val, axom::Array& allVals) const + { + allVals.resize(m_nranks); + int errf = MPI_Allgather(&val, + 1, + mpi_traits::type, + allVals.data(), + 1, + mpi_traits::type, + m_mpiComm); + SLIC_ASSERT(errf == MPI_SUCCESS); + AXOM_UNUSED_VAR(errf); + } + /// Compute bounding box for local part of a mesh. BoxType computeMeshBoundingBox(conduit::Node& xferNode) const { @@ -840,6 +856,21 @@ class DistributedClosestPointExec : public DistributedClosestPointImpl BoxArray allQueryBbs; gatherBoundingBoxes(myQueryBb, allQueryBbs); + // Compute the min of the max distance between m_objectBb and each rank's query box. + // TODO: Move this into a function for readability. + double minMaxSqDist = std::numeric_limits::max(); + for( int i = 0; i < m_nranks; ++i ) + { + auto maxSqDist = maxSqDistBetweenBoxes(m_objectBb, allQueryBbs[i]); + minMaxSqDist = std::min(minMaxSqDist, maxSqDist); + } + + const double sqDistanceThreshold = std::min(m_sqUserDistanceThreshold, minMaxSqDist); + xferNodes[m_rank]->fetch("sqDistanceThreshold") = sqDistanceThreshold; + + axom::Array allSqDistanceThreshold(m_nranks); + gatherPrimitiveValue(sqDistanceThreshold, allSqDistanceThreshold); + { conduit::Node& xferNode = *xferNodes[m_rank]; computeLocalClosestPoints(xferNode); @@ -961,6 +992,7 @@ class DistributedClosestPointExec : public DistributedClosestPointImpl int homeRank = xferNode.fetch_existing("homeRank").value(); BoxType bb; get_bounding_box_from_conduit_node(bb, xferNode.fetch_existing("aabb")); + auto sqDistanceThreshold = xferNode.fetch_existing("sqDistanceThreshold").as_double(); for(int i = 1; i < m_nranks; ++i) { int maybeNextRecip = (m_rank + i) % m_nranks; @@ -1197,6 +1229,47 @@ class DistributedClosestPointExec : public DistributedClosestPointImpl BoxArray m_objectPartitionBbs; std::unique_ptr m_bvh; + + //! @brief Compute maximum squared-distance possible between points in 2 boxes. + AXOM_HOST_DEVICE double maxSqDistBetweenBoxes(const BoxType& a, const BoxType& b) const + { + /* + The following logic is necessary should one box nest inside the + other when projected onto one or more axis directions. + + We look at the distance between each corner of box a and the + opposite corner of b. The max distance is the max among those. + Opposite means that if we choose the lower corner in a, we must + compare with the upper corner in b. A vice versa. + */ + double maxSqDist = 0.0; + int numCorners = 1 << DIM; + for (int i = 0; i separation{aCoords, bCoords}; + double sqDist = separation.squared_norm(); + maxSqDist = std::max(maxSqDist, sqDist); + } + + return maxSqDist; + } + }; // DistributedClosestPointExec } // namespace internal From 1d2c444c566af4d2cfcc4201c5bb770df54df61a Mon Sep 17 00:00:00 2001 From: "Brian T.N. Gunney" Date: Wed, 5 Jun 2024 14:40:28 -0700 Subject: [PATCH 05/17] Use the computed distance threshold to filter communication and searches. --- src/axom/quest/detail/DistributedClosestPointImpl.hpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/axom/quest/detail/DistributedClosestPointImpl.hpp b/src/axom/quest/detail/DistributedClosestPointImpl.hpp index 2618f075bc..8ff79f39e1 100644 --- a/src/axom/quest/detail/DistributedClosestPointImpl.hpp +++ b/src/axom/quest/detail/DistributedClosestPointImpl.hpp @@ -885,7 +885,7 @@ class DistributedClosestPointExec : public DistributedClosestPointImpl const auto& otherQueryBb = allQueryBbs[r]; double sqDistance = axom::primal::squared_distance(otherQueryBb, myObjectBb); - if(sqDistance <= m_sqUserDistanceThreshold) + if(sqDistance <= allSqDistanceThreshold[r]) { ++remainingRecvs; } @@ -1002,7 +1002,7 @@ class DistributedClosestPointExec : public DistributedClosestPointImpl } double sqDistance = primal::squared_distance(bb, m_objectPartitionBbs[maybeNextRecip]); - if(sqDistance <= m_sqUserDistanceThreshold) + if(sqDistance <= sqDistanceThreshold) { return maybeNextRecip; } @@ -1113,7 +1113,7 @@ class DistributedClosestPointExec : public DistributedClosestPointImpl const int rank = m_rank; double* sqDistThresh = axom::allocate(1, m_allocatorID); - *sqDistThresh = m_sqUserDistanceThreshold; + *sqDistThresh = xferNode.fetch("sqDistanceThreshold").as_double(); auto ptCoordsView = m_objectPtCoords.view(); auto ptDomainIdsView = m_objectPtDomainIds.view(); From 3049520bd5ee42a9cd4dec5ee787d7c85b859564 Mon Sep 17 00:00:00 2001 From: "Brian T.N. Gunney" Date: Wed, 5 Jun 2024 17:36:31 -0700 Subject: [PATCH 06/17] Fix bugs in new code. --- src/axom/quest/detail/DistributedClosestPointImpl.hpp | 11 ++++------- 1 file changed, 4 insertions(+), 7 deletions(-) diff --git a/src/axom/quest/detail/DistributedClosestPointImpl.hpp b/src/axom/quest/detail/DistributedClosestPointImpl.hpp index 8ff79f39e1..7bc003c90f 100644 --- a/src/axom/quest/detail/DistributedClosestPointImpl.hpp +++ b/src/axom/quest/detail/DistributedClosestPointImpl.hpp @@ -697,11 +697,8 @@ class DistributedClosestPointExec : public DistributedClosestPointImpl RAJA::forall( RAJA::RangeSegment(0, ptCount), AXOM_LAMBDA(RAJA::Index_type n) { - for (int d=0; d::max(); for( int i = 0; i < m_nranks; ++i ) { - auto maxSqDist = maxSqDistBetweenBoxes(m_objectBb, allQueryBbs[i]); + auto maxSqDist = maxSqDistBetweenBoxes(myQueryBb, m_objectPartitionBbs[i]); minMaxSqDist = std::min(minMaxSqDist, maxSqDist); } From d9b09d720ea18ee46a468c95c6185a7c094c2fc5 Mon Sep 17 00:00:00 2001 From: "Brian T.N. Gunney" Date: Thu, 6 Jun 2024 16:34:43 -0700 Subject: [PATCH 07/17] More robust handling of non-unique strides in ArrayBase. --- src/axom/core/ArrayBase.hpp | 19 +++++++++++++++---- src/axom/core/MDMapping.hpp | 11 +++++++---- 2 files changed, 22 insertions(+), 8 deletions(-) diff --git a/src/axom/core/ArrayBase.hpp b/src/axom/core/ArrayBase.hpp index a24d910af6..4bab2d8fb7 100644 --- a/src/axom/core/ArrayBase.hpp +++ b/src/axom/core/ArrayBase.hpp @@ -387,15 +387,26 @@ class ArrayBase updateStrides(); } - /// \brief Set the shape and stride - AXOM_HOST_DEVICE void setShapeAndStride(const StackArray& shape, - const StackArray& stride) + /*! + \brief Set the shape and stride + + \param [in] shape + \param [in] stride + \param [in] orderPref Preference for resolving non-unique strides. + \a ROW means to advance left index first. \a COLUMN means to advance + right index first. Indexing is correct regardless of the + preference, but the index ordering is dependent on the choice. + */ + AXOM_HOST_DEVICE void setShapeAndStride( + const StackArray& shape, + const StackArray& stride, + axom::ArrayStrideOrder orderPref = axom::ArrayStrideOrder::ROW) { #ifdef AXOM_DEBUG validateShapeAndStride(shape, stride); #endif m_shape = shape; - m_mapping.initializeStrides(stride); + m_mapping.initializeStrides(stride, orderPref); m_minStride = m_mapping.fastestStrideLength(); } diff --git a/src/axom/core/MDMapping.hpp b/src/axom/core/MDMapping.hpp index 74995c8b83..19b303189d 100644 --- a/src/axom/core/MDMapping.hpp +++ b/src/axom/core/MDMapping.hpp @@ -232,10 +232,13 @@ class MDMapping os << strides[d] << ","; } os << strides[DIM - 1] << ")"; - std::cerr << "ERROR: MDMapping: Non-unique strides " << os.str() << ".\n" - << "Likely, multi-dim array shape is 1 in some direction.\n" - << "Impossible to compute index ordering.\n" - << "Please use a different MDMapping initializer.\n"; + std::cerr + << "ERROR: MDMapping: Non-unique strides " << os.str() << ".\n" + << "Caused by multi-dim array shape of 1 in some direction.\n" + << "It is impossible to compute index ordering.\n" + << "Use initializeStrides() with the stride order preference to fix.\n" + << "The resulting slowestDirs() depends on the preference\n" + << "but the array mapping will still be correct.\n"; #endif utilities::processAbort(); } From 81655eb65baef747c85f47df2a2c8a68a4cb303c Mon Sep 17 00:00:00 2001 From: "Brian T.N. Gunney" Date: Thu, 6 Jun 2024 16:35:51 -0700 Subject: [PATCH 08/17] Robust handling of insufficient size specified when creating data fields in MeshViewUtil. --- src/axom/quest/MeshViewUtil.hpp | 15 ++++++++++----- 1 file changed, 10 insertions(+), 5 deletions(-) diff --git a/src/axom/quest/MeshViewUtil.hpp b/src/axom/quest/MeshViewUtil.hpp index 7c3098c6e5..75be596aad 100644 --- a/src/axom/quest/MeshViewUtil.hpp +++ b/src/axom/quest/MeshViewUtil.hpp @@ -124,7 +124,9 @@ static void shapesToStridesAndOffsets( @param offsets [i] Blueprint-style index offsets. @param strides [i] Blueprint-style strides. @param valuesCount [i] Number of values in - ghost-padded data. + ghost-padded data. If this is too small (results in + a negative \a hiPads for the slowest stride + direction), that padding will be bumped to zero. @param paddedShape [o] \a realShape + \a loPads + \a hiPads @param loPads [o] Ghost padding amount on low side. @param hiPads [o] Ghost padding amount ont high side. @@ -165,7 +167,9 @@ static void stridesAndOffsetsToShapes(const axom::StackArray& realSh const int& nextDir = strideOrder[nd + 1]; paddedShape[curDir] = strides[nextDir] / strides[curDir]; } - paddedShape[strideOrder[DIM - 1]] = valuesCount / strides[strideOrder[DIM - 1]]; + const int slowestDir = strideOrder[DIM - 1]; + paddedShape[slowestDir] = std::max(valuesCount / strides[slowestDir], + realShape[slowestDir] + offsets[slowestDir]); for(int d = 0; d < DIM; ++d) { @@ -500,7 +504,7 @@ class MeshViewUtil /*! @brief Return view to a scalar field variable. - WARNING: The view returned has an allocator id determined by + WARNING: The view returned claims an allocator id determined by \a MemSpace, regardless of the memory type. WARNING: Assuming, without checking, that the field contains @@ -647,8 +651,9 @@ class MeshViewUtil @param [in] fieldName @param [in] association "vertex" or "element" - @param [in] dtype Conduit data type to put in the field. Must be at least - big enough for the strides and offsets specified. + @param [in] dtype Conduit data type to put in the field. If this is not + big enough for the strides and offsets specified, a minimally + sufficient size will be used to avoid negative ghost padding. @param [in] strides Data strides. Set to zero for no ghosts and default strides. @param [in] offsets Data index offsets. Set to zero for no ghosts. From 5ae1ab38d0ea44d803560c92fb0ceb9542f651df Mon Sep 17 00:00:00 2001 From: "Brian T.N. Gunney" Date: Thu, 6 Jun 2024 16:37:24 -0700 Subject: [PATCH 09/17] Add cell ranks to mesh output, for debugging. --- ...est_distributed_distance_query_example.cpp | 44 +++++++++++++++++++ 1 file changed, 44 insertions(+) diff --git a/src/axom/quest/examples/quest_distributed_distance_query_example.cpp b/src/axom/quest/examples/quest_distributed_distance_query_example.cpp index 15e798e6b3..77cc2b19d4 100644 --- a/src/axom/quest/examples/quest_distributed_distance_query_example.cpp +++ b/src/axom/quest/examples/quest_distributed_distance_query_example.cpp @@ -339,6 +339,9 @@ struct BlueprintParticleMesh MPI_Allreduce(MPI_IN_PLACE, &m_dimension, 1, MPI_INT, MPI_MAX, MPI_COMM_WORLD); SLIC_ASSERT(m_dimension > 0); + // For debugging, create cell-centered owner-rank field. + m_dimension == 2 ? setCellOwnerRank<2>(mdMesh) : setCellOwnerRank<3>(mdMesh); + if(domCount > 0) { // Put mdMesh into sidre Group. @@ -429,6 +432,43 @@ struct BlueprintParticleMesh } } + /*! + @brief Set the cell-centered owner rank field. + + @param mdMesh A multi-domain structured mesh. + */ + template + void setCellOwnerRank(conduit::Node& mdMesh) + { + std::string fieldName = "owner_rank"; + axom::StackArray zeroPads; + axom::StackArray strideOrder; + for(int d = 0; d < DIM; ++d) + { + zeroPads[d] = 0; + strideOrder[d] = DIM - d - 1; + } + for(conduit::Node& dom : mdMesh.children()) + { + axom::quest::MeshViewUtil domainView(dom, m_topologyName); + domainView.createField(fieldName, + "element", + conduit::DataType::uint32(), + zeroPads, + zeroPads, + strideOrder); + + auto ownerRankView = + domainView.template getFieldView(fieldName); + axom::detail::ArrayOps::fill( + ownerRankView.data(), + 0, + ownerRankView.size(), + ownerRankView.getAllocatorID(), + m_rank); + } + } + template axom::Array> getPoints(int domainIdx) { @@ -1501,12 +1541,16 @@ int main(int argc, char** argv) maxQuery)); } slic::flushStreams(); + SLIC_INFO(axom::fmt::format("Updating closest points.")); + slic::flushStreams(); queryMeshWrapper.update_closest_points(queryMeshNode); int errCount = 0; int localErrCount = 0; if(params.checkResults) { + SLIC_INFO(axom::fmt::format("Checking results.")); + slic::flushStreams(); if(spatialDim == 2) { primal::Point center(params.circleCenter.data()); From 75f94d221e04af1d636676e3b1f5283d0b440a6c Mon Sep 17 00:00:00 2001 From: "Brian T.N. Gunney" Date: Thu, 6 Jun 2024 16:39:49 -0700 Subject: [PATCH 10/17] Count number of times xferNode is locally searched, for analysis and debugging. --- src/axom/quest/detail/DistributedClosestPointImpl.hpp | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/axom/quest/detail/DistributedClosestPointImpl.hpp b/src/axom/quest/detail/DistributedClosestPointImpl.hpp index 7bc003c90f..1b5cee6475 100644 --- a/src/axom/quest/detail/DistributedClosestPointImpl.hpp +++ b/src/axom/quest/detail/DistributedClosestPointImpl.hpp @@ -871,6 +871,7 @@ class DistributedClosestPointExec : public DistributedClosestPointImpl { conduit::Node& xferNode = *xferNodes[m_rank]; computeLocalClosestPoints(xferNode); + xferNode["localSearchCount"].set_int32(1); } const auto& myObjectBb = m_objectPartitionBbs[m_rank]; @@ -951,6 +952,8 @@ class DistributedClosestPointExec : public DistributedClosestPointImpl else { computeLocalClosestPoints(xferNode); + auto tmpCount = xferNode["localSearchCount"].as_int32(); + xferNode["localSearchCount"].set_int32(1 + tmpCount); isendRequests.emplace_back(conduit::relay::mpi::Request()); auto& isendRequest = isendRequests.back(); From 869d456ae7e93e6a14e44adacf3213045128df05 Mon Sep 17 00:00:00 2001 From: "Brian T.N. Gunney" Date: Fri, 7 Jun 2024 08:03:01 -0700 Subject: [PATCH 11/17] Reformat. --- .../detail/DistributedClosestPointImpl.hpp | 40 +++++++++++-------- 1 file changed, 23 insertions(+), 17 deletions(-) diff --git a/src/axom/quest/detail/DistributedClosestPointImpl.hpp b/src/axom/quest/detail/DistributedClosestPointImpl.hpp index 1b5cee6475..7db4e0df61 100644 --- a/src/axom/quest/detail/DistributedClosestPointImpl.hpp +++ b/src/axom/quest/detail/DistributedClosestPointImpl.hpp @@ -690,10 +690,12 @@ class DistributedClosestPointExec : public DistributedClosestPointImpl // Coordinates may be on device but should be compatible with ExecSpace. axom::ArrayView coordsView = m_objectPtCoords.view(); PointType minPt, maxPt; - for (int d=0; d minCoord(std::numeric_limits::max()); - RAJA::ReduceMax maxCoord(std::numeric_limits::min()); + RAJA::ReduceMin minCoord( + std::numeric_limits::max()); + RAJA::ReduceMax maxCoord( + std::numeric_limits::min()); RAJA::forall( RAJA::RangeSegment(0, ptCount), AXOM_LAMBDA(RAJA::Index_type n) { @@ -705,7 +707,9 @@ class DistributedClosestPointExec : public DistributedClosestPointImpl } m_objectBb = BoxType(minPt, maxPt); #else - m_objectBb = axom::primal::BoundingBox{m_objectPtCoords.data(), m_objectPtCoords.size()}; + m_objectBb = + axom::primal::BoundingBox {m_objectPtCoords.data(), + m_objectPtCoords.size()}; #endif gatherBoundingBoxes(m_objectBb, m_objectPartitionBbs); @@ -770,7 +774,7 @@ class DistributedClosestPointExec : public DistributedClosestPointImpl } /// Allgather a primitive value. - template + template void gatherPrimitiveValue(const T& val, axom::Array& allVals) const { allVals.resize(m_nranks); @@ -856,13 +860,14 @@ class DistributedClosestPointExec : public DistributedClosestPointImpl // Compute the min of the max distance between myQueryBb and each rank's object bounding box. // TODO: Move this into a function for readability. double minMaxSqDist = std::numeric_limits::max(); - for( int i = 0; i < m_nranks; ++i ) + for(int i = 0; i < m_nranks; ++i) { auto maxSqDist = maxSqDistBetweenBoxes(myQueryBb, m_objectPartitionBbs[i]); minMaxSqDist = std::min(minMaxSqDist, maxSqDist); } - const double sqDistanceThreshold = std::min(m_sqUserDistanceThreshold, minMaxSqDist); + const double sqDistanceThreshold = + std::min(m_sqUserDistanceThreshold, minMaxSqDist); xferNodes[m_rank]->fetch("sqDistanceThreshold") = sqDistanceThreshold; axom::Array allSqDistanceThreshold(m_nranks); @@ -992,7 +997,8 @@ class DistributedClosestPointExec : public DistributedClosestPointImpl int homeRank = xferNode.fetch_existing("homeRank").value(); BoxType bb; get_bounding_box_from_conduit_node(bb, xferNode.fetch_existing("aabb")); - auto sqDistanceThreshold = xferNode.fetch_existing("sqDistanceThreshold").as_double(); + auto sqDistanceThreshold = + xferNode.fetch_existing("sqDistanceThreshold").as_double(); for(int i = 1; i < m_nranks; ++i) { int maybeNextRecip = (m_rank + i) % m_nranks; @@ -1012,7 +1018,6 @@ class DistributedClosestPointExec : public DistributedClosestPointImpl // Note: following should be private, but nvcc complains about lambdas in private scope public: - void computeLocalClosestPoints(conduit::Node& xferNode) const { using axom::primal::squared_distance; @@ -1231,7 +1236,8 @@ class DistributedClosestPointExec : public DistributedClosestPointImpl std::unique_ptr m_bvh; //! @brief Compute maximum squared-distance possible between points in 2 boxes. - AXOM_HOST_DEVICE double maxSqDistBetweenBoxes(const BoxType& a, const BoxType& b) const + AXOM_HOST_DEVICE double maxSqDistBetweenBoxes(const BoxType& a, + const BoxType& b) const { /* The following logic is necessary should one box nest inside the @@ -1244,25 +1250,25 @@ class DistributedClosestPointExec : public DistributedClosestPointImpl */ double maxSqDist = 0.0; int numCorners = 1 << DIM; - for (int i = 0; i separation{aCoords, bCoords}; + primal::Vector separation {aCoords, bCoords}; double sqDist = separation.squared_norm(); maxSqDist = std::max(maxSqDist, sqDist); } From 43f5b523c08dd6cd044d6c9577ebc8f79aa74d46 Mon Sep 17 00:00:00 2001 From: "Brian T.N. Gunney" Date: Fri, 7 Jun 2024 09:52:30 -0700 Subject: [PATCH 12/17] Special treatment for domain underloading, which causes invalid bounding boxes. --- .../detail/DistributedClosestPointImpl.hpp | 39 +++++++++++++------ 1 file changed, 28 insertions(+), 11 deletions(-) diff --git a/src/axom/quest/detail/DistributedClosestPointImpl.hpp b/src/axom/quest/detail/DistributedClosestPointImpl.hpp index 7db4e0df61..33872e788f 100644 --- a/src/axom/quest/detail/DistributedClosestPointImpl.hpp +++ b/src/axom/quest/detail/DistributedClosestPointImpl.hpp @@ -879,18 +879,28 @@ class DistributedClosestPointExec : public DistributedClosestPointImpl xferNode["localSearchCount"].set_int32(1); } + /* + Count number of remote query partitions to receive, + based on number of those partitions close enough to myObjectBb. + */ const auto& myObjectBb = m_objectPartitionBbs[m_rank]; int remainingRecvs = 0; - for(int r = 0; r < m_nranks; ++r) + if(myObjectBb.isValid()) { - if(r != m_rank) + for(int r = 0; r < m_nranks; ++r) { - const auto& otherQueryBb = allQueryBbs[r]; - double sqDistance = - axom::primal::squared_distance(otherQueryBb, myObjectBb); - if(sqDistance <= allSqDistanceThreshold[r]) + if(r != m_rank) { - ++remainingRecvs; + const auto& otherQueryBb = allQueryBbs[r]; + if(otherQueryBb.isValid()) + { + double sqDistance = + axom::primal::squared_distance(otherQueryBb, myObjectBb); + if(sqDistance <= allSqDistanceThreshold[r]) + { + ++remainingRecvs; + } + } } } } @@ -1235,10 +1245,19 @@ class DistributedClosestPointExec : public DistributedClosestPointImpl std::unique_ptr m_bvh; - //! @brief Compute maximum squared-distance possible between points in 2 boxes. + /*! + @brief Compute maximum squared-distance possible between points in 2 boxes, + or std::numeric_limits::max() if either box is invalid. + */ AXOM_HOST_DEVICE double maxSqDistBetweenBoxes(const BoxType& a, const BoxType& b) const { + if(!a.isValid() || !b.isValid()) + { + return std::numeric_limits::max(); + } + + double maxSqDist = 0.0; /* The following logic is necessary should one box nest inside the other when projected onto one or more axis directions. @@ -1246,9 +1265,8 @@ class DistributedClosestPointExec : public DistributedClosestPointImpl We look at the distance between each corner of box a and the opposite corner of b. The max distance is the max among those. Opposite means that if we choose the lower corner in a, we must - compare with the upper corner in b. A vice versa. + compare with the upper corner in b. And vice versa. */ - double maxSqDist = 0.0; int numCorners = 1 << DIM; for(int i = 0; i < numCorners; ++i) { @@ -1272,7 +1290,6 @@ class DistributedClosestPointExec : public DistributedClosestPointImpl double sqDist = separation.squared_norm(); maxSqDist = std::max(maxSqDist, sqDist); } - return maxSqDist; } From 0a28c48a0e00de4bcbd4080e223fc4227c51e3ce Mon Sep 17 00:00:00 2001 From: "Brian T.N. Gunney" Date: Fri, 7 Jun 2024 11:05:47 -0700 Subject: [PATCH 13/17] Comment and minor changes. --- src/axom/quest/detail/DistributedClosestPointImpl.hpp | 7 +++---- .../examples/quest_distributed_distance_query_example.cpp | 1 + 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/src/axom/quest/detail/DistributedClosestPointImpl.hpp b/src/axom/quest/detail/DistributedClosestPointImpl.hpp index 33872e788f..4e4aabeb98 100644 --- a/src/axom/quest/detail/DistributedClosestPointImpl.hpp +++ b/src/axom/quest/detail/DistributedClosestPointImpl.hpp @@ -12,7 +12,6 @@ #include "axom/primal.hpp" #include "axom/spin.hpp" #include "axom/core/execution/runtime_policy.hpp" -#include "axom/core/WhereMacro.hpp" #include "axom/fmt.hpp" @@ -858,7 +857,6 @@ class DistributedClosestPointExec : public DistributedClosestPointImpl gatherBoundingBoxes(myQueryBb, allQueryBbs); // Compute the min of the max distance between myQueryBb and each rank's object bounding box. - // TODO: Move this into a function for readability. double minMaxSqDist = std::numeric_limits::max(); for(int i = 0; i < m_nranks; ++i) { @@ -881,8 +879,8 @@ class DistributedClosestPointExec : public DistributedClosestPointImpl /* Count number of remote query partitions to receive, - based on number of those partitions close enough to myObjectBb. - */ + i.e., how many partitions close enough to myObjectBb. + */ const auto& myObjectBb = m_objectPartitionBbs[m_rank]; int remainingRecvs = 0; if(myObjectBb.isValid()) @@ -963,6 +961,7 @@ class DistributedClosestPointExec : public DistributedClosestPointImpl if(homeRank == m_rank) { node_copy_xfer_to_query(xferNode, queryMesh, topologyName); + xferNodes.erase(m_rank); } else { diff --git a/src/axom/quest/examples/quest_distributed_distance_query_example.cpp b/src/axom/quest/examples/quest_distributed_distance_query_example.cpp index 77cc2b19d4..0874c1f9bd 100644 --- a/src/axom/quest/examples/quest_distributed_distance_query_example.cpp +++ b/src/axom/quest/examples/quest_distributed_distance_query_example.cpp @@ -1101,6 +1101,7 @@ class QueryMeshWrapper /** * Generates points on a sphere, partitioned into multiple domains. + * The sphere's polar axis is in the z-direction. * Point spacing in the longitudinal direction can be random (default) or uniform. * 3D points cover the given latitude range. */ From 37466d4d2dceb1f72cf40862979b8492e90cfef3 Mon Sep 17 00:00:00 2001 From: "Brian T.N. Gunney" Date: Fri, 7 Jun 2024 12:41:44 -0700 Subject: [PATCH 14/17] Compute analytical distances for non-spherical object mesh. If object "sphere" mesh doesn't extend to the poles, distance assuming a sphere is incorrect. This change gets it correct when the object mesh doesn't have the full latitude range. --- ...est_distributed_distance_query_example.cpp | 33 ++++++++++++++++++- 1 file changed, 32 insertions(+), 1 deletion(-) diff --git a/src/axom/quest/examples/quest_distributed_distance_query_example.cpp b/src/axom/quest/examples/quest_distributed_distance_query_example.cpp index 0874c1f9bd..7c2a91d132 100644 --- a/src/axom/quest/examples/quest_distributed_distance_query_example.cpp +++ b/src/axom/quest/examples/quest_distributed_distance_query_example.cpp @@ -1006,13 +1006,44 @@ class QueryMeshWrapper const double allowableSlack = avgObjectRes / 2; using IndexSet = slam::PositionSet<>; + const PointType circleCenter {params.circleCenter.data(), DIM}; + double zNorth = + params.circleRadius * std::sin(params.latRange[1] * M_PI / 180); + double xyNorth = + params.circleRadius * std::cos(params.latRange[1] * M_PI / 180); + double zSouth = + params.circleRadius * std::sin(params.latRange[0] * M_PI / 180); + double xySouth = + params.circleRadius * std::cos(params.latRange[0] * M_PI / 180); for(auto i : IndexSet(queryPts.size())) { bool errf = false; + // Compute the analytical distance to sphere (or partial sphere + // if the latitude range doesn't go all the way to the poles). const auto& qPt = queryPts[i]; const auto& cpCoord = cpCoords[i]; double analyticalDist = std::fabs(sphere.computeSignedDistance(qPt)); + if(DIM == 3 && params.latRange[0] > -90 && params.latRange[1] < 90) + { + // More complicated analytical distance for partial-sphere object. + axom::primal::Vector cToQ {circleCenter, qPt}; + double z = cToQ[2]; + cToQ[2] = 0.0; + double xy = cToQ.norm(); + double qPtLat = std::atan(z / xy) * 180 / M_PI; + if(qPtLat > params.latRange[1]) + { + analyticalDist = std::sqrt((z - zNorth) * (z - zNorth) + + (xy - xyNorth) * (xy - xyNorth)); + } + else if(qPtLat < params.latRange[0]) + { + analyticalDist = std::sqrt((z - zSouth) * (z - zSouth) + + (xy - xySouth) * (xy - xySouth)); + } + } + const bool closestPointFound = (cpIndices[i] == -1); if(closestPointFound) { @@ -1072,7 +1103,7 @@ class QueryMeshWrapper { errf = true; SLIC_INFO( - axom::fmt::format("***Warning: Closest distance for {} (index " + axom::fmt::format("***Error: Closest distance for {} (index " "{}, cp {}) is {}, off by {}.", qPt, i, From d6b9dccdee2a6042a4032ea0a24a8e32fa4255e9 Mon Sep 17 00:00:00 2001 From: "Brian T.N. Gunney" Date: Fri, 5 Jul 2024 16:23:07 -0700 Subject: [PATCH 15/17] Fix bug in computing object partitions bounding box. --- src/axom/quest/detail/DistributedClosestPointImpl.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/axom/quest/detail/DistributedClosestPointImpl.hpp b/src/axom/quest/detail/DistributedClosestPointImpl.hpp index 4e4aabeb98..0de4122f32 100644 --- a/src/axom/quest/detail/DistributedClosestPointImpl.hpp +++ b/src/axom/quest/detail/DistributedClosestPointImpl.hpp @@ -694,7 +694,7 @@ class DistributedClosestPointExec : public DistributedClosestPointImpl RAJA::ReduceMin minCoord( std::numeric_limits::max()); RAJA::ReduceMax maxCoord( - std::numeric_limits::min()); + -std::numeric_limits::max()); RAJA::forall( RAJA::RangeSegment(0, ptCount), AXOM_LAMBDA(RAJA::Index_type n) { From 662f2ccc3e8b876716a8779e6ddecad470f47011 Mon Sep 17 00:00:00 2001 From: "Brian T.N. Gunney" Date: Fri, 5 Jul 2024 16:45:20 -0700 Subject: [PATCH 16/17] Add performance diagnostics to the example. --- src/axom/quest/DistributedClosestPoint.cpp | 10 +++++ src/axom/quest/DistributedClosestPoint.hpp | 18 ++++++++ .../detail/DistributedClosestPointImpl.hpp | 41 +++++++++++++++++-- ...est_distributed_distance_query_example.cpp | 22 ++++++++++ 4 files changed, 88 insertions(+), 3 deletions(-) diff --git a/src/axom/quest/DistributedClosestPoint.cpp b/src/axom/quest/DistributedClosestPoint.cpp index 9bf6063a19..63553d8c2b 100644 --- a/src/axom/quest/DistributedClosestPoint.cpp +++ b/src/axom/quest/DistributedClosestPoint.cpp @@ -320,5 +320,15 @@ void DistributedClosestPoint::verifyTopologyName(const conduit::Node& meshNode, } } +axom::IndexType DistributedClosestPoint::searchCount() const +{ + return m_impl->searchCount(); +} + +double DistributedClosestPoint::effectiveDistanceThreshold() const +{ + return m_impl->effectiveDistanceThreshold(); +} + } // end namespace quest } // end namespace axom diff --git a/src/axom/quest/DistributedClosestPoint.hpp b/src/axom/quest/DistributedClosestPoint.hpp index dca36f580b..ca91f1b000 100644 --- a/src/axom/quest/DistributedClosestPoint.hpp +++ b/src/axom/quest/DistributedClosestPoint.hpp @@ -8,6 +8,7 @@ #include "axom/config.hpp" #include "axom/core/execution/runtime_policy.hpp" +#include "axom/core/Types.hpp" #include "conduit_node.hpp" @@ -173,6 +174,23 @@ class DistributedClosestPoint void computeClosestPoints(conduit::Node& query_node, const std::string& topology); + /*! + @brief Return the number of searches done on the last query + mesh's local partition. + + This count includes 1 by the owner rank plus however many remote + ranks searched the partition. + */ + axom::IndexType searchCount() const; + /*! + @brief Return the effective distance threshold used for the + last query mesh's local partition. + + Due to optimizations, this may be smaller than the value set + in setDistanceThreshold(). + */ + double effectiveDistanceThreshold() const; + private: /*! @brief Allocate the DistributedClosestPointImpl object, which actually does the work. diff --git a/src/axom/quest/detail/DistributedClosestPointImpl.hpp b/src/axom/quest/detail/DistributedClosestPointImpl.hpp index 0de4122f32..401bfcc1f3 100644 --- a/src/axom/quest/detail/DistributedClosestPointImpl.hpp +++ b/src/axom/quest/detail/DistributedClosestPointImpl.hpp @@ -542,6 +542,23 @@ class DistributedClosestPointImpl virtual void computeClosestPoints(conduit::Node& queryMesh, const std::string& topologyName) const = 0; + /*! + @brief Return the number of searches done on the last query + mesh's local partition. + */ + axom::IndexType searchCount() const + { + return m_searchCount; + } + /*! + @brief Return the effective distance threshold of the last query + partition. + */ + double effectiveDistanceThreshold() const + { + return m_effectiveDistanceThreshold; + } + protected: int m_allocatorID; bool m_isVerbose; @@ -570,6 +587,14 @@ class DistributedClosestPointImpl /// MPI rank of closest element int rank {-1}; }; + + //@{ + //!@name Diagnostic data + //!@brief Number of searches conducted on the last query node. + mutable axom::IndexType m_searchCount = -1; + //!@brief Effective distance threshold of the last query partition. + mutable double m_effectiveDistanceThreshold = 0.0; + //@} }; /*! @@ -856,6 +881,13 @@ class DistributedClosestPointExec : public DistributedClosestPointImpl BoxArray allQueryBbs; gatherBoundingBoxes(myQueryBb, allQueryBbs); + /* + Note: The two m_nranks loops below can be moved to device for the cost of + copying the bounding box arrays to device. Not sure if it's worthwhile, + because in general, we assume that m_nranks is small relative to the + number of points in the partitions. + */ + // Compute the min of the max distance between myQueryBb and each rank's object bounding box. double minMaxSqDist = std::numeric_limits::max(); for(int i = 0; i < m_nranks; ++i) @@ -867,6 +899,7 @@ class DistributedClosestPointExec : public DistributedClosestPointImpl const double sqDistanceThreshold = std::min(m_sqUserDistanceThreshold, minMaxSqDist); xferNodes[m_rank]->fetch("sqDistanceThreshold") = sqDistanceThreshold; + m_effectiveDistanceThreshold = sqrt(sqDistanceThreshold); axom::Array allSqDistanceThreshold(m_nranks); gatherPrimitiveValue(sqDistanceThreshold, allSqDistanceThreshold); @@ -874,7 +907,8 @@ class DistributedClosestPointExec : public DistributedClosestPointImpl { conduit::Node& xferNode = *xferNodes[m_rank]; computeLocalClosestPoints(xferNode); - xferNode["localSearchCount"].set_int32(1); + xferNode["searchCount"].set_int32(1); + m_searchCount = 1; } /* @@ -961,13 +995,14 @@ class DistributedClosestPointExec : public DistributedClosestPointImpl if(homeRank == m_rank) { node_copy_xfer_to_query(xferNode, queryMesh, topologyName); + m_searchCount = xferNode["searchCount"].value(); xferNodes.erase(m_rank); } else { computeLocalClosestPoints(xferNode); - auto tmpCount = xferNode["localSearchCount"].as_int32(); - xferNode["localSearchCount"].set_int32(1 + tmpCount); + auto tmpCount = xferNode["searchCount"].as_int32(); + xferNode["searchCount"].set_int32(1 + tmpCount); isendRequests.emplace_back(conduit::relay::mpi::Request()); auto& isendRequest = isendRequests.back(); diff --git a/src/axom/quest/examples/quest_distributed_distance_query_example.cpp b/src/axom/quest/examples/quest_distributed_distance_query_example.cpp index 7c2a91d132..5e343e0f9b 100644 --- a/src/axom/quest/examples/quest_distributed_distance_query_example.cpp +++ b/src/axom/quest/examples/quest_distributed_distance_query_example.cpp @@ -1550,6 +1550,12 @@ int main(int argc, char** argv) MPI_Allreduce(&inVal, &maxVal, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD); MPI_Allreduce(&inVal, &sumVal, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); }; + auto getIntMinMax = + [](int inVal, int& minVal, int& maxVal, int& sumVal) { + MPI_Allreduce(&inVal, &minVal, 1, MPI_INT, MPI_MIN, MPI_COMM_WORLD); + MPI_Allreduce(&inVal, &maxVal, 1, MPI_INT, MPI_MAX, MPI_COMM_WORLD); + MPI_Allreduce(&inVal, &sumVal, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD); + }; // Output some timing stats { @@ -1559,6 +1565,12 @@ int main(int argc, char** argv) double minQuery, maxQuery, sumQuery; getDoubleMinMax(queryTimer.elapsedTimeInSec(), minQuery, maxQuery, sumQuery); + double minFilterDist, maxFilterDist, sumFilterDist; + getDoubleMinMax(query.effectiveDistanceThreshold(), minFilterDist, maxFilterDist, sumFilterDist); + + int minSearchCount, maxSearchCount, sumSearchCount; + getIntMinMax(int(query.searchCount()), minSearchCount, maxSearchCount, sumSearchCount); + SLIC_INFO(axom::fmt::format( "Initialization with policy {} took {{avg:{}, min:{}, max:{}}} seconds", axom::runtime_policy::s_policyToName.at(params.policy), @@ -1571,6 +1583,16 @@ int main(int argc, char** argv) sumQuery / num_ranks, minQuery, maxQuery)); + SLIC_INFO(axom::fmt::format( + "Search counts {{avg:{}, min:{}, max:{}}}", + sumSearchCount / num_ranks, + minSearchCount, + maxSearchCount)); + SLIC_INFO(axom::fmt::format( + "Effective distance threshold {{avg:{}, min:{}, max:{}}}", + sumFilterDist / num_ranks, + minFilterDist, + maxFilterDist)); } slic::flushStreams(); SLIC_INFO(axom::fmt::format("Updating closest points.")); From 30a4b7293bd6ef818f9468a675fcc885342f5bc4 Mon Sep 17 00:00:00 2001 From: "Brian T.N. Gunney" Date: Wed, 10 Jul 2024 07:20:43 -0700 Subject: [PATCH 17/17] Improve output formatting and fix an error in computing tolerance in 3D. --- .../quest_distributed_distance_query_example.cpp | 14 +++++--------- 1 file changed, 5 insertions(+), 9 deletions(-) diff --git a/src/axom/quest/examples/quest_distributed_distance_query_example.cpp b/src/axom/quest/examples/quest_distributed_distance_query_example.cpp index 5e343e0f9b..9805bd52e0 100644 --- a/src/axom/quest/examples/quest_distributed_distance_query_example.cpp +++ b/src/axom/quest/examples/quest_distributed_distance_query_example.cpp @@ -983,11 +983,6 @@ class QueryMeshWrapper SLIC_ASSERT(queryPts.size() == cpCoords.size()); SLIC_ASSERT(queryPts.size() == cpIndices.size()); - if(params.isVerbose()) - { - SLIC_INFO(axom::fmt::format("Closest points ({}):", cpCoords.size())); - } - /* Allowable slack is half the arclength between 2 adjacent object points. A query point on the object can correctly have that @@ -999,7 +994,7 @@ class QueryMeshWrapper const double longSpacing = 2 * M_PI * params.circleRadius / params.longPointCount; const double latSpacing = params.circleRadius * M_PI / 180 * - (params.latRange[1] - params.latRange[0]) / params.latPointCount; + (params.latRange[1] - params.latRange[0]) / (params.latPointCount - 1); const double avgObjectRes = DIM == 2 ? longSpacing : std::sqrt(longSpacing * longSpacing + latSpacing * latSpacing); @@ -1572,13 +1567,13 @@ int main(int argc, char** argv) getIntMinMax(int(query.searchCount()), minSearchCount, maxSearchCount, sumSearchCount); SLIC_INFO(axom::fmt::format( - "Initialization with policy {} took {{avg:{}, min:{}, max:{}}} seconds", + "Initialization with policy {} took {{avg:{:.5f}, min:{:.5f}, max:{:.5f}}} seconds", axom::runtime_policy::s_policyToName.at(params.policy), sumInit / num_ranks, minInit, maxInit)); SLIC_INFO(axom::fmt::format( - "Query with policy {} took {{avg:{}, min:{}, max:{}}} seconds", + "Query with policy {} took {{avg:{:.5f}, min:{:.5f}, max:{:.5f}}} seconds", axom::runtime_policy::s_policyToName.at(params.policy), sumQuery / num_ranks, minQuery, @@ -1589,7 +1584,7 @@ int main(int argc, char** argv) minSearchCount, maxSearchCount)); SLIC_INFO(axom::fmt::format( - "Effective distance threshold {{avg:{}, min:{}, max:{}}}", + "Effective distance threshold {{avg:{:.5f}, min:{:.5f}, max:{:.5f}}}", sumFilterDist / num_ranks, minFilterDist, maxFilterDist)); @@ -1642,6 +1637,7 @@ int main(int argc, char** argv) queryMeshWrapper.saveMesh(params.distanceFile); + axom::slic::flushStreams(); if(errCount) { SLIC_INFO(axom::fmt::format(" Error exit: {} errors found.", errCount));