diff --git a/src/axom/core/ArrayBase.hpp b/src/axom/core/ArrayBase.hpp index 40a6d65134..611cf6e08b 100644 --- a/src/axom/core/ArrayBase.hpp +++ b/src/axom/core/ArrayBase.hpp @@ -371,16 +371,26 @@ class ArrayBase updateStrides(); } - /// \brief Set the shape and 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_SUPPRESS_HD_WARN AXOM_HOST_DEVICE void setShapeAndStride(const StackArray& shape, - const StackArray& stride) + 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 164dad72b6..0c93695212 100644 --- a/src/axom/core/MDMapping.hpp +++ b/src/axom/core/MDMapping.hpp @@ -232,11 +232,13 @@ class MDMapping } 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"; - utilities::processAbort(); + << "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(); } // 2nd argument doesn't matter because strides are unique. diff --git a/src/axom/quest/DistributedClosestPoint.cpp b/src/axom/quest/DistributedClosestPoint.cpp index caa00be645..e95012c109 100644 --- a/src/axom/quest/DistributedClosestPoint.cpp +++ b/src/axom/quest/DistributedClosestPoint.cpp @@ -206,11 +206,12 @@ bool DistributedClosestPoint::generateBVHTree() void DistributedClosestPoint::computeClosestPoints(conduit::Node& query_node, const std::string& topology) { - SLIC_ASSERT_MSG(m_impl, "Must call 'setObjectMesh' before calling generateBVHTree"); + SLIC_ASSERT_MSG(m_impl, "Must call 'setObjectMesh' before calling computeClosestPoints"); SLIC_ASSERT(this->isValidBlueprint(query_node)); m_impl->setSquaredDistanceThreshold(m_sqDistanceThreshold); + m_impl->setFilterFarPartitions(m_filterFarPartitions); m_impl->setMpiCommunicator(m_mpiComm); m_impl->setOutputSwitches(m_outputRank, m_outputIndex, @@ -301,5 +302,12 @@ 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 af7227d799..b0a37f11b3 100644 --- a/src/axom/quest/DistributedClosestPoint.hpp +++ b/src/axom/quest/DistributedClosestPoint.hpp @@ -9,6 +9,7 @@ #include "axom/config.hpp" #include "axom/core/execution/runtime_policy.hpp" +// #include "axom/core/Types.hpp" #include "axom/slic.hpp" #include "conduit_node.hpp" @@ -104,6 +105,21 @@ class DistributedClosestPoint */ void setDistanceThreshold(double threshold); + /*! + * @brief Sets whether to filter out far partitions. + * \param [in] filterFarPartitions Filter out ranks that are too far + * to give productive results. + * + * Filtering examines the bounding boxes of all the partitions to + * exclude unproductive partition searches. Filtering is enabled by + * default. The benefits are highly configuration-dependent, but + * cost is typically negligible. + */ + void setFilterFarPartitions(bool filterFarPartitions) + { + m_filterFarPartitions = filterFarPartitions; + } + /*! @brief Set what fields to output. @@ -172,6 +188,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. @@ -199,6 +232,7 @@ class DistributedClosestPoint int m_dimension {-1}; bool m_isVerbose {false}; double m_sqDistanceThreshold; + bool m_filterFarPartitions {true}; bool m_outputRank = true; bool m_outputIndex = true; diff --git a/src/axom/quest/MeshViewUtil.hpp b/src/axom/quest/MeshViewUtil.hpp index d488924a6a..7ebc38867f 100644 --- a/src/axom/quest/MeshViewUtil.hpp +++ b/src/axom/quest/MeshViewUtil.hpp @@ -112,7 +112,8 @@ static void shapesToStridesAndOffsets(const axom::StackArray& realSh } /*! - * @brief Convert blueprint-style offsets and strides to shape specifications. + * @brief Convert blueprint-style offsets and strides to + * shape specifications. * * @tparam IType Index type * @tparam DIM Spatial dimension @@ -120,12 +121,17 @@ static void shapesToStridesAndOffsets(const axom::StackArray& realSh * @param realShape [i] * @param offsets [i] Blueprint-style index offsets. * @param strides [i] Blueprint-style strides. - * @param valuesCount [i] Number of values in ghost-padded data. + * @param valuesCount [i] Number of values in + * 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. - * @param minStride [i] Stride of fastest advancing index direction. - * @param strideOrder [i] Fastest-to-slowest advancing index directions. + * @param minStride [i] Stride of fastest advancing + * index direction. + * @param strideOrder [i] Fastest-to-slowest advancing + * index directions. */ template static void stridesAndOffsetsToShapes(const axom::StackArray& realShape, @@ -159,7 +165,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) { @@ -471,10 +479,10 @@ class MeshViewUtil /*! * @brief Return view to a scalar field variable. - * + * * @warning view returned has an allocator id determined by * \a MemSpace, regardless of the memory type. - * + * * @warning Assuming, without checking, that the field contains * data of type \a T. User is responsible for using the correct type. */ @@ -602,8 +610,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. * diff --git a/src/axom/quest/detail/DistributedClosestPointImpl.hpp b/src/axom/quest/detail/DistributedClosestPointImpl.hpp index c3d70854ab..f421a44661 100644 --- a/src/axom/quest/detail/DistributedClosestPointImpl.hpp +++ b/src/axom/quest/detail/DistributedClosestPointImpl.hpp @@ -250,7 +250,7 @@ class DistributedClosestPointImpl , m_mpiComm(MPI_COMM_NULL) , m_rank(-1) , m_nranks(-1) - , m_sqDistanceThreshold(axom::numeric_limits::max()) + , m_sqUserDistanceThreshold(axom::numeric_limits::max()) { } virtual ~DistributedClosestPointImpl() { } @@ -298,7 +298,17 @@ class DistributedClosestPointImpl void setSquaredDistanceThreshold(double sqThreshold) { SLIC_ERROR_IF(sqThreshold < 0.0, "Squared distance-threshold must be non-negative."); - m_sqDistanceThreshold = sqThreshold; + m_sqUserDistanceThreshold = sqThreshold; + } + + /*! + @brief Sets whether to filter partitions to eliminate fruitless searches. + + @param [in] filterFarPartitions Filter out unproductive searches. + */ + void setFilterFarPartitions(bool filterFarPartitions) + { + m_filterFarPartitions = filterFarPartitions; } /*! @@ -374,7 +384,7 @@ class DistributedClosestPointImpl { auto& queryDom = isMultidomain ? queryNode.child(domainNum) : queryNode; conduit::Node& xferDom = xferDoms.child(domainNum); - conduit::Node& fields = queryDom.fetch_existing("fields"); + conduit::Node& fields = queryDom.fetch("fields"); conduit::Node genericHeaders; genericHeaders["association"] = "vertex"; @@ -515,6 +525,17 @@ 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; @@ -523,7 +544,11 @@ class DistributedClosestPointImpl int m_rank; int m_nranks; - double m_sqDistanceThreshold; + //!@brief Distance threshold specified by user. + double m_sqUserDistanceThreshold; + + //!@brief Whether to prefilter domains to search. + bool m_filterFarPartitions = false; bool m_outputRank = true; bool m_outputIndex = true; @@ -542,6 +567,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; + //@} }; /*! @@ -607,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; @@ -637,23 +670,45 @@ class DistributedClosestPointExec : public DistributedClosestPointImpl const int N = internal::extractSize(copySrc); const std::size_t nBytes = sizeof(double) * DIM * N; - axom::copy(coords.data() + copiedCount, copySrc.fetch_existing("x").data_ptr(), nBytes); + axom::copy(m_objectPtCoords.data() + copiedCount, copySrc.fetch_existing("x").data_ptr(), nBytes); tmpValues.reset(); domIds.fill(domainId, N, copiedCount); copiedCount += N; } + + // 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(); + PointType minPt, maxPt; + for(int d = 0; d < DIM; ++d) + { + RAJA::ReduceMin minCoord(std::numeric_limits::max()); + RAJA::ReduceMax maxCoord(-std::numeric_limits::max()); + RAJA::forall( + RAJA::RangeSegment(0, ptCount), + AXOM_LAMBDA(RAJA::Index_type n) { + minCoord.min(coordsView[n][d]); + maxCoord.max(coordsView[n][d]); + }); + minPt[d] = minCoord.get(); + maxPt[d] = maxCoord.get(); + } + m_objectBb = BoxType(minPt, maxPt); +#else + m_objectBb = + axom::primal::BoundingBox {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, @@ -665,16 +720,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. @@ -684,7 +743,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, @@ -705,6 +763,17 @@ 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 { @@ -734,6 +803,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 { @@ -755,22 +837,72 @@ 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. + */ + + axom::Array allSqDistanceThreshold; + + if(m_filterFarPartitions) + { + AXOM_ANNOTATE_BEGIN("FilterFarPartitions"); + // Compute the min of the max distance between myQueryBb and each rank's object bounding box, + // and customize the distance threshold using this min-max. + double minMaxSqDist = std::numeric_limits::max(); + 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); + xferNodes[m_rank]->fetch("sqDistanceThreshold") = sqDistanceThreshold; + m_effectiveDistanceThreshold = sqrt(sqDistanceThreshold); + + allSqDistanceThreshold.resize(m_nranks); + AXOM_ANNOTATE_BEGIN("FilterFarPartitions_gather"); + gatherPrimitiveValue(sqDistanceThreshold, allSqDistanceThreshold); + AXOM_ANNOTATE_END("FilterFarPartitions_gather"); + AXOM_ANNOTATE_END("FilterFarPartitions"); + } + else + { + xferNodes[m_rank]->fetch("sqDistanceThreshold") = m_sqUserDistanceThreshold; + } + { conduit::Node& xferNode = *xferNodes[m_rank]; computeLocalClosestPoints(xferNode); + xferNode["searchCount"].set_int32(1); + m_searchCount = 1; } + /* + Count number of remote query partitions to receive, + i.e., how many 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 <= m_sqDistanceThreshold) + if(r != m_rank) { - ++remainingRecvs; + const auto& otherQueryBb = allQueryBbs[r]; + double sqDistance = axom::primal::squared_distance(otherQueryBb, myObjectBb); + if(sqDistance <= m_sqUserDistanceThreshold) + { + double sqDistanceThreshold = + m_filterFarPartitions ? allSqDistanceThreshold[r] : m_sqUserDistanceThreshold; + if(sqDistance <= sqDistanceThreshold) + { + ++remainingRecvs; + } + } } } } @@ -802,7 +934,9 @@ class DistributedClosestPointExec : public DistributedClosestPointImpl { isendRequests.emplace_back(conduit::relay::mpi::Request()); auto& req = isendRequests.back(); + AXOM_ANNOTATE_BEGIN("computeClosestPoints send"); relay::mpi::isend_using_schema(*xferNodes[m_rank], firstRecipForMyQuery, tag, m_mpiComm, &req); + AXOM_ANNOTATE_END("computeClosestPoints send"); ++remainingRecvs; } } @@ -814,7 +948,9 @@ class DistributedClosestPointExec : public DistributedClosestPointImpl // Receive the next xferNode std::shared_ptr recvXferNodePtr = std::make_shared(); + AXOM_ANNOTATE_BEGIN("computeClosestPoints recv"); conduit::relay::mpi::recv_using_schema(*recvXferNodePtr, MPI_ANY_SOURCE, tag, m_mpiComm); + AXOM_ANNOTATE_END("computeClosestPoints recv"); const int homeRank = recvXferNodePtr->fetch_existing("homeRank").as_int(); --remainingRecvs; @@ -824,28 +960,38 @@ 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["searchCount"].as_int32(); + xferNode["searchCount"].set_int32(1 + tmpCount); isendRequests.emplace_back(conduit::relay::mpi::Request()); auto& isendRequest = isendRequests.back(); int nextRecipient = next_recipient(xferNode); SLIC_ASSERT(nextRecipient != -1); + AXOM_ANNOTATE_BEGIN("computeClosestPoints send"); relay::mpi::isend_using_schema(xferNode, nextRecipient, tag, m_mpiComm, &isendRequest); + AXOM_ANNOTATE_END("computeClosestPoints send"); // Check non-blocking sends to free memory. + AXOM_ANNOTATE_BEGIN("computeClosestPoints check send req"); check_send_requests(isendRequests, false); + AXOM_ANNOTATE_END("computeClosestPoints check send req"); } } // remainingRecvs loop // Complete remaining non-blocking sends. + AXOM_ANNOTATE_BEGIN("computeClosestPoints check send req"); while(!isendRequests.empty()) { check_send_requests(isendRequests, true); } + AXOM_ANNOTATE_END("computeClosestPoints check send req"); MPI_Barrier(m_mpiComm); slic::flushStreams(); @@ -862,6 +1008,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; @@ -870,7 +1017,7 @@ class DistributedClosestPointExec : public DistributedClosestPointImpl return maybeNextRecip; } double sqDistance = primal::squared_distance(bb, m_objectPartitionBbs[maybeNextRecip]); - if(sqDistance <= m_sqDistanceThreshold) + if(sqDistance <= sqDistanceThreshold) { return maybeNextRecip; } @@ -880,31 +1027,9 @@ 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 { + AXOM_ANNOTATE_SCOPE("computeLocalClosestPoints()"); using axom::primal::squared_distance; // Note: There is some additional computation the first time this function @@ -992,7 +1117,7 @@ class DistributedClosestPointExec : public DistributedClosestPointImpl axom::Array sqDistThresh_host(1, 1, axom::execution_space::allocatorID()); - sqDistThresh_host[0] = m_sqDistanceThreshold; + sqDistThresh_host[0] = m_sqUserDistanceThreshold; axom::Array sqDistThresh_device = axom::Array(sqDistThresh_host, m_allocatorID); auto sqDistThresh_device_view = sqDistThresh_device.view(); @@ -1001,7 +1126,7 @@ class DistributedClosestPointExec : public DistributedClosestPointImpl auto ptDomainIdsView = m_objectPtDomainIds.view(); { - AXOM_ANNOTATE_SCOPE("ComputeClosestPoints"); + AXOM_ANNOTATE_SCOPE("computeLocalClosestPoints kernel"); axom::for_all( qPtCount, AXOM_LAMBDA(std::int32_t idx) mutable { @@ -1089,12 +1214,63 @@ 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. */ BoxArray m_objectPartitionBbs; std::unique_ptr m_bvh; + + /*! + @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. + + 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. And vice versa. + */ + int numCorners = 1 << DIM; + for(int i = 0; i < numCorners; ++i) + { + PointType aCoords; // i-th corner of a. + PointType bCoords; // Corner of b opposite from i-th corner of a. + for(int d = 0; d < DIM; ++d) + { + bool upperA_lowerB = i & (1 << d); + if(upperA_lowerB) + { + aCoords[d] = a.getMin()[d]; + bCoords[d] = b.getMax()[d]; + } + else // upperB_lowerA + { + bCoords[d] = b.getMin()[d]; + aCoords[d] = a.getMax()[d]; + } + } + primal::Vector separation {aCoords, bCoords}; + double sqDist = separation.squared_norm(); + maxSqDist = std::max(maxSqDist, sqDist); + } + return maxSqDist; + } + }; // DistributedClosestPointExec } // namespace internal 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 bb39bb93cf..9ca76c75cb 100644 --- a/src/axom/quest/examples/quest_distributed_distance_query_example.cpp +++ b/src/axom/quest/examples/quest_distributed_distance_query_example.cpp @@ -75,16 +75,23 @@ struct Input std::vector latRange {-90.0, 90.0}; int latPointCount {20}; + // Option to swap the 2 meshes for testing purposes. + bool swapMesh {false}; + RuntimePolicy policy {RuntimePolicy::seq}; double distThreshold {axom::numeric_limits::max()}; + bool filterFarPartitions {true}; + bool checkResults {false}; bool randomSpacing {true}; std::vector objDomainCountRange {1, 1}; + std::string annotationMode {"none"}; + private: bool m_verboseOutput {false}; @@ -161,11 +168,19 @@ struct Input ->description("Number of points in the latitudinal direction (3D only)") ->capture_default_str(); + app.add_flag("--swapMesh,!--no-swapMesh", swapMesh) + ->description("Swap the meshes (make spherical mesh into the query mesh)") + ->capture_default_str(); + app.add_option("-d,--dist-threshold", distThreshold) ->check(axom::CLI::NonNegativeNumber) ->description("Distance threshold to search") ->capture_default_str(); + app.add_option("-f,--filter-far-partitions", filterFarPartitions) + ->description("Whether to filter out partitions that are too far for productive searches") + ->capture_default_str(); + app.add_option("-p, --policy", policy) ->description("Set runtime policy for point query method") ->capture_default_str() @@ -175,6 +190,15 @@ struct Input ->description("Enable/disable checking results against analytical solution") ->capture_default_str(); +#ifdef AXOM_USE_CALIPER + app.add_option("--caliper", annotationMode) + ->description( + "caliper annotation mode. Valid options include 'none' and 'report'. " + "Use 'help' to see full list.") + ->capture_default_str() + ->check(axom::utilities::ValidCaliperMode); +#endif + app.get_formatter()->column_width(60); // could throw an exception @@ -325,6 +349,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. @@ -409,6 +436,37 @@ 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::fill(ownerRankView.data(), ownerRankView.size(), m_rank); + } + } + template axom::Array> getPoints(int domainIdx) { @@ -417,12 +475,11 @@ struct BlueprintParticleMesh auto* yView = cGroup->getView("values/y"); auto* zView = DIM >= 3 ? cGroup->getView("values/z") : nullptr; const auto ptCount = xView->getNumElements(); - assert(xView->getStride() == 1); - assert(yView->getStride() == 1); - assert(zView == nullptr || zView->getStride() == 1); - double* xs = xView->getArray(); - double* ys = yView->getArray(); - double* zs = zView ? (double*)(zView->getArray()) : nullptr; + axom::ArrayView xs {xView->getArray(), {xView->getNumElements()}, xView->getStride()}; + axom::ArrayView ys {yView->getArray(), {yView->getNumElements()}, yView->getStride()}; + axom::ArrayView zs = zView + ? axom::ArrayView {zView->getArray(), {zView->getNumElements()}, zView->getStride()} + : axom::ArrayView {}; using PointType = primal::Point; axom::Array pts; @@ -608,10 +665,12 @@ struct BlueprintParticleMesh sidre::View* yv = cvg->getView("y"); sidre::View* zv = ndim == 3 ? cvg->getView("z") : nullptr; axom::IndexType npts = xv->getNumElements(); - double* xp = xv->getData(); - double* yp = yv->getData(); - double* zp = zv ? (double*)(zv->getData()) : nullptr; - double* xyzs[3] {xp, yp, zp}; + axom::ArrayView xp {xv->getArray(), {xv->getNumElements()}, xv->getStride()}; + axom::ArrayView yp {yv->getArray(), {yv->getNumElements()}, yv->getStride()}; + axom::ArrayView zp = zv + ? axom::ArrayView {zv->getArray(), {zv->getNumElements()}, zv->getStride()} + : axom::ArrayView {}; + axom::ArrayView xyzs[3] {xp, yp, zp}; axom::Array rval(npts, ndim); for(int i = 0; i < npts; ++i) { @@ -872,6 +931,12 @@ class QueryMeshWrapper template int checkClosestPoints(const axom::primal::Sphere& sphere, const Input& params) { + if(params.swapMesh) + { + SLIC_INFO( + "Warning: Skipping checkClosestPoints, which doesn't work when meshes are swapped."); + return 0; + } using PointType = axom::primal::Point; m_queryMesh.registerNodalScalarField("error_flag"); @@ -895,11 +960,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 @@ -916,13 +976,38 @@ 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) { @@ -981,7 +1066,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, @@ -1010,6 +1095,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. */ @@ -1277,6 +1363,10 @@ int main(int argc, char** argv) umpire::Allocator umpireAllocator = rm.getAllocator(umpireResourceName); #endif + auto annotation_raii_wrapper = + std::make_unique(params.annotationMode); + AXOM_ANNOTATE_SCOPE("Quest distributed distance query example"); + // Storage for meshes. sidre::DataStore dataStore; @@ -1285,6 +1375,7 @@ int main(int argc, char** argv) // These will be used to query the closest points on the object mesh(es) //--------------------------------------------------------------------------- + AXOM_ANNOTATE_BEGIN("Main: set up example"); QueryMeshWrapper queryMeshWrapper(dataStore.getRoot()->createGroup("queryMesh", true), params.meshFile); @@ -1328,6 +1419,12 @@ int main(int argc, char** argv) } slic::flushStreams(); + if(params.swapMesh) + { + SLIC_INFO(axom::fmt::format("Swaping object and query meshes.")); + std::swap(objectMeshWrapper.getParticleMesh(), queryMeshWrapper.getParticleMesh()); + } + objectMeshWrapper.saveMesh(params.objectFile); slic::flushStreams(); @@ -1370,6 +1467,7 @@ int main(int argc, char** argv) make_coords_contiguous(dom.fetch_existing("coordsets/coords/values")); } } + AXOM_ANNOTATE_END("Main: set up example"); // Create distributed closest point query object and set some parameters quest::DistributedClosestPoint query; @@ -1380,6 +1478,7 @@ int main(int argc, char** argv) query.setMpiCommunicator(MPI_COMM_WORLD, true); query.setVerbosity(params.isVerbose()); query.setDistanceThreshold(params.distThreshold); + query.setFilterFarPartitions(params.filterFarPartitions); // To test support for single-domain format, use single-domain when possible. query.setObjectMesh(objectMeshNode.number_of_children() == 1 ? objectMeshNode[0] : objectMeshNode, objectMeshWrapper.getTopologyName()); @@ -1387,24 +1486,35 @@ int main(int argc, char** argv) // Build the spatial index over the object on each rank SLIC_INFO(init_str); slic::flushStreams(); + MPI_Barrier(MPI_COMM_WORLD); + AXOM_ANNOTATE_BEGIN("query.generateBVHTree()"); initTimer.start(); query.generateBVHTree(); initTimer.stop(); + AXOM_ANNOTATE_END("query.generateBVHTree()"); // Run the distributed closest point query over the nodes of the computational mesh // To test support for single-domain format, use single-domain when possible. slic::flushStreams(); + MPI_Barrier(MPI_COMM_WORLD); + AXOM_ANNOTATE_BEGIN("query.computeClosestPoints()"); queryTimer.start(); query.computeClosestPoints( queryMeshNode.number_of_children() == 1 ? queryMeshNode[0] : queryMeshNode, queryMeshWrapper.getTopologyName()); queryTimer.stop(); + AXOM_ANNOTATE_END("query.computeClosestPoints()"); auto getDoubleMinMax = [](double inVal, double& minVal, double& maxVal, double& sumVal) { MPI_Allreduce(&inVal, &minVal, 1, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD); 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 { @@ -1414,25 +1524,46 @@ 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:{:.5f}, min:{:.5f}, max:{:.5f}}} seconds", + axom::runtime_policy::s_policyToName.at(params.policy), + sumInit / num_ranks, + minInit, + maxInit)); SLIC_INFO( - axom::fmt::format("Initialization with policy {} took {{avg:{}, min:{}, max:{}}} seconds", + axom::fmt::format("Query 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", - axom::runtime_policy::s_policyToName.at(params.policy), - sumQuery / num_ranks, - minQuery, - maxQuery)); + sumQuery / num_ranks, + minQuery, + maxQuery)); + SLIC_INFO(axom::fmt::format("Search counts {{avg:{}, min:{}, max:{}}}", + double(sumSearchCount) / num_ranks, + minSearchCount, + maxSearchCount)); + SLIC_INFO( + axom::fmt::format("Effective distance threshold {{avg:{:.5f}, min:{:.5f}, max:{:.5f}}}", + sumFilterDist / num_ranks, + minFilterDist, + maxFilterDist)); } 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.")); + AXOM_ANNOTATE_SCOPE("Main: check results"); + slic::flushStreams(); if(spatialDim == 2) { primal::Point center(params.circleCenter.data()); @@ -1470,6 +1601,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)); @@ -1480,6 +1612,7 @@ int main(int argc, char** argv) } finalizeLogger(); + annotation_raii_wrapper.reset(); MPI_Finalize(); return errCount != 0;